MODULE config_pem
!-----------------------------------------------------------------------
! NAME
!     config_pem
!
! DESCRIPTION
!     Read and apply parameters for a PEM run from run_pem.def.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE read_rundef(i_myear,nyears_tot)
!-----------------------------------------------------------------------
! NAME
!     read_rundef
!
! DESCRIPTION
!     Read PEM runtime configuration from getin keys and launchPEM.info,
!     then set module-scoped parameters accordingly.
!
! AUTHORS & DATE
!     R. Vandemeulebrouck
!     JB Clement, 2023-2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
#ifdef CPP_IOIPSL
    use IOIPSL,          only: getin
#else
    ! if not using IOIPSL, we still need to use (a local version of) getin
    use ioipsl_getincom, only: getin
#endif
use evolution,        only: year_bp_ini, dt, nyears_max, evol_orbit, var_obl, var_ecc, var_lsp, convert_years
use stopping_crit,    only: h2oice_crit, co2ice_crit, ps_crit
use soil,             only: do_soil, fluxgeo, h2oice_huge_ini, depth_breccia, depth_bedrock, reg_thprop_dependp
use sorption,         only: do_sorption
use glaciers,         only: h2oice_flow, co2ice_flow
use surf_ice,         only: h2oice_cap_threshold
use ice_table,        only: icetable_equilibrium, icetable_dynamic
use layered_deposits, only: layering_algo, d_dust, impose_dust_ratio, dust2ice_ratio
use info_pem,         only: iPCM, iPEM, nPCM, nPCM_ini
use outputs,          only: output_rate

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real, intent(out) :: i_myear, nyears_tot ! Current simulated Martian year and maximum number of Martian years to be simulated

! LOCAL VARIABLES
! ---------------
character(20), parameter :: modname ='read_rundef'
integer                  :: ierr
integer                  :: year_earth_bp_ini ! Initial year (in Earth years) of the simulation of the PEM defined in run.def

! CODE
! ----
! #---------- Martian years parameters from launching script ----------#
open(100,file = 'launchPEM.info',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
if (ierr /= 0) then
    write(*,*) 'Cannot find required file "launchPEM.info"!'
    write(*,*) 'It should be created by the launching script...'
    stop
else
    read(100,*) i_myear, nyears_tot, convert_years, iPCM, iPEM, nPCM, nPCM_ini
endif
close(100)

!#---------- Output parameters ----------#
output_rate = 1 ! Default value: every year
call getin('output_rate',output_rate)

!#---------- Orbital parameters ----------#
evol_orbit = .false.
call getin('evol_orbit',evol_orbit)

year_earth_bp_ini = 0.
call getin('year_earth_bp_ini',year_earth_bp_ini)
year_bp_ini = year_earth_bp_ini/convert_years

var_obl = .true.
call getin('var_obl',var_obl)
write(*,*) 'Does obliquity vary?',var_obl

var_ecc = .true.
call getin('var_ecc',var_ecc)
write(*,*) 'Does eccentricity vary?',var_ecc

var_lsp = .true.
call getin('var_lsp',var_lsp)
write(*,*) 'Does Ls peri vary?',var_lsp

!#---------- Stopping criteria parameters ----------#
nyears_max = 100000000.
call getin('nyears_max',nyears_max)
write(*,*) 'nyears_max =',nyears_max

h2oice_crit = 0.2
call getin('h2oice_crit',h2oice_crit)
write(*,*) 'h2oice_crit =',h2oice_crit

co2ice_crit = 0.2
call getin('co2ice_crit',co2ice_crit)
write(*,*) 'co2ice_crit =',co2ice_crit

ps_crit = 0.15
call getin('ps_crit',ps_crit)
write(*,*) 'ps_crit =',ps_crit

dt = 1.
call getin('dt',dt)

!#---------- Subsurface parameters ----------#
do_soil = .true.
call getin('do_soil',do_soil)

do_sorption = .false.
call getin('do_sorption',do_sorption)

reg_thprop_dependp = .false.
call getin('reg_thprop_dependp',reg_thprop_dependp)
write(*,*) 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp

fluxgeo = 0.
call getin('fluxgeo',fluxgeo)
write(*,*) 'Flux Geothermal is set to',fluxgeo

depth_breccia = 10.
call getin('depth_breccia',depth_breccia)
write(*,*) 'Depth of breccia is set to',depth_breccia

depth_bedrock = 1000.
call getin('depth_bedrock',depth_bedrock)
write(*,*) 'Depth of bedrock is set to',depth_bedrock

icetable_equilibrium = .true.
call getin('icetable_equilibrium',icetable_equilibrium)
write(*,*) 'Is the ice table computed at equilibrium?', icetable_equilibrium

icetable_dynamic = .false.
call getin('icetable_dynamic',icetable_dynamic)
write(*,*) 'Is the ice table computed with the dynamic method?', icetable_dynamic
if ((.not. do_soil) .and. (icetable_equilibrium .or. icetable_dynamic)) then
    write(*,*) 'Ice table (equilibrium or dynamic method) must be used when do_soil = T'
    call abort_physic(modname,"Ice table must be used when do_soil = T",1)
endif
if (icetable_equilibrium .and. icetable_dynamic) then
    write(*,*) 'Ice table is asked to be computed both by the equilibrium and dynamic method.'
    write(*,*) 'The dynamic method is then chosen.'
    icetable_equilibrium = .false.
endif

if ((.not. do_soil) .and. do_sorption) then
    write(*,*) 'Adsorption must be used when do_soil = T'
    call abort_physic(modname,"Adsorption must be used when do_soil = T",1)
endif

if ((.not. do_soil) .and. (fluxgeo > 0.)) then
    write(*,*) 'Soil is not activated but Flux Geo > 0.'
    call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)
endif

if ((.not. do_soil) .and. reg_thprop_dependp) then
    write(*,*) 'Regolith properties vary with Ps only when soil is set to true'
    call abort_physic(modname,'Regolith properties vary with Ps only when soil is set to true',1)
endif

if (evol_orbit .and. abs(year_bp_ini) < 1.e-10) then
    write(*,*) 'You want to follow the file "obl_ecc_lsp.asc" for changing orbital parameters,'
    write(*,*) 'but you did not specify from which year to start.'
    call abort_physic(modname,"evol_orbit=.true. but year_bp_ini = 0",1)
endif

!#---------- Ice management parameters ----------#
h2oice_huge_ini = 9200. ! kg.m-2 (= 10 m)
call getin('h2oice_huge_ini',h2oice_huge_ini)

h2oice_cap_threshold = 460. ! kg.m-2 (= 0.5 m)
call getin('h2oice_cap_threshold',h2oice_cap_threshold)

h2oice_flow = .true.
call getin('h2oice_flow',h2oice_flow)

co2ice_flow = .true.
call getin('co2ice_flow',co2ice_flow)

!#---------- Layering parameters ----------#
layering_algo = .false.
call getin('layering',layering_algo)

d_dust = 1.e-3 ! kg.m-2.y-1
call getin('d_dust',d_dust)

impose_dust_ratio = .false.
call getin('impose_dust_ratio',impose_dust_ratio)

dust2ice_ratio = 0.01
call getin('dust2ice_ratio',dust2ice_ratio)

END SUBROUTINE read_rundef
!=======================================================================

END MODULE config_pem