source: trunk/LMDZ.COMMON/libf/evolution/config.F90 @ 4066

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

! DEPENDENCIES
! ------------
use numerics, only: dp, di, k4, minieps

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

character(7),  parameter :: rundef_name = 'run.def'
character(11), parameter :: runPCMdef_name = 'run_PCM.def'
character(12), parameter :: callphys_name = 'callphys.def'

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

!=======================================================================
SUBROUTINE read_rundef()
!-----------------------------------------------------------------------
! NAME
!     read_rundef
!
! DESCRIPTION
!     Read PEM runtime configuration from getin keys, 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: set_evolution_config
use stopping_crit,    only: set_stopping_crit_config
use soil,             only: set_soil_config
use sorption,         only: set_sorption_config
use glaciers,         only: set_glaciers_config
use surf_ice,         only: set_surf_ice_config
use ice_table,        only: set_ice_table_config
use layered_deposits, only: set_layered_deposits_config
use output,           only: set_output_config
use orbit,            only: set_orbit_config
use atmosphere,       only: set_atmosphere_config
use stoppage,         only: stop_clean
use display,          only: print_msg

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

! LOCAL VARIABLES
! ---------------
logical(k4) :: here, evol_orbit_l, evol_obl_l, evol_ecc_l, evol_lsp_l, do_soil_l, reg_thprop_dependp_l, do_sorption_l, hybrid_alt_coord_l
logical(k4) :: icetable_equilibrium_l, icetable_dynamic_l, h2oice_flow_l, co2ice_flow_l, do_layering_l, impose_dust_ratio_l
integer(di) :: output_rate_l
integer(di) :: pem_ini_earth_date ! Initial year (in Earth years) of the simulation of the PEM defined in run.def
real(dp)    :: dt_l, nmax_yr_run_l, h2oice_crit_l, co2ice_crit_l, ps_crit_l, max_change_obl_l, max_change_ecc_l, max_change_lsp_l
real(dp)    :: flux_geo_l, depth_breccia_l, depth_bedrock_l, threshold_h2oice_cap_l, h2oice_huge_ini_l, d_dust_l, dust2ice_ratio_l

