! moduleInput: Reads JSON configuration file
MODULE moduleInput
  USE json_module
  IMPLICIT NONE

  CONTAINS
    !Main routine to read the input JSON file
    SUBROUTINE readConfig(inputFile)
      USE json_module
      USE moduleErrors
      USE moduleBoundary
      USE moduleInject
      IMPLICIT NONE

      CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
      TYPE(json_file):: config

      !Initialize the json file variable
      CALL config%initialize()

      !Loads the config file
      CALL verboseError('Loading input file...')
      CALL config%load(filename = inputFile)
      CALL checkStatus(config, "load")

      !Reads reference parameters
      CALL verboseError('Reading Reference parameters...')
      CALL readReference(config)
      CALL checkStatus(config, "readReference")

      !Reads output parameters
      CALL verboseError('Reading Output parameters...')
      CALL readOutput(config)
      CALL checkStatus(config, "readOutput")

      !Read species
      CALL verboseError('Reading species information...')
      CALL readSpecies(config)
      CALL checkStatus(config, "readSpecies")

      !Read interactions between species
      CALL verboseError('Reading interaction between species...')
      CALL readInteractions(config)
      CALL checkStatus(config, "readInteractions")

      !Read boundaries
      CALL verboseError('Reading boundary conditions...')
      CALL readBoundary(config)
      CALL checkStatus(config, "readBoundary")

      !Read Geometry
      CALL verboseError('Reading Geometry...')
      CALL readGeometry(config)
      CALL checkStatus(config, "readGeometry")

      !Reads case parameters
      CALL verboseError('Reading Case parameters...')
      CALL readCase(config)
      CALL checkStatus(config, "readCase")

      !Read injection of particles
      CALL verboseError('Reading Interactions between species...')
      CALL readInject(config)
      CALL checkStatus(config, "readInject")

      !Read parallel parameters
      CALL verboseError('Reading Parallel configuration...')
      CALL readParallel(config)
      CALL checkStatus(config, "readParallel")

    END SUBROUTINE readConfig
    
    !Checks the status of the JSON case file and, if failed, exits the execution.
    SUBROUTINE checkStatus(config, step)
      USE moduleErrors
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(LEN=*), INTENT(in):: step

      IF (config%failed()) CALL criticalError("Error reading the JSON input file", TRIM(step))

    END SUBROUTINE checkStatus

    !Reads the reference parameters
    SUBROUTINE readReference(config)
      USE moduleRefParam
      USE moduleConstParam
      USE moduleErrors
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: object

      object = 'reference'
      !Mandatory parameters that define the case and computes derived parameters
      CALL config%get(object // '.density', n_ref, found)
      IF (.NOT. found) CALL criticalError('Reference density not found','readReference')

      CALL config%get(object // '.mass', m_ref, found)
      IF (.NOT. found) CALL criticalError('Reference mass not found','readReference')

      CALL config%get(object // '.temperature', T_ref, found)
      IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')

      !If a reference cross section is given, it is used
      CALL config%get(object // '.crossSection', sigma_ref, found)

      !If not, the reference radius is searched
      IF (.NOT. found) THEN
        CALL config%get(object // '.radius', r_ref, found)
        IF (found) THEN
          sigma_ref = PI*(r_ref+r_ref)**2    !reference cross section

        ELSE
          sigma_ref = 0.D0 !Assume no collisions

        END IF

      END IF
      
      !Derived parameters
      L_ref     = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length
      v_ref     = DSQRT(kb*T_ref/m_ref)          !reference velocity
      ti_ref    = L_ref/v_ref                    !reference time
      Vol_ref   = L_ref**3                       !reference volume
      EF_ref    = qe*n_ref*L_ref/eps_0           !reference electric field
      Volt_ref  = EF_ref*L_ref                   !reference voltage

    END SUBROUTINE readReference

    !Reads the specific case parameters
    SUBROUTINE readCase(config)
      USE moduleRefParam
      USE moduleErrors
      USE moduleCaseParam
      USE moduleSolver
      USE moduleSpecies
      USE moduleCollisions
      USE moduleOutput
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: object
      REAL(8):: time !simulation time in [t]
      CHARACTER(:), ALLOCATABLE:: pusherType, EMType, WSType
      INTEGER:: nTau, nSolver
      INTEGER:: i
      CHARACTER(2):: iString
      CHARACTER(1):: tString

      object = 'case'

      !Time parameters
      CALL config%info(object // '.tau', found, n_children = nTau)
      IF (.NOT. found .OR. nTau == 0) CALL criticalError('Required parameter tau not found','readCase')
      ALLOCATE(tau(1:nSpecies))
      DO i = 1, nTau
        WRITE(iString, '(I2)') i
        CALL config%get(object // '.tau(' // TRIM(iString) // ')', tau(i), found)

      END DO
      IF (nTau < nSpecies) THEN
        CALL warningError('Using minimum time step for some species')
        tau(nTau+1:nSpecies) = MINVAL(tau(1:nTau))

      END IF
      tauMin = MINVAL(tau)

      !Gets the simulation time
      CALL config%get(object // '.time', time, found)
      IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase')
      !Convert simulation time to number of iterations
      tmax = INT(time/tauMin)

      !Gest the pusher for each species
      CALL config%info(object // '.pusher', found, n_children = nSolver)
      IF (.NOT. found .OR. nSolver /= nSpecies) CALL criticalError('Required parameter pusher not found','readCase')
      !Allocates all the pushers for particles
      ALLOCATE(solver%pusher(1:nSpecies))
      !Initialize pushers
      DO i = 1, nSolver
        WRITE(iString, '(I2)') i
        CALL config%get(object // '.pusher(' // TRIM(iString) // ')', pusherType, found)

        !Associate the type of solver
        CALL solver%pusher(i)%init(pusherType, tauMin, tau(i))

      END DO
      
      !Gets the solver for the electromagnetic field
      CALL config%get(object // '.EMSolver', EMType, found)
      CALL solver%initEM(EMType)
      SELECT CASE(EMType)
      CASE("Electrostatic")
        CALL readEMBoundary(config)

      END SELECT

      !Gest the non-analogue scheme
      CALL config%get(object // '.WeightingScheme', WSType, found)
      CALL solver%initWS(WSType)


      !Makes tau(s) non-dimensional
      tau = tau / ti_ref
      tauMin = tauMin / ti_ref

      !Sets the format of output files accordint to iteration number
      iterationDigits = INT(LOG10(REAL(tmax))) + 1
      WRITE(tString, '(I1)') iterationDigits
      iterationFormat = "(I" // tString // "." // tString // ")"

      !Read initial state for species
      CALL verboseError('Reading Initial state...')
      CALL readInitial(config)
      CALL checkStatus(config, "readInitial")

    END SUBROUTINE readCase

    !Reads the initial information for the species
    SUBROUTINE readInitial(config)
      USE moduleSpecies
      USE moduleMesh
      USE moduleOutput
      USE moduleRefParam
      USE moduleRandom
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: object
      INTEGER:: nInitial
      INTEGER:: i, p, e
      CHARACTER(LEN=2):: iString
      CHARACTER(:), ALLOCATABLE:: spName
      INTEGER:: sp
      CHARACTER(:), ALLOCATABLE:: spFile
      INTEGER:: stat
      CHARACTER(100):: dummy
      REAL(8):: density, velocity(1:3), temperature
      INTEGER:: nNewPart = 0.D0
      TYPE(particle), POINTER:: partNew
      REAL(8):: vTh
      TYPE(lNode), POINTER:: partCurr, partNext

      CALL config%info('case.initial', found, n_children = nInitial)

      IF (found) THEN
        !Reads the information from initial species
        DO i = 1, nInitial
          WRITE(iString, '(I2)') i
          object = 'case.initial(' // iString // ')'
          CALL config%get(object // '.speciesName', spName, found)
          sp = speciesName2Index(spName)
          CALL config%get(object // '.initialState', spFile, found)
          OPEN (10, FILE = path // spFile, ACTION = 'READ')
          DO
            READ(10, '(A)', IOSTAT = stat) dummy
            !If EoF, exit reading
            IF (stat /= 0) EXIT
            !If comment, skip
            IF (INDEX(dummy,'#') /= 0) CYCLE
            !Go up one line
            BACKSPACE(10)
            !Read information
            READ(10, *) e, density, velocity, temperature
            !Scale variables
            !Particles in cell volume
            nNewPart = INT(density * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
            !Non-dimensional velocity
            velocity = velocity / v_ref
            !Non-dimensional temperature
            temperature = temperature / T_ref
            !Non-dimensional thermal temperature
            vTh = DSQRT(temperature/species(sp)%obj%m)
            !Allocate new particles
            DO p = 1, nNewPart
              ALLOCATE(partNew)
              partNew%sp     = sp
              partNew%v(1)   = velocity(1) + vTh*randomMaxwellian()
              partNew%v(2)   = velocity(2) + vTh*randomMaxwellian()
              partNew%v(3)   = velocity(3) + vTh*randomMaxwellian()
              partNew%vol    = e
              partNew%r      = mesh%vols(e)%obj%randPos()
              partNew%xi     = mesh%vols(e)%obj%phy2log(partNew%r)
              partNew%n_in   = .TRUE.
              partNew%weight = species(sp)%obj%weight
              !If charged species, add qm to particle
              SELECT TYPE(sp => species(sp)%obj)
              TYPE IS (speciesCharged)
                partNew%qm = sp%qm

              CLASS DEFAULT
                partNew%qm = 0.D0

              END SELECT

              !Assign particle to temporal list of particles
              CALL partInitial%add(partNew)

            END DO

          END DO

        END DO

        !Convert temporal list of particles into initial partOld array
        !Deallocates the list of initial particles
        nNewPart = partInitial%amount
        IF (nNewPart > 0) THEN
          ALLOCATE(partOld(1:nNewPart))
          partCurr => partInitial%head
          DO p = 1, nNewPart
            partNext => partCurr%next
            partOld(p) = partCurr%part
            DEALLOCATE(partCurr)
            partCurr => partNext

          END DO

          IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head)
          IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail)
          partInitial%amount = 0

        END IF

        nPartOld = SIZE(partOld)

      END IF

    END SUBROUTINE readInitial

    !Reads configuration for the output files
    SUBROUTINE readOutput(config)
      USE moduleErrors
      USE moduleOutput
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: object
      CHARACTER(:), ALLOCATABLE:: baseName
      CHARACTER(8)  :: date_now=''
      CHARACTER(10) :: time_now=''

      object = 'output'
      CALL config%get(object // '.path', path, found)
      CALL config%get(object // '.triggerOutput', triggerOutput, found)
      IF (.NOT. found) THEN
        triggerOutput = 100
        CALL warningError('Using default trigger for output file of 100 iterations')

      END IF

      !Gets actual date and time
      CALL DATE_AND_TIME(date_now, time_now)

      !Gets the basename of the folder
      CALL config%get(object // '.folder', baseName, found)
      PRINT *, baseName
      IF (found) THEN
        folder = baseName

      END IF

      !Compose the folder name
      folder = folder // '_' // date_now(1:4) // '-' // date_now(5:6) // '-' // date_now(7:8) // '_' &
                             // time_now(1:2) // '.' // time_now(3:4) // '.' // time_now(5:6)
      
      !Creates the folder
      CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )

      CALL config%get(object // '.cpuTime', timeOutput, found)
      CALL config%get(object // '.numColl', collOutput, found)
      CALL config%get(object // '.EMField', emOutput,   found)

      CALL config%get(object // '.triggerCPUTime', triggerCPUTime, found)
      IF (.NOT. found) THEN
        triggerCPUTime = triggerOutput
        IF (timeOutput) CALL warningError('No triggerCPUTime found, using same vale as triggerOutput')

      END IF

    END SUBROUTINE readOutput

    !Reads information about the case species
    SUBROUTINE readSpecies(config)
      USE moduleSpecies
      USE moduleErrors
      USE moduleRefParam
      USE moduleList
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(2):: iString
      CHARACTER(:), ALLOCATABLE:: object
      CHARACTER(:), ALLOCATABLE:: speciesType
      REAL(8):: mass, charge
      LOGICAL:: found
      INTEGER:: i
      CHARACTER(:), ALLOCATABLE:: linkName
      INTEGER:: linkID

      !Gets the number of species
      CALL config%info('species', found, n_children = nSpecies)
      !Zero species means critical error
      IF (nSpecies == 0) CALL criticalError("No species found", "configRead")

      ALLOCATE(species(1:nSpecies))

      !Reads information of individual species
      DO i = 1, nSpecies
        WRITE(iString, '(I2)') i
        object = 'species(' // TRIM(iString) // ')'
        CALL config%get(object // '.type', speciesType, found)
        CALL config%get(object // '.mass', mass,        found)
        mass = mass/m_ref

        !Allocate species depending on type and assign specific parameters
        SELECT CASE(speciesType)
        CASE ("neutral")
          ALLOCATE(species(i)%obj, source=speciesNeutral())

        CASE ("charged")
          CALL config%get(object // '.charge', charge,      found)
          IF (.NOT. found) CALL criticalError("Required parameter charge not found for species " // object, 'readSpecies')
          ALLOCATE(species(i)%obj, source=speciesCharged(q  = charge, &
                                                         qm = charge/mass))

        CASE DEFAULT
          CALL warningError("Species " // speciesType // " not supported yet")

        END SELECT
        !Assign shared parameters for all species
        CALL config%get(object // '.name',   species(i)%obj%name,   found)
        CALL config%get(object // '.weight', species(i)%obj%weight, found)
        species(i)%obj%sp = i
        species(i)%obj%m  = mass

      END DO

      !Reads relations between species
      DO i = 1, nSpecies
        WRITE(iString, '(I2)') i
        object = 'species(' // TRIM(iString) // ')'
        SELECT TYPE(sp => species(i)%obj)
        TYPE IS (speciesNeutral)
          !Gets species linked ion
          CALL config%get(object // '.ion', linkName, found)
          IF (found) THEN
            linkID =  speciesName2Index(linkName)
            sp%ion => species(linkID)%obj

          END IF

        TYPE IS (speciesCharged)
          !Gets species linked neutral
          CALL config%get(object // '.neutral', linkName, found)
          IF (found) THEN
            linkID     =  speciesName2Index(linkName)
            sp%neutral => species(linkID)%obj

          END IF

          !Gets species linked ion
          CALL config%get(object // '.ion', linkName, found)
          IF (found) THEN
            linkID =  speciesName2Index(linkName)
            sp%ion => species(linkID)%obj

          END IF

        END SELECT

      END DO

      !Set number of particles to 0 for init state
      !TODO: In a future, this should include the particles from init states
      nPartOld = 0

      !Initialize the lock for the non-analogue (NA) list of particles
      CALL OMP_INIT_LOCK(lockWScheme)

    END SUBROUTINE readSpecies

    !Reads information about interactions between species
    SUBROUTINE readInteractions(config)
      USE moduleSpecies
      USE moduleList
      USE moduleCollisions
      USE moduleErrors
      USE OMP_LIB
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(2):: iString, kString
      CHARACTER(:), ALLOCATABLE:: object
      CHARACTER(:), ALLOCATABLE:: species_i, species_j
      CHARACTER(:), ALLOCATABLE:: crossSecFile
      CHARACTER(:), ALLOCATABLE:: crossSecFilePath
      CHARACTER(:), ALLOCATABLE:: cType
      LOGICAL:: found
      INTEGER:: nInteractions, nCollisions
      INTEGER:: i, k, ij
      INTEGER:: pt_i, pt_j
      REAL(8):: energyThreshold, energyBinding
      CHARACTER(:), ALLOCATABLE:: electron

      CALL initInteractionMatrix(interactionMatrix)

      !Path for collision cross-section data files
      CALL config%get('interactions.folderCollisions', pathCollisions, found)

      !Inits lock for list of particles
      CALL OMP_INIT_LOCK(lockCollisions)

      CALL config%info('interactions.collisions', found, n_children = nInteractions)
      DO i = 1, nInteractions
        WRITE(iString, '(I2)') i
        object = 'interactions.collisions(' // TRIM(iString) // ')'
        CALL config%get(object // '.species_i', species_i, found)
        pt_i = speciesName2Index(species_i)
        CALL config%get(object // '.species_j', species_j, found)
        pt_j = speciesName2Index(species_j)
        CALL config%info(object // '.cTypes', found, n_children = nCollisions)
        ij = interactionIndex(pt_i,pt_j)
        !Allocates the required number of collisions per each pair of species ij
        CALL interactionMatrix(ij)%init(nCollisions)

        DO k = 1, nCollisions
          WRITE (kString, '(I2)') k
          object = 'interactions.collisions(' // TRIM(iString) // ').cTypes(' // TRIM(kString) // ')'
          !Reads the cross section file
          CALL config%get(object // '.crossSection', crossSecFile, found)
          crossSecFilePath = pathCollisions // crossSecFile
          IF (.NOT. found) CALL criticalError('crossSection not found for ' // object, 'readInteractions')
          !Reads the type of collision
          CALL config%get(object // '.type', cType, found)
          !Initialize collision type and reads required additional data
          SELECT CASE(cType)
          CASE ('elastic')
            !Elastic collision
            CALL initBinaryElastic(interactionMatrix(ij)%collisions(k)%obj, &
                                   crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)

          CASE ('chargeExchange')
            !Resonant charge exchange
            CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, &
                                          crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)

          CASE ('ionization')
            !Electorn impact ionization
            CALL config%get(object // '.energyThreshold', energyThreshold, found)
            IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
            CALL config%get(object // '.electron', electron, found)
            IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
            CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
                                      crossSecFilePath, energyThreshold, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)

          CASE ('recombination')
            !Electorn impact ionization
            CALL config%get(object // '.energyBinding', energyBinding, found)
            IF (.NOT. found) CALL criticalError('energyThreshold not found for collision' // object, 'readInteractions')
            CALL config%get(object // '.electron', electron, found)
            IF (.NOT. found) CALL criticalError('electron not found for collision' // object, 'readInteractions')
            CALL initBinaryRecombination(interactionMatrix(ij)%collisions(k)%obj, &
                                      crossSecFilePath, energyBinding, species(pt_i)%obj%m, species(pt_j)%obj%m, electron)

          CASE DEFAULT
            CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions')

          END SELECT

        END DO

      END DO

    END SUBROUTINE readInteractions

    !Reads boundary conditions for the mesh
    SUBROUTINE readBoundary(config)
      USE moduleBoundary
      USE moduleErrors
      USE moduleSpecies
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      INTEGER:: i, s
      CHARACTER(2):: istring, sString
      CHARACTER(:), ALLOCATABLE:: object, bType
      REAL(8):: Tw, cw !Wall temperature and specific heat
      LOGICAL:: found
      INTEGER:: nTypes

      CALL config%info('boundary', found, n_children = nBoundary)
      ALLOCATE(boundary(1:nBoundary))
      DO i = 1, nBoundary
        WRITE(istring, '(i2)') i
        object = 'boundary(' // TRIM(istring) // ')'

        boundary(i)%id = i
        CALL config%get(object // '.name',            boundary(i)%name,            found)
        CALL config%get(object // '.physicalSurface', boundary(i)%physicalSurface, found)
        CALL config%info(object // '.bTypes', found, n_children = nTypes)
        IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary')
        ALLOCATE(boundary(i)%bTypes(1:nSpecies))
        DO s = 1, nSpecies
          WRITE(sString,'(i2)') s
          object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
          CALL config%get(object // '.type', bType, found)
          SELECT CASE(bType)
          CASE('reflection')
            ALLOCATE(boundaryReflection:: boundary(i)%bTypes(s)%obj)

          CASE('absorption')
            ALLOCATE(boundaryAbsorption:: boundary(i)%bTypes(s)%obj)

          CASE('transparent')
            ALLOCATE(boundaryTransparent:: boundary(i)%bTypes(s)%obj)

          CASE('wallTemperature')
            CALL config%get(object // '.temperature', Tw, found)
            IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
            CALL config%get(object // '.specificHeat', cw, found)
            IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')

            CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw)

          CASE('axis')
            ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)

          CASE DEFAULT
            CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')

          END SELECT

        END DO

      END DO

    END SUBROUTINE readBoundary

    !Read the geometry (mesh) for the case
    SUBROUTINE readGeometry(config)
      USE moduleMesh
      USE moduleMeshCylRead,    ONLY: meshCylGeneric
      USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric
      USE moduleMesh1DRadRead,  ONLY: mesh1DRadGeneric
      USE moduleErrors
      USE moduleOutput
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: geometryType, meshFormat, meshFile
      CHARACTER(:), ALLOCATABLE:: fullPath

      !Selects the type of geometry.
      CALL config%get('geometry.type', geometryType, found)
      SELECT CASE(geometryType)
      CASE ("2DCyl")
        !Creates a 2D cylindrical mesh
        ALLOCATE(meshCylGeneric:: mesh)

      CASE ("1DCart")
        !Creates a 1D cartesian mesh
        ALLOCATE(mesh1DCartGeneric:: mesh)

      CASE ("1DRad")
        !Creates a 1D cartesian mesh
        ALLOCATE(mesh1DRadGeneric:: mesh)

      CASE DEFAULT
        CALL criticalError("Geometry type " // geometryType // " not supported.", "readGeometry")

      END SELECT

      !Gets the type of mesh
      CALL config%get('geometry.meshType', meshFormat, found)
      CALL mesh%init(meshFormat)
      !Reads the mesh
      CALL config%get('geometry.meshFile', meshFile, found)
      fullpath = path // meshFile
      CALL mesh%readMesh(fullPath)

    END SUBROUTINE readGeometry

    SUBROUTINE readEMBoundary(config)
      USE moduleMesh
      USE moduleOutput
      USE moduleErrors
      USE moduleEM
      USE moduleRefParam
      USE moduleSpecies
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found
      CHARACTER(2):: istring
      INTEGER:: i, e, s
      INTEGER:: info
      EXTERNAL:: dgetrf

      CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
      
      IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
      
      DO i = 1, nBoundaryEM
        WRITE(istring, '(i2)') i
        object = 'boundaryEM(' // trim(istring) // ')'

        CALL config%get(object // '.type', boundEM(i)%typeEM, found)

        SELECT CASE(boundEM(i)%typeEM)
        CASE ("dirichlet")
          CALL config%get(object // '.potential', boundEM(i)%potential, found)
          IF (.NOT. found) &
            CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
          boundEM(i)%potential = boundEM(i)%potential/Volt_ref

          CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
          IF (.NOT. found) &
            CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')

        CASE DEFAULT
          CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')

        END SELECT

      END DO

      ALLOCATE(qSpecies(1:nSpecies))
      DO s = 1, nSpecies
        SELECT TYPE(sp => species(s)%obj)
        TYPE IS (speciesCharged)
          qSpecies(s) = sp%q

        CLASS DEFAULT
          qSpecies(s) = 0.D0

        END SELECT

      END DO

      !TODO: Improve this
      IF (ALLOCATED(boundEM)) THEN
        DO e = 1, mesh%numEdges
          IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
            DO i = 1, nBoundaryEM
              IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
                CALL boundEM(i)%apply(mesh%edges(e)%obj)

              END IF
            END DO
          END IF
        END DO
      END IF

      !Compute the PLU factorization of K once boundary conditions have been read
      CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
      IF (info /= 0) CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')

    END SUBROUTINE readEMBoundary

    !Reads the injection of particles from the boundaries
    SUBROUTINE readInject(config)
      USE moduleSpecies
      USE moduleErrors
      USE moduleInject
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      INTEGER:: i
      CHARACTER(2):: istring
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: speciesName
      CHARACTER(:), ALLOCATABLE:: name
      REAL(8):: v
      REAL(8), ALLOCATABLE:: T(:), normal(:)
      REAL(8):: flow
      CHARACTER(:), ALLOCATABLE:: units
      INTEGER:: physicalSurface
      INTEGER:: sp

      CALL config%info('inject', found, n_children = nInject)
      ALLOCATE(inject(1:nInject))
      nPartInj = 0
      DO i = 1, nInject
        WRITE(istring, '(i2)') i
        object = 'inject(' // trim(istring) // ')'

        !Find species
        CALL config%get(object // '.species', speciesName, found)
        sp = speciesName2Index(speciesName)

        CALL config%get(object // '.name',            name,            found)
        CALL config%get(object // '.v',               v,               found)
        CALL config%get(object // '.T',               T,               found)
        CALL config%get(object // '.n',               normal,          found)
        CALL config%get(object // '.flow',            flow,            found)
        CALL config%get(object // '.units',           units,           found)
        CALL config%get(object // '.physicalSurface', physicalSurface, found)

        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)

        CALL readVelDistr(config, inject(i), object)
        
      END DO

      !Allocate array for injected particles
      IF (nPartInj > 0) THEN
        ALLOCATE(partInj(1:nPartInj))

      END IF
      
    END SUBROUTINE readInject

    !Reads the velocity distribution functions for each inject
    SUBROUTINE readVelDistr(config, inj, object)
      USE moduleErrors
      USE moduleInject
      USE moduleSpecies
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      TYPE(injectGeneric), INTENT(inout):: inj
      CHARACTER(:), ALLOCATABLE, INTENT(in):: object
      INTEGER:: i
      CHARACTER(2):: istring
      CHARACTER(:), ALLOCATABLE:: fvType
      LOGICAL:: found
      REAL(8):: v, T, m

      !Reads species mass
      m = species(inj%sp)%obj%m
      !Reads distribution functions for velocity
      DO i = 1, 3
        WRITE(istring, '(i2)') i
        CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found)
        IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // &
                                            " found for " // object, 'readVelDistr')

        SELECT CASE(fvType)
        CASE ("Maxwellian")
          v = inj%vMod*inj%n(i)
          T = inj%T(i)
          CALL initVelDistMaxwellian(inj%v(i)%obj, v, t, m)

        CASE ("Delta")
          v = inj%vMod*inj%n(i)
          CALL initVelDistDelta(inj%v(i)%obj, v)

        CASE DEFAULT
          CALL criticalError("No velocity distribution type " // fvType // " defined", 'readVelDistr')

        END SELECT

      END DO

    END SUBROUTINE readVelDistr

    !Reads configuration for parallel run
    SUBROUTINE readParallel(config)
      USE moduleParallel
      USE moduleErrors
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found

      !Reads OpenMP parameters
      object = 'parallel.OpenMP'

      CALL config%get(object // '.nThreads', openMP%nThreads, found)

      IF (.NOT. found) THEN
        openMP%nthreads = 8
        CALL warningError('No OpenMP configuration detected, using 8 threads as default.')

      END IF

      CALL initParallel()

    END SUBROUTINE readParallel


END MODULE moduleInput