! 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 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 particle boundaries
      CALL verboseError('Reading particle boundary conditions...')
      CALL readBoundaryParticle(config)
      CALL checkStatus(config, "readBoundaryParticle")

      ! read EM boundaries
      CALL verboseError('Reading EM boundary conditions...')
      CALL readBoundaryEM(config)
      CALL checkStatus(config, "readBoundaryEM")

      ! Read Physical Surfaces
      call verboseError('Reading Physical Surfaces...')
      call readPhysicalSurfaces(config)
      call checkStatus(config, 'readPhysicalSurfaces')

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

      !Read solver parameters
      CALL verboseError('Reading Solver parameters...')
      CALL readSolver(config)
      CALL checkStatus(config, "readSolver")

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

      !Read probes
      CALL verboseError('Reading Probes...')
      CALL readProbes(config)
      CALL checkStatus(config, 'readProbes')

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

      !Read injection of particles
      CALL verboseError('Reading injection of particles from boundaries...')
      CALL readInject(config)
      CALL checkStatus(config, "readInject")

      !Read average scheme
      CALL verboseError('Reading average scheme...')
      CALL readAverage(config)
      CALL checkStatus(config, "readAverage")

      !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 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')

      ! initialize values for cross sections
      sigmaVrel_ref = 0.0d0
      r_ref = 0.0d0

      CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found)
      CALL config%get(object // '.radius', r_ref, found)

      call calculateDerivedRef()

    END SUBROUTINE readReference

    !Reads the specific case parameters
    SUBROUTINE readSolver(config)
      USE moduleRefParam
      USE moduleErrors
      USE moduleCaseParam
      USE moduleSolver
      USE moduleSpecies
      USE moduleCollisions
      USE moduleOutput
      USE moduleMesh
      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, WSType, EMType
      REAL(8):: B(1:3)
      INTEGER:: nTau, nSolver
      INTEGER:: i, n
      CHARACTER(2):: iString
      CHARACTER(1):: tString

      object = 'solver'

      !Time parameters
      CALL config%info(object // '.tau', found, n_children = nTau)
      IF (.NOT. found .OR. nTau == 0) THEN
        CALL criticalError('Required parameter tau not found','readSolver')

      END IF

      ALLOCATE(tau(1:nSpecies))

      DO i = 1, nTau
        WRITE(iString, '(I2)') i
        CALL config%get(object // '.tau(' // TRIM(iString) // ')', tau(i), found)

      END DO

      tauMin = MINVAL(tau(1:nTau))

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

      END IF

      !Gets the simulation final time
      CALL config%get(object // '.finalTime', finalTime, found)
      IF (.NOT. found) THEN
        CALL criticalError('Required parameter finalTime not found','readSolver')

      END IF
      !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) THEN
        tInitial = INT(initialTime / tauMin)

      END IF

      !Gest the pusher for each species
      IF (mesh%dimen > 0) THEN
        CALL config%info(object // '.pusher', found, n_children = nSolver)
        IF (.NOT. found .OR. nSolver /= nSpecies) THEN
          CALL criticalError('Required parameter pusher not found','readSolver')

        END IF

      END IF

      !Allocates all the pushers for particles
      ALLOCATE(solver%pusher(1:nSpecies))
      !Initialize pushers
      DO i = 1, nSpecies
        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

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

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

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

      !Get EM solver
      CALL config%get(object // '.EMSolver', EMType, found)
      IF (found) THEN
        CALL solver%initEM(EMType)
        SELECT CASE(EMType)
        CASE("ConstantB")
          !Read constant magnetic field
          DO i = 1, 3
            WRITE(iString, '(i2)') i
            CALL config%get(object // '.B(' // iString // ')', B(i), found)
            IF (.NOT. found) THEN
              CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver')

            END IF

          END DO
          
          !Non-dimensional magnetic field
          B = B / B_ref

          !Assign it to the nodes
          DO n =1, mesh%numNodes
            mesh%nodes(n)%obj%emData%B(1) = B(1)
            mesh%nodes(n)%obj%emData%B(2) = B(2)
            mesh%nodes(n)%obj%emData%B(3) = B(3)

          END DO

        END SELECT

      END IF

    END SUBROUTINE readSolver

    !Read the initial information for the species
    SUBROUTINE readInitial(config)
      USE moduleSpecies
      USE moduleMesh
      USE moduleOutput
      USE moduleRefParam
      USE moduleRandom
      USE moduleProbe
      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
      CHARACTER(:), ALLOCATABLE:: filename
      REAL(8), ALLOCATABLE, DIMENSION(:):: density, temperature
      REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity
      INTEGER, ALLOCATABLE, DIMENSION(:):: nodes
      INTEGER:: nNodes
      REAL(8), ALLOCATABLE, DIMENSION(:)::   sourceScalar
      REAL(8), ALLOCATABLE, DIMENSION(:,:):: sourceArray
      !Density at the volume centroid
      REAL(8):: densityCen
      !Mean velocity and temperature at particle position
      REAL(8):: velocityXi(1:3), temperatureXi
      INTEGER:: nNewPart = 0
      REAL(8):: weight = 0.D0
      CLASS(meshCell), POINTER:: cell
      TYPE(particle), POINTER:: partNew
      REAL(8):: vTh
      TYPE(lNode), POINTER:: partCurr, partNext

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

      IF (found) THEN
        !Reads the information from initial species
        DO i = 1, nInitial
          WRITE(iString, '(I2)') i
          object = 'solver.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(filename, density, velocity, temperature)
          !Check if initial number of particles is given
          CALL config%get(object // '.particlesPerCell', nNewPart, found)

          !For each volume in the node, create corresponding particles
          DO e = 1, mesh%numCells
            !Scale variables
            !Density at centroid of cell
            nNodes = mesh%cells(e)%obj%nNodes
            nodes = mesh%cells(e)%obj%getNodes(nNodes)
            ALLOCATE(sourceScalar(1:nNodes))
            ALLOCATE(sourceArray(1:nNodes, 1:3))
            sourceScalar = density(nodes)
            densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar)

            !Calculate number of particles
            IF (.NOT. found) THEN
              nNewPart = FLOOR(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)

            END IF
            weight = densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / REAL(nNewPart)

            !Allocate new particles
            DO p = 1, nNewPart
              ALLOCATE(partNew)
              partNew%species => species(sp)%obj
              partNew%r       =  mesh%cells(e)%obj%randPos()
              partNew%Xi      =  mesh%cells(e)%obj%phy2log(partNew%r)
              !Get mean velocity at particle position
              sourceArray(:,1) = velocity((nodes), 1)
              sourceArray(:,2) = velocity((nodes), 2)
              sourceArray(:,3) = velocity((nodes), 3)

              velocityXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceArray)
              velocityXi = velocityXi / v_ref

              !Get temperature at particle position
              sourceScalar = temperature(nodes)
              temperatureXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceScalar)
              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%cell    = e
              IF (doubleMesh) THEN
                partNew%cellColl = findCellBrute(meshColl, partNew%r)

              ELSE
                partNew%cellColl = partNew%cell

              END IF

              partNew%n_in   = .TRUE.

              partNew%weight = weight

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

              !Assign particle to list in volume
              IF (listInCells) THEN
                cell => meshforMCC%cells(partNew%cellColl)%obj
                CALL OMP_SET_LOCK(cell%lock)
                CALL cell%listPart_in(sp)%add(partNew)
                cell%totalWeight(sp) = cell%totalWeight(sp) + partNew%weight
                CALL OMP_UNSET_LOCK(cell%lock)

              END IF

              !Assign particles to probes
              CALL doProbes(partNew)

            END DO

            DEALLOCATE(sourceScalar, sourceArray)

          END DO

        END DO

        !Convert temporal list of particles into initial partOld array
        !Deallocate 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 // '.boundariesParticle', boundaryParticleOutput, found)
      call config%get(object // '.boundariesEM', boundaryEMOutput, found)
      call config%get(object // '.injects', injectOutput, 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
      use moduleMesh, only: qSpecies
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(2):: iString
      CHARACTER(:), ALLOCATABLE:: object
      CHARACTER(:), ALLOCATABLE:: speciesType
      REAL(8):: mass, charge
      LOGICAL:: found
      INTEGER:: s
      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 s = 1, nSpecies
        WRITE(iString, '(I2)') s
        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(s)%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(s)%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(s)%obj%name,   found)
        CALL config%get(object // '.weight', species(s)%obj%weight, found)
        species(s)%obj%n = s
        species(s)%obj%m = mass

      END DO

      ! Allocate the vector with the species charges for calculating the EM field
      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

      !Read relations between species
      DO s = 1, nSpecies
        WRITE(iString, '(I2)') s
        object = 'species(' // TRIM(iString) // ')'
        SELECT TYPE(sp => species(s)%obj)
        TYPE IS (speciesNeutral)
          !Get 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)
          !Get 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(partWScheme%lock)

    END SUBROUTINE readSpecies

    !Reads information about interactions between species
    SUBROUTINE readInteractions(config)
      USE moduleSpecies
      USE moduleList
      USE moduleCollisions
      USE moduleCoulomb
      USE moduleErrors
      USE moduleMesh
      USE moduleCaseParam
      USE moduleRefParam
      use moduleOutput, only: collOutput
      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:: nPairs, nCollisions
      INTEGER:: i, k, ij
      INTEGER:: pt_i, pt_j
      REAL(8):: energyThreshold, energyBinding
      CHARACTER(:), ALLOCATABLE:: electron, electronSecondary
      INTEGER:: e
      CLASS(meshCell), POINTER:: cell

      !Firstly, check if the object 'interactions' exists
      CALL config%info('interactions', found)
      IF (found) THEN
        !Check if MC collisions have been defined
        CALL config%info('interactions.collisions', found)
        IF (found) THEN
          doMCCollisions = .TRUE.
          !Read collision time step
          CALL config%info('interactions.timeStep', found)
          IF (found) THEN
            CALL config%get('interactions.timeStep', tauColl, found)
            tauColl  = tauColl / ti_ref

          ELSE
            tauColl = tauMin

          END IF

          IF (tauColl < tauMin) THEN
            CALL criticalError('Collisional time step below minimum time step.', 'readInteractions')

          END IF

          everyColl = NINT(tauColl / tauMin)

          !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(partCollisions%lock)

          CALL config%info('interactions.collisions', found, n_children = nPairs)
          DO i = 1, nPairs
            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, pt_i, pt_j)

            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)

              CASE ('chargeExchange')
                !Resonant charge exchange
                CALL initBinaryChargeExchange(interactionMatrix(ij)%collisions(k)%obj, crossSecFilePath)

              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 config%get(object // '.electronSecondary', electronSecondary, found)
                IF (found) THEN
                  CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
                                            crossSecFilePath, energyThreshold, electron, electronSecondary)

                ELSE
                  CALL initBinaryIonization(interactionMatrix(ij)%collisions(k)%obj, &
                                            crossSecFilePath, energyThreshold, electron)

                END IF

              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, electron)

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

              END SELECT

            END DO

          END DO

          !Init the required arrays in each volume to account for MCC.
          DO e = 1, meshForMCC%numCells
            cell => meshForMCC%cells(e)%obj

            !Allocate Maximum cross section per collision pair and assign the initial collision rate
            ALLOCATE(cell%sigmaVrelMax(1:nCollPairs))
            cell%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)

            IF (collOutput) THEN
              ALLOCATE(cell%tallyColl(1:nCollPairs))
              DO k = 1, nCollPairs
                ALLOCATE(cell%tallyColl(k)%tally(1:interactionmatrix(k)%amount))
                cell%tallyColl(k)%tally = 0

              END DO

            END IF

          END DO

        END IF

        !Check if Coulomb scattering is defined
        CALL config%info('interactions.Coulomb', found)
        IF (found) THEN

          CALL config%info('interactions.Coulomb', found, n_children = nPairs)
          IF (nPairs > 0) THEN
            nCoulombPairs = nPairs
            doCoulombScattering = .TRUE.
            ALLOCATE(coulombMatrix(1:nPairs))
            DO i = 1, nPairs
              WRITE(iString, '(I2)') i
              object = 'interactions.Coulomb(' // 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 coulombMatrix(i)%init(pt_i, pt_j)

            END DO

          END IF

        END IF

      END IF

      listInCells = doMCCollisions .OR. doCoulombScattering

    END SUBROUTINE readInteractions

    !Reads boundary conditions for the mesh
    SUBROUTINE readBoundaryParticle(config)
      use moduleMesh
      USE moduleErrors
      USE moduleList, ONLY: partSurfaces
      use moduleRefParam, only: m_ref
      use moduleConstParam, only: me
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      integer:: b
      character(2):: bString
      character(:), allocatable:: object
      LOGICAL:: found
      character(:), allocatable:: 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, electronSecondaryID
      character(:), allocatable:: speciesName, crossSection, electronSecondary
      integer:: s_incident

      ! Read models of particles
      object = 'boundaries.particles'
      CALL config%info(object, found, n_children = nBoundariesParticle)
      allocate(boundariesParticle(1:nBoundariesParticle))
      DO b = 1, nBoundariesParticle
        WRITE(bString, '(i2)') b
        object = 'boundaries.particles(' // trim(bString) // ')'

        CALL config%get(object // '.type', bType,     found)
        if (.not. found) then
          call criticalError('Required parameter "type" for particle boundary condition not found', &
                             'initBoundaryParticle')

        end if

        associate(bound => boundariesParticle(b)%obj)
          SELECT CASE(bType)
            CASE('reflection')
              ALLOCATE(boundaryReflection:: bound)

            CASE('absorption')
              ALLOCATE(boundaryAbsorption:: bound)

            CASE('transparent')
              ALLOCATE(boundaryTransparent:: bound)

            CASE('axis')
              ALLOCATE(boundaryAxis:: bound)

            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 wallTemperature_init(bound, Tw, cw)

            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) CALL criticalError("missing parameter 'mass' for neutrals in ionization", 'readBoundary')
              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 // '.electronSecondary', electronSecondary, found)
              IF (found) THEN
                electronSecondaryID = speciesName2Index(electronSecondary)
                CALL ionization_init(bound, me/m_ref, m0, n0, v0, T0, &
                                           speciesID, effTime, crossSection, eThreshold, electronSecondaryID)
              ELSE
                CALL ionization_init(bound, me/m_ref, m0, n0, v0, T0, &
                                           speciesID, effTime, crossSection, eThreshold)

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

            case('quasiNeutrality')
              call config%get(object // '.incident', speciesName, found)
              if (.not. found) call criticalError("Incident species name not found for quasiNeutrality boundary model", 'readBoundary')

              s_incident = speciesName2Index(speciesName)

              call quasiNeutrality_init(bound, s_incident)

            CASE('outflowAdaptive')
              call config%get(object // '.incident', speciesName, found)
              if (.not. found) call criticalError("Incident species name not found for quasiNeutrality boundary model", 'readBoundary')

              s_incident = speciesName2Index(speciesName)

              call outflowAdaptive_init(bound, s_incident)

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

          END SELECT

          bound%n = b

          CALL config%get(object // '.name', bound%name, found)
          if (.not. found) then
            call criticalError('Required parameter "name" for particle boundary condition not found', &
                               'initBoundaryParticle')

          end if

        end associate

      END DO

      !Init the list of particles from surfaces
      CALL OMP_INIT_LOCK(partSurfaces%lock)

    END SUBROUTINE readBoundaryParticle

    SUBROUTINE readBoundaryEM(config)
      USE moduleMesh
      USE moduleOutput
      USE moduleErrors
      USE moduleEM
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found
      INTEGER:: b
      CHARACTER(2):: bString
      character(:), allocatable:: bType

     CALL config%info('boundaries.EM', found, n_children = nBoundariesEM)
      
      IF (found) THEN
        ALLOCATE(boundariesEM(1:nBoundariesEM))

      END IF
      
      do b = 1, nBoundariesEM
        write(bString, '(I2)') b
        object = 'boundaries.EM(' // TRIM(bString) // ')'

        associate(self => boundariesEM(b)%obj)
          call config%get(object // '.type', bType,     found)
          if (.not. found) then
            call criticalError('Required parameter "type" for EM boundary condition not found', &
                               'initBoundaryEM')

          end if

          select case(bType)
          case ("dirichlet")
            allocate(boundaryEMDirichlet:: self)

            CALL initDirichlet(self, config, object)

          case ("dirichlettime")
            allocate(boundaryEMDirichletTime:: self)

            call initDirichletTime(self, config, object)

          case ("neumann")
            allocate(boundaryEMNeumann:: self)

            call initNeumann(self, config, object)

          case ("floating")
            allocate(boundaryEMFloating:: self)

            call initFloating(self, config, object)

          case default
            call criticalError('Boundary type ' // bType // ' not supported', 'readBoundaryEM')

          end select

          self%n = b
          allocate(self%nodes(0))

          call config%get(object // '.name', self%name, found)
          if (.not. found) then
            call criticalError('Required parameter "name" for EM boundary condition not found', &
                               'initBoundaryEM')

          end if

        end associate

      end do

    END SUBROUTINE readBoundaryEM

    subroutine readPhysicalSurfaces(config)
      use json_module
      use moduleMesh
      use moduleErrors
      implicit none
      
      type(json_file), intent(inout):: config
      character(:), allocatable:: object
      logical:: found
      integer:: ps
      character(2):: psString, sSTring
      integer:: nParticleModels
      character(:), allocatable:: particleModel
      character(:), allocatable:: EMModel
      integer:: s, boundaryIndex

      call config%info('physicalSurfaces', found, n_children = nPhysicalSurfaces)
      if (found) then
        allocate(physicalSurfaces(1:nPhysicalSurfaces))

      end if

      do ps = 1, nPhysicalSurfaces
        write(psString, '(I2)') ps
        object = 'physicalSurfaces(' // trim(psString) // ')'

        allocate(physicalSurfaces(ps)%nodes(0))
        allocate(physicalSurfaces(ps)%edges(0))
        call config%get(object // '.index', physicalSurfaces(ps)%index, found) 
        if (.not. found) then
          call criticalError('Physical surface index not found', 'readPhysicalSurfaces')

        end if
        
        ! Link models for particles
        call config%info(object // '.particles', found, n_children = nParticleModels)
        if ((.not. found) .or. &
            (nParticleModels /= nSpecies)) then
          call criticalError('Not enough models for particles provided', 'readPhysicalSurfaces')

        end if
        allocate(physicalSurfaces(ps)%particles(1:nSpecies))
        do s = 1, nSpecies
          write(sString, '(I2)') s
          call config%get(object // '.particles(' // trim(sSTring) // ')', particleModel, found)
          boundaryIndex = boundaryParticleName_to_Index(particleModel)
          physicalSurfaces(ps)%particles(s)%obj => boundariesParticle(boundaryIndex)%obj

        end do

        ! Link electromagnetic boundary condition
        call config%get(object // '.EM', EMModel, found)
        if (found) then
          boundaryIndex = boundaryEMName_to_Index(EMModel)
          physicalSurfaces(ps)%EM => boundariesEM(boundaryIndex)%obj

        end if

      end do

    end subroutine readPhysicalSurfaces

    !Read the geometry (mesh) for the case
    SUBROUTINE readGeometry(config)
      USE moduleMesh
      USE moduleMeshInputGmsh2, ONLY: initGmsh2
      USE moduleMeshInputVTU,   ONLY: initVTU
      USE moduleMeshInput0D,    ONLY: init0D
      USE moduleMeshInputText,  ONLY: initText
      USE moduleMesh3DCart
      USE moduleMesh2DCyl
      USE moduleMesh2DCart
      USE moduleMesh1DRad
      USE moduleMesh1DCart
      USE moduleErrors
      USE moduleOutput
      USE moduleRefParam
      USE moduleSolver
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found
      CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
      REAL(8):: volume
      integer:: b, ps, s
      integer:: e
      integer:: nVolColl

      object = 'geometry'

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

      ELSE
        CALL config%info('interactions', found)
        IF (found) THEN
          !Points the meshForMCC pointer to the Particles Mesh
          meshForMCC => mesh

        END IF
        doubleMesh = .FALSE.

      END IF

      !Get the dimension of the geometry
      CALL config%get(object // '.dimension', mesh%dimen, found)
      IF (.NOT. found) THEN
        CALL criticalError('Required parameter geometry.dimension not found', 'readGeometry')

      END IF
      IF (doubleMesh) THEN
        meshColl%dimen = mesh%dimen

      END IF
        

      SELECT CASE(mesh%dimen)
      CASE (0)
        CALL init0D(mesh)

        !Read the 0D mesh
        CALL mesh%readMesh(pathMeshParticle)

        !Get the volume
        CALL config%get(object // '.volume', volume, found)
        !Rescale the volume
        IF (found) THEN
          mesh%cells(1)%obj%volume = mesh%cells(1)%obj%volume*volume / Vol_ref
          mesh%nodes(1)%obj%v     = mesh%cells(1)%obj%volume

        END IF

      CASE (1:3)
        !Select the type of geometry
        CALL config%get(object // '.type', mesh%geometry, found)
        IF (doubleMesh) THEN
          meshColl%geometry = mesh%geometry

        END IF
        
        !Check the consistency between geometry and dimension
        SELECT CASE(mesh%dimen)
        CASE(3)
          IF (mesh%geometry /= 'Cart') THEN
            CALL criticalError('No available geometry ' // mesh%geometry // ' in 3D', 'readGeometry') 

          END IF

        CASE(2)
          IF (mesh%geometry /= 'Cart' .AND. &
              mesh%geometry /= 'Cyl') THEN
            CALL criticalError('No available geometry ' // mesh%geometry // ' in 2D', 'readGeometry') 

          END IF

        CASE(1)
          IF (mesh%geometry /= 'Cart' .AND. &
              mesh%geometry /= 'Rad') THEN
            CALL criticalError('No available geometry ' // mesh%geometry // ' in 1D', 'readGeometry') 

          END IF

        END SELECT

        !Link the procedure to connect meshes
        SELECT CASE(mesh%dimen)
        CASE(3)
          mesh%connectMesh => connectMesh3DCart

        CASE(2)
          SELECT CASE(mesh%geometry)
          CASE("Cyl")
            mesh%connectMesh => connectMesh2DCyl

          CASE("Cart")
            mesh%connectMesh => connectMesh2DCart

          END SELECT

        CASE(1)
          SELECT CASE(mesh%geometry)
          CASE("Rad")
            mesh%connectMesh => connectMesh1DRad

          CASE("Cart")
            mesh%connectMesh => connectMesh1DCart

          END SELECT

        END SELECT

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

        END IF

        !Get the format of mesh
        CALL config%get(object // '.meshType', meshFormat, found)
        SELECT CASE(meshFormat)
        CASE ("gmsh2")
          CALL initGmsh2(mesh)
          IF (doubleMesh) THEN
            CALL initGmsh2(meshColl)

          END IF

        CASE ("vtu")
          CALL initVTU(mesh)
          IF (doubleMesh) THEN
            CALL initVTU(meshColl)

          END IF

        case ("text")
          !Check if the geometry is right.
          if (mesh%dimen /= 1) then
            call criticalError("Text mesh is only allowed for 1D geometries", 'readGeometry')
          
          end if

          !Read the mesh
          call initText(mesh)
          if (doubleMesh) then
            call initText(meshColl)

          end if

        CASE DEFAULT
          CALL criticalError('Mesh format ' // meshFormat // ' not defined.', 'readGeometry')

        END SELECT

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

        END IF

      CASE DEFAULT
        CALL criticalError("Mesh dimension not recogniced", "readGeometry")

      END SELECT

      !Assign the procedure to find a cell for meshColl
      IF (doubleMesh) THEN
        findCellColl => findCellCollMesh
        ! Link edges with cells in meshColl
        DO e=1, mesh%numEdges
          nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%center())
          IF (nVolColl > 0) THEN
            mesh%edges(e)%obj%eColl => meshColl%cells(nVolColl)%obj

          ELSE
            CALL criticalError("No connection between edge and meshColl", "readGeometry")

          END IF

        end do

      ELSE
        findCellColl => findCellSameMesh
        ! Link edges with cells in meshColl
        DO e=1, mesh%numEdges
          IF (ASSOCIATED(mesh%edges(e)%obj%e1)) THEN
            mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e1

          ELSE
            mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e2

          END IF

        end do

      END IF

      ! If needed, add nodes and edges to boundary models
      ! Particle boundaries
      do b = 1, nBoundariesParticle
        select type(bound => boundariesParticle(b)%obj)
        type is(boundaryQuasiNeutrality)
          ! Loop over all physical surfaces
          do ps = 1, nPhysicalSurfaces
            ! Loop over all species
            do s = 1, nSpecies
              ! If the boundary for the species is linked to the one analysing, add the edges
              if (associated(physicalSurfaces(ps)%particles(s)%obj, bound)) then
                bound%edges = [bound%edges, physicalSurfaces(ps)%edges]

              end if

            end do

          end do

          allocate(bound%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
          bound%alpha = 0.d0

        type is(boundaryOutflowAdaptive)
          ! Loop over all physical surfaces
          do ps = 1, nPhysicalSurfaces
            ! Loop over all species
            do s = 1, nSpecies
              ! If the boundary for the species is linked to the one analysing, add the edges
              if (associated(physicalSurfaces(ps)%particles(s)%obj, bound)) then
                bound%edges = [bound%edges, physicalSurfaces(ps)%edges]

              end if

            end do

          end do

          allocate(bound%velocity_shift(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
          bound%velocity_shift = 0.d0

        end select

      end do

      ! EM Boundaries
      do b = 1, nBoundariesEM
        associate(bound => boundariesEM(b)%obj)
          ! Loop over all physical surfaces
          do ps = 1, nPhysicalSurfaces
            ! If the boundary for the species is linked to the one analysing, add the edges
            if (associated(physicalSurfaces(ps)%EM, bound)) then
              bound%nodes = [bound%nodes, physicalSurfaces(ps)%nodes]
            end if

          end do

          bound%nNodes = size(bound%nodes)

        end associate

      end do

      if (mesh%dimen > 0) then
        ! Builds the K matrix for the Particles mesh
        call mesh%constructGlobalK()

      end if

    END SUBROUTINE readGeometry

    SUBROUTINE readProbes(config)
      USE moduleProbe
      USE moduleErrors
      USE json_module
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      CHARACTER(:), ALLOCATABLE:: object
      LOGICAL:: found
      CHARACTER(2):: iString
      INTEGER:: i
      CHARACTER(:), ALLOCATABLE:: speciesName
      REAL(8), ALLOCATABLE, DIMENSION(:):: r
      REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3
      INTEGER, ALLOCATABLE, DIMENSION(:):: points
      REAL(8):: everyTimeStep

      CALL config%info('output.probes', found, n_children = nProbes)

      IF (found) ALLOCATE(probe(1:nProbes))

      DO i = 1, nProbes
        WRITE(iString, '(I2)') i
        object = 'output.probes(' // trim(iString) // ')'

        CALL config%get(object // '.species', speciesName, found)
        CALL config%get(object // '.position', r, found)
        CALL config%get(object // '.velocity_1', v1, found)
        CALL config%get(object // '.velocity_2', v2, found)
        CALL config%get(object // '.velocity_3', v3, found)
        CALL config%get(object // '.points', points, found)
        CALL config%get(object // '.timeStep', everyTimeStep, found)

        IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes')

        CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep)

      END DO

    END SUBROUTINE readProbes

    !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

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

        CALL initInject(injects(i)%obj, i, object, config)

        CALL readVelDistr(config, injects(i)%obj, object)
        
      END DO

    END SUBROUTINE readInject

    SUBROUTINE readAverage(config)
      USE moduleAverage
      USE moduleCaseParam, ONLY: tauMin
      USE moduleMesh, ONLY: mesh
      USE moduleSpecies, ONLY: nSpecies
      USE moduleRefParam, ONLY: ti_ref
      IMPLICIT NONE

      TYPE(json_file), INTENT(inout):: config
      LOGICAL:: found
      REAL(8):: tStart
      INTEGER:: n

      CALL config%info('average', found)
      IF (found) THEN
        useAverage = .TRUE.
        CALL config%get('average.startTime', tStart, found)

        IF (found) THEN
          tAverageStart = INT(tStart / ti_ref / tauMin)

        ELSE
          tAverageStart = 0
          
        END IF

        ALLOCATE(averageScheme(1:mesh%numNodes))

        DO n = 1, mesh%numNodes
          ALLOCATE(averageScheme(n)%mean%output(1:nSpecies))
          ALLOCATE(averageScheme(n)%deviation%output(1:nSpecies))

        END DO

      END IF

    END SUBROUTINE readAverage

    !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
      class(injectGeneric), INTENT(inout):: inj
      CHARACTER(:), ALLOCATABLE, INTENT(in):: object
      INTEGER:: i
      CHARACTER(2):: iString
      CHARACTER(:), ALLOCATABLE:: fvType
      LOGICAL:: found
      REAL(8):: v, temperature, 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")
          temperature = inj%temperature(i)
          CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m)

        CASE ("Half-Maxwellian")
          temperature = inj%temperature(i)
          CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m)

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

        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

    SUBROUTINE initOutput(inputFile)
      USE moduleMesh,   ONLY: mesh, doubleMesh, pathMeshParticle, pathMeshColl
      USE moduleOutput, ONLY: createOutputFolder, writeReference, copyFileToOutput, writeCommit
      IMPLICIT NONE

      CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
      !If everything is correct, creates the output folder
      call createOutputFolder()
      !Copies input file to output folder
      call copyFileToOutput(inputFile)
      !Copies particle mesh
      IF (mesh%dimen > 0) THEN
        call copyFileToOutput(pathMeshParticle)
        IF (doubleMesh) THEN
          call copyFileToOutput(pathMeshColl)

        END IF

      END IF

      ! Write commit of fpakc
      call writeCommit()

      call writeReference()

    END SUBROUTINE initOutput

END MODULE moduleInput
