!> \file
!> Main model driver for CLASSIC in stand-alone mode using specified boundary
!! conditions and atmospheric forcing.
!! @author D. Verseghy, V. Arora, J. Melton, E. Chan
!! This driver program initializes the run, reads in CLASSIC input files,
!! calls the time-stepping subroutine, calls subroutines
!! that aggregate and write outputs, and closes the run for this grid cell.

module main

  implicit none

  public :: main_driver

contains

  !> \ingroup main_main_driver
  !> @{
  !>
  subroutine main_driver (longitude, latitude, lonIndex, latIndex, lonLocalIndex, latLocalIndex)

    use mainCore,           only : main_core_driver

    use classicParams,      only : nlat, nmos, nmon, monthend, DELT, zbldJobOpt, zrfhJobOpt, zrfmJobOpt
    use landuseChange,      only : initializeLandCover
    use ctemStateVars,      only : vrot, vgat, c_switch, initRowVarsBioGeoChem, resetMonthEnd, &
                                   resetYearEnd, resetMosaicAccum
    use classStateVars,     only : class_rot, resetAccVars, resetClassMon, resetClassYr,initRowVarsPhysics 
    use prepareOutputs,     only : class_monthly_aw, ctem_annual_aw, ctem_monthly_aw, &
                                   ctem_daily_aw, class_annual_aw, class_hh_w, class_daily_aw, convertUnitsCTEM
    use modelStateDrivers,  only : read_initialstate, write_restart, getInput, deallocInput, getMet
    use generalUtils,       only : run_model, findPermafrostVars, initRandomSeed, checksumCalc
    use metDisaggModule,    only : disaggMet
    use outputManager,      only : consecDays

    implicit none

    ! Arguments
    real, intent(in) :: longitude, latitude                 !< Longitude/latitude of grid cell (degrees)
    integer, intent(in) :: lonIndex, latIndex               !< Index of grid cell being run on the input files grid
    integer, intent(in) :: lonLocalIndex, latLocalIndex     !< Index of grid cell being run on the output files grid

    ! Local variables
    integer :: lastDOY             !< Initialized to 365 days, can be overwritten later is leap = true and it is a leap year.
    integer :: metTimeIndex        !< Counter used to move through the meteorological input arrays
    logical :: metDone             !< Logical switch when the end of the stored meteorological array is reached.
    integer :: runyr               !< Year of the model run (counts up starting with readMetStartYear continously, even if metLoop > 1)

    integer :: NCOUNT  !< Counter for daily averaging
    integer :: NDAY    !< Number of short (physics) timesteps in one day. e.g., if physics timestep is 15 min this is 48.
    integer :: IMONTH  !< Month of the year simulation is in.
    integer :: DOM     !< Day of month counter
    integer :: NT      !<
    integer :: IDAY    !< Julian day of the year
    integer :: IYEAR   !< Year of run
    integer :: N       !< Timestep counter

    integer   :: lopcount ! xday, month1, month2,

    !    =================================================================================

    associate( &

    ! This driver is set up to handle one grid cell with any number of mosaic tiles.

    nltest            => nlat,                          & !< integer: Number of grid cells being modelled for this run
    nmtest            => nmos,                          & !< integer: Number of mosaic tiles per grid cell being modelled for this run

    ! Associate names with variables defined in derived types.

    readMetStartYear  => c_switch%readMetStartYear,     & !< integer: First year of meteorological forcing to read in from the met file 

    useTracer         => c_switch%useTracer,            & !< integer: Switch for use of a model tracer. If useTracer is 0 then the tracer code is not used.
                                                          !!          useTracer = 1 turns on a simple tracer that tracks pools and fluxes. The simple 
                                                          !!                        tracer then requires that the tracer values in the init_file and the
                                                          !!                        tracerCO2file are set to meaningful values for the experiment being run.
                                                          !!          useTracer = 2 means the tracer is 14C and will then call a 14C decay scheme.
                                                          !!          useTracer = 3 means the tracer is 13C and will then call a 13C fractionation scheme. 
    metLoop           => c_switch%metLoop,              & !< integer: no. of times the .met file is to be read. this
                                                          !!          option is useful to see how ctem's c pools
                                                          !!          equilibrate when driven with same climate data
                                                          !!          over and over again. 
    jhhstd            => c_switch%jhhstd,               & !< integer: day of the year to start writing the half-hourly output 
    jhhendd           => c_switch%jhhendd,              & !< integer: day of the year to stop writing the half-hourly output 
    jdstd             => c_switch%jdstd,                & !< integer: day of the year to start writing the daily output 
    jdendd            => c_switch%jdendd,               & !< integer: day of the year to stop writing the daily output 
    jhhsty            => c_switch%jhhsty,               & !< integer: simulation year (runyr) to start writing the half-hourly output 
    jhhendy           => c_switch%jhhendy,              & !< integer: simulation year (runyr) to stop writing the half-hourly output 
    jdsty             => c_switch%jdsty,                & !< integer: simulation year (runyr) to start writing the daily output 
    jdendy            => c_switch%jdendy,               & !< integer: simulation year (runyr) to stop writing the daily output 
    jmosty            => c_switch%jmosty,               & !< integer: Year to start writing out the monthly output files. If you want to write monthly outputs right 
    fixedYearLUC      => c_switch%fixedYearLUC,         & !< integer: Set the year to use for land cover if lnduseon is false. If set to -9999,
                                                          !!          we use the PFT distribution found in the initialization file. Any other year
                                                          !!          we search for that year in the LUCFile 
    fixedYearFER      => c_switch%fixedYearFER,         & !< integer: Set the year to use for fertilizer if fertilizeron is false.
    fixedYearDEP      => c_switch%fixedYearDEP,         & !< integer, Set the year to use for deposition if depositionon is false.
    fixedYearOBSWETF  => c_switch%fixedYearOBSWETF,     & !< integer: set the year to use for observed wetland fraction if transientOBSWETF is false. 

    ! Model switches:
    ctem_on           => c_switch%ctem_on,              & !< logical:
    Ncycle_on         => c_switch%Ncycle_on,            & !< logical:
    dofire            => c_switch%dofire,               & !< logical: 
    lnduseon          => c_switch%lnduseon,             & !< logical: 
    fertilizeron      => c_switch%fertilizeron,         & !< logical:
    depositionon      => c_switch%depositionon,         & !< logical: 
    doMethane         => c_switch%doMethane,            & !< logical: 
    transientOBSWETF  => c_switch%transientOBSWETF,     & !< logical: 
    doAnnualOutput    => c_switch%doAnnualOutput,       & !< logical: 
    doMonthOutput     => c_switch%doMonthOutput,        & !< logical: 
    doDayOutput       => c_switch%doDayOutput,          & !< logical: 
    doHhOutput        => c_switch%doHhOutput,           & !< logical: 
    projectedGrid     => c_switch%projectedGrid,        & !< logical: True if you have a projected lon lat grid, false if not. 
                                                          !!          Projected grids can only have regions referenced by the indexes,  
                                                          !!          not coordinates, when running a sub-region 
    ! ROW:
    ipeatlandrow     => vrot%ipeatland,                 & !< integer, dimension(:,:) : This is first set in read_from_ctm. 
    latIndexROW      => class_rot%latIndexROW,          & !< real, dimension(:) : Index of grid cell being run on the input files grid (latitude)
    lonIndexROW      => class_rot%lonIndexROW,          & !< real, dimension(:) : Index of grid cell being run on the input files grid (longitude)
    DLONROW          => class_rot%DLONROW,              & !< real, dimension(:) : 
    DLATROW          => class_rot%DLATROW,              & !< real, dimension(:) : 
    ZBLDROW          => class_rot%ZBLDROW,              & !< real, dimension(:) : 
    ZDHROW           => class_rot%ZDHROW,               & !< real, dimension(:) : 
    ZDMROW           => class_rot%ZDMROW,               & !< real, dimension(:) : 
    ZRFHROW          => class_rot%ZRFHROW,              & !< real, dimension(:) : 
    ZRFMROW          => class_rot%ZRFMROW,              & !< real, dimension(:) : 

    ! GAT:
    grclarea          => vgat%grclarea,                 & !< real, dimension(:) : 
    wetfrac_presgat   => vgat%wetfrac_pres              & !< real, dimension(:) : 
    )

    !    =================================================================================

    ! Variables set in this section are stored via derived types defined in modules, Although they are not used
    ! explicitly in this routine, they need to initialized before being used in other routines.

    DLATROW = latitude
    DLONROW = longitude

    !! Store the index of the cell so that it can be used for a more informative error message
    !! during projected grid simulations. It is much easier to find the problematic cells this way.
    latIndexROW = latIndex
    lonIndexROW = lonIndex

    !> The grid-average height for the momentum diagnostic variables, ZDMROW, and for the
    !! energy diagnostic variables, ZDHROW, are hard-coded to the standard anemometer
    !! height of 10 m and to the screen height of 2 m respectively.
    ZDMROW(:) = 10.0
    ZDHROW(:) = 2.0

    !> ZRFMROW and ZRFHROW, the reference heights at which the momentum variables (wind speed) and energy variables
    !! (temperature and specific humidity) are provided.  In a run using atmospheric model forcing data, these heights
    !! would vary by time step, but since this version of the driver is set up to use field data, ZRFMROW and ZRFHROW
    !! refer to the measurement height of these variables, which is fixed. The value is read in from the job options file.
    ZRFMROW = zrfmJobOpt
    ZRFHROW = zrfhJobOpt

    !> ZBLDROW, the atmospheric blending height.  Technically this variable depends on the length scale of the
    !! patches of roughness elements on the land surface, but this is difficult to ascertain.  Usually it is assigned a value of 50 m.
    !!  The value is read in from the job options file.
    ZBLDROW = zbldJobOpt

    !> The timestep counter N for the run is initialized to 0, the daily
    !! averaging counter NCOUNT is set to 1, and the total number of
    !! timesteps in the day NDAY is calculated as the number of seconds
    !! in a day (86400) divided by the timestep length DELT.
    N = 0
    NCOUNT = 1
    NDAY = 86400/NINT(DELT)
    metTimeIndex = 1    !< Counter used to move through the meteorological input arrays
    metDone = .false.   !< Logical switch when the end of the stored meteorological array is reached.
    run_model = .true.  !< Simple logical switch to either keep run going or finish
    IMONTH = 0          !< Month of the year simulation is in.
    DOM = 1             !< Day of month counter
    lopcount = 1
    lastDOY = 365
    wetfrac_presgat = - 9999. !< If transientOBSWETF or fixedYearOBSWETF != -9999 this variable will be overwritten with
                              !! real wetland fractions. Otherwise the negative is used as a switch so the dynamic
                              !! wetland extent is used instead of the prescribed.

    !> Initialize variables in preparation for the run
    call initRowVarsBioGeoChem
    call initRowVarsPhysics
    call resetAccVars(nlat, nmos)
    call resetMosaicAccum

    !> Read in the model initial state. 
    call read_initialstate(lonIndex, latIndex)
    
    if (ctem_on) then
      !> Read in the inputs for a run with biogeochemical component turned on
      call getInput('CO2') ! CO2 atmospheric concentration
      if (doMethane) call getInput('CH4') ! CH4 atmospheric concentration
      if (useTracer > 0) call getInput('tracerCO2',longitude,latitude) ! tracer atmospheric values
      if (.not. projectedGrid) then
        ! regular lon/lat grid
        if (dofire) call getInput('POPD',longitude,latitude) ! Population density
        if (dofire) call getInput('LGHT',longitude,latitude) ! Cloud-to-ground lightning frequency
        if (doMethane .and. transientOBSWETF .or. fixedYearOBSWETF /= - 9999) call getInput('OBSWETF',longitude,latitude) ! Observed wetland distribution
        if (lnduseon .or. (fixedYearLUC /= - 9999)) call getInput('LUC',longitude,latitude) ! Land use change
      else
        ! Projected grids use the lon and lat indexes, not the actual coordinates
        if (dofire) call getInput('POPD',longitude,latitude,projLonInd = lonIndex,projLatInd = latIndex) ! Population density
        if (dofire) call getInput('LGHT',longitude,latitude,projLonInd = lonIndex,projLatInd = latIndex) ! Cloud-to-ground lightning frequency
        if (doMethane .and. transientOBSWETF .or. fixedYearOBSWETF /= - 9999) &
            call getInput('OBSWETF',longitude,latitude,projLonInd = lonIndex,projLatInd = latIndex) ! Observed wetland distribution
        if (lnduseon .or. (fixedYearLUC /= - 9999)) &
            call getInput('LUC',longitude,latitude,projLonInd = lonIndex,projLatInd = latIndex) ! Land use change
      end if
      !> Regardless of whether lnduseon or not, we need to check the land cover that was read in
      !! and assign the CLASS PFTs as they are not read in when ctem_on.
      call initializeLandCover

      if (Ncycle_on) then
        if (fertilizeron .or. (fixedYearFER .ne. -9999)) call getInput('FER',longitude,latitude) ! Fertilizer
        if (depositionon .or. (fixedYearDEP .ne. -9999)) call getInput('DEP',longitude,latitude) ! Deposition
      end if

    end if

    !> Read in the meteorological forcing data to a suite of arrays
    if (.not. projectedGrid) then
      ! regular lon lat grid
      call getMet(longitude, latitude, nday)
    else
      ! Projected grids use the lon and lat indexes, not the actual coordinates
      call getMet(longitude, latitude, nday, projLonInd = lonIndex, projLatInd = latIndex)
    end if

    !> In preparation for the use of the random number generator by disaggMet,
    !! we need to provide a seed to allow repeatable results.
    call initRandomSeed()

    !> Now disaggregate the meteorological forcing to the right timestep
    !! for this model run (if needed; this is checked for in the subroutine)
    !! NOTE: At present this can't handle leap years as you need to then
    !! change diurnalDistribution to know which years are which 
    call disaggMet(longitude, latitude)

    ! Initialize accumulated array for monthly & yearly outputs
    call resetClassMon(nltest)
    call resetClassYr(nltest)
    if (ctem_on) then
      call resetMonthEnd(nltest, nmtest)
      call resetYearEnd(nltest, nmtest)
    end if

    !     **** LAUNCH RUN. **** ! EC - moved here from main_core_driver 

    !> The do while loop marks the beginning of the time stepping loop
    !! for the actual run.  N is incremented by 1, and the atmospheric forcing
    !! data for the current time step are updated for each grid cell or modelled
    !! area (see the manual section on “Data Requirements”).

    runyr = readMetStartYear  ! Initialize the runyr as the first year of met forcing.

    ! start up the main model loop

    mainModelLoop: do while (run_model)

      call main_core_driver(runyr, IMONTH, DOM, NDAY, metTimeIndex, ncount, & ! in 
                            IYEAR, IDAY, lastDOY, metDone, N) ! out/inout

      if (ncount == nday) then

        DOM = DOM + 1 ! increment the day of month counter

        ! Reset mosaic accumulator arrays.
        if (ctem_on) call resetMosaicAccum

        ! Determine the active layer depth and depth to the frozen water table.
        ! This only occurs once per day since they don't change rapidly.
        call findPermafrostVars(nmtest, nltest, iday)

      end if ! ncount eq nday

      !=======================================================================
              
      ! Half-hourly physics outputs
      if  (doHhOutput .and. &
          (runyr >= jhhsty) .and. &
          (runyr <= jhhendy) .and. &
          (iday >= jhhstd) .and. &
          (iday <= jhhendd) ) call class_hh_w(lonLocalIndex, latLocalIndex, nltest, &
                                              nmtest, ncount, nday, iday, runyr)

      ! Daily physics outputs
      if (doDayOutput .and. &
          (runyr >= jdsty) .and. &
          (runyr <= jdendy) .and. &
          (iday  >= jdstd) .and. &
          (iday  <= jdendd))  call class_daily_aw(lonLocalIndex, latLocalIndex, iday, nltest, nmtest, &
                                                  ncount, nday, lastDOY, runyr)

      do NT = 1,NMON
        if ((IDAY == monthend(NT + 1)) .and. (NCOUNT == NDAY)) then
          IMONTH = NT
          DOM = 1 ! reset the day of month counter
        end if
      end do

      ! Monthly physics outputs
      if (doMonthOutput .and. (runyr >= jmosty)) call class_monthly_aw(lonLocalIndex, &
                                                        latLocalIndex, IDAY, runyr, NCOUNT, &
                                                        NDAY, nltest, nmtest, lastDOY)

      ! Annual physics outputs
      if (doAnnualOutput) call class_annual_aw(lonLocalIndex, latLocalIndex, IDAY, runyr, NCOUNT, NDAY, &
                                              nltest, nmtest, lastDOY)

      if (ctem_on .and. (ncount == nday)) then

        ! Convert units in preparation for output:
        call convertUnitsCTEM(nltest,nmtest)

        ! Daily outputs from biogeochem (CTEM)
        if (doDayOutput .and. &
            (runyr >= jdsty) .and. &
            (runyr <= jdendy) .and. &
            (iday   >= jdstd) .and. &
            (iday   <= jdendd)) call ctem_daily_aw(lonLocalIndex, latLocalIndex, nltest, &
                                                  nmtest, iday, ncount, nday, &
                                                  runyr, grclarea, ipeatlandrow)

        ! Monthly biogeochem outputs
        if (doMonthOutput .and. &
            (runyr >= jmosty)) call ctem_monthly_aw(lonLocalIndex, latLocalIndex, nltest, &
                                                    nmtest, iday, runyr, nday, lastDOY)

        ! Annual biogeochem outputs
        if (doAnnualOutput) call ctem_annual_aw(lonLocalIndex, latLocalIndex, iday, runyr, nltest, nmtest, lastDOY)
        
      end if
        
      ! Increment the consecDays if it is the end of the day. This  
      ! is used to create timestamps for the output files.
      if (NCOUNT == NDAY) consecDays = consecDays + 1.

      ! Check if it is the last timestep of the last day of the year
      if ((IDAY == lastDOY) .and. (NCOUNT == NDAY)) then

        write( * , * )'IYEAR = ',IYEAR,'runyr = ',runyr,'Loop count = ',lopcount,'/',metLoop

        ! Write to the restart file
        call write_restart(lonIndex,latIndex)

        ! Increment the runyr
        runyr = runyr + 1

        ! Reset the month
        imonth = 0

      end if ! last day of year check

      ! Increment the counter for timestep of the day, or reset it to 1.
      NCOUNT = NCOUNT + 1
      if (NCOUNT > NDAY) then
        NCOUNT = 1
      end if

      !> Now check if the met file is done, needs to loop more, or just continues to the next timestep
      if (.not. metDone) then
        !> Increment the metTimeIndex to advance to the next timestep on the next time around
        metTimeIndex = metTimeIndex + 1
      else if (metDone) then
        !> End of met array read-in reached, decide what to do
        if (lopcount == metLoop) then
          !> The lopcount is reached so the run must be over
          run_model = .false.
        else
          !> Loop again so reset the metTimeIndex
          lopcount = lopcount + 1
          metTimeIndex = 1
        end if
      end if


    end do mainModelLoop ! MAIN MODEL LOOP

    ! If we've enabled checksums, now is the time to calculate them.
    if (c_switch%doChecksums) call checksumCalc(lonIndex, latIndex)

    ! deallocate arrays used for input files
    call deallocInput

    ! This is stored in outputManager so that means it retains its values between
    ! grid cells run. It must be reset here to ensure the value doesn't carry over to
    ! the next grid cell !
    consecDays = 1.

    end associate
    return

  end subroutine main_driver
  !! @}
  !> \namespace main
  !> Main model driver for CLASSIC in stand-alone mode using specified boundary
  !! conditions and atmospheric forcing.
  !!
  !! This driver program initializes the run, reads in CLASSIC input files,
  !! calls the time-stepping subroutine (mainCore), calls subroutines
  !! that aggregate and write outputs, and closes the run for this grid cell.

end module main
