MODULE conf_pem_mod

!=======================================================================
!
! Purpose: Read the parameter for a PEM run in run_pem.def
!
! Author: RV, JBC
!=======================================================================

implicit none

!=======================================================================
contains
!=======================================================================

SUBROUTINE conf_pem(i_myear,n_myears)

#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 time_evol_mod,         only: year_bp_ini, dt, h2o_ice_crit, co2_ice_crit, ps_criterion, Max_iter_pem, &
                                 evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years, ecritpem
use comsoil_h_pem,         only: soil_pem, fluxgeo, ini_huge_h2oice, depth_breccia, depth_bedrock, reg_thprop_dependp
use adsorption_mod,        only: adsorption_pem
use glaciers_mod,          only: h2oice_flow, co2ice_flow, inf_h2oice_threshold, metam_co2ice_threshold, metam_h2oice_threshold, metam_h2oice, metam_co2ice
use ice_table_mod,         only: icetable_equilibrium, icetable_dynamic
use layering_mod,          only: layering_algo
use info_PEM_mod,          only: iPCM, iPEM, nPCM, nPCM_ini

implicit none

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

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

!PEM parameters

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

!#---------- Output parameters ----------#
! Frequency of outputs for the PEM
ecritpem = 1 ! Default value: every year
call getin('ecritpem',ecritpem)

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

year_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 ----------#
Max_iter_pem = 100000000
call getin('Max_iter_pem',Max_iter_pem)

h2o_ice_crit = 0.2
call getin('h2o_ice_crit',h2o_ice_crit)

co2_ice_crit = 0.2
call getin('co2_ice_crit',co2_ice_crit)

ps_criterion = 0.15
call getin('ps_criterion',ps_criterion)

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

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

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

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. soil_pem) .and. (icetable_equilibrium .or. icetable_dynamic)) then
    write(*,*) 'Ice table (equilibrium or dynamic method) must be used when soil_pem = T'
    call abort_physic(modname,"Ice table must be used when soil_pem = T",1)
endif
if (icetable_equilibrium .or. 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.'
endif

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

if ((.not. soil_pem) .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. soil_pem) .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_pem .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_pem=.true. but year_bp_ini = 0",1)
endif

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

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

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

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

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

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

metam_co2ice_threshold = 16500. ! kg.m-2 (= 10 m)
call getin('metam_co2ice_threshold',metam_co2ice_threshold)

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

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

END SUBROUTINE conf_pem

END MODULE conf_pem_mod