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

      !Reads reference parameters
      CALL readReference(config)

      !Reads output parameters
      CALL readOutput(config)

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

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

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

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

      !Reads case parameters
      CALL verboseError('Reading Case Parameters...')
      CALL readCase(config)

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

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

    END SUBROUTINE readConfig

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

      CALL config%get(object // '.radius', r_ref, found_r)
      
      !Derived parameters
      v_ref     = DSQRT(kb*T_ref/m_ref)  !reference velocity
      !TODO: Make this solver dependent
      IF (found_r) THEN
        sigma_ref = PI*(r_ref+r_ref)**2    !reference cross section
        L_ref = 1.D0/(sigma_ref*n_ref)     !mean free path

      ELSE
        L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !Debye length
        !TODO: Obtain right sigma_ref for PIC case
        sigma_ref = PI*(4.D-10)**2        !reference cross section

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

      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 non-dimensional
      tau = tau / ti_ref
      tauMin = tauMin / ti_ref

    END SUBROUTINE readCase

    !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(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

      !Creates output folder
      !TODO: Add option for custon name output_folder
      CALL DATE_AND_TIME(date_now, time_now)
      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 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 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

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

      !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 moduleCollisions
      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
      LOGICAL:: found
      INTEGER:: nInteractions, nCollisions
      INTEGER:: i, k, ij
      INTEGER:: pt_i, pt_j

      CALL initInteractionMatrix(interactionMatrix)

      CALL config%get('interactions.folderCollisions', pathCollisions, found)

      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 // '.crossSections', found, n_children = nCollisions)
        ij = interactionIndex(pt_i,pt_j)
        CALL interactionMatrix(ij)%init(nCollisions)

        DO k = 1, nCollisions
          WRITE (kString, '(I2)') k
          CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile, found)
          crossSecFilePath = pathCollisions // crossSecFile
          CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)

        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