source: trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90 @ 3152

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, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
use adsorption_mod,   only: adsorption_pem
use glaciers_mod,     only: co2_ice_flow, h2o_ice_flow
use ice_table_mod,    only: icetable_equilibrium, icetable_dynamic
use time_phylmdz_mod, only: ecritphy

implicit none

integer, 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
endif
close(100)

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

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

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)

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

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

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

ini_h2o_bigreservoir = 1.e4
call getin('ini_h2o_bigreservoir',ini_h2o_bigreservoir)

END SUBROUTINE conf_pem

END MODULE conf_pem_mod