! 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
      USE moduleOutput
      USE moduleMesh
      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, "loading")

      !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 injection of particles from boundaries...')
      CALL readInject(config)
      CALL checkStatus(config, "readInject")

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

      !If everything is correct, creates the output folder
      CALL EXECUTE_COMMAND_LINE('mkdir ' // path // folder )
      !Copies input file to output folder
      CALL EXECUTE_COMMAND_LINE('cp ' // inputFile // ' ' // path // folder)
      !Copies particle mesh
      CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshParticle // ' ' // path // folder)
      IF (doubleMesh) THEN
        CALL EXECUTE_COMMAND_LINE('cp ' // pathMeshColl // ' ' // path // folder)

      END IF

    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
      !simulation final and initial times in [t]
      REAL(8):: finalTime, initialTime
      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 final time
      CALL config%get(object // '.finalTime', finalTime, found)
      IF (.NOT. found) CALL criticalError('Required parameter finalTime not found','readCase')
      !Convert simulation time to number of iterations
      tFinal = INT(finalTime / tauMin)

      !Gets the simulation initial time
      CALL config%get(object // '.initialTime', initialTime, found)
      IF (found) tInitial = INT(initialTime / 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(tFinal))) + 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, j, p, e
      CHARACTER(LEN=2):: iString
      CHARACTER(:), ALLOCATABLE:: spName
      INTEGER:: sp
      CHARACTER(:), ALLOCATABLE:: spFile
      CHARACTER(:), ALLOCATABLE:: filename
      REAL(8), ALLOCATABLE, DIMENSION(:):: density, temperature
      REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity
      INTEGER, ALLOCATABLE, DIMENSION(:):: nodes
      INTEGER:: nNodes
      REAL(8), ALLOCATABLE, DIMENSION(:):: source, fPsi
      !Density at the volume centroid
      REAL(8):: densityCen
      !Mean velocity and temperature at particle position
      REAL(8):: velocityXi(1:3), temperatureXi
      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 // '.species', spName, found)
          sp = speciesName2Index(spName)
          CALL config%get(object // '.file', spFile, found)
          !Reads node values at the nodes
          filename = path // spFile
          CALL mesh%readInitial(sp, filename, density, velocity, temperature)
          !For each volume in the node, create corresponding particles
          DO e = 1, mesh%numVols
            !Scale variables
            !Density at centroid of cell
            nodes = mesh%vols(e)%obj%getNodes()
            nNodes = SIZE(nodes)
            !TODO: Procedure to obtain centroid from element (also for printing Electric Field)
            fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/))
            ALLOCATE(source(1:nNodes))
            DO j = 1, nNodes
              source(j) = density(nodes(j))

            END DO
            densityCen = DOT_PRODUCT(fPsi, source)
            DEALLOCATE(fPsi)

            !Calculate number of particles
            nNewPart = INT(densityCen * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)

            !Allocate new particles
            DO p = 1, nNewPart
              ALLOCATE(partNew)
              partNew%species => species(sp)%obj
              partNew%r       =  mesh%vols(e)%obj%randPos()
              partNew%xi      =  mesh%vols(e)%obj%phy2log(partNew%r)
              !Get mean velocity at particle position
              fPsi            =  mesh%vols(e)%obj%fPsi(partNew%xi)
              DO j = 1, nNodes
                source(j) = velocity(nodes(j), 1)

              END DO
              velocityXi(1) = DOT_PRODUCT(fPsi, source)
              DO j = 1, nNodes
                source(j) = velocity(nodes(j), 2)

              END DO
              velocityXi(2) = DOT_PRODUCT(fPsi, source)
              DO j = 1, nNodes
                source(j) = velocity(nodes(j), 3)

              END DO
              velocityXi(3) = DOT_PRODUCT(fPsi, source)

              velocityXi = velocityXi / v_ref

              !Get temperature at particle position
              DO j = 1, nNodes
                source(j) = temperature(nodes(j))

              END DO
              temperatureXi = DOT_PRODUCT(fPsi, source)
              temperatureXi = temperatureXi / T_ref

              vTh = DSQRT(temperatureXi / species(sp)%obj%m)
              partNew%v(1)    = velocityXi(1) + vTh*randomMaxwellian()
              partNew%v(2)    = velocityXi(2) + vTh*randomMaxwellian()
              partNew%v(3)    = velocityXi(3) + vTh*randomMaxwellian()
              partNew%vol     = e
              IF (ASSOCIATED(meshForMCC, mesh)) THEN
                partNew%volColl = partNew%vol

              ELSE
                partNew%volColl = findCellBrute(meshColl, partNew%r)

              END IF
              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

            DEALLOCATE(source)

          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)
      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)
      
      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%n = 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
      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 moduleMesh
      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

      !Firstly, checks if the object 'interactions' exists
      CALL config%info('interactions', found)
      IF (found) THEN
        !Checks if MC collisions have been defined
        CALL config%info('interactions.collisions', found)
        IF (found) THEN
          !Checks if a mesh for collisions has been defined
          !The mesh will be initialized and reader in readGeometry
          CALL config%info('interactions.meshCollisions', found)
          IF (found) THEN
            !Points meshForMCC to the specific mesh defined
            meshForMCC => meshColl

          ELSE
            !Points the meshForMCC pointer to the Particles Mesh
            meshForMCC => mesh

          END IF

          !Inits the MCC matrix
          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 IF

      END IF

    END SUBROUTINE readInteractions

    !Reads boundary conditions for the mesh
    SUBROUTINE readBoundary(config)
      USE moduleBoundary
      USE moduleErrors
      USE moduleSpecies
      USE moduleRefParam
      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
      !Neutral Properties
      REAL(8):: m0, n0, T0
      REAL(8), DIMENSION(:), ALLOCATABLE:: v0
      REAL(8):: effTime
      REAL(8):: eThreshold !Energy threshold
      INTEGER:: speciesID
      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection
      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)%n = 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('ionization')
            !Neutral parameters
            CALL config%get(object // '.neutral.ion', speciesName, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
            speciesID = speciesName2Index(speciesName)
            CALL config%get(object // '.neutral.mass', m0, found)
            IF (.NOT. found) THEN
              m0 = species(s)%obj%m*m_ref
            END IF
            CALL config%get(object // '.neutral.density', n0, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
            CALL config%get(object // '.neutral.velocity', v0, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
            CALL config%get(object // '.neutral.temperature', T0, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')

            CALL config%get(object // '.effectiveTime', effTime, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')

            CALL config%get(object // '.energyThreshold', eThreshold, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')

            CALL config%get(object // '.crossSection', crossSection, found)
            IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')

            CALL initIonization(boundary(i)%bTypes(s)%obj, species(s)%obj%m, m0, n0, v0, T0, &
                                                           speciesID, effTime, crossSection, eThreshold)

          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 moduleMeshInputGmsh2, ONLY: initGmsh2
      USE moduleMeshInput0D,    ONLY: init0D
      USE moduleMesh3DCart, ONLY: connectMesh3DCart
      USE moduleMesh2DCyl,  ONLY: connectMesh2DCyl
      USE moduleMesh2DCart, ONLY: connectMesh2DCart
      USE moduleMesh1DRad,  ONLY: connectMesh1DRad
      USE moduleMesh1DCart, ONLY: connectMesh1DCart
      USE moduleErrors
      USE moduleOutput
      USE moduleRefParam
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
      REAL(8):: volume

      !Firstly, indicates if a specific mesh for MC collisions is being use
      doubleMesh = ASSOCIATED(meshForMCC, meshColl)

      !Selects the type of geometry.
      CALL config%get('geometry.type', mesh%geometry, found)
      IF (doubleMesh) meshColl%geometry = mesh%geometry

      !Gets the type of mesh
      CALL config%get('geometry.meshType', meshFormat, found)
      SELECT CASE(meshFormat)
      CASE ("gmsh2")
        CALL initGmsh2(mesh)
        IF (doubleMesh) CALL initGmsh2(meshColl)

      CASE ("0D")
        CALL config%get('geometry.meshType', meshFormat, found)
        CALL init0D(mesh)

      CASE DEFAULT
        CALL criticalError("Mesh format " // meshFormat // " not recogniced", "readGeometry")

      END SELECT

      !Reads the mesh file
      CALL config%get('geometry.meshFile', meshFile, found)
      pathMeshParticle = path // meshFile
      CALL mesh%readMesh(pathMeshParticle)
      DEALLOCATE(meshFile)
      IF (doubleMesh) THEN
        !Reads the mesh file for collisions
        CALL config%get('interactions.meshCollisions', meshFile, found)
        pathMeshColl = path // meshFile
        CALL meshColl%readMesh(pathMeshColl)

      END IF

      !Gets the volume for a 0D mesh
      !TODO: Try to constrain this to the inout for 0D
      IF (meshFormat == "0D") THEN
        CALL config%get('geometry.volume', volume, found)
        IF (found) THEN
          mesh%vols(1)%obj%volume = mesh%vols(1)%obj%volume*volume / Vol_ref
          mesh%nodes(1)%obj%v     = mesh%vols(1)%obj%volume

        END IF

      END IF

      !Creates the connectivity between elements
      SELECT CASE(mesh%geometry)
      CASE("3DCart")
        mesh%connectMesh => connectMesh3DCart

      CASE("2DCyl")
        mesh%connectMesh => connectMesh2DCyl

      CASE("2DCart")
        mesh%connectMesh => connectMesh2DCart

      CASE("1DRad")
        mesh%connectMesh => connectMesh1DRad

      CASE("1DCart")
        mesh%connectMesh => connectMesh1DCart

      CASE("0D")
        mesh%connectMesh => NULL()

      END SELECT
      IF (ASSOCIATED(mesh%connectMesh)) CALL mesh%connectMesh

      IF (doubleMesh) THEN
        meshColl%connectMesh => mesh%connectMesh
        CALL meshColl%connectMesh

      END IF

      !Builds the K matrix for the Particles mesh
      CALL mesh%constructGlobalK()

      !Assign the procedure to find a volume for meshColl
      IF (doubleMesh) THEN
        findCellColl => findCellCollMesh

      ELSE
        findCellColl => findCellSameMesh

      END IF

    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 = inj%species%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