! CODE
! ----
inquire(file = rundef_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//rundef_name//'" (PEM)!',1)
call print_msg('> Reading "'//rundef_name//'"')

! Output
! ~~~~~~
output_rate_l = 1_di ! Default value: every year
call getin('output_rate',output_rate_l)

! Orbital parameters
! ~~~~~~~~~~~~~~~~~~
evol_orbit_l = .false.
call getin('evol_orbit',evol_orbit_l)

pem_ini_earth_date = 0_di
call getin('pem_ini_earth_date',pem_ini_earth_date)

evol_obl_l = .true.
call getin('evol_obl',evol_obl_l)

evol_ecc_l = .true.
call getin('evol_ecc',evol_ecc_l)

evol_lsp_l = .true.
call getin('evol_lsp',evol_lsp_l)

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

! Stopping criteria
! ~~~~~~~~~~~~~~~~~
nmax_yr_run_l = 100000000._dp
call getin('nmax_yr_run',nmax_yr_run_l)

h2oice_crit_l = 0.2_dp
call getin('h2oice_crit',h2oice_crit_l)

co2ice_crit_l = 0.2_dp
call getin('co2ice_crit',co2ice_crit_l)

ps_crit_l = 0.15_dp
call getin('ps_crit',ps_crit_l)

max_change_obl_l = 1._dp
call getin('max_change_obl',max_change_obl_l)

max_change_ecc_l = 5.e-3_dp
call getin('max_change_ecc',max_change_ecc_l)

max_change_lsp_l = 20._dp
call getin('max_change_lsp',max_change_lsp_l)

! Subsurface
! ~~~~~~~~~~
do_soil_l = .true.
call getin('do_soil',do_soil_l)

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

reg_thprop_dependp_l = .false.
call getin('reg_thprop_dependp',reg_thprop_dependp_l)

flux_geo_l = 0._dp
call getin('flux_geo',flux_geo_l)

depth_breccia_l = 10._dp
call getin('depth_breccia',depth_breccia_l)

depth_bedrock_l = 1000._dp
call getin('depth_bedrock',depth_bedrock_l)

icetable_equilibrium_l = .true.
call getin('icetable_equilibrium',icetable_equilibrium_l)

icetable_dynamic_l = .false.
call getin('icetable_dynamic',icetable_dynamic_l)

! Ice management
! ~~~~~~~~~~~~~~
h2oice_huge_ini_l = 9200._dp ! [kg/m^2]; Default = 10 m
call getin('h2oice_huge_ini',h2oice_huge_ini_l)

threshold_h2oice_cap_l = 460._dp ! kg.m-2 (= 0.5 m)
call getin('threshold_h2oice_cap',threshold_h2oice_cap_l)

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

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

! Layering
! ~~~~~~~~
do_layering_l = .false.
call getin('do_layering',do_layering_l)

d_dust_l = 5.78e-2_dp ! kg.m-2.y-1 (= 1.e-9 kg.m-2.s-1)
call getin('d_dust',d_dust_l)

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

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

! Setters
call set_output_config(output_rate_l)
call set_orbit_config(evol_orbit_l,evol_obl_l,evol_ecc_l,evol_lsp_l,max_change_obl_l,max_change_ecc_l,max_change_lsp_l)
call set_evolution_config(pem_ini_earth_date,dt_l,nmax_yr_run_l)
call set_stopping_crit_config(h2oice_crit_l,co2ice_crit_l,ps_crit_l)
call set_soil_config(do_soil_l,reg_thprop_dependp_l,flux_geo_l,depth_breccia_l,depth_bedrock_l)
call set_sorption_config(do_sorption_l)
call set_ice_table_config(icetable_equilibrium_l,icetable_dynamic_l)
call set_surf_ice_config(threshold_h2oice_cap_l,h2oice_huge_ini_l)
call set_glaciers_config(h2oice_flow_l,co2ice_flow_l)
call set_layered_deposits_config(do_layering_l,impose_dust_ratio_l,d_dust_l,dust2ice_ratio_l)

! Read "run_PCM.def" parameters
hybrid_alt_coord_l = .true. ! Default setting
call get_hybrid(hybrid_alt_coord_l)
call set_atmosphere_config(hybrid_alt_coord_l)

! Check incompatibilities
call check_config_incompatibility()

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

!=======================================================================
SUBROUTINE check_config_incompatibility()
!-----------------------------------------------------------------------
! NAME
!     check_config_incompatibility
!
! DESCRIPTION
!     Check incompatibilities in the PEM runtime configuration.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage,  only: stop_clean
use soil,      only: do_soil, reg_thprop_dependp, flux_geo
use sorption,  only: do_sorption
use orbit,     only: evol_orbit
use evolution, only: pem_ini_date
use ice_table, only: icetable_equilibrium, icetable_dynamic
use display,   only: print_msg

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

! CODE
! ----
! Warnings (possible incompatibilities)
if (evol_orbit .and. abs(pem_ini_date) < minieps) call print_msg('Warning: evol_orbit = .true. but the initial date of the PEM is set to 0!')

! Errors (true incompatibilities)
if (.not. do_soil) then
    if (icetable_equilibrium .or. icetable_dynamic) call stop_clean(__FILE__,__LINE__,'ice table must be used when do_soil = true!',1)
    if (abs(flux_geo) > minieps) call stop_clean(__FILE__,__LINE__,'soil is not activated but flux_geo /= 0!',1)
    if (reg_thprop_dependp) call stop_clean(__FILE__,__LINE__,'regolith properties vary according to Ps only when soil = true!',1)
    if (do_sorption) call stop_clean(__FILE__,__LINE__,'do_soil must be true when do_sorption = true!',1)
end if

END SUBROUTINE check_config_incompatibility
!=======================================================================

!=======================================================================
SUBROUTINE read_callphys()
!-----------------------------------------------------------------------
! NAME
!     read_callphys
!
! DESCRIPTION
!     Read physics runtime configuration from getin keys, then set
!     module-scoped parameters accordingly.
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!     To work, it needs that "run_PEM.def" hols a line with
!     "INCLUDEDEF=callphys.def".
!-----------------------------------------------------------------------

! 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 atmosphere,      only: CO2cond_ps_PCM
use stoppage,        only: stop_clean
use display,         only: print_msg

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

! LOCAL VARIABLES
! ---------------
logical(k4) :: here

! CODE
! ----
inquire(file = callphys_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//callphys_name//'"!',1)
call print_msg('> Reading "'//callphys_name//'"')

CO2cond_ps_PCM = 1._dp ! Default value
call getin("CO2cond_ps",CO2cond_ps_PCM)
if (CO2cond_ps_PCM < 0._dp .or. CO2cond_ps_PCM > 1._dp) call stop_clean(__FILE__,__LINE__,'Value for ''CO2cond_ps'' is not between 0 and 1. Please, specify a correct value in "'//callphys_name//'"!',1)

END SUBROUTINE read_callphys
!=======================================================================

!=======================================================================
SUBROUTINE get_hybrid(hybrid)
!-----------------------------------------------------------------------
! NAME
!     get_hybrid
!
! DESCRIPTION
!     Get the key definition in "run_PCM.def".
!
! AUTHORS & DATE
!     JB Clement, 12/2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean
use display,  only: print_msg
use utility,  only: bool2str

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

! ARGUMENTS
! ---------
logical(k4), intent(inout) :: hybrid

! LOCAL VARIABLES
! ---------------
integer(di)    :: ierr, funit, eq_pos
logical(k4)    :: here, found, hybrid_in
character(256) :: line
character(128) :: key, res

! CODE
! ----
inquire(file = runPCMdef_name,exist = here)
if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find required file "'//runPCMdef_name//'"!',1)
call print_msg('> Reading "'//runPCMdef_name//'"')
open(newunit = funit,file = runPCMdef_name,status = 'old',form = 'formatted',action = 'read',iostat = ierr)
if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//runPCMdef_name//'"!',ierr)

found = .false.
hybrid_in = hybrid
do
    ! Read the line
    read(funit,'(a)',iostat = ierr) line
    if (ierr /= 0) exit

    ! Skip empty lines and comments
    if (trim(line) == '') cycle
    if (line(1:1) == '#' .or. line(1:1) == '!') cycle

    ! Get the position of equal sign
    eq_pos = index(line,'=')
    if (eq_pos == 0) cycle

    ! Get the key and its value/result
    key = adjustl(trim(line(:eq_pos - 1)))
    res = adjustl(trim(line(eq_pos + 1:)))

    ! Check the key
    if (trim(key) == 'hybrid') then
        read(res,*,iostat = ierr) hybrid
        if (ierr == 0) found = .true.
        exit
    end if
end do

close(funit)

if (.not. found) call print_msg('Warning: key ''hybrid'' not found in the file "'//runPCMdef_name//'"!')
if (hybrid .eqv. hybrid_in) call print_msg('USING DEFAULTS : hybrid = '//bool2str(hybrid_in))

END SUBROUTINE get_hybrid
!=======================================================================

!=======================================================================
SUBROUTINE read_controldata()
!-----------------------------------------------------------------------
! NAME
!     read_controldata
!
! DESCRIPTION
!     Read 'controle' data in "startfi.nc".
!
! AUTHORS & DATE
!     JB Clement, 01/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use io_netcdf, only: open_nc, close_nc, startfi_name, get_dim_nc, get_var_nc
use physics,   only: init_physics
use orbit,     only: init_orbit
use soil,      only: volcapa
use stoppage,  only: stop_clean
use display,   only: print_msg

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

! LOCAL VARIABLES
! ---------------
integer(di)                         :: nindex ! Size of dimension 'index'
real(dp), dimension(:), allocatable :: controle

! CODE
! ----
call print_msg('> Reading control data ("'//startfi_name//'")')
! Open the "startfi.nc" file
call open_nc(startfi_name,'read')

! Get the dimension size of 'index'
call get_dim_nc('index',nindex)

! Get 'controle'
allocate(controle(nindex))
call get_var_nc('controle',controle)

! Close the file
call close_nc(startfi_name)

! Initialize physical data
!     Arguments order: rad, g, mugaz, rcp
call print_msg('    > Initializing physical constants')
call init_physics(controle(5),controle(7),controle(8),controle(9))

! Initialize soil data
call print_msg('    > Initializing soil parameters')
volcapa = controle(35)
if (abs(volcapa) < minieps) call stop_clean(__FILE__,__LINE__,'volcapa is 0 in "'//startfi_name//'"!',1)

! Initialize orbital data
!     Arguments order: Obliquity, Perihelion, Aphelion, Date of perihelion, Year length, Sol length
call print_msg('    > Initializing orbital characteristics of the planet')
call init_orbit(controle(18),controle(15),controle(16),controle(17),controle(14),controle(10))

deallocate(controle)

END SUBROUTINE read_controldata
!=======================================================================

END MODULE config