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_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
                            evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
  use comsoil_h_pem,  only: soil_pem, fluxgeo, water_reservoir_nom, depth_breccia, depth_bedrock, reg_thprop_dependp 
  use adsorption_mod, only: adsorption_pem
  use glaciers_mod,   only: co2glaciersflow, h2oglaciersflow
  use ice_table_mod,  only: icetable_equilibrium, icetable_dynamic 

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

  character(len=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 = 'tmp_PEMyears.txt',status = 'old',form = 'formatted',iostat = ierr)
if (ierr /= 0) then
    write(*,*) 'Cannot find required file "tmp_PEMyears.txt"!'
    write(*,*) 'It should be created by the launching script...'
    stop
else
    read(100,*) i_myear, n_myears, convert_years
endif
close(100)

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

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

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

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

  dt_pem=1
  call getin('dt_pem', dt_pem)

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

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

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

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

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

! #---------- Layering parameters
  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  must be used when soil_pem = T'
       call abort_physic(modname,"Ice table  must be used when soil_pem = T",1)
  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.gt.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. year_bp_ini == 0.) then
     write(*,*)  'You want to follow the file obl_ecc_lsp.asc for changing orb 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 

  water_reservoir_nom = 1e4
  call getin('water_reservoir_nom',water_reservoir_nom)

  END SUBROUTINE conf_pem

END MODULE conf_pem_mod