Ignore:
Timestamp:
Dec 16, 2025, 4:39:24 PM (3 months ago)
Author:
jbclement
Message:

PEM:
Apply documentation template everywhere: standardized headers format with short description, separators between functions/subroutines, normalized code sections, aligned dependencies/arguments/variables declaration.
JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
30 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3989 r3991  
    833833- Simplifies routines and updates calls/interfaces to match the new structure.
    834834This refactor significantly clarifies the codebase and provides a cleaner foundation for forthcoming developments.
     835
     836== 16/12/2025 == JBC
     837Apply documentation template everywhere: standardized headers format with short description, separators between functions/subroutines, normalized code sections, aligned dependencies/arguments/variables declaration.

  • trunk/LMDZ.COMMON/libf/evolution/config_pem.F90

    r3989 r3991  
    11MODULE config_pem
     2!-----------------------------------------------------------------------
     3! NAME
     4!     config_pem
     5!
     6! DESCRIPTION
     7!     Read and apply parameters for a PEM run from run_pem.def.
     8!
     9! AUTHORS & DATE
     10!     R. Vandemeulebrouck
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
     16
     17! DECLARATION
     18! -----------
     19implicit none
     20
     21contains
     22!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    223
    324!=======================================================================
    4 !
    5 ! Purpose: Read the parameter for a PEM run in run_pem.def
    6 !
    7 ! Author: RV, JBC
    8 !=======================================================================
    9 
    10 implicit none
    11 
    12 !=======================================================================
    13 contains
    14 !=======================================================================
    15 
    16 SUBROUTINE read_rundef(i_myear,nyears_tots)
    17 
     25SUBROUTINE read_rundef(i_myear,nyears_tot)
     26!-----------------------------------------------------------------------
     27! NAME
     28!     read_rundef
     29!
     30! DESCRIPTION
     31!     Read PEM runtime configuration from getin keys and launchPEM.info,
     32!     then set module-scoped parameters accordingly.
     33!
     34! AUTHORS & DATE
     35!     R. Vandemeulebrouck
     36!     JB Clement, 2023-2025
     37!
     38! NOTES
     39!
     40!-----------------------------------------------------------------------
     41
     42! DEPENDENCIES
     43! ------------
    1844#ifdef CPP_IOIPSL
    1945    use IOIPSL,          only: getin
    … …  
    2248    use ioipsl_getincom, only: getin
    2349#endif
    24 
    2550use evolution,        only: year_bp_ini, dt, nyears_max, evol_orbit, var_obl, var_ecc, var_lsp, convert_years
    2651use stopping_crit,    only: h2oice_crit, co2ice_crit, ps_crit
    … …  
    3459use outputs,          only: output_rate
    3560
     61! DECLARATION
     62! -----------
    3663implicit none
    3764
    38 real, intent(out) :: i_myear, nyears_tots ! Current simulated Martian year and maximum number of Martian years to be simulated
    39 
     65! ARGUMENTS
     66! ---------
     67real, intent(out) :: i_myear, nyears_tot ! Current simulated Martian year and maximum number of Martian years to be simulated
     68
     69! LOCAL VARIABLES
     70! ---------------
    4071character(20), parameter :: modname ='read_rundef'
    4172integer                  :: ierr
    4273integer                  :: year_earth_bp_ini ! Initial year (in Earth years) of the simulation of the PEM defined in run.def
    4374
    44 !PEM parameters
    45 
    46 ! #---------- Martian years parameters from launching script
     75! CODE
     76! ----
     77! #---------- Martian years parameters from launching script ----------#
    4778open(100,file = 'launchPEM.info',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
    4879if (ierr /= 0) then
    … …  
    5182    stop
    5283else
    53     read(100,*) i_myear, nyears_tots, convert_years, iPCM, iPEM, nPCM, nPCM_ini
     84    read(100,*) i_myear, nyears_tot, convert_years, iPCM, iPEM, nPCM, nPCM_ini
    5485endif
    5586close(100)
    5687
    5788!#---------- Output parameters ----------#
    58 ! Frequency of outputs for the PEM
    5989output_rate = 1 ! Default value: every year
    6090call getin('output_rate',output_rate)
    … …  
    188218
    189219END SUBROUTINE read_rundef
     220!=======================================================================
    190221
    191222END MODULE config_pem

  • trunk/LMDZ.COMMON/libf/evolution/constants_marspem_mod.F90

    r3989 r3991  
    11MODULE constants_marspem_mod
     2!-----------------------------------------------------------------------
     3! NAME
     4!     constants_marspem_mod
     5!
     6! DESCRIPTION
     7!     Physical constants for PEM simulations on Mars
     8!
     9! AUTHORS & DATE
     10!     R. Vandemeulebrouck
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
    216
     17! DECLARATION
     18! -----------
    319implicit none
    420
     21! PARAMETERS
     22! ----------
    523! Duration of a year and day
    624integer, parameter :: sols_per_my = 669    ! Number of sols per year

  • trunk/LMDZ.COMMON/libf/evolution/evolution.F90

    r3989 r3991  
    11MODULE evolution
     2!-----------------------------------------------------------------------
     3! NAME
     4!     evolution
     5!
     6! DESCRIPTION
     7!     Contains global parameters used for the evolution flags.
     8!
     9! AUTHORS & DATE
     10!     R. Vandemeulebrouck
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
    216
     17! DECLARATION
     18! -----------
    319implicit none
    420
     21! PARAMETERS
     22! ----------
    523real    :: year_bp_ini   ! Initial year (in Planetary years) of the simulation of the PEM defined in run.def (in Earth years)
    624real    :: dt            ! Time step used by the PEM in Planetary years
    725real    :: convert_years ! Conversion ratio from Planetary years to Earth years
    8 real    :: nyears_max      ! Maximal number of iterations when converging to a steady state, read in evol.def
     26real    :: nyears_max    ! Maximal number of iterations when converging to a steady state, read in evol.def
    927logical :: evol_orbit    ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
    1028logical :: var_obl       ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity

  • trunk/LMDZ.COMMON/libf/evolution/glaciers.F90

    r3989 r3991  
    11MODULE glaciers
    2 
    3 implicit none
    4 
    5 ! Flags for ice flow
    6 logical :: h2oice_flow  ! True by default, flag to compute H2O ice flow
    7 logical :: co2ice_flow  ! True by default, flag to compute CO2 ice flow
    8 
    9 !=======================================================================
     2!-----------------------------------------------------------------------
     3! NAME
     4!     glaciers
     5!
     6! DESCRIPTION
     7!     Compute flow and transfer of CO2 and H2O ice glaciers on slopes
     8!     based on maximum ice thickness and basal drag conditions.
     9!
     10! AUTHORS & DATE
     11!     L. Lange
     12!     JB Clement, 2023-2025
     13!
     14! NOTES
     15!
     16!-----------------------------------------------------------------------
     17
     18! DECLARATION
     19! -----------
     20implicit none
     21
     22! PARAMETERS
     23! ----------
     24logical :: h2oice_flow ! True by default, flag to compute H2O ice flow
     25logical :: co2ice_flow ! True by default, flag to compute CO2 ice flow
     26
    1027contains
    11 !=======================================================================
    12 
     28!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     29
     30!=======================================================================
    1331SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,co2ice,flag_co2flow)
    14 
    15 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    16 !!!
    17 !!! Purpose: Main for CO2 glaciers evolution: compute maximum thickness, and do
    18 !!!          the ice transfer
    19 !!!
    20 !!!
    21 !!! Author: LL
    22 !!!
    23 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    24 
    25 implicit none
    26 
    27 ! Inputs
    28 !-------
    29 integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
    30 real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Grid points x Slopes: Distribution of the subgrid slopes
    31 real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Grid points: values of the sub grid slope angles
    32 real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Grid points x Time field : VMR of co2 in the first layer [mol/mol]
    33 real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Grid points x Time field: surface pressure given by the PCM [Pa]
    34 real,                              intent(in) :: ps_avg_global_ini             ! Global averaged surface pressure at the beginning [Pa]
    35 real,                              intent(in) :: ps_avg_global                 ! Global averaged surface pressure during the PEM iteration [Pa]
    36 ! Ouputs
    37 !-------
    38 real, dimension(ngrid,nslope), intent(inout) :: co2ice ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    39 integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_co2flow ! Flag to see if there is flow on the subgrid slopes
    40 ! Local
    41 !------
    42 real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
    43 real, dimension(ngrid,nslope) :: hmax  ! Grid points x Slope field: maximum thickness for co2  glacier before flow
    44 
     32!-----------------------------------------------------------------------
     33! NAME
     34!     flow_co2glaciers
     35!
     36! DESCRIPTION
     37!     Compute maximum thickness and transfer CO2 ice between subslopes.
     38!
     39! AUTHORS & DATE
     40!     L. Lange
     41!     JB Clement, 2023-2025
     42!
     43! NOTES
     44!
     45!-----------------------------------------------------------------------
     46
     47! DECLARATION
     48! -----------
     49implicit none
     50
     51! ARGUMENTS
     52! ---------
     53integer,                                  intent(in)    :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     54real, dimension(ngrid,nslope),            intent(in)    :: subslope_dist                 ! Grid points x Slopes: Distribution of the subgrid slopes
     55real, dimension(ngrid),                   intent(in)    :: def_slope_mean                ! Grid points: values of the sub grid slope angles
     56real, dimension(ngrid,timelen),           intent(in)    :: vmr_co2_PEM                   ! Grid points x Time field : VMR of CO2 in the first layer [mol/mol]
     57real, dimension(ngrid,timelen),           intent(in)    :: ps_PCM                        ! Grid points x Time field: surface pressure given by the PCM [Pa]
     58real,                                     intent(in)    :: ps_avg_global_ini             ! Global averaged surface pressure at the beginning [Pa]
     59real,                                     intent(in)    :: ps_avg_global                 ! Global averaged surface pressure during the PEM iteration [Pa]
     60real, dimension(ngrid,nslope),            intent(inout) :: co2ice                        ! Grid points x Slope field: CO2 ice on the subgrid slopes [kg/m^2]
     61integer(kind=1), dimension(ngrid,nslope), intent(out)   :: flag_co2flow                  ! Flag to see if there is flow on the subgrid slopes
     62
     63! LOCAL VARIABLES
     64! ---------------
     65real, dimension(ngrid,nslope) :: Tcond ! CO2 condensation temperature [K]
     66real, dimension(ngrid,nslope) :: hmax  ! Maximum thickness before flow
     67
     68! CODE
     69! ----
    4570write(*,*) "> Flow of CO2 glaciers"
    4671call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
    … …  
    4974
    5075END SUBROUTINE flow_co2glaciers
    51 
    52 !=======================================================================
    53 
     76!=======================================================================
     77
     78!=======================================================================
    5479SUBROUTINE flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow)
    55 
    56 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    57 !!!
    58 !!! Purpose: Main for H2O glaciers evolution: compute maximum thickness, and do
    59 !!!          the ice transfer
    60 !!!
    61 !!!
    62 !!! Author: LL
    63 !!!
    64 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    65 
    66 implicit none
    67 
    68 ! arguments
    69 ! ---------
    70 
    71 ! Inputs
    72 ! ------
    73 integer,                       intent(in) :: ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
    74 real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Grid points x Slopes : Distribution of the subgrid slopes
    75 real, dimension(ngrid),        intent(in) :: def_slope_mean ! Slopes: values of the sub grid slope angles
    76 real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice Temperature [K]
    77 ! Outputs
    78 !--------
    79 real, dimension(ngrid,nslope), intent(inout) :: h2oice ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    80 integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_h2oflow ! Flag to see if there is flow on the subgrid slopes
    81 ! Local
    82 ! -----
    83 real, dimension(ngrid,nslope) :: hmax ! Grid points x Slope field: maximum thickness for co2  glacier before flow
    84 
     80!-----------------------------------------------------------------------
     81! NAME
     82!     flow_h2oglaciers
     83!
     84! DESCRIPTION
     85!     Compute maximum thickness and transfer H2O ice between subslopes.
     86!
     87! AUTHORS & DATE
     88!     L. Lange
     89!     JB Clement, 2023-2025
     90!
     91! NOTES
     92!
     93!-----------------------------------------------------------------------
     94
     95! DECLARATION
     96! -----------
     97implicit none
     98
     99! ARGUMENTS
     100! ---------
     101integer,                                  intent(in)    ::  ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     102real, dimension(ngrid,nslope),            intent(in)    ::  subslope_dist        ! Distribution of the subgrid slopes
     103real, dimension(ngrid),                   intent(in)    ::  def_slope_mean       ! Values of the sub grid slope angles
     104real, dimension(ngrid,nslope),            intent(in)    ::  Tice                 ! Ice Temperature [K]
     105real, dimension(ngrid,nslope),            intent(inout) ::  h2oice               ! H2O ice on the subgrid slopes [kg/m^2]
     106integer(kind=1), dimension(ngrid,nslope), intent(out)   ::  flag_h2oflow         ! Flow flag on subgrid slopes
     107
     108! LOCAL VARIABLES
     109! ---------------
     110real, dimension(ngrid,nslope) :: hmax ! Maximum thickness before flow
     111
     112! CODE
     113! ----
    85114write(*,*) "> Flow of H2O glaciers"
    86115call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
    … …  
    88117
    89118END SUBROUTINE flow_h2oglaciers
    90 
    91 !=======================================================================
    92 
     119!=======================================================================
     120
     121!=======================================================================
    93122SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
    94 
    95 use ice_table,  only: rho_ice
    96 use stop_pem,  only: stop_clean
     123!-----------------------------------------------------------------------
     124! NAME
     125!     compute_hmaxglaciers
     126!
     127! DESCRIPTION
     128!     Compute the maximum thickness of CO2 and H2O glaciers given a
     129!     slope angle before initiating flow.
     130!
     131! AUTHORS & DATE
     132!     L. Lange
     133!     JB Clement, 2023-2025
     134!
     135! NOTES
     136!     Based on work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
     137!
     138!-----------------------------------------------------------------------
     139
     140! DEPENDENCIES
     141! ------------
     142use ice_table,      only: rho_ice
     143use stop_pem,       only: stop_clean
    97144use phys_constants, only: pi, g
    98145
    99 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    100 !!!
    101 !!! Purpose: Compute the maximum thickness of CO2 and H2O glaciers given a slope angle before initating flow
    102 !!!
    103 !!! Author: LL, based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
    104 !!!
    105 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    106 
    107 implicit none
    108 
    109 ! Inputs
    110 ! ------
    111 integer,                       intent(in) :: ngrid, nslope  ! # of grid points and subslopes
    112 integer,                       intent(in) :: iflat          ! index of the flat subslope
    113 real, dimension(nslope),       intent(in) :: def_slope_mean ! Slope field: Values of the subgrid slope angles [deg]
    114 real, dimension(ngrid,nslope), intent(in) :: Tice           ! Physical field:  ice temperature [K]
    115 character(3),                  intent(in) :: name_ice       ! Nature of ice
    116 ! Outputs
    117 ! -------
    118 real, dimension(ngrid,nslope), intent(out) :: hmax ! Physical grid x Slope field: maximum  thickness before flaw [m]
    119 ! Local
    120 ! -----
    121 real    :: tau_d      ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
    122 integer :: ig, islope ! loop variables
     146! DECLARATION
     147! -----------
     148implicit none
     149
     150! ARGUMENTS
     151! ---------
     152integer,                       intent(in)  :: ngrid, nslope  ! # of grid points and subslopes
     153integer,                       intent(in)  :: iflat          ! index of the flat subslope
     154real, dimension(nslope),       intent(in)  :: def_slope_mean ! Values of the subgrid slope angles [deg]
     155real, dimension(ngrid,nslope), intent(in)  :: Tice           ! Ice temperature [K]
     156character(3),                  intent(in)  :: name_ice       ! Nature of ice
     157real, dimension(ngrid,nslope), intent(out) :: hmax           ! Maximum thickness before flow [m]
     158
     159! LOCAL VARIABLES
     160! ---------------
     161real    :: tau_d ! characteristic basal drag, understood as the stress that an ice mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
     162integer :: ig, islope
    123163real    :: slo_angle
    124164
     165! CODE
     166! ----
    125167select case (trim(adjustl(name_ice)))
    126168    case('h2o')
    … …  
    144186
    145187END SUBROUTINE compute_hmaxglaciers
    146 
    147 !=======================================================================
    148 
     188!=======================================================================
     189
     190!=======================================================================
    149191SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,qice,flag_flow)
    150 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    151 !!!
    152 !!! Purpose: Transfer the excess of ice from one subslope to another
    153 !!!          No transfer between mesh at the time
    154 !!! Author: LL
    155 !!!
    156 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    157 
     192!-----------------------------------------------------------------------
     193! NAME
     194!     transfer_ice_duringflow
     195!
     196! DESCRIPTION
     197!     Transfer the excess of ice from one subslope to another.
     198!     No transfer between mesh at the time.
     199!
     200! AUTHORS & DATE
     201!     L. Lange
     202!     JB Clement, 2023-2025
     203!
     204! NOTES
     205!
     206!-----------------------------------------------------------------------
     207
     208! DEPENDENCIES
     209! ------------
    158210use ice_table,      only: rho_ice
    159211use stop_pem,       only: stop_clean
    160212use phys_constants, only: pi
    161213
    162 implicit none
    163 
    164 ! Inputs
    165 !-------
    166 integer,                       intent(in) :: ngrid, nslope  ! # of physical points and subslope
    167 integer,                       intent(in) :: iflat          ! index of the flat subslope
    168 real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Distribution of the subgrid slopes within the mesh
    169 real, dimension(nslope),       intent(in) :: def_slope_mean ! values of the subgrid slopes
    170 real, dimension(ngrid,nslope), intent(in) :: hmax           ! maximum height of the  glaciers before initiating flow [m]
    171 real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice temperature[K]
    172 ! Outputs
    173 !--------
    174 real, dimension(ngrid,nslope), intent(inout) :: qice ! CO2 in the subslope [kg/m^2]
    175 integer(kind=1), dimension(ngrid,nslope), intent(out) :: flag_flow ! Flag to see if there is flow on the subgrid slopes
    176 ! Local
    177 !------
    178 integer :: ig, islope ! loop
    179 integer :: iaval      ! ice will be transfered here
    180 
     214! DECLARATION
     215! -----------
     216implicit none
     217
     218! ARGUMENTS
     219! ---------
     220integer,                                  intent(in)    :: ngrid, nslope  ! # of physical points and subslope
     221integer,                                  intent(in)    :: iflat          ! Index of the flat subslope
     222real, dimension(ngrid,nslope),            intent(in)    :: subslope_dist  ! Distribution of the subgrid slopes
     223real, dimension(nslope),                  intent(in)    :: def_slope_mean ! Values of the subgrid slopes
     224real, dimension(ngrid,nslope),            intent(in)    :: hmax           ! Maximum height before initiating flow [m]
     225real, dimension(ngrid,nslope),            intent(in)    :: Tice           ! Ice temperature [K]
     226real, dimension(ngrid,nslope),            intent(inout) :: qice           ! Ice in the subslope [kg/m^2]
     227integer(kind=1), dimension(ngrid,nslope), intent(out)   :: flag_flow      ! Flow flag on subgrid slopes
     228
     229! LOCAL VARIABLES
     230! ---------------
     231integer :: ig, islope
     232integer :: iaval ! index where ice is transferred
     233
     234! CODE
     235! ----
    181236flag_flow = 0
    182237
    … …  
    188243! Second: determine the flatest slopes possible:
    189244                if (islope > iflat) then
    190                     iaval=islope-1
     245                    iaval = islope - 1
    191246                else
    192247                    iaval = islope + 1
    … …  
    212267
    213268!=======================================================================
    214 
    215269SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
    216 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    217 !!!
    218 !!! Purpose: Compute CO2 condensation temperature
    219 !!!
    220 !!! Author: LL
    221 !!!
    222 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    223 
     270!-----------------------------------------------------------------------
     271! NAME
     272!     computeTcondCO2
     273!
     274! DESCRIPTION
     275!     Compute CO2 condensation temperature.
     276!
     277! AUTHORS & DATE
     278!     L. Lange
     279!     JB Clement, 2023-2025
     280!
     281! NOTES
     282!
     283!-----------------------------------------------------------------------
     284
     285! DEPENDENCIES
     286! ------------
    224287use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2
    225288
    226 implicit none
    227 
    228 ! arguments:
    229 ! ----------
    230 
    231 ! INPUT
    232 integer,                        intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes
    233 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Grid points x times field: VMR of CO2 in the first layer [mol/mol]
    234 real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Grid points x times field: surface pressure in the PCM [Pa]
    235 real,                           intent(in) :: ps_avg_global_ini      ! Global averaged surfacepressure in the PCM [Pa]
    236 real,                           intent(in) :: ps_avg_global          ! Global averaged surface pressure computed during the PEM iteration
    237 
    238 ! OUTPUT
    239 real, dimension(ngrid,nslope), intent(out) :: Tcond ! Grid points: condensation temperature of CO2, yearly averaged
    240 
    241 ! LOCAL
    242 integer :: ig, it ! For loop
    243 
     289! DECLARATION
     290! -----------
     291implicit none
     292
     293! ARGUMENTS
     294! ---------
     295integer,                        intent(in)  :: timelen, ngrid, nslope
     296real, dimension(ngrid,timelen), intent(in)  :: vmr_co2_PEM       ! VMR of CO2 in the first layer [mol/mol]
     297real, dimension(ngrid,timelen), intent(in)  :: ps_PCM            ! Surface pressure in the PCM [Pa]
     298real,                           intent(in)  :: ps_avg_global_ini ! Global averaged surfacepressure in the PCM [Pa]
     299real,                           intent(in)  :: ps_avg_global     ! Global averaged surface pressure computed during the PEM iteration
     300real, dimension(ngrid,nslope),  intent(out) :: Tcond             ! Condensation temperature of CO2, yearly averaged
     301
     302! LOCAL VARIABLES
     303! ---------------
     304integer :: ig, it
     305
     306! CODE
     307! ----
    244308do ig = 1,ngrid
    245309    Tcond(ig,:) = 0
    … …  
    251315
    252316END SUBROUTINE computeTcondCO2
     317!=======================================================================
    253318
    254319END MODULE glaciers

  • trunk/LMDZ.COMMON/libf/evolution/grid_conversion.F90

    r3980 r3991  
    11MODULE grid_conversion
     2!-----------------------------------------------------------------------
     3! NAME
     4!     grid_conversion
     5!
     6! DESCRIPTION
     7!     Provides tools to convert data between lon x lat grid format
     8!     and vector format.
     9!
     10! AUTHORS & DATE
     11!     JB Clement, 12/2025
     12!
     13! NOTES
     14!     Handles pole duplication differences between the grid format
     15!     and vector format.
     16!-----------------------------------------------------------------------
    217
     18! DECLARATION
     19! -----------
    320implicit none
    421
    5 !=======================================================================
    622contains
     23!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     24
    725!=======================================================================
    826SUBROUTINE lonlat2vect(nlon,nlat,ngrid,v_ll,v_vect)
    9 ! To convert data from lon x lat grid (where values at the poles are duplicated) into a vector
    10 ! /!\ The longitudes -180 and +180 are not duplicated like in the PCM dynamics
     27!-----------------------------------------------------------------------
     28! NAME
     29!     lonlat2vect
     30!
     31! DESCRIPTION
     32!     Convert data from lon x lat grid (where values at the poles are
     33!     duplicated) into a vector.
     34!
     35! AUTHORS & DATE
     36!     JB Clement, 12/2025
     37!
     38! NOTES
     39!     The longitudes -180 and +180 are not duplicated like in the PCM
     40!     dynamics.
     41!
     42!-----------------------------------------------------------------------
    1143
     44! DECLARATION
     45! -----------
    1246implicit none
    1347
    1448! Arguments
    1549!----------
    16 integer,                    intent(in) :: nlon, nlat, ngrid
    17 real, dimension(nlon,nlat), intent(in) :: v_ll
    18 real, dimension(ngrid), intent(out) :: v_vect
     50integer,                    intent(in)  :: nlon, nlat, ngrid
     51real, dimension(nlon,nlat), intent(in)  :: v_ll
     52real, dimension(ngrid),     intent(out) :: v_vect
    1953
    2054! Local variables
    … …  
    4882
    4983END SUBROUTINE lonlat2vect
     84!=======================================================================
    5085
    5186!=======================================================================
    5287SUBROUTINE vect2lonlat(nlon,nlat,ngrid,v_vect,v_ll)
    53 ! To convert data from a vector into lon x lat grid (where values at the poles are duplicated)
    54 ! /!\ The longitudes -180 and +180 are not duplicated like in the PCM dynamics
     88!-----------------------------------------------------------------------
     89! NAME
     90!     vect2lonlat
     91!
     92! DESCRIPTION
     93!     Convert data from a vector into lon x lat grid (where values
     94!     at the poles are duplicated).
     95!
     96! AUTHORS & DATE
     97!     JB Clement, 12/2025
     98!
     99! NOTES
     100!     The longitudes -180 and +180 are not duplicated like in the PCM
     101!     dynamics.
     102!-----------------------------------------------------------------------
    55103
     104! DECLARATION
     105! -----------
    56106implicit none
    57107
    58108! Arguments
    59109!----------
    60 integer,                intent(in) :: nlon, nlat, ngrid
    61 real, dimension(ngrid), intent(in) :: v_vect
     110integer,                    intent(in) :: nlon, nlat, ngrid
     111real, dimension(ngrid),     intent(in) :: v_vect
    62112real, dimension(nlon,nlat), intent(out) :: v_ll
    63113
    … …  
    92142
    93143END SUBROUTINE vect2lonlat
     144!=======================================================================
    94145
    95146END MODULE grid_conversion

  • trunk/LMDZ.COMMON/libf/evolution/ice_table.F90

    r3989 r3991  
    11MODULE ice_table
    2 
    3 implicit none
    4 
    5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    6 !!!
    7 !!! Purpose: Ice table (pore-filling) variables and methods to compute it (dynamic and static)
    8 !!! Author:  LL, 02/2023
    9 !!!
    10 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    11 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     ice_table
     5!
     6! DESCRIPTION
     7!     Ice table variables and methods to compute it (dynamic and static).
     8!
     9! AUTHORS & DATE
     10!     L. Lange, 02/2023
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
     16
     17! DECLARATION
     18! -----------
     19implicit none
     20
     21! MODULE VARIABLES
     22! ----------------
    1223logical                             :: icetable_equilibrium ! Boolean to say if the PEM needs to recompute the icetable depth when at equilibrium
    1324logical                             :: icetable_dynamic     ! Boolean to say if the PEM needs to recompute the icetable depth with the dynamic method
    14 real, allocatable, dimension(:,:)   :: icetable_depth       ! ngrid x nslope: Depth of the ice table [m]
    15 real, allocatable, dimension(:,:)   :: icetable_thickness   ! ngrid x nslope: Thickness of the ice table [m]
    16 real, allocatable, dimension(:,:,:) :: ice_porefilling      ! the amout of porefilling in each layer in each grid [m^3/m^3]
    17 
    18 !-----------------------------------------------------------------------
     25real, allocatable, dimension(:,:)   :: icetable_depth       ! Depth of the ice table [m]
     26real, allocatable, dimension(:,:)   :: icetable_thickness   ! Thickness of the ice table [m]
     27real, allocatable, dimension(:,:,:) :: ice_porefilling      ! Amount of porefilling in each layer in each grid [m^3/m^3]
     28
    1929contains
    20 
    21 !-----------------------------------------------------------------------
     30!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     31
     32!=======================================================================
    2233SUBROUTINE ini_ice_table(ngrid,nslope,nsoil)
    23 
    24 implicit none
    25 
    26 integer, intent(in) :: ngrid  ! number of atmospheric columns
    27 integer, intent(in) :: nslope ! number of slope within a mesh
    28 integer, intent(in) :: nsoil  ! number of soil layers
    29 
     34!-----------------------------------------------------------------------
     35! NAME
     36!     ini_ice_table
     37!
     38! DESCRIPTION
     39!     Allocate module arrays for ice table fields.
     40!
     41! AUTHORS & DATE
     42!     L. Lange
     43!     JB Clement, 2023-2025
     44!
     45! NOTES
     46!
     47!-----------------------------------------------------------------------
     48
     49! DECLARATION
     50! -----------
     51implicit none
     52
     53! ARGUMENTS
     54! ---------
     55integer, intent(in) :: ngrid  ! Number of atmospheric columns
     56integer, intent(in) :: nslope ! Number of slope within a mesh
     57integer, intent(in) :: nsoil  ! Number of soil layers
     58
     59! CODE
     60! ----
    3061allocate(icetable_depth(ngrid,nslope))
    3162allocate(icetable_thickness(ngrid,nslope))
    … …  
    3364
    3465END SUBROUTINE ini_ice_table
    35 
    36 !-----------------------------------------------------------------------
     66!=======================================================================
     67
     68!=======================================================================
    3769SUBROUTINE end_ice_table
    38 
    39 implicit none
    40 
     70!-----------------------------------------------------------------------
     71! NAME
     72!     end_ice_table
     73!
     74! DESCRIPTION
     75!     Deallocate ice table arrays.
     76!
     77! AUTHORS & DATE
     78!     L. Lange
     79!     JB Clement, 2023-2025
     80!
     81! NOTES
     82!
     83!-----------------------------------------------------------------------
     84
     85! DECLARATION
     86! -----------
     87implicit none
     88
     89! CODE
     90! ----
    4191if (allocated(icetable_depth)) deallocate(icetable_depth)
    4292if (allocated(icetable_thickness)) deallocate(icetable_thickness)
    … …  
    4494
    4595END SUBROUTINE end_ice_table
    46 
    47 !------------------------------------------------------------------------------------------------------
     96!=======================================================================
     97
     98!=======================================================================
    4899SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_depth,ice_table_thickness)
    49 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    50 !!!
    51 !!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water
    52 !!! density at the surface and at depth.
    53 !!! Computations are made following the methods in Schorgofer et al., 2005
    54 !!! This subroutine only gives the ice table at equilibrium and does not consider exchange with the atmosphere
    55 !!!
    56 !!! Author: LL
    57 !!!
    58 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     100!-----------------------------------------------------------------------
     101! NAME
     102!     computeice_table_equilibrium
     103!
     104! DESCRIPTION
     105!     Compute the ice table depth knowing the yearly average water
     106!     density at the surface and at depth. Computations are made following
     107!     the methods in Schorghofer et al., 2005.
     108!
     109! AUTHORS & DATE
     110!     L. Lange
     111!     JB Clement, 12/12/2025
     112!
     113! NOTES
     114!     This subroutine only gives the ice table at equilibrium and does
     115!     not consider exchange with the atmosphere.
     116!-----------------------------------------------------------------------
     117
     118! DEPENDENCIES
     119! ------------
    59120use math_toolkit,         only: findroot
    60121use soil,                 only: mlayer_PEM, layer_PEM ! Depth of the vertical grid
    61122use soil_thermal_inertia, only: get_ice_TI
    62123
    63 implicit none
    64 
    65 ! Inputs
    66 ! ------
    67 integer,                             intent(in) :: ngrid, nslope, nsoil ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    68 logical, dimension(ngrid),           intent(in) :: watercaptag          ! Boolean to check the presence of a perennial glacier
    69 real, dimension(ngrid,nslope),       intent(in) :: rhowatersurf_ave     ! Water density at the surface, yearly averaged [kg/m^3]
    70 real, dimension(ngrid,nsoil,nslope), intent(in) :: rhowatersoil_ave     ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3]
    71 real, dimension(ngrid,nslope),       intent(in) :: regolith_inertia     ! Thermal inertia of the regolith layer [SI]
    72 ! Ouputs
    73 ! ------
    74 real, dimension(ngrid,nslope), intent(out) :: ice_table_depth     ! ice table depth [m]
    75 real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m]
    76 ! Local variables
     124! DECLARATION
     125! -----------
     126implicit none
     127
     128! ARGUMENTS
     129! ---------
     130integer,                                intent(in)  :: ngrid, nslope, nsoil ! Size of the physical grid, number of subslope, number of soil layer in the PEM
     131logical, dimension(ngrid),              intent(in)  :: watercaptag          ! Boolean to check the presence of a perennial glacier
     132real,    dimension(ngrid,nslope),       intent(in)  :: rhowatersurf_ave     ! Water density at the surface, yearly averaged [kg/m^3]
     133real,    dimension(ngrid,nsoil,nslope), intent(in)  :: rhowatersoil_ave     ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged [kg/m^3]
     134real,    dimension(ngrid,nslope),       intent(in)  :: regolith_inertia     ! Thermal inertia of the regolith layer [SI]
     135real,    dimension(ngrid,nslope),       intent(out) :: ice_table_depth      ! Ice table depth [m]
     136real,    dimension(ngrid,nslope),       intent(out) :: ice_table_thickness  ! Ice table thickness [m]
     137
     138! LOCAL VARIABLES
    77139! ---------------
    78 integer                       :: ig, islope, isoil, isoilend ! loop variables
    79 real, dimension(nsoil)        :: diff_rho                    ! difference of water vapor density between the surface and at depth [kg/m^3]
    80 real                          :: ice_table_end               ! depth of the end of the ice table  [m]
     140integer                       :: ig, islope, isoil, isoilend
     141real, dimension(nsoil)        :: diff_rho                    ! Difference of water vapor density between the surface and at depth [kg/m^3]
     142real                          :: ice_table_end               ! Depth of the end of the ice table  [m]
    81143real, dimension(ngrid,nslope) :: previous_icetable_depth     ! Ice table computed at previous ice depth [m]
    82 real                          :: stretch                     ! stretch factor to improve the convergence of the ice table
     144real                          :: stretch                     ! Stretch factor to improve the convergence of the ice table
    83145real                          :: wice_inertia                ! Water Ice thermal Inertia [USI]
    84 ! Code
     146
     147! CODE
    85148! ----
    86149previous_icetable_depth = ice_table_depth
    … …  
    139202
    140203END SUBROUTINE computeice_table_equilibrium
    141 
    142 !-----------------------------------------------------------------------
     204!=======================================================================
     205
     206!=======================================================================
    143207SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil,icetable_thickness_old,ice_porefilling_old,tsurf,tsoil,delta_m_h2o)
    144 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    145 !!!
    146 !!! Purpose: Compute the mass of H2O that has sublimated from the ice table / condensed
    147 !!!
    148 !!! Author: LL
    149 !!!
    150 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    151 use soil,         only: mlayer_PEM
     208!-----------------------------------------------------------------------
     209! NAME
     210!     compute_massh2o_exchange_ssi
     211!
     212! DESCRIPTION
     213!     Compute the mass of H2O that has sublimated from the ice table /
     214!     condensed.
     215!
     216! AUTHORS & DATE
     217!     L. Lange
     218!     JB Clement, 2023-2025
     219!
     220! NOTES
     221!
     222!-----------------------------------------------------------------------
     223
     224! DEPENDENCIES
     225! ------------
     226use soil,                  only: mlayer_PEM
    152227use comslope_mod,          only: subslope_dist, def_slope_mean
    153228use constants_marspem_mod, only: porosity
    154229use phys_constants,        only: pi
    155230
    156 implicit none
    157 
    158 ! Inputs
    159 ! ------
    160 integer,                             intent(in) :: ngrid, nslope, nsoil   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    161 real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness_old ! Ice table thickness at the previous iteration [m]
    162 real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling_old    ! Ice pore filling at the previous iteration [m]
    163 real, dimension(ngrid,nslope),       intent(in) :: tsurf                  ! Surface temperature [K]
    164 real, dimension(ngrid,nsoil,nslope), intent(in) :: tsoil                  ! Soil temperature [K]
    165 ! Outputs
    166 ! -------
    167 real, dimension(ngrid), intent(out) :: delta_m_h2o ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]
    168 ! Locals
    169 !-------
    170 integer :: ig, islope, isoil ! loop index
    171 real    :: rho               ! density of water ice [kg/m^3]
    172 real    :: Tice              ! ice temperature [k]
    173 ! Code
    174 ! ----
    175 
     231! DECLARATION
     232! -----------
     233implicit none
     234
     235
     236! ARGUMENTS
     237! ---------
     238integer,                                intent(in)  :: ngrid, nslope, nsoil
     239real,    dimension(ngrid,nslope),       intent(in)  :: icetable_thickness_old ! Ice table thickness at the previous iteration [m]
     240real,    dimension(ngrid,nsoil,nslope), intent(in)  :: ice_porefilling_old    ! Ice pore filling at the previous iteration [m]
     241real,    dimension(ngrid,nslope),       intent(in)  :: tsurf                  ! Surface temperature [K]
     242real,    dimension(ngrid,nsoil,nslope), intent(in)  :: tsoil                  ! Soil temperature [K]
     243real,    dimension(ngrid),              intent(out) :: delta_m_h2o            ! Mass of H2O ice that has been condensed on the ice table / sublimates from the ice table [kg/m^2]
     244
     245! LOCAL VARIABLES
     246! ---------------
     247integer :: ig, islope, isoil
     248real    :: rho  ! Density of water ice [kg/m^3]
     249real    :: Tice ! Ice temperature [k]
     250
     251! CODE
     252! ----
    176253! Let's compute the amount of ice that has sublimated in each subslope
    177254if (icetable_equilibrium) then
    … …  
    198275
    199276END SUBROUTINE compute_massh2o_exchange_ssi
    200 
    201 !-----------------------------------------------------------------------
     277!=======================================================================
     278
     279!=======================================================================
    202280SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
    203 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    204 !!!
    205 !!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM
    206 !!!
    207 !!! Author: LL
    208 !!!
    209 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    210 use soil, only: nsoilmx_PEM, mlayer_PEM
    211 use math_toolkit,  only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
    212 
    213 implicit none
    214 
    215 ! Inputs
    216 ! ------
    217 real, dimension(nsoilmx_PEM), intent(in) :: porefill    ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
    218 real, dimension(nsoilmx_PEM), intent(in) :: psat_soil   ! Soil water pressure at saturation, yearly averaged [Pa]
    219 real,                         intent(in) :: psat_surf   ! surface water pressure at saturation, yearly averaged [Pa]
    220 real,                         intent(in) :: pwat_surf   ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
    221 real,                         intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
    222 real,                         intent(in) :: B           ! constant (Eq. 8 from  Schorgofer, Icarus (2010).)
    223 integer,                      intent(in) :: index_IS    ! index of the soil layer where the ice sheet begins [1]
    224 real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
    225 ! Outputs
    226 ! -------
    227 integer,                      intent(out) :: index_filling    ! index where the pore filling begins [1]
    228 integer,                      intent(out) :: index_geothermal ! index where the ice table stops [1]
    229 real,                         intent(out) :: depth_geothermal ! depth where the ice table stops [m]
    230 real, dimension(nsoilmx_PEM), intent(out) :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
    231 ! Locals
    232 !-------
    233 real, dimension(nsoilmx_PEM) :: eta                          ! constriction
    234 real, dimension(nsoilmx_PEM) :: dz_psat                      ! first derivative of the vapor pressure at saturation
     281!-----------------------------------------------------------------------
     282! NAME
     283!     find_layering_icetable
     284!
     285! DESCRIPTION
     286!     Compute layering between dry soil, pore filling ice and ice
     287!     sheet based on Schorghofer, Icarus (2010). Adapted from NS MSIM.
     288!
     289! AUTHORS & DATE
     290!     L. Lange
     291!     JB Clement, 2023-2025
     292!
     293! NOTES
     294!
     295!-----------------------------------------------------------------------
     296
     297! DEPENDENCIES
     298! ------------
     299use soil,         only: nsoilmx_PEM, mlayer_PEM
     300use math_toolkit, only: deriv1, deriv1_onesided, colint, findroot, deriv2_simple
     301
     302! DECLARATION
     303! -----------
     304implicit none
     305
     306! ARGUMENTS
     307! ---------
     308real, dimension(nsoilmx_PEM), intent(in)    :: porefill         ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
     309real, dimension(nsoilmx_PEM), intent(in)    :: psat_soil        ! Soil water pressure at saturation, yearly averaged [Pa]
     310real,                         intent(in)    :: psat_surf        ! Surface water pressure at saturation, yearly averaged [Pa]
     311real,                         intent(in)    :: pwat_surf        ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
     312real,                         intent(in)    :: psat_bottom      ! Boundary conditions for soil vapor pressure [Pa]
     313real,                         intent(in)    :: B                ! Constant (Eq. 8 from  Schorgofer, Icarus (2010).)
     314integer,                      intent(in)    :: index_IS         ! Index of the soil layer where the ice sheet begins [1]
     315real,                         intent(inout) :: depth_filling    ! Depth where pore filling begins [m]
     316integer,                      intent(out)   :: index_filling    ! Index where the pore filling begins [1]
     317integer,                      intent(out)   :: index_geothermal ! Index where the ice table stops [1]
     318real,                         intent(out)   :: depth_geothermal ! Depth where the ice table stops [m]
     319real, dimension(nsoilmx_PEM), intent(out)   :: dz_etadz_rho     ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
     320
     321! LOCAL VARIABLES
     322! ---------------
     323real, dimension(nsoilmx_PEM) :: eta                          ! Constriction
     324real, dimension(nsoilmx_PEM) :: dz_psat                      ! First derivative of the vapor pressure at saturation
    235325real, dimension(nsoilmx_PEM) :: dz_eta                       ! \partial z \eta
    236326real, dimension(nsoilmx_PEM) :: dzz_psat                     ! \partial \partial psat
    237 integer                      :: ilay, index_tmp              ! index for loop
    238 real                         :: old_depth_filling            ! depth_filling saved
    239 real                         :: Jdry                         ! flux trought the dry layer
    240 real                         :: Jsat                         ! flux trought the ice layer
    241 real                         :: Jdry_prevlay, Jsat_prevlay   ! same but for the previous ice layer
    242 integer                      :: index_firstice               ! first index where ice appears (i.e., f > 0)
    243 real                         :: massfillabove, massfillafter ! h2O mass above and after index_geothermal
    244 ! Constant
    245 !---------
    246 real, parameter :: pvap2rho = 18.e-3/8.314
    247 ! Code
    248 !-----
     327integer                      :: ilay, index_tmp              ! Index for loop
     328real                         :: old_depth_filling            ! Depth_filling saved
     329real                         :: Jdry                         ! Flux trought the dry layer
     330real                         :: Jsat                         ! Flux trought the ice layer
     331real                         :: Jdry_prevlay, Jsat_prevlay   ! Same but for the previous ice layer
     332integer                      :: index_firstice               ! First index where ice appears (i.e., f > 0)
     333real                         :: massfillabove, massfillafter ! H2O mass above and after index_geothermal
     334real, parameter              :: pvap2rho = 18.e-3/8.314
     335
     336! CODE
     337! ----
    249338! 0. Compute constriction over the layer
    250339! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
    … …  
    342431
    343432END SUBROUTINE find_layering_icetable
    344 
    345 !-----------------------------------------------------------------------
     433!=======================================================================
     434
     435!=======================================================================
    346436FUNCTION constriction(porefill) RESULT(eta)
    347 
    348 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    349 !!!
    350 !!! Purpose: Compute the constriction of vapor flux by pore ice
    351 !!!
    352 !!! Author: LL
    353 !!!
    354 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    355 
    356 implicit none
    357 
     437!-----------------------------------------------------------------------
     438! NAME
     439!     constriction
     440!
     441! DESCRIPTION
     442!     Compute the constriction of vapor flux by pore ice.
     443!
     444! AUTHORS & DATE
     445!     L. Lange
     446!     JB Clement, 2023-2025
     447!
     448! NOTES
     449!
     450!-----------------------------------------------------------------------
     451
     452! DECLARATION
     453! -----------
     454implicit none
     455
     456! ARGUMENTS
     457! ---------
    358458real, intent(in) :: porefill ! pore filling fraction
     459
     460! LOCAL VARIABLES
     461! ---------------
    359462real :: eta ! constriction
    360463
     464! CODE
     465! ----
    361466if (porefill <= 0.) then
    362467    eta = 1.
    … …  
    368473
    369474END FUNCTION constriction
    370 
    371 !-----------------------------------------------------------------------
     475!=======================================================================
     476
     477!=======================================================================
    372478SUBROUTINE compute_Tice_pem(nsoil,ptsoil,ptsurf,ice_depth,Tice)
    373 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    374 !!!
    375 !!! Purpose: Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells.
    376 !!!
    377 !!! Author: LL
    378 !!!
    379 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    380 
    381 use soil,      only: mlayer_PEM
     479!-----------------------------------------------------------------------
     480! NAME
     481!     compute_Tice_pem
     482!
     483! DESCRIPTION
     484!     Compute subsurface ice temperature by interpolating the
     485!     temperatures between the two adjacent cells.
     486!
     487! AUTHORS & DATE
     488!     L. Lange
     489!     JB Clement, 12/12/2025
     490!
     491! NOTES
     492!
     493!-----------------------------------------------------------------------
     494
     495! DEPENDENCIES
     496! ------------
     497use soil,     only: mlayer_PEM
    382498use stop_pem, only: stop_clean
    383499
    384 implicit none
    385 
    386 ! Inputs
    387 ! ------
    388 integer,                intent(in) :: nsoil     ! Number of soil layers
    389 real, dimension(nsoil), intent(in) :: ptsoil    ! Soil temperature (K)
    390 real,                   intent(in) :: ptsurf    ! Surface temperature (K)
    391 real,                   intent(in) :: ice_depth ! Ice depth (m)
    392 
    393 ! Outputs
    394 ! -------
    395 real, intent(out) :: Tice ! Ice temperatures (K)
    396 
    397 ! Local variables
     500! DECLARATION
     501! -----------
     502implicit none
     503
     504! ARGUMENTS
     505! ---------
     506integer,                intent(in)  :: nsoil     ! Number of soil layers
     507real, dimension(nsoil), intent(in)  :: ptsoil    ! Soil temperature (K)
     508real,                   intent(in)  :: ptsurf    ! Surface temperature (K)
     509real,                   intent(in)  :: ice_depth ! Ice depth (m)
     510real,                   intent(out) :: Tice      ! Ice temperatures (K)
     511
     512! LOCAL VARIABLES
    398513! ---------------
    399514integer :: ik       ! Loop variables
    400515integer :: indexice ! Index of the ice
    401516
    402 ! Code
    403 !-----
     517! CODE
     518! ----
    404519indexice = -1
    405520if (ice_depth >= mlayer_PEM(nsoil - 1)) then
    … …  
    428543
    429544END SUBROUTINE compute_Tice_pem
    430 
    431 !-----------------------------------------------------------------------
     545!=======================================================================
     546
     547!=======================================================================
    432548FUNCTION rho_ice(T,ice) RESULT(rho)
    433 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    434 !!!
    435 !!! Purpose: Compute subsurface ice density
    436 !!!
    437 !!! Author: JBC
    438 !!!
    439 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    440 
     549!-----------------------------------------------------------------------
     550! NAME
     551!     rho_ice
     552!
     553! DESCRIPTION
     554!     Compute subsurface ice density.
     555!
     556! AUTHORS & DATE
     557!     JB Clement, 12/12/2025
     558!
     559! NOTES
     560!
     561!-----------------------------------------------------------------------
     562
     563! DEPENDENCIES
     564! ------------
    441565use stop_pem, only: stop_clean
    442566
    443 implicit none
    444 
     567! DECLARATION
     568! -----------
     569implicit none
     570
     571! ARGUMENTS
     572! ---------
    445573real,         intent(in) :: T
    446574character(3), intent(in) :: ice
     575
     576! LOCAL VARIABLES
     577! ---------------
    447578real :: rho
    448579
     580! CODE
     581! ----
    449582select case (trim(adjustl(ice)))
    450583    case('h2o')
    … …  
    457590
    458591END FUNCTION rho_ice
     592!=======================================================================
    459593
    460594END MODULE ice_table

  • trunk/LMDZ.COMMON/libf/evolution/info_pem.F90

    r3989 r3991  
    11MODULE info_pem
     2!-----------------------------------------------------------------------
     3! NAME
     4!     info_pem
     5!
     6! DESCRIPTION
     7!     Handles counters for PCM/PEM coupled runs and simulation metadata.
     8!
     9! AUTHORS & DATE
     10!     R. Vandemeulebrouck
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
    216
     17! DECLARATION
     18! -----------
    319implicit none
    420
     21! MODULE VARIABLES
     22! ----------------
    523integer :: iPCM, iPEM, nPCM, nPCM_ini ! Data about the chained simulation of PCM/PEM runs
    624
    7 !=======================================================================
    825contains
    9 !=======================================================================
    10 
    11 SUBROUTINE update_info(i_myear_leg,stopPEM,i_myear,nyears_tot)
     26!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    1227
    1328!=======================================================================
     29SUBROUTINE update_info(i_myear_leg,stopPEM,i_myear,nyears_tot)
     30!-----------------------------------------------------------------------
     31! NAME
     32!    update_info
    1433!
    15 ! Purpose: Update the first line of "launchPEM.info" to count the number of simulated Martian years
    16 !          Write in "launchPEM.info" the reason why the PEM stopped and the number of simulated years
     34! DESCRIPTION
     35!    Update the first line of "launchPEM.info" to count the number of
     36!    simulated Martian years. Write in "launchPEM.info" the reason why
     37!    the PEM stopped and the number of simulated years.
    1738!
    18 ! Author: RV, JBC
    19 !=======================================================================
     39! AUTHORS & DATE
     40!    R. Vandemeulebrouck
     41!    JB Clement, 2023-2025
     42!
     43! NOTES
     44!
     45!-----------------------------------------------------------------------
    2046
     47! DEPENDENCIES
     48! ------------
    2149use evolution, only: convert_years, year_bp_ini
    2250
     51! DECLARATION
     52! -----------
    2353implicit none
    2454
    25 !----- Arguments
    26 integer, intent(in) :: stopPEM          ! Reason to stop
    27 real,    intent(in) :: i_myear_leg      ! # of years
    28 real,    intent(in) :: i_myear, nyears_tot ! Current simulated Martian year and maximum number of Martian years to be simulated
     55! ARGUMENTS
     56! ---------
     57integer, intent(in) :: stopPEM     ! Reason to stop
     58real,    intent(in) :: i_myear_leg ! # of years
     59real,    intent(in) :: i_myear     ! Current simulated Martian year
     60real,    intent(in) :: nyears_tot  ! Maximum number of Martian years to be simulated
    2961
    30 !----- Local variables
     62! LOCAL VARIABLES
     63! ---------------
    3164logical       :: ok
    3265integer       :: cstat
    3366character(20) :: fch1, fch2, fch3
    3467
    35 !----- Code
     68! CODE
     69! ----
    3670inquire(file = 'launchPEM.info',exist = ok)
    3771if (ok) then
    … …  
    5589
    5690END SUBROUTINE update_info
     91!=======================================================================
    5792
    5893!=======================================================================
     94FUNCTION int2str(i) RESULT(str)
     95!-----------------------------------------------------------------------
     96! NAME
     97!     int2str
     98!
     99! DESCRIPTION
     100!     Convert an integer into a string.
     101!
     102! AUTHORS & DATE
     103!     JB Clement, 2023-2025
     104!
     105! NOTES
     106!
     107!-----------------------------------------------------------------------
    59108
    60 FUNCTION int2str(i) RESULT(str)
    61 ! Function to convert an integer into a string
     109! DECLARATION
     110! -----------
     111implicit none
    62112
    63 integer, intent(in)       :: i
     113! ARGUMENTS
     114! ---------
     115integer, intent(in) :: i
     116
     117! LOCAL VARIABLES
     118! ---------------
    64119character(20)             :: str_tmp
    65120character(:), allocatable :: str
    66121
     122! CODE
     123! ----
    67124if (nb_digits(real(i)) > len(str_tmp)) error stop 'int2str [info_pem]: invalid integer for conversion!'
    68125write(str_tmp,'(i0)') i
    … …  
    70127
    71128END FUNCTION int2str
     129!=======================================================================
    72130
    73131!=======================================================================
     132FUNCTION nb_digits(x) RESULT(idigits)
     133!-----------------------------------------------------------------------
     134! NAME
     135!     nb_digits
     136!
     137! DESCRIPTION
     138!     Give the number of digits for the integer part of a real number.
     139!
     140! AUTHORS & DATE
     141!     JB Clement, 2023-2025
     142!
     143! NOTES
     144!
     145!-----------------------------------------------------------------------
    74146
    75 FUNCTION nb_digits(x) RESULT(idigits)
    76 ! Function to give the number of digits for the integer part of a real number
     147! DECLARATION
     148! -----------
     149implicit none
    77150
     151! ARGUMENTS
     152! ---------
    78153real, intent(in) :: x
    79 integer          :: idigits
    80154
     155! LOCAL VARIABLES
     156! ---------------
     157integer :: idigits
     158
     159! CODE
     160! ----
    81161idigits = 1
    82162! If x /= 0 then:
    … …  
    84164
    85165END FUNCTION nb_digits
     166!=======================================================================
    86167
    87168END MODULE info_pem

  • trunk/LMDZ.COMMON/libf/evolution/iostart_pem.F90

    r3989 r3991  
    11MODULE iostart_pem
    2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    3 !!!
    4 !!! Purpose: Read and write specific netcdf for the PEM
    5 !!!
    6 !!!
    7 !!! Author: LL, inspired by iostart from the PCM
    8 !!!
    9 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    10 
    11     IMPLICIT NONE
     2!-----------------------------------------------------------------------
     3! NAME
     4!     iostart_pem
     5!
     6! DESCRIPTION
     7!     Read and write specific netcdf files for the PEM (start and restart)
     8!     Provides interfaces for field and variable I/O operations
     9!
     10! AUTHORS & DATE
     11!     L. Lange
     12!     JB Clement, 2023-2025
     13!
     14! NOTES
     15!     Inspired by iostart from the PCM.
     16!
     17!-----------------------------------------------------------------------
     18
     19    ! DECLARATION
     20! -----------
     21implicit none
     22
     23! MODULE VARIABLES
     24! ----------------
    1225PRIVATE
    1326    INTEGER,SAVE :: nid_start ! NetCDF file identifier for startfi.nc file
    … …  
    5164
    5265CONTAINS
    53 
     66!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     67
     68!=======================================================================
    5469  SUBROUTINE open_startphy(filename)
    55   USE netcdf, only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, nf90_strerror
     70!-----------------------------------------------------------------------
     71! NAME
     72!     open_startphy
     73!
     74! DESCRIPTION
     75!     Open the startphy netCDF file for reading.
     76
     77! AUTHORS & DATE
     78!     L. Lange
     79!     JB Clement, 2023-2025
     80
     81! NOTES
     82!
     83!-----------------------------------------------------------------------
     84
     85! DEPENDENCIES
     86! ------------
     87  USE netcdf,             only: NF90_OPEN, NF90_NOERR, NF90_NOWRITE, nf90_strerror
    5688  USE mod_phys_lmdz_para, only: is_master, bcast
    57   IMPLICIT NONE
     89
     90! DECLARATION
     91! -----------
     92implicit none
     93
     94! ARGUMENTS
     95! ---------
    5896    CHARACTER(LEN=*) :: filename
    59     INTEGER          :: ierr
    60 
     97
     98! LOCAL VARIABLES
     99! ---------------
     100    INTEGER :: ierr
     101
     102! CODE
     103! ----
    61104    IF (is_master) THEN
    62105      ierr = NF90_OPEN (filename, NF90_NOWRITE, nid_start)
    … …  
    71114
    72115  END SUBROUTINE open_startphy
    73 
     116!=======================================================================
     117
     118!=======================================================================
    74119  SUBROUTINE close_startphy
    75   USE netcdf, only: NF90_CLOSE
     120!-----------------------------------------------------------------------
     121! NAME
     122!     close_startphy
     123!
     124! DESCRIPTION
     125!     Close the startphy netCDF file.
     126
     127! AUTHORS & DATE
     128!     L. Lange
     129!     JB Clement, 2023-2025
     130
     131! NOTES
     132!
     133!-----------------------------------------------------------------------
     134
     135! DEPENDENCIES
     136! ------------
     137  USE netcdf,             only: NF90_CLOSE
    76138  USE mod_phys_lmdz_para, only: is_master
    77   IMPLICIT NONE
    78     INTEGER          :: ierr
    79 
     139
     140! DECLARATION
     141! -----------
     142implicit none
     143
     144! LOCAL VARIABLES
     145! ---------------
     146    INTEGER :: ierr
     147
     148! CODE
     149! ----
    80150    IF (is_master) THEN
    81151        ierr = NF90_CLOSE (nid_start)
    … …  
    83153
    84154  END SUBROUTINE close_startphy
    85 
    86 
     155!=======================================================================
     156
     157!=======================================================================
    87158  FUNCTION inquire_field(Field_name)
    88   ! check if a given field is present in the input file
    89   USE netcdf, only: NF90_INQ_VARID, NF90_NOERR
     159!-----------------------------------------------------------------------
     160! NAME
     161!     inquire_field
     162!
     163! DESCRIPTION
     164!     Check if a given field is present in the input file.
     165
     166! AUTHORS & DATE
     167!     L. Lange
     168!     JB Clement, 2023-2025
     169
     170! NOTES
     171!
     172!-----------------------------------------------------------------------
     173
     174! DEPENDENCIES
     175! ------------
     176  USE netcdf,             only: NF90_INQ_VARID, NF90_NOERR
    90177  USE mod_phys_lmdz_para, only: is_master, bcast
    91   IMPLICIT NONE
    92     CHARACTER(LEN=*),INTENT(IN) :: Field_name
     178
     179! DECLARATION
     180! -----------
     181implicit none
     182
     183! ARGUMENTS
     184! ---------
     185    CHARACTER(LEN=*), INTENT(IN) :: Field_name
     186
     187! LOCAL VARIABLES
     188! ---------------
    93189    LOGICAL :: inquire_field
    94190    INTEGER :: varid
    95191    INTEGER :: ierr
    96192
     193! CODE
     194! ----
    97195    IF (is_master) THEN
    98196      ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
    … …  
    107205
    108206  END FUNCTION inquire_field
    109 
    110 
     207!=======================================================================
     208
     209!=======================================================================
    111210  FUNCTION inquire_field_ndims(Field_name)
    112   ! give the number of dimensions of "Field_name" stored in the input file
    113   USE netcdf, only: nf90_inq_varid, nf90_inquire_variable, &
    114                     NF90_NOERR, nf90_strerror
     211!-----------------------------------------------------------------------
     212! NAME
     213!     inquire_field_ndims
     214!
     215! DESCRIPTION
     216!     Give the number of dimensions of "Field_name" stored in the input file.
     217
     218! AUTHORS & DATE
     219!     L. Lange
     220!     JB Clement, 2023-2025
     221
     222! NOTES
     223!
     224!-----------------------------------------------------------------------
     225
     226! DEPENDENCIES
     227! ------------
     228  USE netcdf,             only: nf90_inq_varid, nf90_inquire_variable, &
     229                                NF90_NOERR, nf90_strerror
    115230  USE mod_phys_lmdz_para, only: is_master, bcast
    116   IMPLICIT NONE
    117     CHARACTER(LEN=*),INTENT(IN) :: Field_name
     231
     232! DECLARATION
     233! -----------
     234implicit none
     235
     236! ARGUMENTS
     237! ---------
     238    CHARACTER(LEN=*), INTENT(IN) :: Field_name
     239
     240! LOCAL VARIABLES
     241! ---------------
    118242    INTEGER :: inquire_field_ndims
    119243    INTEGER :: varid
    120244    INTEGER :: ierr
    121245
     246! CODE
     247! ----
    122248    IF (is_master) THEN
    123249      ierr=nf90_inq_varid(nid_start,Field_name,varid)
    … …  
    135261
    136262  END FUNCTION inquire_field_ndims
    137 
    138 
     263!=======================================================================
     264
     265!=======================================================================
    139266  FUNCTION inquire_dimension(Field_name)
    140   ! check if a given dimension is present in the input file
     267!-----------------------------------------------------------------------
     268! NAME
     269!     inquire_dimension
     270!
     271! DESCRIPTION
     272!     Check if a given dimension is present in the input file.
     273
     274! AUTHORS & DATE
     275!     L. Lange
     276!     JB Clement, 2023-2025
     277
     278! NOTES
     279!
     280!-----------------------------------------------------------------------
     281
     282! DEPENDENCIES
     283! ------------
    141284  USE netcdf, only: nf90_inq_dimid, NF90_NOERR
    142285  USE mod_phys_lmdz_para, only: is_master, bcast
    143   IMPLICIT NONE
    144     CHARACTER(LEN=*),INTENT(IN) :: Field_name
     286
     287! DECLARATION
     288! -----------
     289implicit none
     290
     291! ARGUMENTS
     292! ---------
     293    CHARACTER(LEN=*), INTENT(IN) :: Field_name
     294
     295! LOCAL VARIABLES
     296! ---------------
    145297    LOGICAL :: inquire_dimension
    146298    INTEGER :: varid
    147299    INTEGER :: ierr
    148300
     301! CODE
     302! ----
    149303    IF (is_master) THEN
    150304      ierr=NF90_INQ_DIMID(nid_start,Field_name,varid)
    … …  
    159313
    160314  END FUNCTION inquire_dimension
    161 
     315!=======================================================================
     316
     317!=======================================================================
    162318  FUNCTION inquire_dimension_length(Field_name)
    163   ! give the length of the "Field_name" dimension stored in the input file
     319!-----------------------------------------------------------------------
     320! NAME
     321!     inquire_dimension_length
     322!
     323! DESCRIPTION
     324!     Give the length of the "Field_name" dimension stored in the input file.
     325
     326! AUTHORS & DATE
     327!     L. Lange
     328!     JB Clement, 2023-2025
     329
     330! NOTES
     331!
     332!-----------------------------------------------------------------------
     333
     334! DEPENDENCIES
     335! ------------
    164336  USE netcdf, only: nf90_inquire_dimension, nf90_inq_dimid, &
    165337                    NF90_NOERR, nf90_strerror
    166338  USE mod_phys_lmdz_para, only: is_master, bcast
    167   IMPLICIT NONE
    168     CHARACTER(LEN=*),INTENT(IN) :: Field_name
     339
     340! DECLARATION
     341! -----------
     342implicit none
     343
     344! ARGUMENTS
     345! ---------
     346    CHARACTER(LEN=*), INTENT(IN) :: Field_name
     347
     348! LOCAL VARIABLES
     349! ---------------
    169350    INTEGER :: inquire_dimension_length
    170351    INTEGER :: varid
    171352    INTEGER :: ierr
    172353
     354! CODE
     355! ----
    173356    IF (is_master) THEN
    174357      ierr=nf90_inq_dimid(nid_start,Field_name,varid)
    … …  
    186369
    187370  END FUNCTION inquire_dimension_length
    188 
    189 
    190 
     371!=======================================================================
     372
     373!=======================================================================
    191374  SUBROUTINE Get_Field_r1(field_name,field,found,timeindex)
    192   ! For a surface field
    193   use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
    194   IMPLICIT NONE
    195     CHARACTER(LEN=*),INTENT(IN)    :: Field_name
    196     REAL,INTENT(INOUT)               :: Field(:)
    197     LOGICAL,INTENT(OUT),OPTIONAL   :: found
    198     INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
    199 
     375!-----------------------------------------------------------------------
     376! NAME
     377!     Get_Field_r1
     378!
     379! DESCRIPTION
     380!     Read a surface field (1D) from the start file.
     381
     382! AUTHORS & DATE
     383!     L. Lange
     384!     JB Clement, 2023-2025
     385
     386! NOTES
     387!
     388!-----------------------------------------------------------------------
     389
     390! DEPENDENCIES
     391! ------------
     392  use mod_grid_phy_lmdz, only: klon_glo
     393
     394! DECLARATION
     395! -----------
     396implicit none
     397
     398! ARGUMENTS
     399! ---------
     400  CHARACTER(LEN=*), INTENT(IN)            :: field_name
     401  REAL,             INTENT(INOUT)         :: field(:)
     402  LOGICAL,          INTENT(OUT), OPTIONAL :: found
     403  INTEGER,          INTENT(IN),  OPTIONAL :: timeindex
     404
     405! LOCAL VARIABLES
     406! ---------------
    200407    integer :: edges(4), corners(4)
    201408
     409! CODE
     410! ----
    202411    edges(1)=klon_glo
    203412    edges(2:4)=1
    … …  
    214423
    215424  END SUBROUTINE Get_Field_r1
    216 
     425!=======================================================================
     426
     427!=======================================================================
    217428  SUBROUTINE Get_Field_r2(field_name,field,found,timeindex)
    218   ! For a "3D" horizontal-vertical field
    219   use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
    220   IMPLICIT NONE
    221     CHARACTER(LEN=*),INTENT(IN)    :: Field_name
    222     REAL,INTENT(INOUT)               :: Field(:,:)
    223     LOGICAL,INTENT(OUT),OPTIONAL   :: found
    224     INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
    225 
     429!-----------------------------------------------------------------------
     430! NAME
     431!     Get_Field_r2
     432!
     433! DESCRIPTION
     434!     Read a "3D" horizontal-vertical field (2D) from the start file.
     435
     436! AUTHORS & DATE
     437!     L. Lange
     438!     JB Clement, 2023-2025
     439
     440! NOTES
     441!
     442!-----------------------------------------------------------------------
     443
     444! DEPENDENCIES
     445! ------------
     446  use mod_grid_phy_lmdz, only: klon_glo
     447
     448! DECLARATION
     449! -----------
     450implicit none
     451
     452! ARGUMENTS
     453! ---------
     454  CHARACTER(LEN=*), INTENT(IN)            :: field_name
     455  REAL,             INTENT(INOUT)         :: field(:,:)
     456  LOGICAL,          INTENT(OUT), OPTIONAL :: found
     457  INTEGER,          INTENT(IN),  OPTIONAL :: timeindex
     458
     459! LOCAL VARIABLES
     460! ---------------
    226461    integer :: edges(4), corners(4)
    227462
     463! CODE
     464! ----
    228465    edges(1)=klon_glo
    229466    edges(2)=size(field,2)
    … …  
    244481
    245482  END SUBROUTINE Get_Field_r2
    246 
     483!=======================================================================
     484
     485!=======================================================================
    247486  SUBROUTINE Get_Field_r3(field_name,field,found,timeindex)
    248   ! for a "4D" field surf/alt/??
    249   use mod_grid_phy_lmdz, only: klon_glo ! number of atmospheric columns (full grid)
    250   IMPLICIT NONE
    251     CHARACTER(LEN=*),INTENT(IN)    :: Field_name
    252     REAL,INTENT(INOUT)               :: Field(:,:,:)
    253     LOGICAL,INTENT(OUT),OPTIONAL   :: found
    254     INTEGER,INTENT(IN),OPTIONAL    :: timeindex ! time index of sought data
    255 
     487!-----------------------------------------------------------------------
     488! NAME
     489!     Get_Field_r3
     490!
     491! DESCRIPTION
     492!     Read a "4D" field surf/alt/?? (3D) from the start file.
     493
     494! AUTHORS & DATE
     495!     L. Lange
     496!     JB Clement, 2023-2025
     497
     498! NOTES
     499!
     500!-----------------------------------------------------------------------
     501
     502! DEPENDENCIES
     503! ------------
     504  use mod_grid_phy_lmdz, only: klon_glo
     505
     506! DECLARATION
     507! -----------
     508implicit none
     509
     510! ARGUMENTS
     511! ---------
     512  CHARACTER(LEN=*), INTENT(IN)            :: field_name
     513  REAL,             INTENT(INOUT)         :: field(:,:,:)
     514  LOGICAL,          INTENT(OUT), OPTIONAL :: found
     515  INTEGER,          INTENT(IN),  OPTIONAL :: timeindex
     516
     517! LOCAL VARIABLES
     518! ---------------
    256519    integer :: edges(4), corners(4)
    257520
     521! CODE
     522! ----
    258523    edges(1)=klon_glo
    259524    edges(2)=size(field,2)
    … …  
    274539
    275540  END SUBROUTINE Get_Field_r3
    276 
     541!=======================================================================
     542
     543!=======================================================================
    277544  SUBROUTINE Get_field_rgen(field_name,field,field_size, &
    278545                            corners,edges,found)
     546!-----------------------------------------------------------------------
     547! NAME
     548!     Get_field_rgen
     549!
     550! DESCRIPTION
     551!     Generic subroutine to read fields from start file with scatter to processes.
     552
     553! AUTHORS & DATE
     554!     L. Lange
     555!     JB Clement, 2023-2025
     556
     557! NOTES
     558!
     559!-----------------------------------------------------------------------
     560
     561! DEPENDENCIES
     562! ------------
    279563  USE netcdf
    280564  USE dimphy
    281565  USE mod_grid_phy_lmdz
    282566  USE mod_phys_lmdz_para
    283   IMPLICIT NONE
    284     CHARACTER(LEN=*) :: Field_name
    285     INTEGER          :: field_size
    286     REAL             :: field(klon,field_size)
    287     INTEGER,INTENT(IN) :: corners(4)
    288     INTEGER,INTENT(IN) :: edges(4)
    289     LOGICAL,OPTIONAL :: found
    290 
     567
     568! DECLARATION
     569! -----------
     570implicit none
     571
     572! ARGUMENTS
     573! ---------
     574    CHARACTER(LEN=*), INTENT(IN)            :: field_name
     575    INTEGER,          INTENT(IN)            :: field_size
     576    REAL,             INTENT(OUT)           :: field(klon,field_size)
     577    INTEGER,          INTENT(IN)            :: corners(4)
     578    INTEGER,          INTENT(IN)            :: edges(4)
     579    LOGICAL,          INTENT(OUT), OPTIONAL :: found
     580
     581! LOCAL VARIABLES
     582! ---------------
    291583    REAL    :: field_glo(klon_glo,field_size)
    292584    LOGICAL :: tmp_found
    … …  
    294586    INTEGER :: ierr
    295587
     588! CODE
     589! ----
    296590    IF (is_master) THEN
    297591
    298       ierr=NF90_INQ_VARID(nid_start,Field_name,varid)
     592      ierr=NF90_INQ_VARID(nid_start,field_name,varid)
    299593
    300594      IF (ierr==NF90_NOERR) THEN
    … …  
    323617
    324618
    325     CONTAINS
    326 
    327      SUBROUTINE body(field_glo)
    328        REAL :: field_glo(klon_glo*field_size)
    329          ierr=NF90_GET_VAR(nid_start,varid,field_glo,corners,edges)
    330          IF (ierr/=NF90_NOERR) THEN
    331            ! La variable exist dans le fichier mais la lecture a echouee.
    332            write(*,*) 'get_field_rgen: Failed reading <'//field_name//'>'
     619        CONTAINS
     620    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     621
     622    !=======================================================================
     623      SUBROUTINE body(field_glo)
     624    !-----------------------------------------------------------------------
     625    ! NAME
     626    !     body
     627    !
     628    ! DESCRIPTION
     629    !     Read the requested field from the start file into the global array.
     630    !
     631    ! AUTHORS & DATE
     632    !     L. Lange
     633    !     JB Clement, 2023-2025
     634    !
     635    ! NOTES
     636    !     Helper contained in Get_field_rgen.
     637    !-----------------------------------------------------------------------
     638
     639    ! DECLARATION
     640    ! -----------
     641      ! DECLARATION
     642! -----------
     643implicit none
     644
     645    ! ARGUMENTS
     646    ! ---------
     647      REAL, INTENT(OUT) :: field_glo(klon_glo*field_size) ! Buffer for global field
     648
     649    ! CODE
     650    ! ----
     651          ierr=NF90_GET_VAR(nid_start,varid,field_glo,corners,edges)
     652          IF (ierr/=NF90_NOERR) THEN
     653         ! La variable exist dans le fichier mais la lecture a echouee.
     654         write(*,*) 'get_field_rgen: Failed reading <'//field_name//'>'
    333655
    334656!           IF (field_name=='CLWCON' .OR. field_name=='RNEBCON' .OR. field_name=='RATQS') THEN
    … …  
    350672
    351673  END SUBROUTINE Get_field_rgen
    352 
    353 
     674!=======================================================================
     675
     676!=======================================================================
    354677  SUBROUTINE get_var_r0(var_name,var,found)
    355   ! Get a scalar from input file
    356   IMPLICIT NONE
    357     CHARACTER(LEN=*),INTENT(IN)  :: var_name
    358     REAL,INTENT(INOUT)             :: var
    359     LOGICAL,OPTIONAL,INTENT(OUT) :: found
    360 
    361     REAL                         :: varout(1)
    362 
     678!-----------------------------------------------------------------------
     679! NAME
     680!     get_var_r0
     681!
     682! DESCRIPTION
     683!     Get a scalar variable from input file.
     684!
     685! AUTHORS & DATE
     686!     L. Lange
     687!     JB Clement, 2023-2025
     688!
     689! NOTES
     690!
     691!-----------------------------------------------------------------------
     692
     693! DECLARATION
     694! -----------
     695implicit none
     696
     697! ARGUMENTS
     698! ---------
     699  CHARACTER(LEN=*), INTENT(IN)            :: var_name
     700  REAL,             INTENT(INOUT)         :: var
     701  LOGICAL,          INTENT(OUT), OPTIONAL :: found
     702
     703! LOCAL VARIABLES
     704! ---------------
     705    REAL :: varout(1)
     706
     707! CODE
     708! ----
    363709    IF (PRESENT(found)) THEN
    364710      CALL Get_var_rgen(var_name,varout,size(varout),found)
    … …  
    369715
    370716  END SUBROUTINE get_var_r0
    371 
     717!=======================================================================
     718
     719!=======================================================================
    372720  SUBROUTINE get_var_r1(var_name,var,found)
    373   ! Get a vector from input file
    374   IMPLICIT NONE
    375     CHARACTER(LEN=*),INTENT(IN)  :: var_name
    376     REAL,INTENT(INOUT)             :: var(:)
    377     LOGICAL,OPTIONAL,INTENT(OUT) :: found
    378 
     721!-----------------------------------------------------------------------
     722! NAME
     723!     get_var_r1
     724!
     725! DESCRIPTION
     726!     Get a vector (1D) variable from input file.
     727!
     728! AUTHORS & DATE
     729!     L. Lange
     730!     JB Clement, 2023-2025
     731!
     732! NOTES
     733!
     734!-----------------------------------------------------------------------
     735
     736! DECLARATION
     737! -----------
     738implicit none
     739
     740! ARGUMENTS
     741! ---------
     742  CHARACTER(LEN=*), INTENT(IN)            :: var_name
     743  REAL,             INTENT(INOUT)         :: var(:)
     744  LOGICAL,          INTENT(OUT), OPTIONAL :: found
     745
     746! CODE
     747! ----
    379748    IF (PRESENT(found)) THEN
    380749      CALL Get_var_rgen(var_name,var,size(var),found)
    … …  
    384753
    385754  END SUBROUTINE get_var_r1
    386 
     755!=======================================================================
     756
     757!=======================================================================
    387758  SUBROUTINE get_var_r2(var_name,var,found)
    388   ! Get a 2D field from input file
    389   IMPLICIT NONE
    390     CHARACTER(LEN=*),INTENT(IN)  :: var_name
    391     REAL,INTENT(OUT)             :: var(:,:)
    392     LOGICAL,OPTIONAL,INTENT(OUT) :: found
    393 
     759!-----------------------------------------------------------------------
     760! NAME
     761!     get_var_r2
     762!
     763! DESCRIPTION
     764!     Get a 2D field variable from input file.
     765!
     766! AUTHORS & DATE
     767!     L. Lange
     768!     JB Clement, 2023-2025
     769!
     770! NOTES
     771!
     772!-----------------------------------------------------------------------
     773
     774! DECLARATION
     775! -----------
     776implicit none
     777
     778! ARGUMENTS
     779! ---------
     780  CHARACTER(LEN=*), INTENT(IN)            :: var_name
     781  REAL,             INTENT(OUT)           :: var(:,:)
     782  LOGICAL,          INTENT(OUT), OPTIONAL :: found
     783
     784! CODE
     785! ----
    394786    IF (PRESENT(found)) THEN
    395787      CALL Get_var_rgen(var_name,var,size(var),found)
    … …  
    399791
    400792  END SUBROUTINE get_var_r2
    401 
     793!=======================================================================
     794
     795!=======================================================================
    402796  SUBROUTINE get_var_r3(var_name,var,found)
    403   ! Get a 3D field frominput file
    404   IMPLICIT NONE
    405     CHARACTER(LEN=*),INTENT(IN)  :: var_name
    406     REAL,INTENT(INOUT)             :: var(:,:,:)
    407     LOGICAL,OPTIONAL,INTENT(OUT) :: found
    408 
     797!-----------------------------------------------------------------------
     798! NAME
     799!     get_var_r3
     800!
     801! DESCRIPTION
     802!     Get a 3D field variable from input file.
     803!
     804! AUTHORS & DATE
     805!     L. Lange
     806!     JB Clement, 2023-2025
     807!
     808! NOTES
     809!
     810!-----------------------------------------------------------------------
     811
     812! DECLARATION
     813! -----------
     814implicit none
     815
     816! ARGUMENTS
     817! ---------
     818  CHARACTER(LEN=*), INTENT(IN)            :: var_name
     819  REAL,             INTENT(INOUT)         :: var(:,:,:)
     820  LOGICAL,          INTENT(OUT), OPTIONAL :: found
     821
     822! CODE
     823! ----
    409824    IF (PRESENT(found)) THEN
    410825      CALL Get_var_rgen(var_name,var,size(var),found)
    … …  
    414829
    415830  END SUBROUTINE get_var_r3
    416 
     831!=======================================================================
     832
     833!=======================================================================
    417834  SUBROUTINE Get_var_rgen(var_name,var,var_size,found)
     835!-----------------------------------------------------------------------
     836! NAME
     837!     Get_var_rgen
     838!
     839! DESCRIPTION
     840!     Generic subroutine to read variables from input file.
     841!
     842! AUTHORS & DATE
     843!     L. Lange
     844!     JB Clement, 2023-2025
     845!
     846! NOTES
     847!
     848!-----------------------------------------------------------------------
     849
     850! DEPENDENCIES
     851! ------------
    418852  USE netcdf
    419853  USE dimphy
    420854  USE mod_grid_phy_lmdz
    421855  USE mod_phys_lmdz_para
    422   IMPLICIT NONE
    423     CHARACTER(LEN=*) :: var_name
    424     INTEGER          :: var_size
    425     REAL             :: var(var_size)
    426     LOGICAL,OPTIONAL :: found
    427 
     856
     857! DECLARATION
     858! -----------
     859implicit none
     860
     861! ARGUMENTS
     862! ---------
     863    CHARACTER(LEN=*), INTENT(IN)            :: var_name
     864    INTEGER,          INTENT(IN)            :: var_size
     865    REAL,             INTENT(OUT)           :: var(var_size)
     866    LOGICAL,          INTENT(OUT), OPTIONAL :: found
     867
     868! LOCAL VARIABLES
     869! ---------------
    428870    LOGICAL :: tmp_found
    429871    INTEGER :: varid
    430872    INTEGER :: ierr
    431873
     874! CODE
     875! ----
    432876    IF (is_mpi_root .AND. is_omp_root) THEN
    433877
    … …  
    463907
    464908  END SUBROUTINE Get_var_rgen
    465 
    466 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    467 
     909!=======================================================================
     910
     911!=======================================================================
    468912  SUBROUTINE open_restartphy(filename)
     913!-----------------------------------------------------------------------
     914! NAME
     915!     open_restartphy
     916!
     917! DESCRIPTION
     918!     Open (or create) the restart netCDF file for writing.
     919!     Define dimensions and global attributes on first call.
     920!
     921! AUTHORS & DATE
     922!     L. Lange
     923!     JB Clement, 2023-2025
     924!
     925! NOTES
     926!
     927!-----------------------------------------------------------------------
     928
     929! DEPENDENCIES
     930! ------------
    469931  USE netcdf, only: NF90_CREATE, NF90_CLOBBER, NF90_64BIT_OFFSET, &
    470932                    NF90_NOERR, nf90_strerror, &
    … …  
    473935                    NF90_WRITE, NF90_OPEN
    474936  USE mod_phys_lmdz_para, only: is_master
    475   USE mod_grid_phy_lmdz, only: klon_glo
    476   USE dimphy, only: klev, klevp1
    477   USE tracer_mod, only: nqmx
    478   USE soil, only:  nsoilmx_PEM
    479   USE comslope_mod, only: nslope
    480   use layered_deposits, only: nb_str_max
    481   IMPLICIT NONE
    482     CHARACTER(LEN=*),INTENT(IN) :: filename
    483     INTEGER                     :: ierr
    484     LOGICAL,SAVE :: already_created=.false.
    485 
     937  USE mod_grid_phy_lmdz,  only: klon_glo
     938  USE dimphy,             only: klev, klevp1
     939  USE tracer_mod,         only: nqmx
     940  USE soil,               only:  nsoilmx_PEM
     941  USE comslope_mod,       only: nslope
     942  use layered_deposits,   only: nb_str_max
     943
     944! DECLARATION
     945! -----------
     946implicit none
     947
     948! ARGUMENTS
     949! ---------
     950  CHARACTER(LEN=*), INTENT(IN) :: filename
     951
     952! LOCAL VARIABLES
     953! ---------------
     954    INTEGER           :: ierr
     955    LOGICAL,SAVE      :: already_created=.false.
     956
     957! CODE
     958! ----
    486959    IF (is_master) THEN
    487960      if (.not.already_created) then
    … …  
    5971070
    5981071  END SUBROUTINE open_restartphy
    599 
     1072!=======================================================================
     1073
     1074!=======================================================================
    6001075  SUBROUTINE close_restartphy
    601   USE netcdf, only: NF90_CLOSE
     1076!-----------------------------------------------------------------------
     1077! NAME
     1078!     close_restartphy
     1079!
     1080! DESCRIPTION
     1081!     Close the restart netCDF file.
     1082!
     1083! AUTHORS & DATE
     1084!     L. Lange
     1085!     JB Clement, 2023-2025
     1086!
     1087! NOTES
     1088!
     1089!-----------------------------------------------------------------------
     1090
     1091! DEPENDENCIES
     1092! ------------
     1093  USE netcdf,             only: NF90_CLOSE
    6021094  USE mod_phys_lmdz_para, only: is_master
    603   IMPLICIT NONE
    604     INTEGER          :: ierr
    605 
     1095
     1096! DECLARATION
     1097! -----------
     1098implicit none
     1099
     1100! LOCAL VARIABLES
     1101! ---------------
     1102    INTEGER :: ierr
     1103
     1104! CODE
     1105! ----
    6061106    IF (is_master) THEN
    6071107      ierr = NF90_CLOSE (nid_restart)
    … …  
    6091109
    6101110  END SUBROUTINE close_restartphy
    611 
     1111!=======================================================================
     1112
     1113!=======================================================================
    6121114  SUBROUTINE put_field_r1(field_name,title,field,time)
    613   ! For a surface field
    614   IMPLICIT NONE
    615   CHARACTER(LEN=*),INTENT(IN)    :: field_name
    616   CHARACTER(LEN=*),INTENT(IN)    :: title
    617   REAL,INTENT(IN)                :: field(:)
    618   REAL,OPTIONAL,INTENT(IN)       :: time
    619 
     1115!-----------------------------------------------------------------------
     1116! NAME
     1117!     put_field_r1
     1118!
     1119! DESCRIPTION
     1120!     Write a surface field to the restart file.
     1121!
     1122! AUTHORS & DATE
     1123!     L. Lange
     1124!     JB Clement, 2023-2025
     1125!
     1126! NOTES
     1127!
     1128!-----------------------------------------------------------------------
     1129
     1130! DECLARATION
     1131! -----------
     1132implicit none
     1133
     1134
     1135! ARGUMENTS
     1136! ---------
     1137  CHARACTER(LEN=*), INTENT(IN)           :: field_name
     1138  CHARACTER(LEN=*), INTENT(IN)           :: title
     1139  REAL,             INTENT(IN)           :: field(:)
     1140  REAL,             INTENT(IN), OPTIONAL :: time
     1141
     1142! CODE
     1143! ----
    6201144  IF (present(time)) THEN
    6211145    ! if timeindex is present, it is a time-dependent variable
    … …  
    6261150
    6271151  END SUBROUTINE put_field_r1
    628 
     1152!=======================================================================
     1153
     1154!=======================================================================
    6291155  SUBROUTINE put_field_r2(field_name,title,field,time)
    630   ! For a "3D" horizontal-vertical field
    631   IMPLICIT NONE
    632   CHARACTER(LEN=*),INTENT(IN)    :: field_name
    633   CHARACTER(LEN=*),INTENT(IN)    :: title
    634   REAL,INTENT(IN)                :: field(:,:)
    635   REAL,OPTIONAL,INTENT(IN)       :: time
    636 
     1156!-----------------------------------------------------------------------
     1157! NAME
     1158!     put_field_r2
     1159!
     1160! DESCRIPTION
     1161!     Write a "3D" horizontal-vertical field (2D) to the restart file.
     1162
     1163! AUTHORS & DATE
     1164!     L. Lange
     1165!     JB Clement, 2023-2025
     1166
     1167! NOTES
     1168!
     1169!-----------------------------------------------------------------------
     1170
     1171! DECLARATION
     1172! -----------
     1173implicit none
     1174
     1175! ARGUMENTS
     1176! ---------
     1177  CHARACTER(LEN=*), INTENT(IN)           :: field_name
     1178  CHARACTER(LEN=*), INTENT(IN)           :: title
     1179  REAL,             INTENT(IN)           :: field(:,:)
     1180  REAL,             INTENT(IN), OPTIONAL :: time
     1181
     1182! CODE
     1183! ----
    6371184  IF (present(time)) THEN
    6381185    ! if timeindex is present, it is a time-dependent variable
    … …  
    6431190
    6441191  END SUBROUTINE put_field_r2
    645 
     1192!=======================================================================
     1193
     1194!=======================================================================
    6461195  SUBROUTINE put_field_r3(field_name,title,field,time)
    647   ! For a "4D" field surf/alt/??
    648   IMPLICIT NONE
    649   CHARACTER(LEN=*),INTENT(IN)    :: field_name
    650   CHARACTER(LEN=*),INTENT(IN)    :: title
    651   REAL,INTENT(IN)                :: field(:,:,:)
    652   REAL,OPTIONAL,INTENT(IN)       :: time
    653 
     1196!-----------------------------------------------------------------------
     1197! NAME
     1198!     put_field_r3
     1199!
     1200! DESCRIPTION
     1201!     Write a "4D" field surf/alt/?? (3D) to the restart file.
     1202
     1203! AUTHORS & DATE
     1204!     L. Lange
     1205!     JB Clement, 2023-2025
     1206
     1207! NOTES
     1208!
     1209!-----------------------------------------------------------------------
     1210
     1211! DECLARATION
     1212! -----------
     1213implicit none
     1214
     1215! ARGUMENTS
     1216! ---------
     1217  CHARACTER(LEN=*), INTENT(IN)           :: field_name
     1218  CHARACTER(LEN=*), INTENT(IN)           :: title
     1219  REAL,             INTENT(IN)           :: field(:,:,:)
     1220  REAL,             INTENT(IN), OPTIONAL :: time
     1221
     1222! CODE
     1223! ----
    6541224  IF (present(time)) THEN
    6551225    ! if timeindex is present, it is a time-dependent variable
    … …  
    6611231
    6621232  END SUBROUTINE put_field_r3
    663 
     1233!=======================================================================
     1234
     1235!=======================================================================
    6641236  SUBROUTINE put_field_rgen(field_name,title,field,field_size,time)
     1237!-----------------------------------------------------------------------
     1238! NAME
     1239!     put_field_rgen
     1240!
     1241! DESCRIPTION
     1242!     Generic subroutine to write fields to restart file with gather from processes
     1243
     1244! AUTHORS & DATE
     1245!     L. Lange
     1246!     JB Clement, 2023-2025
     1247
     1248! NOTES
     1249!
     1250!-----------------------------------------------------------------------
     1251
     1252! DEPENDENCIES
     1253! ------------
    6651254  USE netcdf
    6661255  USE dimphy
    667   USE soil, only: nsoilmx_PEM
     1256  USE soil,         only: nsoilmx_PEM
    6681257  USE comslope_mod, ONLY: nslope
    6691258  USE mod_grid_phy_lmdz
    6701259  USE mod_phys_lmdz_para
    671   IMPLICIT NONE
     1260
     1261! DECLARATION
     1262! -----------
     1263implicit none
     1264
     1265! ARGUMENTS
     1266! ---------
    6721267  CHARACTER(LEN=*),INTENT(IN)    :: field_name
    6731268  CHARACTER(LEN=*),INTENT(IN)    :: title
    … …  
    6761271  REAL,OPTIONAL,INTENT(IN)       :: time
    6771272
     1273! LOCAL VARIABLES
     1274! ---------------
    6781275  REAL                           :: field_glo(klon_glo,field_size)
    6791276  INTEGER                        :: ierr
    6801277  INTEGER                        :: nvarid
    6811278
     1279! CODE
     1280! ----
    6821281    CALL gather(field,field_glo)
    6831282
    … …  
    9571556
    9581557  END SUBROUTINE put_field_rgen
    959 
     1558!=======================================================================
     1559
     1560!=======================================================================
    9601561  SUBROUTINE put_var_r0(var_name,title,var)
    961   ! Put a scalar in file
    962    IMPLICIT NONE
    963      CHARACTER(LEN=*),INTENT(IN) :: var_name
    964      CHARACTER(LEN=*),INTENT(IN) :: title
    965      REAL,INTENT(IN)             :: var
    966      REAL                        :: varin(1)
    967 
     1562!-----------------------------------------------------------------------
     1563! NAME
     1564!     put_var_r0
     1565!
     1566! DESCRIPTION
     1567!     Write a scalar variable to the restart file
     1568
     1569! AUTHORS & DATE
     1570!     L. Lange
     1571!     JB Clement, 2023-2025
     1572
     1573! NOTES
     1574!
     1575!-----------------------------------------------------------------------
     1576
     1577! DECLARATION
     1578! -----------
     1579implicit none
     1580
     1581! ARGUMENTS
     1582! ---------
     1583  CHARACTER(LEN=*), INTENT(IN) :: var_name
     1584  CHARACTER(LEN=*), INTENT(IN) :: title
     1585  REAL,             INTENT(IN) :: var
     1586
     1587! LOCAL VARIABLES
     1588! ---------------
     1589     REAL :: varin(1)
     1590
     1591! CODE
     1592! ----
    9681593     varin(1)=var
    9691594
    … …  
    9711596
    9721597  END SUBROUTINE put_var_r0
    973 
    974 
     1598!=======================================================================
     1599
     1600!=======================================================================
    9751601  SUBROUTINE put_var_r1(var_name,title,var)
    976   ! Put a vector in file
    977    IMPLICIT NONE
    978      CHARACTER(LEN=*),INTENT(IN) :: var_name
    979      CHARACTER(LEN=*),INTENT(IN) :: title
    980      REAL,INTENT(IN)             :: var(:)
    981 
    982      CALL put_var_rgen(var_name,title,var,size(var))
     1602!-----------------------------------------------------------------------
     1603! NAME
     1604!    put_var_r1
     1605!
     1606! DESCRIPTION
     1607!     Write a vector (1D) variable to the restart file
     1608
     1609! AUTHORS & DATE
     1610!     L. Lange
     1611!     JB Clement, 2023-2025
     1612
     1613! NOTES
     1614!
     1615!-----------------------------------------------------------------------
     1616
     1617! DECLARATION
     1618! -----------
     1619implicit none
     1620
     1621! ARGUMENTS
     1622! ---------
     1623  CHARACTER(LEN=*), INTENT(IN) :: var_name
     1624  CHARACTER(LEN=*), INTENT(IN) :: title
     1625  REAL,             INTENT(IN) :: var(:)
     1626
     1627! CODE
     1628! ----
     1629  CALL put_var_rgen(var_name,title,var,size(var))
    9831630
    9841631  END SUBROUTINE put_var_r1
    985 
     1632!=======================================================================
     1633
     1634!=======================================================================
    9861635  SUBROUTINE put_var_r2(var_name,title,var)
    987   ! Put a 2D field in file
    988    IMPLICIT NONE
    989      CHARACTER(LEN=*),INTENT(IN) :: var_name
    990      CHARACTER(LEN=*),INTENT(IN) :: title
    991      REAL,INTENT(IN)             :: var(:,:)
    992 
    993      CALL put_var_rgen(var_name,title,var,size(var))
     1636!-----------------------------------------------------------------------
     1637! NAME
     1638!     put_var_r2
     1639!
     1640! DESCRIPTION
     1641!     Write a 3D field variable to the restart file
     1642
     1643! AUTHORS & DATE
     1644!     L. Lange
     1645!     JB Clement, 2023-2025
     1646
     1647! NOTES
     1648!
     1649!-----------------------------------------------------------------------
     1650
     1651! DECLARATION
     1652! -----------
     1653implicit none
     1654
     1655! ARGUMENTS
     1656! ---------
     1657  CHARACTER(LEN=*), INTENT(IN) :: var_name
     1658  CHARACTER(LEN=*), INTENT(IN) :: title
     1659  REAL,             INTENT(IN) :: var(:,:)
     1660
     1661! CODE
     1662! ----
     1663  CALL put_var_rgen(var_name,title,var,size(var))
    9941664
    9951665  END SUBROUTINE put_var_r2
    996 
     1666!=======================================================================
     1667
     1668!=======================================================================
    9971669  SUBROUTINE put_var_r3(var_name,title,var)
    998   ! Put a 3D field in file
    999    IMPLICIT NONE
    1000      CHARACTER(LEN=*),INTENT(IN) :: var_name
    1001      CHARACTER(LEN=*),INTENT(IN) :: title
    1002      REAL,INTENT(IN)             :: var(:,:,:)
    1003 
    1004      CALL put_var_rgen(var_name,title,var,size(var))
     1670!-----------------------------------------------------------------------
     1671! NAME
     1672!     put_var_r3
     1673!
     1674! DESCRIPTION
     1675!     Write a 3D field variable to the restart file
     1676!
     1677! AUTHORS & DATE
     1678!     L. Lange
     1679!     JB Clement, 2023-2025
     1680!
     1681! NOTES
     1682!
     1683!-----------------------------------------------------------------------
     1684
     1685! DECLARATION
     1686! -----------
     1687implicit none
     1688
     1689! ARGUMENTS
     1690! ---------
     1691  CHARACTER(LEN=*), INTENT(IN) :: var_name
     1692  CHARACTER(LEN=*), INTENT(IN) :: title
     1693  REAL,             INTENT(IN) :: var(:,:,:)
     1694
     1695! CODE
     1696! ----
     1697  CALL put_var_rgen(var_name,title,var,size(var))
    10051698
    10061699  END SUBROUTINE put_var_r3
    1007 
     1700!=======================================================================
     1701
     1702!=======================================================================
    10081703  SUBROUTINE put_var_rgen(var_name,title,var,var_size)
     1704!-----------------------------------------------------------------------
     1705! NAME
     1706!     put_var_rgen
     1707!
     1708! DESCRIPTION
     1709!     Generic subroutine to write variables to the restart file
     1710
     1711! AUTHORS & DATE
     1712!     L. Lange
     1713!     JB Clement, 2023-2025
     1714
     1715! NOTES
     1716!
     1717!-----------------------------------------------------------------------
     1718
     1719! DEPENDENCIES
     1720! ------------
    10091721  USE netcdf, only: NF90_REDEF, NF90_DEF_VAR, NF90_ENDDEF, NF90_PUT_VAR, &
    10101722                    NF90_FLOAT, NF90_DOUBLE, &
    10111723                    NF90_PUT_ATT, NF90_NOERR, nf90_strerror, &
    10121724                    nf90_inq_dimid, nf90_inquire_dimension, NF90_INQ_VARID
    1013   USE soil, only:  nsoilmx_PEM
    1014   USE comslope_mod, only: nslope
     1725  USE soil,               only:  nsoilmx_PEM
     1726  USE comslope_mod,       only: nslope
    10151727  USE mod_phys_lmdz_para, only: is_master
    1016   IMPLICIT NONE
    1017      CHARACTER(LEN=*),INTENT(IN) :: var_name
    1018      CHARACTER(LEN=*),INTENT(IN) :: title
    1019      INTEGER,INTENT(IN)          :: var_size
    1020      REAL,INTENT(IN)             :: var(var_size)
    1021 
    1022      INTEGER :: ierr
    1023      INTEGER :: nvarid
    1024      INTEGER :: idim1d
    1025      logical,save :: firsttime=.true.
    1026 
     1728
     1729! DECLARATION
     1730! -----------
     1731implicit none
     1732
     1733! ARGUMENTS
     1734! ---------
     1735  CHARACTER(LEN=*), INTENT(IN) :: var_name
     1736  CHARACTER(LEN=*), INTENT(IN) :: title
     1737  INTEGER,          INTENT(IN) :: var_size
     1738  REAL,             INTENT(IN) :: var(var_size)
     1739
     1740! LOCAL VARIABLES
     1741! ---------------
     1742  INTEGER       :: ierr
     1743  INTEGER       :: nvarid
     1744  INTEGER       :: idim1d
     1745  LOGICAL, SAVE :: firsttime = .true.
     1746
     1747! CODE
     1748! ----
    10271749    IF (is_master) THEN
    10281750
    … …  
    10991821
    11001822  END SUBROUTINE put_var_rgen
     1823!=======================================================================
    11011824
    11021825END MODULE iostart_pem

  • trunk/LMDZ.COMMON/libf/evolution/layered_deposits.F90

    r3989 r3991  
    11MODULE layered_deposits
    2 
    3 !=======================================================================
    4 ! Purpose: Manage layered layerings_map in the PEM with a linked list data structure
    5 !
    6 ! Author: JBC
    7 ! Date: 08/03/2024
    8 !=======================================================================
    9 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     layered_deposits
     5!
     6! DESCRIPTION
     7!     Manage layered deposits in the PEM with a linked list data structure.
     8!     Provides data types and subroutines to manipulate layers of ice and dust.
     9!
     10! AUTHORS & DATE
     11!     JB Clement, 2024
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
     16
     17! DEPENDENCIES
     18! ------------
    1019use surf_ice, only: rho_co2ice, rho_h2oice
    1120
    12 implicit none
    13 
    14 !---- Declarations
     21! DECLARATION
     22! -----------
     23implicit none
     24
     25! MODULE VARIABLES
     26! ----------------
    1527! Numerical parameter
    1628real, parameter :: eps = epsilon(1.), tol = 1.e2*eps
    … …  
    3749real, parameter :: fred_sublco2 = 0.1 ! Reduction factor of sublimation rate
    3850
     51! DATA STRUCTURES
     52! ---------------
    3953! Stratum = node of the linked list
    4054type :: stratum
    … …  
    6478end type ptrarray
    6579
    66 !=======================================================================
    6780contains
    68 ! Procedures to manipulate the data structure:
    69 !     > print_stratum
    70 !     > add_stratum
    71 !     > insert_stratum
    72 !     > del_stratum
    73 !     > get_stratum
    74 !     > ini_layering
    75 !     > del_layering
    76 !     > print_layering
    77 !     > get_nb_str_max
    78 !     > map2array
    79 !     > array2map
    80 !     > print_layerings_map
    81 ! Procedures to get information about a stratum:
    82 !     > thickness
    83 !     > is_co2ice_str
    84 !     > is_h2oice_str
    85 !     > is_dust_str
    86 !     > is_dust_lag
    87 !     > subsurface_ice_layering
    88 ! Procedures for the algorithm to evolve the layering:
    89 !     > make_layering
    90 !     > lift_dust
    91 !     > subl_ice_str
    92 !=======================================================================
    93 ! To display a stratum
     81!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     82
     83!=======================================================================
    9484SUBROUTINE print_stratum(this)
    95 
    96 implicit none
    97 
    98 !---- Arguments
     85!-----------------------------------------------------------------------
     86! NAME
     87!     print_stratum
     88!
     89! DESCRIPTION
     90!     Print a stratum.
     91!
     92! AUTHORS & DATE
     93!     JB Clement, 2024
     94!
     95! NOTES
     96!
     97!-----------------------------------------------------------------------
     98
     99! DECLARATION
     100! -----------
     101implicit none
     102
     103! ARGUMENTS
     104! ---------
    99105type(stratum), intent(in) :: this
    100106
    101 !---- Code
     107! CODE
     108! ----
    102109write(*,'(a,es13.6)') 'Top elevation       [m]: ', this%top_elevation
    103110write(*,'(a,es13.6)') 'Height of CO2 ice   [m]: ', this%h_co2ice
    … …  
    108115
    109116END SUBROUTINE print_stratum
    110 
    111 !=======================================================================
    112 ! To add a stratum to the top of a layering
     117!=======================================================================
     118
     119!=======================================================================
    113120SUBROUTINE add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice)
    114 
    115 implicit none
    116 
    117 !---- Arguments
    118 type(layering), intent(inout)  :: this
    119 real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice
    120 
    121 !---- Local variables
     121!-----------------------------------------------------------------------
     122! NAME
     123!     add_stratum
     124!
     125! DESCRIPTION
     126!     Add a stratum to the top of a layering.
     127!
     128! AUTHORS & DATE
     129!     JB Clement, 2024
     130!
     131! NOTES
     132!
     133!-----------------------------------------------------------------------
     134
     135! DECLARATION
     136! -----------
     137implicit none
     138
     139! ARGUMENTS
     140! ---------
     141type(layering), intent(inout)        :: this
     142real,           intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice
     143
     144! LOCAL VARIABLES
     145! ---------------
    122146type(stratum), pointer :: str
    123147
    124 !---- Code
     148! CODE
     149! ----
    125150! Creation of the stratum
    126151allocate(str)
    … …  
    153178
    154179END SUBROUTINE add_stratum
    155 
    156 !=======================================================================
    157 ! To insert a stratum after another one in a layering
     180!=======================================================================
     181
     182!=======================================================================
    158183SUBROUTINE insert_stratum(this,str_prev,top_elevation,co2ice,h2oice,dust,pore,poreice)
    159 
    160 implicit none
    161 
    162 !---- Arguments
    163 type(layering),         intent(inout) :: this
    164 type(stratum), pointer, intent(inout) :: str_prev
    165 real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice
    166 
    167 !---- Local variables
     184!-----------------------------------------------------------------------
     185! NAME
     186!     insert_stratum
     187!
     188! DESCRIPTION
     189!     Insert a stratum after a given previous stratum in the layering.
     190!
     191! AUTHORS & DATE
     192!     JB Clement, 2024
     193!
     194! NOTES
     195!
     196!-----------------------------------------------------------------------
     197
     198! DECLARATION
     199! -----------
     200implicit none
     201
     202! ARGUMENTS
     203! ---------
     204type(layering),         intent(inout)        ::  this
     205type(stratum), pointer, intent(inout)        ::  str_prev
     206real,                   intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice
     207
     208! LOCAL VARIABLES
     209! ---------------
    168210type(stratum), pointer :: str, current
    169211
    170 !---- Code
     212! CODE
     213! ----
    171214if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering
    172215    call add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice)
    … …  
    206249
    207250END SUBROUTINE insert_stratum
    208 
    209 !=======================================================================
    210 ! To remove a stratum in a layering
     251!=======================================================================
     252
     253!=======================================================================
    211254SUBROUTINE del_stratum(this,str)
    212 
    213 implicit none
    214 
    215 !---- Arguments
     255!-----------------------------------------------------------------------
     256! NAME
     257!     del_stratum
     258!
     259! DESCRIPTION
     260!     Delete a stratum from the layering and fix links.
     261!
     262! AUTHORS & DATE
     263!     JB Clement, 2024
     264!
     265! NOTES
     266!
     267!-----------------------------------------------------------------------
     268
     269! DECLARATION
     270! -----------
     271implicit none
     272
     273! ARGUMENTS
     274! ---------
    216275type(layering),         intent(inout) :: this
    217276type(stratum), pointer, intent(inout) :: str
    218277
    219 !---- Local variables
    220 type(stratum), pointer :: new_str ! To make the argument point towards something
    221 
    222 !---- Code
     278! LOCAL VARIABLES
     279! ---------------
     280type(stratum), pointer :: new_str
     281
     282! CODE
     283! ----
    223284! Protect the sub-surface strata from deletion
    224285if (str%top_elevation <= 0.) return
    … …  
    247308
    248309END SUBROUTINE del_stratum
     310!=======================================================================
    249311
    250312!=======================================================================
    251313! To get a specific stratum in a layering
    252314FUNCTION get_stratum(lay,i) RESULT(str)
    253 
    254 implicit none
    255 
    256 !---- Arguments
     315!-----------------------------------------------------------------------
     316! NAME
     317!     get_stratum
     318!
     319! DESCRIPTION
     320!     Return the i-th stratum in the layering.
     321!
     322! AUTHORS & DATE
     323!     JB Clement, 2024
     324!
     325! NOTES
     326!
     327!-----------------------------------------------------------------------
     328
     329! DECLARATION
     330! -----------
     331implicit none
     332
     333! ARGUMENTS
     334! ---------
    257335type(layering), intent(in) :: lay
    258336integer,        intent(in) :: i
    259337type(stratum), pointer     :: str
    260338
    261 !---- Local variables
     339! LOCAL VARIABLES
     340! ---------------
    262341integer :: k
    263342
    264 !---- Code
     343! CODE
     344! ----
    265345if (i < 1 .or. i > lay%nb_str) error stop 'get_stratum: requested number is beyond the number of strata of the layering!'
    266346k = 1
    … …  
    272352
    273353END FUNCTION get_stratum
     354!=======================================================================
    274355
    275356!=======================================================================
    276357! To initialize the layering
    277358SUBROUTINE ini_layering(this)
    278 
     359!-----------------------------------------------------------------------
     360! NAME
     361!     del_layering
     362!
     363! DESCRIPTION
     364!     Delete all strata in a layering and reset counters.
     365!
     366! AUTHORS & DATE
     367!     JB Clement, 2024
     368!
     369! NOTES
     370!
     371!-----------------------------------------------------------------------
     372
     373! DEPENDENCIES
     374! ------------
    279375use soil, only: do_soil, nsoilmx_PEM, layer_PEM, index_breccia, index_bedrock
    280376
    281 implicit none
    282 
    283 !---- Arguments
     377! DECLARATION
     378! -----------
     379implicit none
     380
     381! ARGUMENTS
     382! ---------
    284383type(layering), intent(inout) :: this
    285384
    286 !---- Local variables
     385! LOCAL VARIABLES
     386! ---------------
    287387integer :: i
    288388real    :: h_soil, h_pore, h_tot
    289389
    290 !---- Code
     390! CODE
     391! ----
    291392! Creation of strata at the bottom of the layering to describe the sub-surface
    292393if (do_soil) then
    … …  
    313414
    314415END SUBROUTINE ini_layering
    315 
    316 !=======================================================================
    317 ! To delete the entire layering
     416!=======================================================================
     417
     418!=======================================================================
    318419SUBROUTINE del_layering(this)
    319 
    320 implicit none
    321 
    322 !---- Arguments
     420!-----------------------------------------------------------------------
     421! NAME
     422!     del_layering
     423!
     424! DESCRIPTION
     425!     Delete all strata in a layering and reset counters.
     426!
     427! AUTHORS & DATE
     428!     JB Clement, 2024
     429!
     430! NOTES
     431!
     432!-----------------------------------------------------------------------
     433
     434! DECLARATION
     435! -----------
     436implicit none
     437
     438! ARGUMENTS
     439! ---------
    323440type(layering), intent(inout) :: this
    324441
    325 !---- Local variables
     442! LOCAL VARIABLES
     443! ---------------
    326444type(stratum), pointer :: current, next
    327445
    328 !---- Code
     446! CODE
     447! ----
    329448current => this%bottom
    330449do while (associated(current))
    … …  
    337456
    338457END SUBROUTINE del_layering
    339 
    340 !=======================================================================
    341 ! To display a layering
     458!=======================================================================
     459
     460!=======================================================================
    342461SUBROUTINE print_layering(this)
    343 
    344 implicit none
    345 
    346 !---- Arguments
     462!-----------------------------------------------------------------------
     463! NAME
     464!     print_layering
     465!
     466! DESCRIPTION
     467!     Display all strata from top to bottom for a layering.
     468!
     469! AUTHORS & DATE
     470!     JB Clement, 2024
     471!
     472! NOTES
     473!
     474!-----------------------------------------------------------------------
     475
     476! DECLARATION
     477! -----------
     478implicit none
     479
     480! ARGUMENTS
     481! ---------
    347482type(layering), intent(in) :: this
    348483
    349 !---- Local variables
     484! LOCAL VARIABLES
     485! ---------------
    350486type(stratum), pointer :: current
    351487integer                :: i
    352488
    353 !---- Code
     489! CODE
     490! ----
    354491current => this%top
    355492
    … …  
    364501
    365502END SUBROUTINE print_layering
    366 
    367 !=======================================================================
    368 ! To get the maximum number of "stratum" in the stratification (layerings)
     503!=======================================================================
     504
     505!=======================================================================
    369506FUNCTION get_nb_str_max(layerings_map,ngrid,nslope) RESULT(nb_str_max)
    370 
    371 implicit none
    372 
    373 !---- Arguments
     507!-----------------------------------------------------------------------
     508! NAME
     509!     get_nb_str_max
     510!
     511! DESCRIPTION
     512!     Return the maximum number of strata across all grid/slope layerings.
     513!
     514! AUTHORS & DATE
     515!     JB Clement, 2024
     516!
     517! NOTES
     518!
     519!-----------------------------------------------------------------------
     520
     521! DECLARATION
     522! -----------
     523implicit none
     524
     525! ARGUMENTS
     526! ---------
    374527type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map
    375528integer,                                 intent(in) :: ngrid, nslope
    376 integer :: nb_str_max
    377 
    378 !---- Local variables
     529integer                                             :: nb_str_max
     530
     531! LOCAL VARIABLES
     532! ---------------
    379533integer :: ig, islope
    380534
    381 !---- Code
     535! CODE
     536! ----
    382537nb_str_max = 0
    383538do islope = 1,nslope
    … …  
    388543
    389544END FUNCTION get_nb_str_max
    390 
    391 !=======================================================================
    392 ! To convert the layerings map into an array able to be outputted
     545!=======================================================================
     546
     547!=======================================================================
    393548SUBROUTINE map2array(layerings_map,ngrid,nslope,layerings_array)
    394 
    395 implicit none
    396 
    397 !---- Arguments
    398 integer,                                 intent(in) :: ngrid, nslope
    399 type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map
    400 real, dimension(:,:,:,:), allocatable, intent(inout) :: layerings_array
    401 
    402 !---- Local variables
     549!-----------------------------------------------------------------------
     550! NAME
     551!     map2array
     552!
     553! DESCRIPTION
     554!     Convert the layerings map into an allocatable array for output.
     555!
     556! AUTHORS & DATE
     557!     JB Clement, 2024
     558!
     559! NOTES
     560!
     561!-----------------------------------------------------------------------
     562
     563! DECLARATION
     564! -----------
     565implicit none
     566
     567! ARGUMENTS
     568! ---------
     569integer,                                 intent(in)    :: ngrid, nslope
     570type(layering), dimension(ngrid,nslope), intent(in)    :: layerings_map
     571real, dimension(:,:,:,:), allocatable,   intent(inout) :: layerings_array
     572
     573! LOCAL VARIABLES
     574! ---------------
    403575type(stratum), pointer :: current
    404576integer                :: ig, islope, k
    405577
    406 !---- Code
     578! CODE
     579! ----
    407580layerings_array = 0.
    408581do islope = 1,nslope
    … …  
    424597
    425598END SUBROUTINE map2array
    426 
    427 !=======================================================================
    428 ! To convert the stratification array into the layerings map
     599!=======================================================================
     600
     601!=======================================================================
    429602SUBROUTINE array2map(layerings_array,ngrid,nslope,layerings_map)
    430 
    431 implicit none
    432 
    433 !---- Arguments
    434 integer,                               intent(in) :: ngrid, nslope
    435 real, dimension(:,:,:,:), allocatable, intent(in) :: layerings_array
     603!-----------------------------------------------------------------------
     604! NAME
     605!     array2map
     606!
     607! DESCRIPTION
     608!     Convert the stratification array back into the layerings map.
     609!
     610! AUTHORS & DATE
     611!     JB Clement, 2024
     612!
     613! NOTES
     614!
     615!-----------------------------------------------------------------------
     616
     617! DECLARATION
     618! -----------
     619implicit none
     620
     621! ARGUMENTS
     622! ---------
     623integer,                                 intent(in)    :: ngrid, nslope
     624real, dimension(:,:,:,:), allocatable,   intent(in)    :: layerings_array
    436625type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map
    437626
    438 !---- Local variables
     627! LOCAL VARIABLES
     628! ---------------
    439629integer :: ig, islope, k
    440630
    441 !---- Code
     631! CODE
     632! ----
    442633do islope = 1,nslope
    443634    do ig = 1,ngrid
    … …  
    449640
    450641END SUBROUTINE array2map
    451 
    452 !=======================================================================
    453 ! To display the map of layerings
     642!=======================================================================
     643
     644!=======================================================================
    454645SUBROUTINE print_layerings_map(this,ngrid,nslope)
    455 
    456 implicit none
    457 
    458 !---- Arguments
     646!-----------------------------------------------------------------------
     647! NAME
     648!     print_layerings_map
     649!
     650! DESCRIPTION
     651!     Display the entire layerings map for all grids and slopes.
     652!
     653! AUTHORS & DATE
     654!     JB Clement, 2024
     655!
     656! NOTES
     657!
     658!-----------------------------------------------------------------------
     659
     660! DECLARATION
     661! -----------
     662implicit none
     663
     664! ARGUMENTS
     665! ---------
    459666integer,                                 intent(in) :: ngrid, nslope
    460667type(layering), dimension(ngrid,nslope), intent(in) :: this
    461668
    462 !---- Local variables
     669! LOCAL VARIABLES
     670! ---------------
    463671integer :: ig, islope
    464672
    465 !---- Code
     673! CODE
     674! ----
    466675do ig = 1,ngrid
    467676    write(*,'(a10,i4)') 'Grid cell ', ig
    … …  
    475684
    476685END SUBROUTINE print_layerings_map
    477 
    478 !=======================================================================
    479 !=======================================================================
    480 ! To compute the thickness of a stratum
     686!=======================================================================
     687
     688!=======================================================================
    481689FUNCTION thickness(this) RESULT(h_str)
    482 
    483 implicit none
    484 
    485 !---- Arguments
     690!-----------------------------------------------------------------------
     691! NAME
     692!     thickness
     693!
     694! DESCRIPTION
     695!     Compute the total thickness of a stratum.
     696!
     697! AUTHORS & DATE
     698!     JB Clement, 2024
     699!
     700! NOTES
     701!
     702!-----------------------------------------------------------------------
     703
     704! DECLARATION
     705! -----------
     706implicit none
     707
     708! ARGUMENTS
     709! ---------
    486710type(stratum), intent(in) :: this
     711
     712! LOCAL VARIABLES
     713! ---------------
    487714real :: h_str
    488715
    489 !---- Code
     716! CODE
     717! ----
    490718h_str = this%h_co2ice + this%h_h2oice + this%h_dust + this%h_pore
    491719
    492720END FUNCTION thickness
    493 
    494 !=======================================================================
    495 ! To know whether the stratum can be considered as a CO2 ice layer
     721!=======================================================================
     722
     723!=======================================================================
    496724FUNCTION is_co2ice_str(str) RESULT(co2ice_str)
    497 
    498 implicit none
    499 
    500 !---- Arguments
     725!-----------------------------------------------------------------------
     726! NAME
     727!     is_co2ice_str
     728!
     729! DESCRIPTION
     730!     Check if a stratum can be considered as a CO2 ice layer.
     731!
     732! AUTHORS & DATE
     733!     JB Clement, 2024
     734!
     735! NOTES
     736!
     737!-----------------------------------------------------------------------
     738
     739! DECLARATION
     740! -----------
     741implicit none
     742
     743! ARGUMENTS
     744! ---------
    501745type(stratum), intent(in) :: str
     746
     747! LOCAL VARIABLES
     748! ---------------
    502749logical :: co2ice_str
    503750
    504 !---- Code
     751! CODE
     752! ----
    505753co2ice_str = .false.
    506754if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) co2ice_str = .true.
    507755
    508756END FUNCTION is_co2ice_str
    509 
    510 !=======================================================================
    511 ! To know whether the stratum can be considered as a H2O ice layer
     757!=======================================================================
     758
     759!=======================================================================
    512760FUNCTION is_h2oice_str(str) RESULT(h2oice_str)
    513 
    514 implicit none
    515 
    516 !---- Arguments
     761!-----------------------------------------------------------------------
     762! NAME
     763!     is_h2oice_str
     764!
     765! DESCRIPTION
     766!     Check if a stratum can be considered as a H2O ice layer.
     767!
     768! AUTHORS & DATE
     769!     JB Clement, 2024
     770!
     771! NOTES
     772!
     773!-----------------------------------------------------------------------
     774
     775! DECLARATION
     776! -----------
     777implicit none
     778
     779! ARGUMENTS
     780! ---------
    517781type(stratum), intent(in) :: str
     782
     783! LOCAL VARIABLES
     784! ---------------
    518785logical :: h2oice_str
    519786
    520 !---- Code
     787! CODE
     788! ----
    521789h2oice_str = .false.
    522790if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) h2oice_str = .true.
    523791
    524792END FUNCTION is_h2oice_str
    525 
    526 !=======================================================================
    527 ! To know whether the stratum can be considered as a dust layer
     793!=======================================================================
     794
     795!=======================================================================
    528796FUNCTION is_dust_str(str) RESULT(dust_str)
    529 
    530 implicit none
    531 
    532 !---- Arguments
     797!-----------------------------------------------------------------------
     798! NAME
     799!     is_dust_str
     800!
     801! DESCRIPTION
     802!     Check if a stratum can be considered as a dust layer.
     803!
     804! AUTHORS & DATE
     805!     JB Clement, 2024
     806!
     807! NOTES
     808!
     809!-----------------------------------------------------------------------
     810
     811! DECLARATION
     812! -----------
     813implicit none
     814
     815! ARGUMENTS
     816! ---------
    533817type(stratum), intent(in) :: str
     818
     819! LOCAL VARIABLES
     820! ---------------
    534821logical :: dust_str
    535822
    536 !---- Code
     823! CODE
     824! ----
    537825dust_str = .false.
    538826if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) dust_str = .true.
    539827
    540828END FUNCTION is_dust_str
    541 
    542 !=======================================================================
    543 ! To know whether the stratum can be considered as a dust lag
     829!=======================================================================
     830
     831!=======================================================================
    544832FUNCTION is_dust_lag(str) RESULT(dust_lag)
    545 
    546 implicit none
    547 
    548 !---- Arguments
     833!-----------------------------------------------------------------------
     834! NAME
     835!     is_dust_lag
     836!
     837! DESCRIPTION
     838!     Check if a stratum is a dust lag (no ice content).
     839!
     840! AUTHORS & DATE
     841!     JB Clement, 2024
     842!
     843! NOTES
     844!
     845!-----------------------------------------------------------------------
     846
     847! DECLARATION
     848! -----------
     849implicit none
     850
     851! ARGUMENTS
     852! ---------
    549853type(stratum), intent(in) :: str
     854
     855! LOCAL VARIABLES
     856! ---------------
    550857logical :: dust_lag
    551858
    552 !---- Code
     859! CODE
     860! ----
    553861dust_lag = .false.
    554862if (str%h_co2ice < tol .and. str%h_h2oice < tol .and. str%poreice_volfrac < tol) dust_lag = .true.
    555863
    556864END FUNCTION is_dust_lag
    557 
    558 !=======================================================================
    559 ! To get data about possible subsurface ice in a layering
     865!=======================================================================
     866
     867!=======================================================================
    560868SUBROUTINE subsurface_ice_layering(this,h_toplag,h2o_ice,co2_ice)
    561 
    562 implicit none
    563 
    564 !---- Arguments
    565 type(layering), intent(in) :: this
    566 real, intent(out) :: h2o_ice, co2_ice, h_toplag
    567 
    568 !---- Local variables
     869!-----------------------------------------------------------------------
     870! NAME
     871!     subsurface_ice_layering
     872!
     873! DESCRIPTION
     874!     Get data about possible subsurface ice in a layering.
     875!
     876! AUTHORS & DATE
     877!     JB Clement, 2024
     878!
     879! NOTES
     880!
     881!-----------------------------------------------------------------------
     882
     883! DECLARATION
     884! -----------
     885implicit none
     886
     887! ARGUMENTS
     888! ---------
     889type(layering), intent(in)  :: this
     890real,           intent(out) :: h2o_ice, co2_ice, h_toplag
     891
     892! LOCAL VARIABLES
     893! ---------------
    569894type(stratum), pointer :: current
    570895
    571 !---- Code
     896! CODE
     897! ----
    572898h_toplag = 0.
    573899h2o_ice = 0.
    … …  
    602928
    603929END SUBROUTINE subsurface_ice_layering
    604 
    605 !=======================================================================
    606 ! To lift dust in a stratum
     930!=======================================================================
     931
     932!=======================================================================
    607933SUBROUTINE lift_dust(this,str,h2lift)
    608 
    609 implicit none
    610 
    611 !---- Arguments
     934!-----------------------------------------------------------------------
     935! NAME
     936!     lift_dust
     937!
     938! DESCRIPTION
     939!     Lift dust in a stratum by wind erosion or other processes.
     940!
     941! AUTHORS & DATE
     942!     JB Clement, 2024
     943!
     944! NOTES
     945!
     946!-----------------------------------------------------------------------
     947
     948! DECLARATION
     949! -----------
     950implicit none
     951
     952! ARGUMENTS
     953! ---------
    612954type(layering),         intent(inout) :: this
    613955type(stratum), pointer, intent(inout) :: str
    614 real, intent(in) :: h2lift
    615 
    616 !---- Code
     956real,                   intent(in)    :: h2lift
     957
     958! CODE
     959! ----
    617960! Update of properties in the eroding dust lag
    618961str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust)
    … …  
    624967
    625968END SUBROUTINE lift_dust
    626 
    627 !=======================================================================
    628 ! To sublimate ice in a stratum
     969!=======================================================================
     970
     971!=======================================================================
    629972SUBROUTINE sublimate_ice(this,str,h2subl_co2ice,h2subl_h2oice,new_lag,zlag)
    630 
    631 implicit none
    632 
    633 !---- Arguments
     973!-----------------------------------------------------------------------
     974! NAME
     975!     sublimate_ice
     976!
     977! DESCRIPTION
     978!     Handle ice sublimation in a stratum and create lag layer.
     979!
     980! AUTHORS & DATE
     981!     JB Clement, 2024
     982!
     983! NOTES
     984!
     985!-----------------------------------------------------------------------
     986
     987! DECLARATION
     988! -----------
     989implicit none
     990
     991! ARGUMENTS
     992! ---------
    634993type(layering),         intent(inout) :: this
    635994type(stratum), pointer, intent(inout) :: str
    636995logical,                intent(inout) :: new_lag
    637996real,                   intent(inout) :: zlag
    638 real, intent(in) :: h2subl_co2ice, h2subl_h2oice
    639 
    640 !---- Local variables
     997real,                   intent(in)    :: h2subl_co2ice, h2subl_h2oice
     998
     999! LOCAL VARIABLES
     1000! ---------------
    6411001real                   :: h2subl, h_ice, h_pore, h_pore_new, h_tot
    6421002real                   :: hlag_dust, hlag_h2oice
    6431003type(stratum), pointer :: current
    6441004
    645 !---- Code
     1005! CODE
     1006! ----
    6461007h_ice = str%h_co2ice + str%h_h2oice
    6471008h_pore = str%h_pore
    … …  
    6921053
    6931054END SUBROUTINE sublimate_ice
    694 
    695 !=======================================================================
    696 ! Layering algorithm
     1055!=======================================================================
     1056
     1057!=======================================================================
    6971058SUBROUTINE make_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current)
    698 
    699 implicit none
    700 
    701 !---- Arguments
    702 ! Inputs
    703 
    704 ! Outputs
    705 type(layering),         intent(inout) :: this
    706 type(stratum), pointer, intent(inout) :: current
    707 real,                   intent(inout) :: d_co2ice, d_h2oice
    708 logical,                intent(inout) :: new_str, new_lag
    709 real, intent(out) :: zshift_surf, zlag
    710 
    711 !---- Local variables
     1059!-----------------------------------------------------------------------
     1060! NAME
     1061!     make_layering
     1062!
     1063! DESCRIPTION
     1064!     Main algorithm for layering evolution. Manages ice condensation/sublimation
     1065!     and dust sedimentation/lifting. Handles multiple scenarios: building new
     1066!     strata from ice/dust, retreating strata from sublimation/lifting, and mixed
     1067!     scenarios with CO2 ice sublimation combined with H2O ice condensation.
     1068!
     1069! AUTHORS & DATE
     1070!     JB Clement, 2024
     1071!
     1072! NOTES
     1073!     Creates dust lag layers when ice sublimation occurs.
     1074!-----------------------------------------------------------------------
     1075
     1076! DECLARATION
     1077! -----------
     1078implicit none
     1079
     1080! ARGUMENTS
     1081! ---------
     1082type(layering),         intent(inout) ::  this
     1083type(stratum), pointer, intent(inout) ::  current
     1084real,                   intent(inout) ::  d_co2ice, d_h2oice
     1085logical,                intent(inout) ::  new_str, new_lag
     1086real,                   intent(out)   ::  zshift_surf, zlag
     1087
     1088! LOCAL VARIABLES
     1089! ---------------
    7121090real                   :: dh_co2ice, dh_h2oice, dh_dust, dh_dust_co2, dh_dust_h2o
    7131091real                   :: h_co2ice_old, h_h2oice_old, h_dust_old
    … …  
    7171095type(stratum), pointer :: currentloc
    7181096
    719 !---- Code
     1097! CODE
     1098! ----
    7201099dh_co2ice = d_co2ice/rho_co2ice
    7211100dh_h2oice = d_h2oice/rho_h2oice
    … …  
    7361115unexpected = .false.
    7371116
    738 !-----------------------------------------------------------------------
     1117!~~~~~~~~~~
    7391118! NOTHING HAPPENS
    7401119if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then
    7411120    ! So we do nothing
    742 !-----------------------------------------------------------------------
     1121!~~~~~~~~~~
    7431122! BUILDING A STRATUM (ice condensation and/or dust desimentation)
    7441123else if (dh_co2ice >= 0. .and. dh_h2oice >= 0. .and. dh_dust >= 0.) then
    … …  
    7651144        this%top%h_pore        = this%top%h_pore        + h_pore
    7661145    endif
    767 !-----------------------------------------------------------------------
     1146!~~~~~~~~~~
    7681147 ! RETREATING A STRATUM (ice sublimation and/or dust lifting)
    7691148else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. dh_dust <= 0.) then
    … …  
    8371216        if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0
    8381217        if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
    839 !-----------------------------------------------------------------------
     1218!~~~~~~~~~~
    8401219! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION
    8411220else if (dh_co2ice <= 0. .and. dh_h2oice >= 0.) then
    … …  
    10001379        if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0
    10011380    endif
    1002 !-----------------------------------------------------------------------
     1381!~~~~~~~~~~
    10031382! NOT EXPECTED SITUATION
    10041383else
    … …  
    10171396
    10181397END SUBROUTINE make_layering
     1398!=======================================================================
    10191399
    10201400END MODULE layered_deposits

  • trunk/LMDZ.COMMON/libf/evolution/math_toolkit.F90

    r3989 r3991  
    11module math_toolkit
    2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    3 !!!
    4 !!! Purpose: The module contains all the mathematical SUBROUTINE used in the PEM
    5 !!!
    6 !!! Author: Adapted from Schorgofer MSIM (N.S, Icarus 2010), impletented here by LL
    7 !!!
    8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    9 
    10 implicit none
    11 
    12 !=======================================================================
    13   contains
    14 !=======================================================================
    15 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     math_toolkit
     5!
     6! DESCRIPTION
     7!     The module contains all the mathematical subroutines used in the PEM.
     8!
     9! AUTHORS & DATE
     10!     L. Lange
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!     Adapted from Schorghofer MSIM (N.S, Icarus 2010)
     15!-----------------------------------------------------------------------
     16
     17! DECLARATION
     18! -----------
     19implicit none
     20
     21contains
     22!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     23
     24!=======================================================================
    1625SUBROUTINE deriv1(z,nz,y,y0,ybot,dzY)
    17 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    18 !!!
    19 !!! Purpose: Compute the first derivative of a function y(z) on an irregular grid
    20 !!!
    21 !!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
    22 !!!
    23 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    24 ! first derivative of a function y(z) on irregular grid
    25 ! upper boundary conditions: y(0) = y0
    26 ! lower boundary condition.: yp = ybottom
    27 
    28 implicit none
    29 
    30 ! Inputs
    31 !-------
    32 integer,             intent(in) :: nz       ! number of layer
    33 real, dimension(nz), intent(in) :: z        ! depth layer
    34 real, dimension(nz), intent(in) :: y        ! function which needs to be derived
    35 real,                intent(in) :: y0, ybot ! boundary conditions
    36 ! Outputs
    37 !--------
    38 real, dimension(nz), intent(out) :: dzY ! derivative of y w.r.t depth
    39 ! Local variables
    40 !----------------
     26!-----------------------------------------------------------------------
     27! NAME
     28!     deriv1
     29!
     30! DESCRIPTION
     31!     Compute the first derivative of a function y(z) on an irregular grid.
     32!
     33! AUTHORS & DATE
     34!     N.S (Icarus 2010)
     35!     L. Lange
     36!     JB Clement, 2023-2025
     37!
     38! NOTES
     39!     Upper boundary conditions: y(0) = y0
     40!     Lower boundary condition: yp = ybottom
     41!-----------------------------------------------------------------------
     42
     43! DECLARATION
     44! -----------
     45implicit none
     46
     47! ARGUMENTS
     48! ---------
     49integer,             intent(in)  :: nz       ! number of layer
     50real, dimension(nz), intent(in)  :: z        ! depth layer
     51real, dimension(nz), intent(in)  :: y        ! function which needs to be derived
     52real,                intent(in)  :: y0, ybot ! boundary conditions
     53real, dimension(nz), intent(out) :: dzY      ! derivative of y w.r.t depth
     54
     55! LOCAL VARIABLES
     56! ---------------
    4157integer :: j
    4258real    :: hm, hp, c1, c2, c3
    4359
     60! CODE
     61! ----
    4462hp = z(2) - z(1)
    4563c1 = z(1)/(hp*z(2))
    … …  
    5876
    5977END SUBROUTINE deriv1
    60 
    61 !=======================================================================
    62 
     78!=======================================================================
     79
     80!=======================================================================
    6381SUBROUTINE deriv2_simple(z,nz,y,y0,yNp1,yp2)
    64 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    65 !!!
    66 !!! Purpose: Compute the second derivative of a function y(z) on an irregular grid
    67 !!!
    68 !!! Author: N.S (raw copy/paste from MSIM)
    69 !!!
    70 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    71 ! second derivative y_zz on irregular grid
    72 ! boundary conditions: y(0) = y0, y(nz + 1) = yNp1
    73 
    74 implicit none
    75 
    76 ! Inputs
    77 !-------
    78 integer,             intent(in) :: nz
    79 real, dimension(nz), intent(in) :: z, y
    80 real,                intent(in) :: y0, yNp1
    81 ! Outputs
    82 !--------
     82!-----------------------------------------------------------------------
     83! NAME
     84!     deriv2_simple
     85!
     86! DESCRIPTION
     87!     Compute the second derivative of a function y(z) on an irregular grid.
     88!
     89! AUTHORS & DATE
     90!     N.S (Icarus 2010)
     91!     JB Clement, 2023-2025
     92!
     93! NOTES
     94!     Boundary conditions: y(0) = y0, y(nz + 1) = yNp1
     95!-----------------------------------------------------------------------
     96
     97! DECLARATION
     98! -----------
     99implicit none
     100
     101! ARGUMENTS
     102! ---------
     103integer,             intent(in)  :: nz
     104real, dimension(nz), intent(in)  :: z, y
     105real,                intent(in)  :: y0, yNp1
    83106real, dimension(nz), intent(out) :: yp2
    84 ! Local variables
    85 !----------------
     107
     108! LOCAL VARIABLES
     109! ---------------
    86110integer :: j
    87111real    :: hm, hp, c1, c2, c3
    88112
     113! CODE
     114! ----
    89115c1 = +2./((z(2) - z(1))*z(2))
    90116c2 = -2./((z(2) - z(1))*z(1))
    … …  
    102128
    103129END SUBROUTINE deriv2_simple
    104 
    105 !=======================================================================
    106 
     130!=======================================================================
     131
     132!=======================================================================
    107133SUBROUTINE  deriv1_onesided(j,z,nz,y,dy_zj)
    108 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    109 !!!
    110 !!! Purpose: First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
    111 !!!
    112 !!! Author: N.S (raw copy/paste from MSIM)
    113 !!!
    114 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    115 
    116 implicit none
    117 
    118 ! Inputs
    119 !-------
    120 integer,             intent(in) :: nz, j
    121 real, dimension(nz), intent(in) :: z, y
    122 ! Outputs
    123 !--------
    124 real, intent(out) :: dy_zj
    125 ! Local viariables
    126 !-----------------
     134!-----------------------------------------------------------------------
     135! NAME
     136!     deriv1_onesided
     137!
     138! DESCRIPTION
     139!     First derivative of function y(z) at z(j) one-sided derivative
     140!     on irregular grid.
     141!
     142! AUTHORS & DATE
     143!     N.S (Icarus 2010)
     144!     JB Clement, 2023-2025
     145!
     146! NOTES
     147!
     148!-----------------------------------------------------------------------
     149
     150! DECLARATION
     151! -----------
     152implicit none
     153
     154! ARGUMENTS
     155! ---------
     156integer,             intent(in)  :: nz, j
     157real, dimension(nz), intent(in)  :: z, y
     158real,                intent(out) :: dy_zj
     159
     160! LOCAL VARIABLES
     161! ---------------
    127162real :: h1, h2, c1, c2, c3
    128163
     164! CODE
     165! ----
    129166if (j < 1 .or. j > nz - 2) then
    130167    dy_zj = -1.
    … …  
    139176
    140177END SUBROUTINE deriv1_onesided
    141 
    142 !=======================================================================
    143 
     178!=======================================================================
     179
     180!=======================================================================
    144181PURE SUBROUTINE colint(y,z,nz,i1,i2,integral)
    145 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    146 !!!
    147 !!! Purpose:  Column integrates y on irregular grid
    148 !!!
    149 !!! Author: N.S (raw copy/paste from MSIM)
    150 !!!
    151 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    152 
    153 implicit none
    154 
    155 ! Inputs
    156 !-------
    157 integer,             intent(in) :: nz, i1, i2
    158 real, dimension(nz), intent(in) :: y, z
    159 ! Outputs
    160 !--------
    161 real, intent(out) :: integral
    162 ! Local viariables
    163 !-----------------
     182!-----------------------------------------------------------------------
     183! NAME
     184!     colint
     185!
     186! DESCRIPTION
     187!     Column integrates y on irregular grid.
     188!
     189! AUTHORS & DATE
     190!     N.S (Icarus 2010)
     191!     JB Clement, 2023-2025
     192!
     193! NOTES
     194!
     195!-----------------------------------------------------------------------
     196
     197! DECLARATION
     198! -----------
     199implicit none
     200
     201! ARGUMENTS
     202! ---------
     203integer,             intent(in)  :: nz, i1, i2
     204real, dimension(nz), intent(in)  :: y, z
     205real,                intent(out) :: integral
     206
     207! LOCAL VARIABLES
     208! ---------------
    164209integer             :: i
    165210real, dimension(nz) :: dz
    166211
     212! CODE
     213! ----
    167214dz(1) = (z(2) - 0.)/2
    168215do i = 2,nz - 1
    … …  
    173220
    174221END SUBROUTINE colint
    175 
    176 !=======================================================================
    177 
     222!=======================================================================
     223
     224!=======================================================================
    178225SUBROUTINE findroot(y1,y2,z1,z2,zr)
    179 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    180 !!!
    181 !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
    182 !!!
    183 !!! Author: LL
    184 !!!
    185 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    186 
    187 implicit none
    188 
    189 ! Inputs
    190 !-------
    191 real, intent(in) :: y1, y2 ! difference between surface water density and at depth  [kg/m^3]
    192 real, intent(in) :: z1, z2 ! depth [m]
    193 ! Outputs
    194 !--------
    195 real, intent(out) :: zr ! depth at which we have zero
    196 
     226!-----------------------------------------------------------------------
     227! NAME
     228!     findroot
     229!
     230! DESCRIPTION
     231!     Compute the root zr, between two values y1 and y2 at depth z1,z2.
     232!
     233! AUTHORS & DATE
     234!     L. Lange
     235!     JB Clement, 2023-2025
     236!
     237! NOTES
     238!
     239!-----------------------------------------------------------------------
     240
     241! DECLARATION
     242! -----------
     243implicit none
     244
     245! ARGUMENTS
     246! ---------
     247real, intent(in)  :: y1, y2
     248real, intent(in)  :: z1, z2
     249real, intent(out) :: zr
     250
     251! CODE
     252! ----
    197253zr = (y1*z2 - y2*z1)/(y1 - y2)
    198254
    199255END SUBROUTINE findroot
    200 
    201 !=======================================================================
    202 
     256!=======================================================================
     257
     258!=======================================================================
    203259SUBROUTINE solve_tridiag(a,b,c,d,n,x,error)
    204 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    205 !!!
    206 !!! Purpose: Solve a tridiagonal system Ax = d using the Thomas' algorithm
    207 !!!          a: sub-diagonal
    208 !!!          b: main diagonal
    209 !!!          c: super-diagonal
    210 !!!          d: right-hand side
    211 !!!          x: solution
    212 !!!
    213 !!! Author: JBC
    214 !!!
    215 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    216 
    217 implicit none
    218 
    219 ! Inputs
    220 !-------
    221 integer,                intent(in) :: n
    222 real, dimension(n),     intent(in) :: b, d
    223 real, dimension(n - 1), intent(in) :: a, c
    224 ! Outputs
    225 !--------
    226 real, dimension(n), intent(out) :: x
    227 integer,            intent(out) :: error
    228 ! Local viariables
    229 !-----------------
     260!-----------------------------------------------------------------------
     261! NAME
     262!     solve_tridiag
     263!
     264! DESCRIPTION
     265!     Solve a tridiagonal system Ax = d using the Thomas' algorithm
     266!     a: sub-diagonal
     267!     b: main diagonal
     268!     c: super-diagonal
     269!     d: right-hand side
     270!     x: solution
     271!
     272! AUTHORS & DATE
     273!     JB Clement, 2025
     274!
     275! NOTES
     276!
     277!-----------------------------------------------------------------------
     278
     279! DECLARATION
     280! -----------
     281implicit none
     282
     283! ARGUMENTS
     284! ---------
     285integer,                intent(in)  :: n
     286real, dimension(n),     intent(in)  :: b, d
     287real, dimension(n - 1), intent(in)  :: a, c
     288real, dimension(n),     intent(out) :: x
     289integer,                intent(out) :: error
     290
     291! LOCAL VARIABLES
     292! ---------------
    230293integer            :: i
    231294real               :: m
    232295real, dimension(n) :: cp, dp
    233296
     297! CODE
     298! ----
    234299! Check stability: diagonally dominant condition
    235300error = 0
    … …  
    269334
    270335END SUBROUTINE solve_tridiag
    271 
    272 !=======================================================================
    273 
     336!=======================================================================
     337
     338!=======================================================================
    274339SUBROUTINE solve_steady_heat(n,z,mz,kappa,mkappa,T_left,q_right,T)
    275 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    276 !!!
    277 !!! Purpose: Solve 1D steady-state heat equation with space-dependent diffusivity
    278 !!!          Left boudary condition : prescribed temperature T_left
    279 !!!          Right boudary condition: prescribed thermal flux q_right
    280 !!!
    281 !!!          z     : grid points
    282 !!!          mz    : mid-grid points
    283 !!!          kappa : thermal diffusivity at grid points
    284 !!!          mkappa: thermal diffusivity at mid-grid points
    285 !!!
    286 !!! Author: JBC
    287 !!!
    288 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    289 
     340!-----------------------------------------------------------------------
     341! NAME
     342!     solve_steady_heat
     343!
     344! DESCRIPTION
     345!     Solve 1D steady-state heat equation with space-dependent thermal diffusivity.
     346!     Uses Thomas algorithm to solve tridiagonal system.
     347!     Left boundary: prescribed temperature (Dirichlet).
     348!     Right boundary: prescribed thermal flux (Neumann).
     349!
     350! AUTHORS & DATE
     351!     JB Clement, 2025
     352!
     353! NOTES
     354!     Grid points at z, mid-grid points at mz.
     355!     Thermal diffusivity specified at both grids (kappa, mkappa).
     356!-----------------------------------------------------------------------
     357
     358! DEPENDENCIES
     359! ------------
    290360use stop_pem, only: stop_clean
    291361
    292 implicit none
    293 
    294 ! Inputs
    295 !-------
    296 integer,                intent(in) :: n
    297 real, dimension(n),     intent(in) :: z, kappa
    298 real, dimension(n - 1), intent(in) :: mz, mkappa
    299 real,                   intent(in) :: T_left, q_right
    300 ! Outputs
    301 !--------
    302 real, dimension(n), intent(out) :: T
    303 ! Local viariables
    304 !-----------------
     362! DECLARATION
     363! -----------
     364implicit none
     365
     366! ARGUMENTS
     367! ---------
     368integer,                intent(in)  :: n
     369real, dimension(n),     intent(in)  :: z, kappa
     370real, dimension(n - 1), intent(in)  :: mz, mkappa
     371real,                   intent(in)  :: T_left, q_right
     372real, dimension(n),     intent(out) :: T
     373
     374! LOCAL VARIABLES
     375! ---------------
    305376integer                :: i, error
    306377real, dimension(n)     :: b, d
    307378real, dimension(n - 1) :: a, c
    308379
     380! CODE
     381! ----
    309382! Initialization
    310383a = 0.; b = 0.; c = 0.; d = 0.
    … …  
    331404
    332405END SUBROUTINE solve_steady_heat
     406!=======================================================================
    333407
    334408end module math_toolkit

  • trunk/LMDZ.COMMON/libf/evolution/metamorphism.F90

    r3989 r3991  
    11MODULE metamorphism
    2 
    3 implicit none
    4 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     metamorphism
     5!
     6! DESCRIPTION
     7!     Module for managing frost variables.
     8!
     9! AUTHORS & DATE
     10!     JB Clement, 12/2025
     11!
     12! NOTES
     13!
     14!-----------------------------------------------------------------------
     15
     16! DECLARATION
     17! -----------
     18implicit none
     19
     20! MODULE VARIABLES
     21! ----------------
    522! Different types of frost retained by the PEM to give back to the PCM at the end
    623real, dimension(:,:), allocatable :: h2o_frost4PCM
    … …  
    1027integer :: iPCM_h2ofrost, iPCM_co2frost
    1128
    12 !=======================================================================
    1329contains
    14 !=======================================================================
    15 
     30!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     31
     32!=======================================================================
    1633SUBROUTINE ini_frost_id(nqtot,noms)
    17 
    18 implicit none
    19 
    20 ! Arguments
    21 !----------
     34!-----------------------------------------------------------------------
     35! NAME
     36!     ini_frost_id
     37!
     38! DESCRIPTION
     39!     Initialize frost indices from PCM variable names.
     40!
     41! AUTHORS & DATE
     42!     JB Clement, 12/2025
     43!
     44! NOTES
     45!
     46!-----------------------------------------------------------------------
     47
     48! DECLARATION
     49! -----------
     50implicit none
     51
     52! ARGUMENTS
     53! ---------
    2254integer,                        intent(in) :: nqtot
    2355character(*), dimension(nqtot), intent(in) :: noms
    2456
    25 ! Local variables
    26 !----------------
     57! LOCAL VARIABLES
     58! ---------------
    2759integer :: i
    2860
    29 ! Code
    30 !-----
     61! CODE
     62! ----
    3163! Initialization
    3264iPCM_h2ofrost = -1
    … …  
    4678!=======================================================================
    4779
     80!=======================================================================
    4881SUBROUTINE compute_frost(ngrid,nslope,h2ofrost_PCM,min_h2ofrost,co2frost_PCM,min_co2frost)
    49 
    50 implicit none
    51 
    52 ! Arguments
    53 !----------
     82!-----------------------------------------------------------------------
     83! NAME
     84!     compute_frost
     85!
     86! DESCRIPTION
     87!     Compute the frost to give back to the PCM.
     88!
     89! AUTHORS & DATE
     90!     JB Clement, 12/2025
     91!
     92! NOTES
     93!     Frost for the PEM is the extra part of the PCM frost above the
     94!     yearly minimum.
     95!-----------------------------------------------------------------------
     96
     97! DECLARATION
     98! -----------
     99implicit none
     100
     101! ARGUMENTS
     102! ---------
    54103integer,                       intent(in) :: ngrid, nslope
    55104real, dimension(ngrid,nslope), intent(in) :: h2ofrost_PCM, min_h2ofrost, co2frost_PCM, min_co2frost
    56105
    57 ! Local variables
    58 !----------------
    59 
    60 ! Code
    61 !-----
     106! CODE
     107! ----
    62108write(*,*) '> Computing frost to give back to the PCM'
    63109
    … …  
    76122!=======================================================================
    77123
     124!=======================================================================
    78125SUBROUTINE set_frost4PCM(PCMfrost)
    79 
    80 implicit none
    81 
    82 ! Arguments
    83 !----------
     126!-----------------------------------------------------------------------
     127! NAME
     128!     set_frost4PCM
     129!
     130! DESCRIPTION
     131!     Reconstruct frost for the PCM from PEM computations.
     132!
     133! AUTHORS & DATE
     134!     JB Clement, 12/2025
     135!
     136! NOTES
     137!
     138!-----------------------------------------------------------------------
     139
     140! DECLARATION
     141! -----------
     142implicit none
     143
     144! ARGUMENTS
     145! ---------
    84146real, dimension(:,:,:), intent(inout) :: PCMfrost
    85147
    86 ! Local variables
    87 !----------------
    88 
    89 ! Code
    90 !-----
     148! CODE
     149! ----
    91150write(*,*) '> Reconstructing frost for the PCM'
    92151PCMfrost(:,iPCM_h2ofrost,:) = h2o_frost4PCM(:,:)
    … …  
    99158!=======================================================================
    100159
     160!=======================================================================
    101161SUBROUTINE ini_frost(ngrid,nslope)
    102 
    103 implicit none
    104 
    105 ! Arguments
    106 !----------
     162!-----------------------------------------------------------------------
     163! NAME
     164!     ini_frost
     165!
     166! DESCRIPTION
     167!     Initialize frost arrays.
     168!
     169! AUTHORS & DATE
     170!     JB Clement, 12/2025
     171!
     172! NOTES
     173!
     174!-----------------------------------------------------------------------
     175
     176! DECLARATION
     177! -----------
     178implicit none
     179
     180! ARGUMENTS
     181! ---------
    107182integer, intent(in) :: ngrid, nslope
    108183
    109 ! Local variables
    110 !----------------
    111 
    112 ! Code
    113 !-----
     184! CODE
     185! ----
    114186if (.not. allocated(h2o_frost4PCM)) allocate(h2o_frost4PCM(ngrid,nslope))
    115187if (.not. allocated(co2_frost4PCM)) allocate(co2_frost4PCM(ngrid,nslope))
    … …  
    118190!=======================================================================
    119191
     192!=======================================================================
    120193SUBROUTINE end_frost()
    121 
    122 implicit none
    123 
    124 ! Arguments
    125 !----------
    126 
    127 ! Local variables
    128 !----------------
    129 
    130 ! Code
    131 !-----
     194!-----------------------------------------------------------------------
     195! NAME
     196!     end_frost
     197!
     198! DESCRIPTION
     199!     Deallocate frost arrays.
     200!
     201! AUTHORS & DATE
     202!     JB Clement, 12/2025
     203!
     204! NOTES
     205!
     206!-----------------------------------------------------------------------
     207
     208! DECLARATION
     209! -----------
     210implicit none
     211
     212! CODE
     213! ----
    132214if (allocated(h2o_frost4PCM)) deallocate(h2o_frost4PCM)
    133215if (allocated(co2_frost4PCM)) deallocate(co2_frost4PCM)
    134216
    135217END SUBROUTINE end_frost
     218!=======================================================================
    136219
    137220END MODULE metamorphism

  • trunk/LMDZ.COMMON/libf/evolution/orbit.F90

    r3989 r3991  
    11MODULE orbit
    2 
    3 !=======================================================================
    4 !
    5 ! Purpose: Compute the maximum number of iteration for PEM based on the
    6 ! stopping criterion given by the orbital parameters
    7 !
    8 ! Author: RV, JBC
    9 !=======================================================================
    10 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     orbit
     5! DESCRIPTION
     6!     Compute orbital-parameter-based stopping criteria and update
     7!     planetary orbit parameters.
     8! AUTHORS & DATE
     9!     R. Vandemeulebrouck
     10!     JB Clement, 2023-2025
     11! NOTES
     12!
     13!-----------------------------------------------------------------------
     14
     15! DECLARATION
     16! -----------
    1117implicit none
    1218
     19! MODULE VARIABLES
     20! ----------------
    1321real, dimension(:), allocatable :: year_lask  ! year before present from Laskar et al.
    14 real, dimension(:), allocatable :: obl_lask   ! obliquity    [deg]
     22real, dimension(:), allocatable :: obl_lask   ! obliquity [deg]
    1523real, dimension(:), allocatable :: ecc_lask   ! eccentricity
    1624real, dimension(:), allocatable :: lsp_lask   ! ls perihelie [deg]
    1725integer                         :: iyear_lask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
    1826
    19 !=======================================================================
    2027contains
    21 !=======================================================================
    22 
     28!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     29
     30!=======================================================================
    2331SUBROUTINE read_orbitpm(Year)
    24 
     32!-----------------------------------------------------------------------
     33! NAME
     34!     read_orbitpm
     35!
     36! DESCRIPTION
     37!     Read Laskar orbital file and prepare arrays; locate index for
     38!     closest lower year relative to current simulation year.
     39!
     40! AUTHORS & DATE
     41!     R. Vandemeulebrouck
     42!     JB Clement, 2023-2025
     43!
     44! NOTES
     45!     Converts kilo Earth years to Martian years using 'convert_years'.
     46!-----------------------------------------------------------------------
     47
     48! DEPENDENCIES
     49! ------------
    2550use evolution, only: convert_years
    2651
     52! DECLARATION
     53! -----------
    2754implicit none
    2855
     56! ARGUMENTS
     57! ---------
    2958real, intent(in) :: Year ! Martian year of the simulation
    3059
    31 integer :: n ! number of lines in Laskar files
     60! LOCAL VARIABLES
     61! ---------------
     62integer :: n ! Number of lines in Laskar files
    3263integer :: ierr, i
    3364
     65! CODE
     66! ----
    3467n = 0
    3568open(1,file = 'obl_ecc_lsp.asc',status = 'old',action = 'read')
    … …  
    5689
    5790END SUBROUTINE read_orbitpm
    58 
    59 !=======================================================================
    60 
     91!=======================================================================
     92
     93!=======================================================================
    6194SUBROUTINE end_orbit
    62 
     95!-----------------------------------------------------------------------
     96! NAME
     97!     end_orbit
     98!
     99! DESCRIPTION
     100!     Deallocate orbital data arrays.
     101!
     102! AUTHORS & DATE
     103!     R. Vandemeulebrouck
     104!     JB Clement, 2023-2025
     105!
     106! NOTES
     107!
     108!-----------------------------------------------------------------------
     109
     110! DECLARATION
     111! -----------
    63112implicit none
    64113
     114! CODE
     115! ----
    65116if (allocated(year_lask)) deallocate(year_lask)
    66 if (allocated(obl_lask)) deallocate(obl_lask)
    67 if (allocated(ecc_lask)) deallocate(ecc_lask)
    68 if (allocated(lsp_lask)) deallocate(lsp_lask)
     117if (allocated(obl_lask))  deallocate(obl_lask)
     118if (allocated(ecc_lask))  deallocate(ecc_lask)
     119if (allocated(lsp_lask))  deallocate(lsp_lask)
    69120
    70121END SUBROUTINE end_orbit
    71122!=======================================================================
    72123
     124!=======================================================================
    73125SUBROUTINE orbit_param_criterion(i_myear,nyears_maxorb)
    74 
     126!-----------------------------------------------------------------------
     127! NAME
     128!     orbit_param_criterion
     129!
     130! DESCRIPTION
     131!     Determine maximum PEM iterations before orbital params change
     132!     beyond admissible thresholds.
     133!
     134! AUTHORS & DATE
     135!     R. Vandemeulebrouck
     136!     JB Clement, 2023-2025
     137!
     138! NOTES
     139!
     140!-----------------------------------------------------------------------
    75141#ifdef CPP_IOIPSL
    76142    use IOIPSL,          only: getin
    … …  
    83149use evolution,      only: year_bp_ini, var_obl, var_ecc, var_lsp, convert_years
    84150
     151! DECLARATION
     152! -----------
    85153implicit none
    86154
    87 !--------------------------------------------------------
    88 ! Input Variables
    89 !--------------------------------------------------------
    90 real, intent(in) :: i_myear ! Martian year of the simulation
    91 
    92 !--------------------------------------------------------
    93 ! Output Variables
    94 !--------------------------------------------------------
     155! ARGUMENTS
     156! ---------
     157real, intent(in)  :: i_myear       ! Martian year of the simulation
    95158real, intent(out) :: nyears_maxorb ! Maximum number of iteration of the PEM before orb changes too much
    96159
    97 !--------------------------------------------------------
    98 ! Local variables
    99 !--------------------------------------------------------
     160! LOCAL VARIABLES
     161! ---------------
    100162real                            :: Year                                           ! Year of the simulation
    101163real                            :: Lsp                                            ! Ls perihelion
    … …  
    103165real                            :: max_obl, max_ecc, max_lsp                      ! Maximum value of orbit param given the admissible change
    104166real                            :: min_obl, min_ecc, min_lsp                      ! Minimum value of orbit param given the admissible change
    105 real                            :: nyears_maxobl, nyears_maxecc, nyears_maxlsp       ! Maximum year iteration before reaching an inadmissible value of orbit param
    106 integer                         :: i                                          ! Loop variable
    107 real                            :: xa, xb, ya, yb                                 ! Variables for interpolation
    108 logical                         :: found_obl, found_ecc, found_lsp                ! Flag variables for orbital parameters
    109 real                            :: lsp_corr                                       ! Lsp corrected if the "modulo is crossed"
    110 real                            :: delta                                          ! Intermediate variable
    111 real, dimension(:), allocatable :: lsplask_unwrap                                 ! Unwrapped sequence of Lsp from Laskar's file
    112 
     167real                            :: nyears_maxobl, nyears_maxecc, nyears_maxlsp    ! Maximum year iteration before reaching an inadmissible value of orbit param
     168real                            :: xa, xb, ya, yb                  ! Variables for interpolation
     169logical                         :: found_obl, found_ecc, found_lsp ! Flag variables for orbital parameters
     170real                            :: lsp_corr                        ! Lsp corrected if the "modulo is crossed"
     171real                            :: delta                           ! Intermediate variable
     172real, dimension(:), allocatable :: lsplask_unwrap                  ! Unwrapped sequence of Lsp from Laskar's file
     173integer                         :: i
     174
     175! CODE
     176! ----
    113177! **********************************************************************
    114178! 0. Initializations
    … …  
    246310
    247311END SUBROUTINE orbit_param_criterion
    248 
    249 !***********************************************************************
    250 ! Purpose: Recompute orbit parameters based on values by Laskar et al.
     312!=======================================================================
     313
     314!=======================================================================
    251315SUBROUTINE set_new_orbitpm(i_myear,i_myear_leg)
     316!-----------------------------------------------------------------------
     317! NAME
     318!     set_new_orbitpm
     319!
     320! DESCRIPTION
     321!     Recompute orbit parameters from Laskar values at a new year and
     322!     update planetary constants; handles Lsp modulo crossings.
     323!
     324! AUTHORS & DATE
     325!     R. Vandemeulebrouck
     326!     JB Clement, 2023-2025
     327!
     328! NOTES
     329!
     330!-----------------------------------------------------------------------
    252331
    253332use evolution,        only: year_bp_ini, var_obl, var_ecc, var_lsp
    … …  
    256335use call_dayperi_mod, only: call_dayperi
    257336
     337! DECLARATION
     338! -----------
    258339implicit none
    259340
    260 !--------------------------------------------------------
    261 ! Input Variables
    262 !--------------------------------------------------------
     341! ARGUMENTS
     342! ---------
    263343real, intent(in) :: i_myear     ! Number of simulated Martian years
    264344real, intent(in) :: i_myear_leg ! Number of iterations of the PEM
    265345
    266 !--------------------------------------------------------
    267 ! Output Variables
    268 !--------------------------------------------------------
    269 
    270 !--------------------------------------------------------
    271 ! Local variables
    272 !--------------------------------------------------------
     346! LOCAL VARIABLES
     347! ---------------
    273348real    :: Year           ! Year of the simulation
    274349real    :: Lsp            ! Ls perihelion
    275 integer :: i          ! Loop variable
    276350real    :: halfaxe        ! Million of km
     351integer :: i
    277352real    :: unitastr
    278353real    :: xa, xb, ya, yb ! Variables for interpolation
    279354logical :: found_year     ! Flag variable
    280355
     356! CODE
     357! ----
    281358! **********************************************************************
    282359! 0. Initializations
    … …  
    353430
    354431END SUBROUTINE set_new_orbitpm
     432!=======================================================================
    355433
    356434END MODULE orbit

  • trunk/LMDZ.COMMON/libf/evolution/outputs.F90

    r3989 r3991  
    11MODULE outputs
    2 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     outputs
     5!
     6! DESCRIPTION
     7!     Tools to write PEM diagnostic outputs.
     8!
     9! AUTHORS & DATE
     10!     LMDZ team
     11!     E. Leconte, 2010
     12!     F. Forget, 2011
     13!     JB Clement, 2023–2025
     14!
     15! NOTES
     16!     Uses NetCDF low-level API and supports parallel runs.
     17!-----------------------------------------------------------------------
     18
     19! DECLARATION
     20! -----------
    321implicit none
    422
     23! MODULE VARIABLES
     24! ----------------
    525integer :: output_rate ! Output rate
    626
    7 !=======================================================================
    827contains
    9 !=======================================================================
     28!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    1029
    1130SUBROUTINE write_diagpem(ngrid,nom,titre,unite,dim,px)
    … …  
    4867!      dim : dimension de px : 0, 1, 2, ou 3 dimensions
    4968!
    50 !=================================================================
     69!-----------------------------------------------------------------------
     70! NAME
     71!     write_diagpem
     72!
     73! DESCRIPTION
     74!     Write selected diagnostic variables to NetCDF file 'diagpem.nc'.
     75!     Supports 3D, 2D, 1D (column), and 0D (scalar) variables.
     76!
     77! AUTHORS & DATE
     78!     E. Leconte, 2010
     79!     F. Forget, 2011
     80!     JB Clement, 2023–2025
     81!
     82! NOTES
     83!     Output cadence is controlled by 'output_rate' from run.def. Parallel-safe
     84!     with master task handling NetCDF I/O. Can use 'diagpem.def' to select vars.
     85!-----------------------------------------------------------------------
     86
     87! DEPENDENCIES
     88! ------------
    5189use surfdat_h,          only: phisfi
    5290use geometry_mod,       only: cell_area
    … …  
    5492use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat, nbp_lev, grid_type, unstructured
    5593
     94! DECLARATION
     95! -----------
    5696implicit none
    5797
    … …  
    5999include "netcdf.inc"
    60100
    61 ! Arguments on input:
     101! ARGUMENTS
     102! ---------
    62103integer,                        intent(in) :: ngrid
    63104character(len=*),               intent(in) :: nom, titre, unite
    … …  
    65106real, dimension(ngrid,nbp_lev), intent(in) :: px
    66107
    67 ! Local variables:
     108! LOCAL VARIABLES
     109! ---------------
    68110real*4, dimension(nbp_lon + 1,nbp_lat,nbp_lev) :: dx3    ! to store a 3D data set
    69111real*4, dimension(nbp_lon + 1,nbp_lat)         :: dx2    ! to store a 2D (surface) data set
    … …  
    117159#endif
    118160
     161! CODE
     162! ----
    119163if (grid_type == unstructured) return
    120164
    … …  
    580624
    581625END SUBROUTINE write_diagpem
    582 
    583626!=================================================================
    584627
     628!=================================================================
    585629SUBROUTINE write_diagsoilpem(ngrid,name,title,units,dimpx,px)
     630!-----------------------------------------------------------------------
     631! NAME
     632!     write_diagsoilpem
     633!
     634! DESCRIPTION
     635!     Write soil-related diagnostic variables to 'diagsoilpem.nc'.
     636!     Supports 3D (lon,lat,depth), 2D (lon,lat), and 0D scalars.
     637!
     638! AUTHORS & DATE
     639!     E. Leconte, 2010
     640!     JB Clement, 2023–2025
     641!
     642! NOTES
     643!     Output cadence uses 'output_rate'. Only lon-lat (or 1D) grids supported.
     644!-----------------------------------------------------------------------
    586645
    587646! Write variable 'name' to NetCDF file 'diagsoilpem.nc'.
    … …  
    596655! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
    597656
     657! DEPENDENCIES
     658! ------------
    598659use soil,               only: mlayer_PEM, nsoilmx_PEM, inertiedat_PEM
    599660use geometry_mod,       only: cell_area
    … …  
    603664use iniwritesoil_mod,   only: iniwritesoil
    604665
     666! DECLARATION
     667! -----------
    605668implicit none
    606669
    607670include"netcdf.inc"
    608671
    609 ! Arguments:
    610 integer,intent(in) :: ngrid ! number of (horizontal) points of physics grid
     672! ARGUMENTS
     673! ---------
     674integer,                            intent(in) :: ngrid ! number of (horizontal) points of physics grid
    611675! i.e. ngrid = 2+(jjm-1)*iim - 1/jjm
    612 character(len=*),intent(in) :: name ! 'name' of the variable
    613 character(len=*),intent(in) :: title ! 'long_name' attribute of the variable
    614 character(len=*),intent(in) :: units ! 'units' attribute of the variable
    615 integer,intent(in) :: dimpx ! dimension of the variable (3,2 or 0)
    616 real,dimension(ngrid,nsoilmx_PEM),intent(in) :: px ! variable
    617 
    618 ! Local variables:
     676character(len=*),                   intent(in) :: name  ! 'name' of the variable
     677character(len=*),                   intent(in) :: title ! 'long_name' attribute of the variable
     678character(len=*),                   intent(in) :: units ! 'units' attribute of the variable
     679integer,                            intent(in) :: dimpx ! dimension of the variable (3,2 or 0)
     680real, dimension(ngrid,nsoilmx_PEM), intent(in) :: px    ! variable
     681
     682! LOCAL VARIABLES
     683! ---------------
    619684real*4,dimension(nbp_lon+1,nbp_lat,nsoilmx_PEM) :: data3 ! to store 3D data
    620685real*4,dimension(nbp_lon+1,nbp_lat) :: data2 ! to store 2D data
    … …  
    657722#endif
    658723
     724! CODE
     725! ----
    659726! 0. Do we ouput a diagsoilpem.nc file? If not just bail out now.
    660727
    … …  
    9891056
    9901057END SUBROUTINE write_diagsoilpem
     1058!=================================================================
    9911059
    9921060END MODULE outputs

  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3989 r3991  
    2626
    2727PROGRAM pem
    28 
     28!-----------------------------------------------------------------------
     29! NAME
     30!     pem
     31!
     32! DESCRIPTION
     33!     Main entry point for the Planetary Evolution Model (PEM).
     34!
     35! AUTHORS & DATE
     36!     R. Vandemeulebrouck, 2022
     37!     L. Lange, 2022
     38!     JB Clement, 2023-2025
     39!
     40! NOTES
     41!     Drives initialization, run loop, and outputs for PEM/PCM coupling.
     42!-----------------------------------------------------------------------
     43
     44! DEPENDENCIES
     45! ------------
    2946use phyetat0_mod,          only: phyetat0
    3047use phyredem,              only: physdem0, physdem1
    … …  
    99116#endif
    100117
     118! DECLARATION
     119! -----------
    101120implicit none
    102121
    … …  
    106125include "iniprint.h"
    107126
     127! LOCAL VARIABLES
     128! ---------------
    108129! Same variable names as in the PCM
    109130integer, parameter :: ngridmx = 2 + (jjm - 1)*iim - 1/jjm
    … …  
    256277#endif
    257278
    258 ! Loop variables
    259279integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
    260280real    :: totmass_ini
    261281logical :: num_str
    262282
    263 ! Code
    264 !-----
     283! CODE
     284! ----
    265285! Elapsed time with system clock
    266286call system_clock(count_rate = cr)
    … …  
    682702        call stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopCrit)
    683703        call balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
    684    
     704
    685705        do islope = 1,nslope
    686706            do ig = 1,ngrid

  • trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90

    r3989 r3991  
    11MODULE pemetat0_mod
    2 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     pemetat0_mod
     5!
     6! DESCRIPTION
     7!     Initialize PEM state from start files and PCM/XIOS data.
     8!
     9! AUTHORS & DATE
     10!     L. Lange
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
     16
     17! DECLARATION
     18! -----------
    319implicit none
    420
     21contains
     22!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     23
    524!=======================================================================
    6 contains
    7 !=======================================================================
    8 
    925SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    1026                    tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,d_h2oice,d_co2ice,co2_ice,h2o_ice,                         &
    1127                    watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,m_h2o_regolith_phys,deltam_h2o_regolith_phys,layerings_map)
    12 
     28!-----------------------------------------------------------------------
     29! NAME
     30!     pemetat0
     31!
     32! DESCRIPTION
     33!     Read PEM start file (if present) and initialize model state.
     34!     Reconstructs when missing.
     35!
     36! AUTHORS & DATE
     37!     L. Lange
     38!     JB Clement, 2023-2025
     39!
     40! NOTES
     41!     Order of operations (must be respected):
     42!       0) Previous-year context and ice fields
     43!       1) Thermal inertia
     44!       2) Soil temperature
     45!       3) Ice table (equilibrium or dynamic)
     46!       4) Regolith adsorption (CO2 & H2O)
     47!-----------------------------------------------------------------------
     48
     49! DEPENDENCIES
     50! ------------
    1351use iostart_pem,           only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
    1452use soil,                  only: do_soil, layer_PEM, fluxgeo, inertiedat_PEM, h2oice_huge_ini, depth_breccia, depth_bedrock, index_breccia, index_bedrock
    … …  
    2866use tracers,               only: mmol, iPCM_qh2o
    2967
     68! DECLARATION
     69! -----------
    3070implicit none
    3171
    32 character(*),                   intent(in) :: filename      ! name of the startfi_PEM.nc
    33 integer,                        intent(in) :: ngrid         ! # of physical grid points
    34 integer,                        intent(in) :: nsoil_PCM     ! # of vertical grid points in the PCM
    35 integer,                        intent(in) :: nsoil_PEM     ! # of vertical grid points in the PEM
    36 integer,                        intent(in) :: nslope        ! # of sub-grid slopes
    37 integer,                        intent(in) :: timelen       ! # time samples
    38 real,                           intent(in) :: timestep      ! time step [s]
    39 real,                           intent(in) :: ps_avg_global ! mean average pressure on the planet [Pa]
    40 real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1 ! surface temperature at the first year of PCM call [K]
    41 real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2 ! surface temperature at the second year of PCM call [K]
    42 real, dimension(ngrid,timelen), intent(in) :: q_co2         ! MMR tracer co2 [kg/kg]
    43 real, dimension(ngrid,timelen), intent(in) :: q_h2o         ! MMR tracer h2o [kg/kg]
    44 real, dimension(ngrid,timelen), intent(in) :: ps_timeseries ! surface pressure [Pa]
    45 real, dimension(ngrid,nslope),  intent(in) :: d_h2oice      ! tendencies on h2o ice
    46 real, dimension(ngrid,nslope),  intent(in) :: d_co2ice      ! tendencies on co2 ice
    47 real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg ! surface water ice density, yearly averaged (kg/m^3)
    48 ! outputs
    49 real, dimension(ngrid),                  intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
    50 real, dimension(ngrid),                  intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
    51 real, dimension(ngrid,nslope),           intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
    52 real, dimension(ngrid,nslope),           intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
     72! ARGUMENTS
     73! ---------
     74character(*),                            intent(in)    :: filename      ! Name of the startfi_PEM.nc
     75integer,                                 intent(in)    :: ngrid         ! # of physical grid points
     76integer,                                 intent(in)    :: nsoil_PCM     ! # of vertical grid points in the PCM
     77integer,                                 intent(in)    :: nsoil_PEM     ! # of vertical grid points in the PEM
     78integer,                                 intent(in)    :: nslope        ! # of sub-grid slopes
     79integer,                                 intent(in)    :: timelen       ! # time samples
     80real,                                    intent(in)    :: timestep      ! Time step [s]
     81real,                                    intent(in)    :: ps_avg_global ! Mean average pressure on the planet [Pa]
     82real, dimension(ngrid,nslope),           intent(in)    :: tsurf_avg_yr1 ! Surface temperature at the first year of PCM call [K]
     83real, dimension(ngrid,nslope),           intent(in)    :: tsurf_avg_yr2 ! Surface temperature at the second year of PCM call [K]
     84real, dimension(ngrid,timelen),          intent(in)    :: q_co2         ! MMR tracer co2 [kg/kg]
     85real, dimension(ngrid,timelen),          intent(in)    :: q_h2o         ! MMR tracer h2o [kg/kg]
     86real, dimension(ngrid,timelen),          intent(in)    :: ps_timeseries  ! Surface pressure [Pa]
     87real, dimension(ngrid,nslope),           intent(in)    :: d_h2oice      ! Tendencies on h2o ice
     88real, dimension(ngrid,nslope),           intent(in)    :: d_co2ice      ! Tendencies on co2 ice
     89real, dimension(ngrid,nslope),           intent(in)    :: watersurf_avg ! Surface water ice density, yearly averaged (kg/m^3)
     90real, dimension(ngrid),                  intent(out)   :: deltam_co2_regolith_phys ! Mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
     91real, dimension(ngrid),                  intent(out)   :: deltam_h2o_regolith_phys ! Mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
     92real, dimension(ngrid,nslope),           intent(out)   :: h2o_ice                  ! H2O ice amount [kg/m^2]
     93real, dimension(ngrid,nslope),           intent(out)   :: co2_ice                  ! CO2 ice amount [kg/m^2]
    5394type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map       ! Layerings
    54 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
    55 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
     95real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM              ! Soil (mid-layer) thermal inertia in the PEM grid [SI]
     96real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM           ! Soil (mid-layer) temperature [K]
    5697real, dimension(ngrid,nslope),           intent(inout) :: icetable_depth      ! Ice table depth [m]
    5798real, dimension(ngrid,nslope),           intent(inout) :: icetable_thickness  ! Ice table thickness [m]
    5899real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling     ! Subsurface ice pore filling [m3/m3]
    59 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
    60 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
    61 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
    62 ! local
    63 real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem    ! soil temperature saved in the start [K]
    64 real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM       ! soil thermal inertia saved in the start [SI]
    65 logical                                 :: found, found2     ! check if variables are found in the start
    66 integer                                 :: iloop, ig, islope ! index for loops
     100real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! Mass of co2 adsorbed [kg/m^2]
     101real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! Mass of h2o adsorbed [kg/m^2]
     102real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg       ! Surface water ice density, yearly averaged (kg/m^3)
     103
     104! LOCAL VARIABLES
     105! ---------------
     106real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem    ! Soil temperature saved in the start [K]
     107real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM       ! Soil thermal inertia saved in the start [SI]
     108logical                                 :: found, found2     ! Check if variables are found in the start
     109integer                                 :: iloop, ig, islope ! Index for loops
    67110real                                    :: delta             ! Depth of the interface regolith-breccia, breccia -bedrock [m]
    68 character(2)                            :: num               ! intermediate string to read PEM start sloped variables
    69 logical                                 :: startpem_file     ! boolean to check if we read the startfile or not
    70 real, dimension(:,:,:,:), allocatable   :: layerings_array     ! Array for layerings
    71 
    72 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    73 !!!
    74 !!! Purpose: read start_pem. Need a specific iostart
    75 !!!
    76 !!! Order: 0. Previous year of the PEM run
    77 !!!           Ice initialization
    78 !!!        1. Thermal Inertia
    79 !!!        2. Soil Temperature
    80 !!!        3. Ice table
    81 !!!        4. Mass of CO2 & H2O adsorbed
    82 !!!
    83 !!! /!\ This order must be respected!
    84 !!! Author: LL
    85 !!!
    86 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    87 
     111character(2)                            :: num               ! Intermediate string to read PEM start sloped variables
     112logical                                 :: startpem_file     ! Boolean to check if we read the startfile or not
     113real, dimension(:,:,:,:), allocatable   :: layerings_array   ! Array for layerings
     114
     115! CODE
     116! ----
    88117!0.1 Check if the start_PEM exist.
    89118
    … …  
    516545
    517546END SUBROUTINE pemetat0
     547!========================================================================
    518548
    519549END MODULE pemetat0_mod

  • trunk/LMDZ.COMMON/libf/evolution/pemredem.F90

    r3989 r3991  
    11MODULE pemredem
     2!-----------------------------------------------------------------------
     3! NAME
     4!     pemredem
     5!
     6! DESCRIPTION
     7!     Write PEM-specific NetCDF restart files.
     8!
     9! AUTHORS & DATE
     10!     L. Lange
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!     Inspired by phyredem from the PCM. Handles time-independent and
     15!     time-dependent variables for restart functionality.
     16!-----------------------------------------------------------------------
    217
    3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4 !!!
    5 !!! Purpose: Write specific netcdf restart for the PEM
    6 !!!
    7 !!!
    8 !!! Author: LL, inspired by phyredem from the PCM
    9 !!!
    10 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    11 
     18! DECLARATION
     19! -----------
    1220implicit none
    1321
     22contains
     23!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     24
    1425!=======================================================================
    15 contains
    16 !=======================================================================
     26SUBROUTINE pemdem0(filename,lonfi,latfi,cell_area,ngrid,nslope,def_slope,subslope_dist)
     27!-----------------------------------------------------------------------
     28! NAME
     29!     pemdem0
     30!
     31! DESCRIPTION
     32!     Create physics restart file and write time-independent variables.
     33!
     34! AUTHORS & DATE
     35!     L. Lange
     36!     JB Clement, 2023-2025
     37!
     38! NOTES
     39!
     40!-----------------------------------------------------------------------
    1741
    18 SUBROUTINE pemdem0(filename,lonfi,latfi,cell_area,ngrid,nslope,def_slope,subslope_dist)
    19 
    20 ! create physics restart file and write time-independent variables
    21 use soil,    only: mlayer_PEM
     42! DEPENDENCIES
     43! ------------
     44use soil,        only: mlayer_PEM
    2245use iostart_pem, only: open_restartphy, close_restartphy, put_var, put_field, length
    2346
     47! DECLARATION
     48! -----------
    2449implicit none
    2550
     51! ARGUMENTS
     52! ---------
    2653character(*),                  intent(in) :: filename
    2754integer,                       intent(in) :: ngrid, nslope
    … …  
    3158real, dimension(ngrid,nslope), intent(in) :: subslope_dist ! undermesh statistics
    3259
    33 ! Create physics start file
     60! CODE
     61! ----
    3462call open_restartphy(filename)
    3563
    … …  
    5482
    5583END SUBROUTINE pemdem0
     84!=======================================================================
    5685
    5786!=======================================================================
    58 
    5987SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, &
    6088                   icetable_depth,icetable_thickness,ice_porefilling,m_co2_regolith,m_h2o_regolith,h2o_ice,co2_ice,layerings_map)
     89!-----------------------------------------------------------------------
     90! NAME
     91!     pemdem1
     92!
     93! DESCRIPTION
     94!     Write time-dependent variables to restart file.
     95!
     96! AUTHORS & DATE
     97!     L. Lange
     98!     JB Clement, 2023-2025
     99!
     100! NOTES
     101!
     102!-----------------------------------------------------------------------
    61103
    62 ! write time-dependent variable to restart file
     104! DEPENDENCIES
     105! ------------
    63106use iostart_pem,      only: open_restartphy, close_restartphy, put_var, put_field
    64107use soil,             only: inertiedat_PEM, do_soil
    … …  
    66109use layered_deposits, only: layering, nb_str_max, map2array, print_layering, layering_algo
    67110
     111! DECLARATION
     112! -----------
    68113implicit none
    69114
     115! ARGUMENTS
     116! ---------
    70117character(*),                            intent(in) :: filename
    71118integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope
    … …  
    80127type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map ! Layerings
    81128
     129! LOCAL VARIABLES
     130! ---------------
    82131integer                               :: islope
    83132character(2)                          :: num
    … …  
    85134real, dimension(:,:,:,:), allocatable :: layerings_array ! Array for stratification (layerings)
    86135
     136! CODE
     137! ----
    87138! Open file
    88139call open_restartphy(filename)

  • trunk/LMDZ.COMMON/libf/evolution/phys_constants.F90

    r3989 r3991  
    11MODULE phys_constants
     2!-----------------------------------------------------------------------
     3! NAME
     4!     phys_constants
     5!
     6! DESCRIPTION
     7!     Physical constants used across PEM modules.
     8!
     9! AUTHORS & DATE
     10!     JB Clement, 12/2025
     11!
     12! NOTES
     13!     Constants are initialized from NetCDF 'startfi.nc' via read_constants.
     14!-----------------------------------------------------------------------
    215
     16! DECLARATION
     17! -----------
    318implicit none
    419
     20! MODULE PARAMETERS
     21! -----------------
    522real, parameter :: pi = 4.*atan(1.) ! PI = 3.14159...
    623
    … …  
    1229real :: rcp   ! r/cpp
    1330
     31contains
     32!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     33
    1434!=======================================================================
    15 contains
    16 !=======================================================================
     35SUBROUTINE read_constants(filename)
     36!-----------------------------------------------------------------------
     37! NAME
     38!     read_constants
     39!
     40! DESCRIPTION
     41!     Read physical constants from NetCDF file 'startfi.nc' in variable
     42!     'controle' and initialize module-level constants.
     43!
     44! AUTHORS & DATE
     45!     JB Clement, 12/2025
     46!
     47! NOTES
     48!     Reads controle(5,7,8,9) for rad, g, mugaz, rcp and computes r.
     49!-----------------------------------------------------------------------
    1750
    18 SUBROUTINE read_constants(filename)
    19 ! To read the necessary constants in the variable 'controle' of the file "startfi.nc"
    20 
     51! DEPENDENCIES
     52! ------------
    2153use netcdf
    2254
     55! DECLARATION
     56! -----------
    2357implicit none
    2458
    25 ! Arguments
    26 !----------
     59! ARGUMENTS
     60! ---------
    2761character(*), intent(in) :: filename
    2862
    29 ! Local variables
    30 !----------------
     63! LOCAL VARIABLES
     64! ---------------
    3165real, dimension(:), allocatable :: controle
    3266integer :: ncid           ! File ID
    … …  
    3670integer :: ierr           ! Return codes
    3771
    38 ! Code
    39 !-----
     72! CODE
     73! ----
    4074! Open the NetCDF file
    4175ierr = nf90_open(trim(filename),NF90_NOWRITE,ncid)
    … …  
    87121
    88122END SUBROUTINE read_constants
     123!=======================================================================
    89124
    90125END MODULE phys_constants

  • trunk/LMDZ.COMMON/libf/evolution/soil.F90

    r3989 r3991  
    11MODULE soil
    2 
    3 implicit none
    4 
    5 integer, parameter                  :: nsoilmx_PEM = 68   ! number of layers in the PEM
     2!-----------------------------------------------------------------------
     3! NAME
     4!     soil
     5!
     6! DESCRIPTION
     7!     Soil state and properties for PEM.
     8!
     9! AUTHORS & DATE
     10!     L. Lange, 04/2023
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
     16
     17! DECLARATION
     18! -----------
     19implicit none
     20
     21! MODULE PARAMETERS
     22! -----------------
     23integer, parameter :: nsoilmx_PEM = 68 ! number of layers in the PEM
     24
     25! MODULE VARIABLES
     26! -----------------
    627real, allocatable, dimension(:)     :: layer_PEM          ! soil layer depths [m]
    728real, allocatable, dimension(:)     :: mlayer_PEM         ! soil mid-layer depths [m]
    … …  
    2546logical                             :: reg_thprop_dependp ! thermal properites of the regolith vary with the surface pressure
    2647
    27 !=======================================================================
    2848contains
    29 !=======================================================================
    30 
     49!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     50
     51!=======================================================================
    3152SUBROUTINE ini_soil(ngrid,nslope)
    32 
    33 implicit none
    34 
     53!-----------------------------------------------------------------------
     54! NAME
     55!     ini_soil
     56!
     57! DESCRIPTION
     58!     Allocate soil arrays based on grid dimensions.
     59!
     60! AUTHORS & DATE
     61!     L. Lange, 2023
     62!     JB Clement, 2023-2025
     63!
     64! NOTES
     65!
     66!-----------------------------------------------------------------------
     67
     68! DECLARATION
     69! -----------
     70implicit none
     71
     72! ARGUMENTS
     73! ---------
    3574integer, intent(in) :: ngrid  ! number of atmospheric columns
    3675integer, intent(in) :: nslope ! number of slope within a mesh
    3776
     77! CODE
     78! ----
    3879allocate(layer_PEM(nsoilmx_PEM))
    3980allocate(mlayer_PEM(0:nsoilmx_PEM - 1))
    … …  
    4990
    5091END SUBROUTINE ini_soil
    51 
    52 !=======================================================================
    53 
     92!=======================================================================
     93
     94!=======================================================================
    5495SUBROUTINE end_soil
    55 
    56 implicit none
    57 
     96!-----------------------------------------------------------------------
     97! NAME
     98!     end_soil
     99!
     100! DESCRIPTION
     101!     Deallocate all soil arrays.
     102!
     103! AUTHORS & DATE
     104!     L. Lange, 2023
     105!     JB Clement, 2023-2025
     106!
     107! NOTES
     108!
     109!-----------------------------------------------------------------------
     110
     111! DECLARATION
     112! -----------
     113implicit none
     114
     115! CODE
     116! ----
    58117if (allocated(layer_PEM)) deallocate(layer_PEM)
    59118if (allocated(mlayer_PEM)) deallocate(mlayer_PEM)
    … …  
    69128
    70129END SUBROUTINE end_soil
     130!=======================================================================
    71131
    72132!=======================================================================
    73133SUBROUTINE set_soil(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PCM,TI_PEM)
    74 
     134!-----------------------------------------------------------------------
     135! NAME
     136!     set_soil
     137!
     138! DESCRIPTION
     139!     Initialize soil depths and properties. Builds layer depths using
     140!     exponential then power-law distribution; interpolates thermal inertia
     141!     from PCM to PEM grid.
     142!
     143! AUTHORS & DATE
     144!     E. Millour, 07/2006
     145!     L. Lange, 2023
     146!     JB Clement, 2023-2025
     147!
     148! NOTES
     149!     Modifications:
     150!          Aug.2010 EM: use NetCDF90 to load variables (enables using
     151!                       r4 or r8 restarts independently of having compiled
     152!                       the GCM in r4 or r8)
     153!          June 2013 TN: Possibility to read files with a time axis
     154!     The various actions and variable read/initialized are:
     155!     1. Read/build layer (and midlayer) depth
     156!     2. Interpolate thermal inertia and temperature on the new grid, if necessary
     157!-----------------------------------------------------------------------
     158
     159! DEPENDENCIES
     160! ------------
    75161use iostart, only: inquire_field_ndims, get_var, get_field, inquire_field, inquire_dimension_length
    76162
    77 implicit none
    78 
    79 !=======================================================================
    80 !  Author: LL, based on work by Ehouarn Millour (07/2006)
    81 !
    82 !  Purpose: Read and/or initialise soil depths and properties
    83 !
    84 ! Modifications: Aug.2010 EM: use NetCDF90 to load variables (enables using
    85 !                             r4 or r8 restarts independently of having compiled
    86 !                             the GCM in r4 or r8)
    87 !                June 2013 TN: Possibility to read files with a time axis
    88 !
    89 !  The various actions and variable read/initialized are:
    90 !  1. Read/build layer (and midlayer) depth
    91 !  2. Interpolate thermal inertia and temperature on the new grid, if necessary
    92 !=======================================================================
    93 
    94 !=======================================================================
    95 !  arguments
    96 !  ---------
    97 !  inputs:
    98 integer,                                 intent(in) :: ngrid     ! # of horizontal grid points
    99 integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
    100 integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
    101 integer,                                 intent(in) :: nsoil_PCM ! # of soil layers in the PCM
    102 real, dimension(ngrid,nsoil_PCM,nslope), intent(in) :: TI_PCM    ! Thermal inertia  in the PCM [SI]
    103 
    104 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Thermal inertia   in the PEM [SI]
    105 !=======================================================================
    106 ! local variables:
     163! DECLARATION
     164! -----------
     165implicit none
     166
     167! ARGUMENTS
     168! ---------
     169integer,                                 intent(in)    :: ngrid     ! # of horizontal grid points
     170integer,                                 intent(in)    :: nslope    ! # of subslope wihtin the mesh
     171integer,                                 intent(in)    :: nsoil_PEM ! # of soil layers in the PEM
     172integer,                                 intent(in)    :: nsoil_PCM ! # of soil layers in the PCM
     173real, dimension(ngrid,nsoil_PCM,nslope), intent(in)    :: TI_PCM    ! Thermal inertia in the PCM [SI]
     174real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM    ! Thermal inertia in the PEM [SI]
     175
     176! LOCAL VARIABLES
     177! ---------------
    107178integer :: ig, iloop, islope ! loop counters
    108179real    :: alpha, lay1       ! coefficients for building layers
    109180real    :: index_powerlaw    ! coefficient to match the power low grid with the exponential one
    110 !=======================================================================
     181
     182! CODE
     183! ----
    111184! 1. Depth coordinate
    112185! -------------------
    … …  
    173246
    174247END SUBROUTINE set_soil
     248!=======================================================================
    175249
    176250END MODULE soil

  • trunk/LMDZ.COMMON/libf/evolution/soil_temp.F90

    r3989 r3991  
    11MODULE soil_temp
    2 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     soil_temp
     5!
     6! DESCRIPTION
     7!     Routines to compute soil temperature evolution and initialization.
     8!
     9! AUTHORS & DATE
     10!     L. Lange, 2023
     11!     JB Clement, 2023-2025
     12!
     13! NOTES
     14!
     15!-----------------------------------------------------------------------
     16
     17! DECLARATION
     18! -----------
    319implicit none
    4 !-----------------------------------------------------------------------
    5 !  Author: LL
    6 !  Purpose: This module gathers the different routines used in the PEM to compute the soil temperature evolution and initialisation.
    7 !
    8 !  Note: depths of layers and mid-layers, soil thermal inertia and
    9 !        heat capacity are commons in comsoil_PEM.h
    10 !-----------------------------------------------------------------------
     20
    1121contains
    12 !=======================================================================
     22!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    1323
    1424!=======================================================================
    1525SUBROUTINE compute_tsoil(ngrid,nsoil,firstcall,therm_i,timestep,tsurf,tsoil)
    16 
    17 use soil, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
    18 use comsoil_h,     only: volcapa
    19 
     26!-----------------------------------------------------------------------
     27! NAME
     28!     compute_tsoil
     29!
     30! DESCRIPTION
     31!     Compute soil temperature using an implicit 1st order scheme.
     32!
     33! AUTHORS & DATE
     34!     L. Lange, 2023
     35!     JB Clement, 2023-2025
     36!
     37! NOTES
     38!
     39!-----------------------------------------------------------------------
     40
     41! DEPENDENCIES
     42! ------------
     43use soil,      only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
     44use comsoil_h, only: volcapa
     45
     46! DECLARATION
     47! -----------
    2048implicit none
    2149
    22 !-----------------------------------------------------------------------
    23 !  Author: LL
    24 !  Purpose: Compute soil temperature using an implict 1st order scheme
    25 !
    26 !  Note: depths of layers and mid-layers, soil thermal inertia and
    27 !        heat capacity are commons in comsoil_PEM.h
    28 !-----------------------------------------------------------------------
    29 
    3050#include "dimensions.h"
    3151
    32 ! Inputs:
    33 ! -------
    34 integer,                      intent(in) :: ngrid     ! number of (horizontal) grid-points
    35 integer,                      intent(in) :: nsoil     ! number of soil layers
    36 logical,                      intent(in) :: firstcall ! identifier for initialization call
    37 real, dimension(ngrid,nsoil), intent(in) :: therm_i   ! thermal inertia [SI]
    38 real,                         intent(in) :: timestep  ! time step [s]
    39 real, dimension(ngrid),       intent(in) :: tsurf     ! surface temperature [K]
    40 ! Outputs:
    41 !---------
    42 real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
    43 ! Local:
    44 !-------
     52! ARGUMENTS
     53! ---------
     54integer,                      intent(in)    :: ngrid     ! number of (horizontal) grid-points
     55integer,                      intent(in)    :: nsoil     ! number of soil layers
     56logical,                      intent(in)    :: firstcall ! identifier for initialization call
     57real, dimension(ngrid,nsoil), intent(in)    :: therm_i   ! thermal inertia [SI]
     58real,                         intent(in)    :: timestep  ! time step [s]
     59real, dimension(ngrid),       intent(in)    :: tsurf     ! surface temperature [K]
     60real, dimension(ngrid,nsoil), intent(inout) :: tsoil     ! soil (mid-layer) temperature [K]
     61
     62! LOCAL VARIABLES
     63! ---------------
    4564integer :: ig, ik
    4665
     66! CODE
     67! ----
    4768! 0. Initialisations and preprocessing step
    4869 if (firstcall) then
    … …  
    121142!=======================================================================
    122143SUBROUTINE ini_tsoil_pem(ngrid,nsoil,therm_i,tsurf,tsoil)
    123 
    124 use soil, only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
    125 use comsoil_h,     only: volcapa
    126 
     144!-----------------------------------------------------------------------
     145! NAME
     146!     ini_tsoil_pem
     147!
     148! DESCRIPTION
     149!     Initialize soil with the solution of the stationary heat conduction problem.
     150!     Boundary conditions: Tsurf averaged from PCM; geothermal flux at bottom.
     151!
     152! AUTHORS & DATE
     153!     L. Lang, 2023
     154!     JB Clement, 2023-2025
     155!
     156! NOTES
     157!
     158!-----------------------------------------------------------------------
     159
     160! DEPENDENCIES
     161! ------------
     162use soil,      only: layer_PEM, mlayer_PEM, mthermdiff_PEM, thermdiff_PEM, coefq_PEM, coefd_PEM, mu_PEM, alph_PEM, beta_PEM, fluxgeo
     163use comsoil_h, only: volcapa
     164
     165! DECLARATION
     166! -----------
    127167implicit none
    128168
    129 !-----------------------------------------------------------------------
    130 !  Author: LL
    131 !  Purpose: Initialize the soil with the solution of the stationnary problem of Heat Conduction. Boundarry conditions: Tsurf averaged from the PCM; Geothermal flux at the bottom layer
    132 !
    133 !  Note: depths of layers and mid-layers, soil thermal inertia and
    134 !        heat capacity are commons in comsoil_PEM.h
    135 !-----------------------------------------------------------------------
    136 
    137169#include "dimensions.h"
    138170
    139 ! Inputs:
    140 !--------
    141 integer,                      intent(in) :: ngrid   ! number of (horizontal) grid-points
    142 integer,                      intent(in) :: nsoil   ! number of soil layers
    143 real, dimension(ngrid,nsoil), intent(in) :: therm_i ! thermal inertia [SI]
    144 real, dimension(ngrid),       intent(in) :: tsurf   ! surface temperature [K]
    145 ! Outputs:
    146 !---------
    147 real, dimension(ngrid,nsoil), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
    148 ! Local:
    149 !-------
     171! ARGUMENTS
     172! ---------
     173integer,                      intent(in)    :: ngrid   ! number of (horizontal) grid-points
     174integer,                      intent(in)    :: nsoil   ! number of soil layers
     175real, dimension(ngrid,nsoil), intent(in)    :: therm_i ! thermal inertia [SI]
     176real, dimension(ngrid),       intent(in)    :: tsurf   ! surface temperature [K]
     177real, dimension(ngrid,nsoil), intent(inout) :: tsoil   ! soil (mid-layer) temperature [K]
     178
     179! LOCAL VARIABLES
     180! ---------------
    150181integer :: ig, ik, iloop
    151182
     183! CODE
     184! ----
    152185! 0. Initialisations and preprocessing step
    153186! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities
    … …  
    220253!=======================================================================
    221254SUBROUTINE shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf,tsoil)
    222 
    223 use soil, only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM
    224 use math_toolkit,  only: solve_steady_heat
    225 
     255!-----------------------------------------------------------------------
     256! NAME
     257!     shift_tsoil2surf
     258!
     259! DESCRIPTION
     260!     Shift soil temperature profile to follow surface evolution due to ice
     261!     condensation/sublimation.
     262!
     263! AUTHORS & DATE
     264!     JB Clement, 2025
     265!
     266! NOTES
     267!
     268!-----------------------------------------------------------------------
     269
     270! DEPENDENCIES
     271! ------------
     272use soil,         only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM
     273use math_toolkit, only: solve_steady_heat
     274
     275! DECLARATION
     276! -----------
    226277implicit none
    227278
    228 !-----------------------------------------------------------------------
    229 !  Author: JBC
    230 !  Purpose: Shifting the soil temperature profile to follow the surface evolution due to ice condensation/sublimation
    231 !-----------------------------------------------------------------------
    232 ! Inputs:
    233 ! -------
    234 integer,                       intent(in) :: ngrid       ! number of (horizontal) grid-points
    235 integer,                       intent(in) :: nsoil       ! number of soil layers
    236 integer,                       intent(in) :: nslope      ! number of sub-slopes
    237 real, dimension(ngrid,nslope), intent(in) :: zshift_surf ! elevation shift for the surface [m]
    238 real, dimension(ngrid,nslope), intent(in) :: zlag        ! newly built lag thickness [m]
    239 real, dimension(ngrid,nslope), intent(in) :: tsurf       ! surface temperature [K]
    240 ! Outputs:
    241 ! --------
    242 real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
    243 ! Local:
    244 ! ------
     279! ARGUMENTS
     280! ---------
     281integer,                             intent(in)    :: ngrid       ! number of (horizontal) grid-points
     282integer,                             intent(in)    :: nsoil       ! number of soil layers
     283integer,                             intent(in)    :: nslope      ! number of sub-slopes
     284real, dimension(ngrid,nslope),       intent(in)    :: zshift_surf ! elevation shift for the surface [m]
     285real, dimension(ngrid,nslope),       intent(in)    :: zlag        ! newly built lag thickness [m]
     286real, dimension(ngrid,nslope),       intent(in)    :: tsurf       ! surface temperature [K]
     287real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil       ! soil (mid-layer) temperature [K]
     288
     289! LOCAL VARIABLES
     290! ---------------
    245291integer                             :: ig, isoil, islope, ishift, ilag, iz
    246292real                                :: a, z, zshift_surfloc, tsoil_minus, mlayer_minus
    247293real, dimension(ngrid,nsoil,nslope) :: tsoil_old
    248294
     295! CODE
     296! ----
    249297write(*,*) "> Shifting soil temperature profile to match surface evolution"
    250298tsoil_old = tsoil
    … …  
    317365!=======================================================================
    318366FUNCTION itp_tsoil(tsoil,tsurf,z) RESULT(tsoil_z)
    319 
     367!-----------------------------------------------------------------------
     368! NAME
     369!     itp_tsoil
     370!
     371! DESCRIPTION
     372!     Interpolate soil temperature profile.
     373!
     374! AUTHORS & DATE
     375!     JB Clement, 2025
     376!
     377! NOTES
     378!
     379!-----------------------------------------------------------------------
     380
     381! DEPENDENCIES
     382! ------------
    320383use soil, only: mlayer_PEM
    321384
     385! DECLARATION
     386! -----------
    322387implicit none
    323388
     389! ARGUMENTS
     390! ---------
    324391real, dimension(:), intent(in) :: tsoil
    325392real,               intent(in) :: z, tsurf
    326393
     394! LOCAL VARIABLES
     395! ---------------
    327396real    :: tsoil_z, tsoil_minus, mlayer_minus, a
    328397integer :: iz
    329398
     399! CODE
     400! ----
    330401! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs
    331402iz = 0

  • trunk/LMDZ.COMMON/libf/evolution/soil_thermal_inertia.F90

    r3989 r3991  
    11MODULE soil_thermal_inertia
    2 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     soil_thermal_inertia
     5!
     6! DESCRIPTION
     7!     Compute and update soil thermal properties based on ice content,
     8!     pressure, and cementation state.
     9!
     10! AUTHORS & DATE
     11!     L. Lange, 2023
     12!     JB Clement, 2025
     13!
     14! NOTES
     15!
     16!-----------------------------------------------------------------------
     17
     18! DECLARATION
     19! -----------
    320implicit none
    421
    5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    6 !!!
    7 !!! Purpose: Compute the soil thermal properties
    8 !!! Author:  LL, 04/2023
    9 !!!
    10 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    11 
    12 ! Constants:
     22! MODULE PARAMETERS
     23! -----------------
    1324real, parameter :: reg_inertie_thresold = 370. ! Above this thermal inertia, the regolith has too much cementation to be dependant on the pressure [J/m^2/K/s^1/2]
    1425real, parameter :: reg_inertie_minvalue = 50.  ! Minimum value of the Thermal Inertia at low pressure (Piqueux & Christensen 2009) [J/m^2/K/s^1/2]
    … …  
    1930real, parameter :: Pa2Torr = 1./133.3          ! Conversion from Pa to tor [Pa/Torr]
    2031
    21 !=======================================================================
    2232contains
    23 !=======================================================================
    24 
     33!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     34
     35!=======================================================================
    2536SUBROUTINE get_ice_TI(ispureice,pore_filling,surf_thermalinertia,ice_thermalinertia)
    26 !=======================================================================
    27 !   subject: Compute ice thermal properties
    28 !   --------
    29 !
    30 !   author: LL, 04/2023
    31 !   -------
    32 !
    33 !=======================================================================
    34 
     37!-----------------------------------------------------------------------
     38! NAME
     39!     get_ice_TI
     40!
     41! DESCRIPTION
     42!     Compute ice thermal properties based on pore filling and purity.
     43!
     44! AUTHORS & DATE
     45!     L. Lange, 2023
     46!     JB Clement, 2025
     47!
     48! NOTES
     49!     Uses Siegler et al. (2012) formula for pore-filling ice;
     50!     uses Mellon et al. (2004) value for pure water ice.
     51!-----------------------------------------------------------------------
     52
     53! DEPENDENCIES
     54! ------------
    3555use constants_marspem_mod, only: porosity
    3656
     57! DECLARATION
     58! -----------
    3759implicit none
    3860
    39 !-----------------------------------------------------------------------
    40 !=======================================================================
    41 !    Declarations :
    42 !=======================================================================
    43 !
    44 !    Input/Output
    45 !    ------------
    46 logical, intent(in) :: ispureice           ! Boolean to check if ice is massive or just pore filling
    47 real,    intent(in) :: pore_filling        ! ice pore filling in each layer (1)
    48 real,    intent(in) :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)
    49 real,    intent(out) :: ice_thermalinertia ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)
    50 
    51 !    Local Variables
    52 !    --------------
    53 REAL :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004)
    54 !=======================================================================
    55 !   Beginning of the code
    56 !=======================================================================
    57 
     61! ARGUMENTS
     62! ---------
     63logical, intent(in)  :: ispureice           ! Boolean to check if ice is massive or just pore filling
     64real,    intent(in)  :: pore_filling        ! ice pore filling in each layer (1)
     65real,    intent(in)  :: surf_thermalinertia ! surface thermal inertia (J/m^2/K/s^1/2)
     66real,    intent(out) :: ice_thermalinertia  ! Thermal inertia of ice when present in the pore (J/m^2/K/s^1/2)
     67
     68! LOCAL VARIABLES
     69! ---------------
     70real :: inertie_purewaterice = 2100 ! 2050 is better according to my computations with the formula from Siegler et al., 2012, but let's follow Mellon et al. (2004)
     71
     72! CODE
     73! ----
    5874if (ispureice) then
    5975    ice_thermalinertia = inertie_purewaterice
    … …  
    6581!=======================================================================
    6682
     83!=======================================================================
    6784SUBROUTINE update_soil_TI(ngrid,nslope,nsoil,tendencies_waterice,waterice,p_avg_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
    68 
     85!-----------------------------------------------------------------------
     86! NAME
     87!     update_soil_TI
     88!
     89! DESCRIPTION
     90!     Update soil thermal inertia based on ice table, regolith properties,
     91!     and pressure-dependent cementation.
     92!
     93! AUTHORS & DATE
     94!     L. Lange, 2023
     95!     JB Clement, 2025
     96!
     97! NOTES
     98!
     99!-----------------------------------------------------------------------
     100
     101! DEPENDENCIES
     102! ------------
    69103use comsoil_h,             only: volcapa
    70104use soil,                  only: layer_PEM, inertiedat_PEM, depth_breccia, depth_bedrock, index_breccia, index_bedrock, reg_thprop_dependp
    71105use constants_marspem_mod, only: TI_breccia, TI_bedrock, TI_regolith_avg
    72106
     107! DECLARATION
     108! -----------
    73109implicit none
    74110
    75 ! Input:
    76 integer,                             intent(in) :: ngrid, nslope, nsoil ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
    77 real,                                intent(in) :: p_avg_new            ! Global average surface pressure [Pa]
    78 real, dimension(ngrid,nslope),       intent(in) :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
    79 real, dimension(ngrid,nslope),       intent(in) :: waterice             ! Surface Water ice [kg/m^2]
    80 real, dimension(ngrid,nslope),       intent(in) :: icetable_depth       ! Ice table depth [m]
    81 real, dimension(ngrid,nslope),       intent(in) :: icetable_thickness   ! Ice table thickness [m]
    82 real, dimension(ngrid,nsoil,nslope), intent(in) :: ice_porefilling      ! Ice porefilling [m^3/m^3]
    83 logical,                             intent(in) :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
    84 
    85 ! Outputs:
    86 real, dimension(ngrid,nsoil,nslope), intent(inout) :: TI_PEM ! Soil Thermal Inertia [J/m^2/K/s^1/2]
    87 
    88 ! Local variables:
    89 integer :: ig, islope, isoil, iref, iend          ! Loop variables
    90 real, dimension(ngrid,nslope) :: regolith_inertia ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
    91 real    :: delta                                  ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
    92 real    :: ice_bottom_depth                       ! Bottom depth of the subsurface ice [m]
    93 real    :: d_part                                 ! Regolith particle size [micrometer]
    94 real    :: ice_inertia                            ! Inertia of water ice [SI]
    95 
     111! ARGUMENTS
     112! ---------
     113integer,                             intent(in)    :: ngrid, nslope, nsoil ! Shape of the arrays: physical grid, number of sub-grid slopes, number of layer in the soil
     114real,                                intent(in)    :: p_avg_new            ! Global average surface pressure [Pa]
     115real, dimension(ngrid,nslope),       intent(in)    :: tendencies_waterice  ! Tendencies on the water ice [kg/m^2/year]
     116real, dimension(ngrid,nslope),       intent(in)    :: waterice             ! Surface Water ice [kg/m^2]
     117real, dimension(ngrid,nslope),       intent(in)    :: icetable_depth       ! Ice table depth [m]
     118real, dimension(ngrid,nslope),       intent(in)    :: icetable_thickness   ! Ice table thickness [m]
     119real, dimension(ngrid,nsoil,nslope), intent(in)    :: ice_porefilling      ! Ice porefilling [m^3/m^3]
     120logical,                             intent(in)    :: icetable_equilibrium, icetable_dynamic ! Computing method for ice table
     121real, dimension(ngrid,nsoil,nslope), intent(inout) :: TI_PEM               ! Soil Thermal Inertia [J/m^2/K/s^1/2]
     122
     123! LOCAL VARIABLES
     124! ---------------
     125integer                       :: ig, islope, isoil, iref, iend ! Loop variables
     126real, dimension(ngrid,nslope) :: regolith_inertia              ! Thermal inertia of the regolith (first layer in the GCM) [J/m^2/K/s^1/2]
     127real                          :: delta                         ! Difference of depth to compute the  thermal inertia in a mixed layer [m]
     128real                          :: ice_bottom_depth              ! Bottom depth of the subsurface ice [m]
     129real                          :: d_part                        ! Regolith particle size [micrometer]
     130real                          :: ice_inertia                   ! Inertia of water ice [SI]
     131
     132! CODE
     133! ----
    96134write(*,*) "> Updating soil properties"
    97135
    … …  
    222260
    223261END SUBROUTINE update_soil_TI
     262!=======================================================================
    224263
    225264END MODULE soil_thermal_inertia

  • trunk/LMDZ.COMMON/libf/evolution/sorption.F90

    r3989 r3991  
    11MODULE sorption
    2 
    3 implicit none
    4 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     sorption
     5!
     6! DESCRIPTION
     7!     CO2 and H2O adsorption/desorption in regolith following Zent & Quinn
     8!     (1995) and Jackosky et al. (1997).
     9!
     10! AUTHORS & DATE
     11!     L. Lange, 2023
     12!     JB Clement, 2025
     13!
     14! NOTES
     15!
     16!-----------------------------------------------------------------------
     17
     18! DECLARATION
     19! -----------
     20implicit none
     21
     22! MODULE VARIABLES
     23! ----------------
    524logical                             :: do_sorption       ! Flag to compute adsorption/desorption
    625real, allocatable, dimension(:,:,:) :: co2_adsorbed_phys ! co2 that is in the regolith [kg/m^2]
    726real, allocatable, dimension(:,:,:) :: h2o_adsorbed_phys ! h2o that is in the regolith [kg/m^2]
    827
    9 !=======================================================================
    1028contains
    11 !=======================================================================
    12 
    13 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    14 !!!
    15 !!! Purpose: Compute CO2 and H2O adsorption, following the methods from Zent & Quinn 1995, Jackosky et al., 1997
    16 !!!
    17 !!! Author: LL, 01/2023
    18 !!!
    19 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     29!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     30
     31!=======================================================================
    2032SUBROUTINE ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
    21 
    22 implicit none
    23 
     33!-----------------------------------------------------------------------
     34! NAME
     35!     ini_adsorption_h_PEM
     36!
     37! DESCRIPTION
     38!     Allocate CO2 and H2O adsorption arrays.
     39!
     40! AUTHORS & DATE
     41!     L. Lange, 2023
     42!     JB Clement, 2025
     43!
     44! NOTES
     45!
     46!-----------------------------------------------------------------------
     47
     48! DECLARATION
     49! -----------
     50implicit none
     51
     52! ARGUMENTS
     53! ---------
    2454integer, intent(in) :: ngrid       ! number of atmospheric columns
    2555integer, intent(in) :: nslope      ! number of slope within a mesh
    2656integer, intent(in) :: nsoilmx_PEM ! number of soil layer in the PEM
    2757
     58! CODE
     59! ----
    2860allocate(co2_adsorbed_phys(ngrid,nsoilmx_PEM,nslope))
    2961allocate(h2o_adsorbed_phys(ngrid,nsoilmx_PEM,nslope))
    3062
    3163END SUBROUTINE ini_adsorption_h_PEM
     64!=======================================================================
    3265
    3366!=======================================================================
    3467SUBROUTINE end_adsorption_h_PEM
    35 
     68!-----------------------------------------------------------------------
     69! NAME
     70!     end_adsorption_h_PEM
     71!
     72! DESCRIPTION
     73!     Deallocate adsorption arrays.
     74!
     75! AUTHORS & DATE
     76!     L. Lange, 2023
     77!     JB Clement, 2025
     78!
     79! NOTES
     80!
     81!-----------------------------------------------------------------------
     82
     83! DECLARATION
     84! -----------
     85implicit none
     86
     87! CODE
     88! ----
    3689if (allocated(co2_adsorbed_phys)) deallocate(co2_adsorbed_phys)
    3790if (allocated(h2o_adsorbed_phys)) deallocate(h2o_adsorbed_phys)
    3891
    3992END SUBROUTINE end_adsorption_h_PEM
     93!=======================================================================
    4094
    4195!=======================================================================
    4296SUBROUTINE regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps,q_co2,q_h2o, &
    4397                               m_h2o_completesoil,delta_mh2oreg, m_co2_completesoil,delta_mco2reg)
    44 
    45 implicit none
    46 
    47 ! Inputs
    48 integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM, timelen  ! size dimension: physics x subslope x soil x timeseries
    49 real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! tendancies on the glaciers [1]
    50 real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! water ice at the surface [kg/m^2]
    51 real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! co2 ice at the surface [kg/m^2]
    52 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Soil Thermal inertia (J/K/^2/s^1/2)
    53 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Soil temperature (K)
    54 real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Average surface pressure [Pa]
    55 real,    dimension(ngrid,timelen),          intent(in) :: q_co2                              ! Mass mixing ratio of co2 in the first layer (kg/kg)
    56 real,    dimension(ngrid,timelen),          intent(in) :: q_h2o                              ! Mass mixing ratio of H2o in the first layer (kg/kg)
    57 
    58 ! Outputs
    59 real, dimension(ngrid), intent(out) :: delta_mh2oreg ! Difference density of h2o adsorbed (kg/m^3)
    60 real, dimension(ngrid), intent(out) :: delta_mco2reg ! Difference density of co2 adsorbed (kg/m^3)
    61 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil ! Density of co2 adsorbed (kg/m^3)
    62 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil ! Density of h2o adsorbed (kg/m^3)
    63 
    64 ! Local variables
     98!-----------------------------------------------------------------------
     99! NAME
     100!     regolith_adsorption
     101!
     102! DESCRIPTION
     103!     Compute both H2O and CO2 adsorption in regolith.
     104!
     105! AUTHORS & DATE
     106!     L. Lange, 2023
     107!     JB Clement, 2025
     108!
     109! NOTES
     110!
     111!-----------------------------------------------------------------------
     112
     113! DECLARATION
     114! -----------
     115implicit none
     116
     117! ARGUMENTS
     118! ---------
     119integer,                                    intent(in)    :: ngrid, nslope, nsoil_PEM, timelen  ! size dimension: physics x subslope x soil x timeseries
     120real,    dimension(ngrid,nslope),           intent(in)    :: d_h2oglaciers, d_co2glaciers ! tendancies on the glaciers [1]
     121real,    dimension(ngrid,nslope),           intent(in)    :: waterice                           ! water ice at the surface [kg/m^2]
     122real,    dimension(ngrid,nslope),           intent(in)    :: co2ice                             ! co2 ice at the surface [kg/m^2]
     123real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: TI_PEM                             ! Soil Thermal inertia (J/K/^2/s^1/2)
     124real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: tsoil_PEM                          ! Soil temperature (K)
     125real,    dimension(ngrid,timelen),          intent(in)    :: ps                                 ! Average surface pressure [Pa]
     126real,    dimension(ngrid,timelen),          intent(in)    :: q_co2                              ! Mass mixing ratio of co2 in the first layer (kg/kg)
     127real,    dimension(ngrid,timelen),          intent(in)    :: q_h2o                              ! Mass mixing ratio of H2o in the first layer (kg/kg)
     128real,    dimension(ngrid),                  intent(out)   :: delta_mh2oreg                      ! Difference density of h2o adsorbed (kg/m^3)
     129real,    dimension(ngrid),                  intent(out)   :: delta_mco2reg                      ! Difference density of co2 adsorbed (kg/m^3)
     130real,    dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil                 ! Density of co2 adsorbed (kg/m^3)
     131real,    dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil                 ! Density of h2o adsorbed (kg/m^3)
     132
     133! LOCAL VARIABLES
     134! ---------------
    65135real, dimension(ngrid,nsoil_PEM,nslope) :: theta_h2o_adsorbed ! Fraction of the pores occupied by H2O molecules
    66 ! -------------
    67 ! Compute H2O adsorption, then CO2 adsorption
     136
     137! CODE
     138! ----
    68139call regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
    69140                            theta_h2o_adsorbed,m_h2o_completesoil,delta_mh2oreg)
    … …  
    72143
    73144END SUBROUTINE regolith_adsorption
     145!=======================================================================
    74146
    75147!=======================================================================
    76148SUBROUTINE regolith_h2oadsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM, &
    77149                                  theta_h2o_adsorbed,m_h2o_completesoil,delta_mreg)
    78 
    79 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    80 !!!
    81 !!! Purpose: Compute H2O adsorption, following the methods from Jackosky et al., 1997
    82 !!!
    83 !!! Author: LL, 01/2023
    84 !!!
    85 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    86 
    87 use soil,         only: layer_PEM, index_breccia
     150!-----------------------------------------------------------------------
     151! NAME
     152!     regolith_h2oadsorption
     153!
     154! DESCRIPTION
     155!     Compute H2O adsorption in regolith following Jackosky et al. (1997).
     156!
     157! AUTHORS & DATE
     158!     L. Lange, 2023
     159!     JB Clement, 2025
     160!
     161! NOTES
     162!
     163!-----------------------------------------------------------------------
     164
     165! DEPENDENCIES
     166! ------------
     167use soil,                  only: layer_PEM, index_breccia
    88168use comslope_mod,          only: subslope_dist, def_slope_mean
    89169use vertical_layers_mod,   only: ap, bp
    … …  
    91171use phys_constants,        only: pi
    92172
    93 implicit none
    94 
    95 ! Inputs
    96 integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM,timelen   ! Size dimension
    97 real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! Tendencies on the glaciers ()
    98 real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
    99 real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! CO2 ice at the surface [kg/m^2]
    100 real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Surface pressure (Pa)
    101 real,    dimension(ngrid,timelen),          intent(in) :: q_co2                              ! Mass mixing ratio of co2 in the first layer (kg/kg)
    102 real,    dimension(ngrid,timelen),          intent(in) :: q_h2o                              ! Mass mixing ratio of H2o in the first layer (kg/kg)
    103 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Soil Thermal inertia (J/K/^2/s^1/2)
    104 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Soil temperature (K)
    105 
    106 ! Outputs
    107 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)
    108 real, dimension(ngrid,nsoil_PEM,nslope), intent(out) :: theta_h2o_adsorbed ! Fraction of the pores occupied by H2O molecules
    109 real, dimension(ngrid),                  intent(out) :: delta_mreg          ! Difference density of h2o adsorbed (kg/m^3)
    110 
    111 ! Constants
     173! DECLARATION
     174! -----------
     175implicit none
     176
     177! ARGUMENTS
     178! ---------
     179integer,                                    intent(in)    :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension
     180real,    dimension(ngrid,nslope),           intent(in)    :: d_h2oglaciers, d_co2glaciers     ! Tendencies on the glaciers ()
     181real,    dimension(ngrid,nslope),           intent(in)    :: waterice                         ! Water ice at the surface [kg/m^2]
     182real,    dimension(ngrid,nslope),           intent(in)    :: co2ice                           ! CO2 ice at the surface [kg/m^2]
     183real,    dimension(ngrid,timelen),          intent(in)    :: ps                               ! Surface pressure (Pa)
     184real,    dimension(ngrid,timelen),          intent(in)    :: q_co2                            ! Mass mixing ratio of co2 in the first layer (kg/kg)
     185real,    dimension(ngrid,timelen),          intent(in)    :: q_h2o                            ! Mass mixing ratio of H2o in the first layer (kg/kg)
     186real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: TI_PEM                           ! Soil Thermal inertia (J/K/^2/s^1/2)
     187real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: tsoil_PEM                        ! Soil temperature (K)
     188real,    dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_completesoil               ! Density of h2o adsorbed (kg/m^3)(ngrid,nsoil_PEM,nslope)
     189real,    dimension(ngrid,nsoil_PEM,nslope), intent(out)   :: theta_h2o_adsorbed               ! Fraction of the pores occupied by H2O molecules
     190real,    dimension(ngrid),                  intent(out)   :: delta_mreg                       ! Difference density of h2o adsorbed (kg/m^3)
     191
     192! LOCAL VARIABLES
     193! ---------------
     194! Constants:
    112195real :: Ko =  1.57e-8            ! Jackosky et al. 1997
    113196real :: e = 2573.9               ! Jackosky et al. 1997
    … …  
    117200real :: as = 9.48e4              ! Specific area, Zent
    118201real ::  inertie_thresold = 800. ! TI > 800 means cementation
    119 
    120 ! Local variables
    121202real,    dimension(ngrid,index_breccia,nslope) :: deltam_reg_complete     ! Difference in the mass per slope and soil layer (kg/m^3)
    122203real                                           :: K                       ! Used to compute theta
    … …  
    133214real,    dimension(:)  , allocatable           :: pvapor_avg              ! Yearly averaged
    134215
     216! CODE
     217! ----
    135218! 0. Some initializations
    136219allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pvapor(ngrid,timelen),pvapor_avg(ngrid))
    … …  
    207290
    208291END SUBROUTINE regolith_h2oadsorption
    209 
    210 !=======================================================================
    211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    212 !!!
    213 !!! Purpose: Compute CO2  following the methods from Zent & Quinn 1995
    214 !!!
    215 !!! Author: LL, 01/2023
    216 !!!
    217 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     292!=======================================================================
     293
     294!=======================================================================
    218295SUBROUTINE regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oglaciers,d_co2glaciers,waterice,co2ice,ps,q_co2,q_h2o,tsoil_PEM,TI_PEM,m_co2_completesoil,delta_mreg)
    219296
    220 use soil,                  only: layer_PEM, index_breccia, index_breccia
     297!-----------------------------------------------------------------------
     298! NAME
     299!     regolith_co2adsorption
     300!
     301! DESCRIPTION
     302!     Compute CO2 adsorption in regolith following Zent & Quinn (1995).
     303!
     304! AUTHORS & DATE
     305!     L. Lange, 2023
     306!     JB Clement, 2025
     307!
     308! NOTES
     309!
     310!-----------------------------------------------------------------------
     311
     312! DEPENDENCIES
     313! ------------
     314use soil,                  only: layer_PEM, index_breccia
    221315use comslope_mod,          only: subslope_dist, def_slope_mean
    222316use vertical_layers_mod,   only: ap, bp
    … …  
    224318use phys_constants,        only: pi
    225319
    226 implicit none
    227 
    228 ! Inputs:
    229 integer,                                    intent(in) :: ngrid, nslope, nsoil_PEM,timelen   ! Size dimension
    230 real,    dimension(ngrid,timelen),          intent(in) :: ps                                 ! Average surface pressure [Pa]
    231 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_PEM                          ! Mean Soil Temperature [K]
    232 real,    dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM                             ! Mean Thermal Inertia [USI]
    233 real,    dimension(ngrid,nslope),           intent(in) :: d_h2oglaciers, d_co2glaciers ! Tendencies on the glaciers ()
    234 real,    dimension(ngrid,timelen),          intent(in) :: q_co2, q_h2o                       ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
    235 real,    dimension(ngrid,nslope),           intent(in) :: waterice                           ! Water ice at the surface [kg/m^2]
    236 real,    dimension(ngrid,nslope),           intent(in) :: co2ice                             ! CO2 ice at the surface [kg/m^2]
    237 
    238 ! Outputs:
    239 real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil ! Density of co2 adsorbed (kg/m^3)
    240 real, dimension(ngrid), intent(out) :: delta_mreg ! Difference density of co2 adsorbed (kg/m^3)
    241 
     320! DECLARATION
     321! -----------
     322implicit none
     323
     324! ARGUMENTS
     325! ---------
     326integer,                                    intent(in)    :: ngrid, nslope, nsoil_PEM,timelen ! Size dimension
     327real,    dimension(ngrid,timelen),          intent(in)    :: ps                               ! Average surface pressure [Pa]
     328real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: tsoil_PEM                        ! Mean Soil Temperature [K]
     329real,    dimension(ngrid,nsoil_PEM,nslope), intent(in)    :: TI_PEM                           ! Mean Thermal Inertia [USI]
     330real,    dimension(ngrid,nslope),           intent(in)    :: d_h2oglaciers, d_co2glaciers     ! Tendencies on the glaciers ()
     331real,    dimension(ngrid,timelen),          intent(in)    :: q_co2, q_h2o                     ! Mass mixing ratio of co2 and h2o in the first layer (kg/kg)
     332real,    dimension(ngrid,nslope),           intent(in)    :: waterice                         ! Water ice at the surface [kg/m^2]
     333real,    dimension(ngrid,nslope),           intent(in)    :: co2ice                           ! CO2 ice at the surface [kg/m^2]
     334real,    dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_completesoil               ! Density of co2 adsorbed (kg/m^3)
     335real,    dimension(ngrid),                  intent(out)   :: delta_mreg                       ! Difference density of co2 adsorbed (kg/m^3)
     336
     337! LOCAL VARIABLES
     338! ---------------
    242339! Constants:
    243340real :: alpha = 7.512e-6        ! Zent & Quinn 1995
    … …  
    247344! real :: as = 18.9e3             ! Specific area, Buhler & Piqueux 2021
    248345real :: as = 9.48e4             ! Same as previous but from zent
    249 
    250 ! Local
    251346real                                           :: A, B                    ! Used to compute the mean mass above the surface
    252347integer                                        :: ig, islope, iloop, it   ! For loops
    … …  
    265360real,    dimension(:),   allocatable           :: pco2_avg                ! Yearly averaged
    266361
     362! CODE
     363! ----
    267364! 0. Some initializations
    268365allocate(mass_mean(ngrid,timelen),zplev_mean(ngrid,timelen),pco2(ngrid,timelen),pco2_avg(ngrid))
    … …  
    348445
    349446END SUBROUTINE regolith_co2adsorption
     447!=======================================================================
    350448
    351449END MODULE sorption

  • trunk/LMDZ.COMMON/libf/evolution/stop_pem.F90

    r3989 r3991  
    11MODULE stop_pem
     2!-----------------------------------------------------------------------
     3! NAME
     4!     stop_pem
     5!
     6! DESCRIPTION
     7!     Clean stopping utilities for PEM: close outputs and report reason.
     8!
     9! AUTHORS & DATE
     10!     JB Clement, 2025
     11!
     12! NOTES
     13!
     14!-----------------------------------------------------------------------
    215
     16! DECLARATION
     17! -----------
    318implicit none
    419
    5 !=======================================================================
    620contains
    7 !=======================================================================
     21!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    822
    923!=======================================================================
    1024SUBROUTINE stop_clean(modname,message,ierr)
    11 ! Stop the simulation cleanly, closing files and printing various comments
    12 ! Inputs: modname = name of calling program
    13 !         message = stuff to print
    14 !         ierr    = severity of situation (= 0 normal)
     25!-----------------------------------------------------------------------
     26! NAME
     27!     stop_clean
     28!
     29! DESCRIPTION
     30!     Stop simulation cleanly, closing files and printing diagnostics.
     31!
     32! AUTHORS & DATE
     33!     JB Clement, 2025
     34!
     35! NOTES
     36!     Taken from Mars PCM.
     37!-----------------------------------------------------------------------
    1538
     39! DEPENDENCIES
     40! ------------
    1641#ifdef CPP_IOIPSL
    1742    use IOIPSL
    … …  
    2045    use ioipsl_getincom
    2146#endif
    22 
    2347#ifdef CPP_XIOS
    2448    use wxios ! For XIOS outputs
    2549#endif
    2650
     51! DECLARATION
     52! -----------
    2753implicit none
    2854
    2955#include "iniprint.h"
    3056
    31 ! Arguments
    32 !----------
    33 character(*), intent(in) :: modname
    34 integer,      intent(in) :: ierr
    35 character(*), intent(in) :: message
     57! ARGUMENTS
     58! ---------
     59character(*), intent(in) :: modname ! name of calling program
     60integer,      intent(in) :: ierr    ! severity of situation (= 0 normal)
     61character(*), intent(in) :: message ! stuff to print
    3662
    37 ! Code
    38 !-----
     63! CODE
     64! ----
    3965#ifdef CPP_XIOS
    4066    CALL wxios_close() ! Closing XIOS properly

  • trunk/LMDZ.COMMON/libf/evolution/stopping_crit.F90

    r3989 r3991  
    11MODULE stopping_crit
    2 
    3 implicit none
    4 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     stopping_crit
     5!
     6! DESCRIPTION
     7!     Stopping criteria for PEM.
     8!
     9! AUTHORS & DATE
     10!     JB Clement, 2025
     11!
     12! NOTES
     13!
     14!-----------------------------------------------------------------------
     15
     16! DECLARATION
     17! -----------
     18implicit none
     19
     20! MODULE VARIABLES
     21! ----------------
    522real :: h2oice_crit ! Percentage of change of the surface of h2o ice sublimating before stopping the PEM
    623real :: co2ice_crit ! Percentage of change of the surface of co2 ice sublimating before stopping the PEM
    7 real :: ps_crit ! Percentage of change of averaged surface pressure before stopping the PEM
     24real :: ps_crit     ! Percentage of change of averaged surface pressure before stopping the PEM
    825
    926type :: stopFlags
    … …  
    2239end type
    2340
    24 !=======================================================================
    2541contains
    26 !=======================================================================
     42!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    2743
    2844!=======================================================================
    2945FUNCTION is_any_set(this) RESULT(stopPEM)
    30 ! To detect if any flag is true
    31 
    32 class(stopFlags), intent(in) :: this
    33 logical :: stopPEM
    34 
     46!-----------------------------------------------------------------------
     47! NAME
     48!     is_any_set
     49!
     50! DESCRIPTION
     51!     Detect if any stopping flag is set to true.
     52!
     53! AUTHORS & DATE
     54!     JB Clement, 2025
     55!
     56! NOTES
     57!
     58!-----------------------------------------------------------------------
     59
     60! DECLARATION
     61! -----------
     62implicit none
     63
     64! ARGUMENTS
     65! ---------
     66class(stopFlags), intent(in) :: this    ! Stopping flags object
     67logical                      :: stopPEM ! True if any flag is set
     68
     69! CODE
     70! ----
    3571stopPEM = this%surface_h2oice_change .or. &
    3672          this%zero_dh2oice          .or. &
    … …  
    4682
    4783!=======================================================================
    48 ! To give the digit code corresponding to the stopping flag
    4984FUNCTION stop_code(this) result(code)
    50    
    51 class(StopFlags), intent(in) :: this
    52 integer :: code
    53 
     85!-----------------------------------------------------------------------
     86! NAME
     87!     stop_code
     88!
     89! DESCRIPTION
     90!     Return digit code corresponding to the active stopping flag.
     91!
     92! AUTHORS & DATE
     93!     JB Clement, 2025
     94!
     95! NOTES
     96!
     97!-----------------------------------------------------------------------
     98
     99! DECLARATION
     100! -----------
     101implicit none
     102
     103! ARGUMENTS
     104! ---------
     105class(StopFlags), intent(in) :: this ! Stopping flags object
     106integer                      :: code ! Code corresponding to active flag (0 if none)
     107
     108! CODE
     109! ----
    54110code = 0
    55111if     (this%surface_h2oice_change) then; code = 1
    … …  
    67123
    68124!=======================================================================
    69 ! To give the message corresponding to the stopping flag
    70125FUNCTION stop_message(this) result(msg)
    71    
    72 class(StopFlags), intent(in) :: this
    73 character(:), allocatable :: msg
    74 integer :: stopCode
    75 
     126!-----------------------------------------------------------------------
     127! NAME
     128!     stop_message
     129!
     130! DESCRIPTION
     131!     Return descriptive message corresponding to the active stopping flag.
     132!
     133! AUTHORS & DATE
     134!     JB Clement, 2025
     135!
     136! NOTES
     137!
     138!-----------------------------------------------------------------------
     139
     140! DECLARATION
     141! -----------
     142implicit none
     143
     144! ARGUMENTS
     145! ---------
     146class(StopFlags), intent(in) :: this ! Stopping flags object
     147character(:), allocatable    :: msg  ! Descriptive message for active flag
     148
     149! CODE
     150! ----
    76151if     (this%surface_h2oice_change) then; msg = " **** STOPPING because surface of h2o ice has changed too much (see message above)."
    77152elseif (this%zero_dh2oice)          then; msg = " **** STOPPING because tendencies on h2o ice = 0 (see message above)."
    … …  
    87162!=======================================================================
    88163
    89 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    90 !!!
    91 !!! Purpose: Criterions to check if the PEM needs to call the PCM
    92 !!! Author: RV & LL, 02/2023
    93 !!!
    94 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    95 
     164!=======================================================================
    96165SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopCrit,ngrid)
    97166
     167!-----------------------------------------------------------------------
     168! NAME
     169!     stopping_crit_h2o_ice
     170!
     171! DESCRIPTION
     172!     Check if H2O ice surface change criterion to stop PEM is reached.
     173!
     174! AUTHORS & DATE
     175!     R. Vandemeulebrouck
     176!     L. Lange, 2023
     177!     JB Clement, 2023-2025
     178!
     179! NOTES
     180!
     181!-----------------------------------------------------------------------
     182
     183! DEPENDENCIES
     184! ------------
    98185use comslope_mod,  only: subslope_dist, nslope
    99186
    100 implicit none
    101 
    102 !=======================================================================
    103 !
    104 ! Routine to check if the h2o ice criterion to stop the PEM is reached
    105 !
    106 !=======================================================================
    107 ! Inputs
    108 !-------
    109 integer,                          intent(in) :: ngrid                ! # of physical grid points
    110 real,    dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
    111 real,    dimension(ngrid,nslope), intent(in) :: h2o_ice              ! Actual density of h2o ice
    112 real,                             intent(in) :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
    113 logical, dimension(ngrid,nslope), intent(in) :: is_h2oice_sublim_ini ! Grid points where h2o ice was initially sublimating
    114 ! Outputs
    115 !--------
    116 type(stopFlags), intent(inout) :: stopCrit ! Stopping criterion
    117 ! Local variables
     187! DECLARATION
     188! -----------
     189implicit none
     190
     191! ARGUMENTS
     192! ---------
     193integer,                          intent(in)    :: ngrid                ! # of physical grid points
     194real,    dimension(ngrid),        intent(in)    :: cell_area            ! Area of the cells
     195real,    dimension(ngrid,nslope), intent(in)    :: h2o_ice              ! Actual density of h2o ice
     196real,                             intent(in)    :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
     197logical, dimension(ngrid,nslope), intent(in)    :: is_h2oice_sublim_ini ! Grid points where h2o ice was initially sublimating
     198type(stopFlags),                  intent(inout) :: stopCrit             ! Stopping criterion
     199
     200! LOCAL VARIABLES
    118201! ---------------
    119202integer :: i, islope       ! Loop
    120203real    :: h2oice_now_surf ! Current surface of h2o ice
    121204
    122 !=======================================================================
     205! CODE
     206! ----
    123207if (stopCrit%stop_code() > 0) return
    124208
    … …  
    153237!=======================================================================
    154238SUBROUTINE stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopCrit,ngrid,ps_avg_global_ini,ps_avg_global,nslope)
    155 
     239!-----------------------------------------------------------------------
     240! NAME
     241!     stopping_crit_co2
     242!
     243! DESCRIPTION
     244!     Check if CO2 ice surface and pressure change criteria are reached.
     245!
     246! AUTHORS & DATE
     247!     R. Vandemeulebrouck
     248!     L. Lange, 2023
     249!     JB Clement, 2023-2025
     250!
     251! NOTES
     252!
     253!-----------------------------------------------------------------------
     254
     255! DEPENDENCIES
     256! ------------
    156257use comslope_mod,  only: subslope_dist
    157258
    158 implicit none
    159 
    160 !=======================================================================
    161 !
    162 ! Routine to check if the co2 and pressure criteria to stop the PEM are reached
    163 !
    164 !=======================================================================
    165 ! Inputs
    166 !-------
    167 integer,                          intent(in) :: ngrid, nslope          ! # of grid physical grid points
    168 real,    dimension(ngrid),        intent(in) :: cell_area              ! Area of the cells
    169 real,    dimension(ngrid,nslope), intent(in) :: co2_ice                ! Actual density of co2 ice
    170 real,                             intent(in) :: co2ice_sublim_surf_ini ! Initial surface of sublimatingco2 ice
    171 logical, dimension(ngrid,nslope), intent(in) :: is_co2ice_sublim_ini   ! Grid points where co2 ice was initially sublimating
    172 real,                             intent(in) :: ps_avg_global_ini      ! Planet average pressure from the PCM start files
    173 real,                             intent(in) :: ps_avg_global          ! Planet average pressure from the PEM computations
    174 ! Outputs
    175 !--------
    176 type(stopFlags), intent(inout) :: stopCrit ! Stopping criterion
    177 
    178 ! Local variables
     259! DECLARATION
     260! -----------
     261implicit none
     262
     263! ARGUMENTS
     264! ---------
     265integer,                          intent(in)    :: ngrid, nslope          ! # of grid physical grid points
     266real,    dimension(ngrid),        intent(in)    :: cell_area              ! Area of the cells
     267real,    dimension(ngrid,nslope), intent(in)    :: co2_ice                ! Actual density of co2 ice
     268real,                             intent(in)    :: co2ice_sublim_surf_ini ! Initial surface of sublimatingco2 ice
     269logical, dimension(ngrid,nslope), intent(in)    :: is_co2ice_sublim_ini   ! Grid points where co2 ice was initially sublimating
     270real,                             intent(in)    :: ps_avg_global_ini      ! Planet average pressure from the PCM start files
     271real,                             intent(in)    :: ps_avg_global          ! Planet average pressure from the PEM computations
     272type(stopFlags),                  intent(inout) :: stopCrit               ! Stopping criterion
     273
     274! LOCAL VARIABLES
    179275! ---------------
    180276integer :: i, islope       ! Loop
    181277real    :: co2ice_now_surf ! Current surface of co2 ice
    182278
    183 !=======================================================================
     279! CODE
     280! ----
    184281if (stopCrit%stop_code() > 0) return
    185282
    … …  
    230327!=======================================================================
    231328SUBROUTINE stopping_crit_h2o(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopCrit)
    232 
     329!-----------------------------------------------------------------------
     330! NAME
     331!     stopping_crit_h2o
     332!
     333! DESCRIPTION
     334!     Check if PEM oscillates infinitely with H2O (no net exchange).
     335!
     336! AUTHORS & DATE
     337!     L. Lange, 2024
     338!     JB Clement, 2025
     339!
     340! NOTES
     341!
     342!-----------------------------------------------------------------------
     343
     344! DEPENDENCIES
     345! ------------
    233346use evolution,      only: dt
    234347use comslope_mod,   only: subslope_dist, def_slope_mean
    235348use phys_constants, only: pi
    236349
    237 implicit none
    238 
    239 !=======================================================================
    240 !
    241 ! Routine to check if the h2o is only exchanged between grid points
    242 !
    243 !=======================================================================
    244 ! Inputs
    245 !-------
    246 integer,                       intent(in) :: ngrid, nslope            ! # of grid points, # of subslopes
    247 real, dimension(ngrid),        intent(in) :: cell_area                ! Area of each mesh grid (m^2)
    248 real, dimension(ngrid),        intent(in) :: delta_h2o_adsorbed       ! Mass of H2O adsorbed/desorbded in the soil (kg/m^2)
    249 real, dimension(ngrid),        intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that condensed/sublimated at the ice table
    250 real, dimension(ngrid,nslope), intent(in) :: h2o_ice                  ! H2O ice (kg/m^2)
    251 real, dimension(ngrid,nslope), intent(in) :: d_h2oice                 ! Tendency of H2O ice (kg/m^2/year)
    252 ! Outputs
    253 !--------
    254 type(stopFlags), intent(inout) :: stopCrit ! Stopping criterion
    255 real, intent(out) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to conserve H2O
    256 ! Local variables
     350! DECLARATION
     351! -----------
     352implicit none
     353
     354! ARGUMENTS
     355! ---------
     356integer,                       intent(in)    :: ngrid, nslope            ! # of grid points, # of subslopes
     357real, dimension(ngrid),        intent(in)    :: cell_area                ! Area of each mesh grid (m^2)
     358real, dimension(ngrid),        intent(in)    :: delta_h2o_adsorbed       ! Mass of H2O adsorbed/desorbded in the soil (kg/m^2)
     359real, dimension(ngrid),        intent(in)    :: delta_h2o_icetablesublim ! Mass of H2O that condensed/sublimated at the ice table
     360real, dimension(ngrid,nslope), intent(in)    :: h2o_ice                  ! H2O ice (kg/m^2)
     361real, dimension(ngrid,nslope), intent(in)    :: d_h2oice                 ! Tendency of H2O ice (kg/m^2/year)
     362real,                          intent(out)   :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to conserve H2O
     363type(stopFlags),               intent(inout) :: stopCrit                 ! Stopping criterion
     364
     365! LOCAL VARIABLES
    257366! ---------------
    258367integer :: i, islope ! Loop
    259 !=======================================================================
     368
     369! CODE
     370! ----
    260371! Checks to escape the subroutine (adjusment not relevant in 1D or if the PEM should stop already)
    261372if (ngrid == 1 .or. stopCrit%stop_code() > 0) then
    … …  
    307418!=======================================================================
    308419
    309 END MODULE
     420END MODULE stopping_crit

  • trunk/LMDZ.COMMON/libf/evolution/subsurface_ice.F90

    r3989 r3991  
    11MODULE subsurface_ice
    2 
    3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4 !!!
    5 !!! Purpose:  Retreat and growth of subsurface ice on Mars
    6 !!!           orbital elements remain constant
    7 !!!
    8 !!!
    9 !!! Author: EV, updated NS MSIM dynamical program for the PEM
    10 !!!         Taken from Norbert Schorgofer's code
    11 !!!
    12 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    13 
    14 implicit none
    15 
    16   ! parameters for thermal model
    17   ! they are only used in the subroutines below
    18   real(8), parameter :: dt = 0.02  ! in units of Mars solar days
    19   !real(8), parameter :: Fgeotherm = 0.
    20   real(8), parameter :: Fgeotherm = 0 !0.028  ! [W/m^2]
    21   real(8), parameter :: Lco2frost=6.0e5, co2albedo=0.60, co2emiss=1.
    22   real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
    23   integer, parameter :: EQUILTIME =15 ! [Mars years]
    24 
    25   integer, parameter :: NMAX = 1000
     2!-----------------------------------------------------------------------
     3! NAME
     4!     subsurface_ice
     5!
     6! DESCRIPTION
     7!     Retreat and growth of subsurface ice on Mars with constant orbital
     8!     elements.
     9!
     10! AUTHORS & DATE
     11!     N. Schorghofer (MSIM dynamical program)
     12!     E. Vos, 2025
     13!     JB Clement, 2025
     14!
     15! NOTES
     16!    Based on Norbert Schorgofer's code. Updated for PEM.
     17!-----------------------------------------------------------------------
     18
     19! DECLARATION
     20! -----------
     21implicit none
     22
     23! MODULE PARAMETERS
     24! -----------------
     25real(8), parameter :: dt = 0.02  ! in units of Mars solar days
     26!real(8), parameter :: Fgeotherm = 0.
     27real(8), parameter :: Fgeotherm = 0 !0.028  ! [W/m^2]
     28real(8), parameter :: Lco2frost=6.0e5, co2albedo=0.60, co2emiss=1.
     29real(8), parameter :: emiss0 = 1.  ! emissivity of dry surface
     30integer, parameter :: EQUILTIME =15 ! [Mars years]
     31
     32integer, parameter :: NMAX = 1000
    2633
    2734contains
     35!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    2836
    2937
    … …  
    4149  use constants_marspem_mod, only: sec_per_sol
    4250  use phys_constants,        only: pi
    43   implicit none
     51  ! DECLARATION
     52! -----------
     53implicit none
    4454  integer, parameter :: NP=1   ! # of sites
    4555  integer nz, i, k, iloop
    … …  
    212222!     output are newti and newrhoc
    213223!***********************************************************************
    214   implicit none
     224  ! DECLARATION
     225! -----------
     226implicit none
    215227  integer, intent(IN) :: layertype
    216228  real(8), intent(IN) :: porosity, fill, rhocobs, tiobs
    … …  
    282294!     input is partial pressure [Pascal]
    283295!     output is temperature [Kelvin]
    284       implicit none
     296      ! DECLARATION
     297! -----------
     298implicit none
    285299      real*8 p
    286300
    … …  
    326340  use ice_table,      only: rho_ice
    327341  !use omp_lib
    328   implicit none
     342  ! DECLARATION
     343! -----------
     344implicit none
    329345  integer, intent(IN) :: nz, NP
    330346  real(8), intent(IN) :: bigstep
    … …  
    442458!***********************************************************************
    443459  use constants_marspem_mod, only: sols_per_my, sec_per_sol
    444   implicit none
     460  ! DECLARATION
     461! -----------
     462implicit none
    445463  integer, intent(IN) :: nz, typeT
    446464!  real(8), intent(IN) :: latitude  ! in radians
    … …  
    595613!     the inverse of function psvco2
    596614!     input is pressure in Pascal, output is temperature in Kelvin
     615! DECLARATION
     616! -----------
    597617implicit none
    598618real*8 p
    … …  
    605625!     saturation vapor pressure of H2O ice [Pascal]
    606626!     input is temperature [Kelvin]
    607       implicit none
     627      ! DECLARATION
     628! -----------
     629implicit none
    608630      real*8 T
    609631
    … …  
    642664pure function zint(y1,y2,z1,z2)
    643665  ! interpolate between two values, y1*y2<0
    644   implicit none
     666  ! DECLARATION
     667! -----------
     668implicit none
    645669  real(8), intent(IN) :: y1,y2,z1,z2
    646670  real(8) zint
    … …  
    655679!  this is not the true (final) equilibrium depth
    656680!***********************************************************************
    657   implicit none
     681  ! DECLARATION
     682! -----------
     683implicit none
    658684  integer, intent(IN) :: nz
    659685  real(8), intent(IN) :: z(nz), rhosatav(nz)
    … …  
    691717  use math_toolkit, only: deriv2_simple, deriv1_onesided, deriv1, colint
    692718  use ice_table,    only: constriction
    693   implicit none
     719  ! DECLARATION
     720! -----------
     721implicit none
    694722  integer, intent(IN) :: nz, typeT
    695723  real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill
    … …  
    822850  use math_toolkit, only: colint
    823851  use ice_table,    only: rho_ice
    824   implicit none
     852  ! DECLARATION
     853! -----------
     854implicit none
    825855  integer, intent(IN) :: nz, typeF, typeG
    826856  real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP

  • trunk/LMDZ.COMMON/libf/evolution/surf_ice.F90

    r3989 r3991  
    11MODULE surf_ice
    2 
    3 implicit none
    4 
    5 ! Different types of ice
    6 real, dimension(:,:), allocatable :: h2o_ice
    7 real, dimension(:,:), allocatable :: co2_ice
    8 
    9 ! Threshold to consider H2O ice as watercap (infinite reservoir) [kg.m-2]
    10 real :: h2oice_cap_threshold
    11 
    12 ! Ice properties
     2!-----------------------------------------------------------------------
     3! NAME
     4!     surf_ice
     5!
     6! DESCRIPTION
     7!     Surface ice management.
     8!
     9! AUTHORS & DATE
     10!     JB Clement, 12/2025
     11!
     12! NOTES
     13!
     14!-----------------------------------------------------------------------
     15
     16! DECLARATION
     17! -----------
     18implicit none
     19
     20! MODULE VARIABLES
     21! ----------------
     22real, dimension(:,:), allocatable :: h2o_ice ! H2O ice surface [kg.m-2]
     23real, dimension(:,:), allocatable :: co2_ice ! CO2 ice surface [kg.m-2]
     24real :: h2oice_cap_threshold                 ! Threshold to consider H2O ice as infinite reservoir [kg.m-2]
     25
     26! MODULE PARAMETERS
     27! -----------------
    1328real, parameter :: rho_co2ice = 1650. ! Density of CO2 ice [kg.m-3]
    1429real, parameter :: rho_h2oice = 920.  ! Density of H2O ice [kg.m-3]
    1530
    16 !=======================================================================
    1731contains
    18 !=======================================================================
     32!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    1933
    2034!=======================================================================
    2135SUBROUTINE set_perice4PCM(ngrid,nslope,PCMfrost,is_h2o_perice,h2oice_PCM,co2ice_PCM)
    22 
     36!-----------------------------------------------------------------------
     37! NAME
     38!     set_perice4PCM
     39!
     40! DESCRIPTION
     41!     Reconstruct perennial ice for PCM from PEM ice computations.
     42!
     43! AUTHORS & DATE
     44!     JB Clement, 12/2025
     45!
     46! NOTES
     47!
     48!-----------------------------------------------------------------------
     49
     50! DEPENDENCIES
     51! ------------
    2352use metamorphism,   only: iPCM_h2ofrost
    2453use comslope_mod,   only: subslope_dist, def_slope_mean
    2554use phys_constants, only: pi
    2655
    27 
    28 implicit none
    29 
    30 ! Arguments
    31 !----------
    32 integer, intent(in) :: ngrid, nslope
    33 real, dimension(:,:,:), intent(inout) :: PCMfrost
    34 logical, dimension(ngrid),        intent(out) :: is_h2o_perice
    35 real,    dimension(ngrid,nslope), intent(out) :: co2ice_PCM, h2oice_PCM
    36 
    37 ! Local variables
    38 !----------------
     56! DECLARATION
     57! -----------
     58implicit none
     59
     60! ARGUMENTS
     61! ---------
     62integer,                       intent(in)    :: ngrid, nslope
     63real, dimension(:,:,:),        intent(inout) :: PCMfrost               ! PCM frost
     64logical, dimension(ngrid),     intent(out)   :: is_h2o_perice          ! H2O perennial ice flag
     65real, dimension(ngrid,nslope), intent(out)   :: h2oice_PCM, co2ice_PCM ! Ice for PCM
     66
     67! LOCAL VARIABLES
     68! ---------------
    3969integer :: i
    4070
    41 ! Code
    42 !-----
     71! CODE
     72! ----
    4373write(*,*) '> Reconstructing perennial ic for the PCM'
    4474co2ice_PCM(:,:) = co2_ice(:,:)
    … …  
    6292!=======================================================================
    6393SUBROUTINE ini_surf_ice(ngrid,nslope)
    64 
    65 implicit none
    66 
    67 ! Arguments
    68 !----------
     94!-----------------------------------------------------------------------
     95! NAME
     96!     ini_surf_ice
     97!
     98! DESCRIPTION
     99!     Initialize surface ice arrays.
     100!
     101! AUTHORS & DATE
     102!     JB Clement 12/2025
     103!
     104! NOTES
     105!
     106!-----------------------------------------------------------------------
     107
     108! DECLARATION
     109! -----------
     110implicit none
     111
     112! ARGUMENTS
     113! ---------
    69114integer, intent(in) :: ngrid, nslope
    70115
    71 ! Local variables
    72 !----------------
    73 
    74 ! Code
    75 !-----
     116! CODE
     117! ----
    76118if (.not. allocated(h2o_ice)) allocate(h2o_ice(ngrid,nslope))
    77119if (.not. allocated(co2_ice)) allocate(co2_ice(ngrid,nslope))
    … …  
    82124!=======================================================================
    83125SUBROUTINE end_surf_ice()
    84 
    85 implicit none
    86 
    87 ! Arguments
    88 !----------
    89 
    90 ! Local variables
    91 !----------------
    92 
    93 ! Code
    94 !-----
     126!-----------------------------------------------------------------------
     127! NAME
     128!     end_surf_ice
     129!
     130! DESCRIPTION
     131!     Deallocate surface ice arrays.
     132!
     133! AUTHORS & DATE
     134!     JB Clement, 12/2025
     135!
     136! NOTES
     137!
     138!-----------------------------------------------------------------------
     139
     140! DECLARATION
     141! -----------
     142implicit none
     143
     144! CODE
     145! ----
    95146if (allocated(h2o_ice)) deallocate(h2o_ice)
    96147if (allocated(co2_ice)) deallocate(co2_ice)
    … …  
    101152!=======================================================================
    102153SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
    103 ! Routine to compute the evolution of CO2 ice
    104 
     154!-----------------------------------------------------------------------
     155! NAME
     156!     evol_co2_ice
     157!
     158! DESCRIPTION
     159!     Compute the evolution of CO2 ice over one year.
     160!
     161! AUTHORS & DATE
     162!     R. Vandemeulebrouck
     163!     L. Lange
     164!     JB Clement, 2023-2025
     165!
     166! NOTES
     167!
     168!-----------------------------------------------------------------------
     169
     170! DEPENDENCIES
     171! ------------
    105172use evolution, only: dt
    106173
    107 implicit none
    108 
    109 ! Arguments
    110 ! ---------
    111 integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes
    112 real, dimension(ngrid,nslope), intent(inout) :: co2_ice     ! Previous and actual density of CO2 ice
    113 real, dimension(ngrid,nslope), intent(inout) :: d_co2ice    ! Evolution of perennial ice over one year
    114 real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
    115 
    116 ! Local variables
     174! DECLARATION
     175! -----------
     176implicit none
     177
     178! ARGUMENTS
     179! ---------
     180integer,                       intent(in)    :: ngrid, nslope
     181real, dimension(ngrid,nslope), intent(inout) :: co2_ice     ! CO2 ice
     182real, dimension(ngrid,nslope), intent(inout) :: d_co2ice    ! Evolution of CO2 ice over one year
     183real, dimension(ngrid,nslope), intent(out)   :: zshift_surf ! Elevation shift for the surface [m]
     184
     185! LOCAL VARIABLES
    117186! ---------------
    118187real, dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice
    119188
    120 ! Code
     189! CODE
    121190! ----
    122191! Evolution of CO2 ice for each physical point
    … …  
    137206!=======================================================================
    138207SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopCrit)
    139 ! Routine to compute the evolution of h2o ice
    140 
    141 use evolution, only: dt
     208!-----------------------------------------------------------------------
     209! NAME
     210!     evol_h2o_ice
     211!
     212! DESCRIPTION
     213!     Compute the evolution of H2O ice over one year.
     214!
     215! AUTHORS & DATE
     216!     R. Vandemeulebrouck
     217!     L. Lange
     218!     JB Clement, 2023-2025
     219!
     220! NOTES
     221!
     222!-----------------------------------------------------------------------
     223
     224! DEPENDENCIES
     225! ------------
     226use evolution,     only: dt
    142227use stopping_crit, only: stopping_crit_h2o, stopFlags
    143228
    144 implicit none
    145 
    146 ! Arguments
    147 ! ---------
    148 integer,                intent(in) :: ngrid, nslope            ! # of grid points, # of subslopes
    149 real, dimension(ngrid), intent(in) :: cell_area                ! Area of each mesh grid (m^2)
    150 real, dimension(ngrid), intent(in) :: delta_h2o_adsorbed       ! Mass of H2O adsorbed/desorbded in the soil (kg/m^2)
    151 real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2)
    152 real, dimension(ngrid,nslope), intent(inout) :: h2o_ice  ! H2O ice (kg/m^2)
    153 real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Tendency of H2O ice (kg/m^2/year)
    154 type(stopFlags),               intent(inout) :: stopCrit ! Stopping criterion
    155 real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
    156 
    157 ! Local variables
     229! DECLARATION
     230! -----------
     231implicit none
     232
     233! ARGUMENTS
     234! ---------
     235integer,                       intent(in)    :: ngrid, nslope
     236real, dimension(ngrid),        intent(in)    :: cell_area                ! Area of each mesh grid (m^2)
     237real, dimension(ngrid),        intent(in)    :: delta_h2o_adsorbed       ! Mass of H2O adsorbed/desorbded in soil (kg/m^2)
     238real, dimension(ngrid),        intent(in)    :: delta_h2o_icetablesublim ! Mass condensed/sublimated at ice table (kg/m^2)
     239real, dimension(ngrid,nslope), intent(inout) :: h2o_ice                  ! H2O ice (kg/m^2)
     240real, dimension(ngrid,nslope), intent(inout) :: d_h2oice                 ! Tendency of H2O ice (kg/m^2/year)
     241type(stopFlags),               intent(inout) :: stopCrit                 ! Stopping criterion
     242real, dimension(ngrid,nslope), intent(out)   :: zshift_surf              ! Elevation shift for the surface [m]
     243
     244! LOCAL VARIABLES
    158245! ---------------
    159 integer                       :: i, islope                                                ! Loop variables
    160 real                          :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to balance H2O
     246integer                       :: i, islope
     247real                          :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! H2O balance variables
    161248real, dimension(ngrid,nslope) :: d_h2oice_new                                             ! Tendencies computed to keep balance
    162249
    163 ! Code
     250! CODE
    164251! ----
    165252write(*,*) '> Evolution of H2O ice'
    … …  
    182269!=======================================================================
    183270SUBROUTINE balance_h2oice_reservoirs(ngrid,nslope,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new)
    184 ! Routine to balance the H2O flux from and into the atmosphere accross reservoirs
    185 
     271!-----------------------------------------------------------------------
     272! NAME
     273!     balance_h2oice_reservoirs
     274!
     275! DESCRIPTION
     276!     Balance H2O flux from and into atmosphere across reservoirs.
     277!
     278! AUTHORS & DATE
     279!     JB Clement, 2025
     280!
     281! NOTES
     282!
     283!-----------------------------------------------------------------------
     284
     285! DEPENDENCIES
     286! ------------
    186287use evolution, only: dt
    187288
    188 implicit none
    189 
    190 ! Arguments
    191 ! ---------
    192 integer,                       intent(in) :: ngrid, nslope ! # of grid points, # of subslopes
    193 real,                          intent(in) :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! ! Variables to conserve H2O
    194 real, dimension(ngrid,nslope), intent(in) :: h2o_ice ! H2O ice (kg/m^2)
    195 real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Tendency of H2O ice (kg/m^2/year)
    196 real, dimension(ngrid,nslope), intent(out) :: d_h2oice_new ! Adjusted tendencies to keep balance between donor and recipient reservoirs
    197 
    198 ! Local variables
     289! DECLARATION
     290! -----------
     291implicit none
     292
     293! ARGUMENTS
     294! ---------
     295integer,                       intent(in)    :: ngrid, nslope ! # of grid points, # of subslopes
     296real,                          intent(in)    :: S_atm_2_h2o, S_h2o_2_atm, S_atm_2_h2oice, S_h2oice_2_atm ! Variables to conserve H2O
     297real, dimension(ngrid,nslope), intent(in)    :: h2o_ice       ! H2O ice (kg/m^2)
     298real, dimension(ngrid,nslope), intent(inout) :: d_h2oice      ! Tendency of H2O ice (kg/m^2/year)
     299real, dimension(ngrid,nslope), intent(out)   :: d_h2oice_new  ! Adjusted tendencies to keep balance between donor and recipient reservoirs
     300
     301! LOCAL VARIABLES
    199302! ---------------
    200303integer :: i, islope
    201 real    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice, d_target
    202 
    203 ! Code
     304real    :: S_target, S_target_subl_h2oice, S_target_cond_h2oice, S_ghostice, d_target ! Balance variables
     305
     306! CODE
    204307! ----
    205308S_target = (S_atm_2_h2o + S_h2o_2_atm)/2.

  • trunk/LMDZ.COMMON/libf/evolution/surf_temp.F90

    r3989 r3991  
    11MODULE surf_temp
    2 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     surf_temp
     5!
     6! DESCRIPTION
     7!     Surface temperature management.
     8!
     9! AUTHORS & DATE
     10!     JB Clement, 2025
     11!
     12! NOTES
     13!
     14!-----------------------------------------------------------------------
     15
     16! DECLARATION
     17! -----------
    318implicit none
    419
    5 !=======================================================================
    620contains
    7 !=======================================================================
     21!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    822
    923!=======================================================================
    1024SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,nlon,nlat,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
    11 
     25!-----------------------------------------------------------------------
     26! NAME
     27!     update_tsurf_nearest_baresoil
     28!
     29! DESCRIPTION
     30!     Update surface temperature where CO2 ice disappeared using nearby
     31!     bare soil temperature.
     32!
     33! AUTHORS & DATE
     34!     JB Clement, 2025
     35!
     36! NOTES
     37!
     38!-----------------------------------------------------------------------
     39
     40! DEPENDENCIES
     41! ------------
    1242use grid_conversion, only: vect2lonlat, lonlat2vect
    1343
     44! DECLARATION
     45! -----------
    1446implicit none
    1547
    16 ! Inputs:
    17 integer,                          intent(in) :: nlon, nlat, nslope, ngrid
    18 real,    dimension(ngrid,nslope), intent(in) :: co2_ice
    19 real,    dimension(ngrid),        intent(in) :: latitude
    20 logical, dimension(ngrid,nslope), intent(in) :: is_co2ice_ini
    21 ! Outputs:
    22 real,    dimension(ngrid,nslope), intent(inout) :: tsurf_avg
    23 logical, dimension(ngrid,nslope), intent(inout) :: co2ice_disappeared
    24 ! Local variables:
    25 real, parameter                   :: eps = 1.e-10
    26 integer                           :: islope, i, j, k, radius, rmax, di, dj, ii, jj
    27 logical                           :: found
    28 real, dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
    29 real, dimension(nlon,nlat)        :: latitude_ll
    30 real, dimension(ngrid)            :: tmp
    31 integer, dimension(nslope - 1)    :: priority
    32 
     48! ARGUMENTS
     49! ---------
     50integer,                          intent(in)    :: nlon, nlat, nslope, ngrid ! Grid dimensions
     51real,    dimension(ngrid,nslope), intent(in)    :: co2_ice                   ! CO2 ice density
     52real,    dimension(ngrid),        intent(in)    :: latitude                  ! Latitude
     53logical, dimension(ngrid,nslope), intent(in)    :: is_co2ice_ini             ! Initial CO2 ice flag
     54real,    dimension(ngrid,nslope), intent(inout) :: tsurf_avg                 ! Average surface temperature
     55logical, dimension(ngrid,nslope), intent(inout) :: co2ice_disappeared        ! Ice disappeared flag
     56
     57! LOCAL VARIABLES
     58! ---------------
     59real, parameter                      :: eps = 1.e-10
     60integer                              :: islope, i, j, k, radius, rmax, di, dj, ii, jj
     61logical                              :: found
     62real,    dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
     63real,    dimension(nlon,nlat)        :: latitude_ll
     64real,    dimension(ngrid)            :: tmp
     65integer, dimension(nslope - 1)       :: priority
     66
     67! CODE
     68! ----
    3369! Check to escape the subroutine (not relevant in 1D)
    3470if (ngrid == 1) return
    … …  
    110146!=======================================================================
    111147SUBROUTINE get_slope_priority(lat,nslope,islope,priority)
    112 ! Priority given to equator-ward slope which are most likely to hold no ice
    113 
     148!-----------------------------------------------------------------------
     149! NAME
     150!     get_slope_priority
     151!
     152! DESCRIPTION
     153!     Determine slope priority based on latitude (equator-ward favored).
     154!
     155! AUTHORS & DATE
     156!     JB Clement, 2025
     157!
     158! NOTES
     159!     Equator-ward slopes are most likely to hold no ice.
     160!-----------------------------------------------------------------------
     161
     162! DECLARATION
     163! -----------
    114164implicit none
    115165
    116 ! Inputs:
    117 real,    intent(in) :: lat
    118 integer, intent(in) :: nslope, islope
    119 ! Outputs:
    120 integer, dimension(nslope - 1), intent(out) :: priority
    121 ! Locals:
     166! ARGUMENTS
     167! ---------
     168real,                           intent(in)  :: lat      ! Latitude [degrees]
     169integer,                        intent(in)  :: nslope, islope
     170integer, dimension(nslope - 1), intent(out) :: priority ! Priority ordering of slopes
     171
     172! LOCAL VARIABLES
     173! ---------------
    122174integer :: i, k
    123175
    124 ! Code
    125 !-----
     176! CODE
     177! ----
    126178k = 1
    127179

  • trunk/LMDZ.COMMON/libf/evolution/tendencies.F90

    r3989 r3991  
    11MODULE tendencies
    2 
    3 implicit none
    4 
    5 !=======================================================================
     2!-----------------------------------------------------------------------
     3! NAME
     4!     tendencies
     5!
     6! DESCRIPTION
     7!     Computation and update of PEM ice evolution tendencies.
     8!
     9! AUTHORS & DATE
     10!     R. Vandemeulebrouck
     11!     L. Lange
     12!     JB Clement, 2023-2025
     13!
     14! NOTES
     15!
     16!-----------------------------------------------------------------------
     17
     18! DECLARATION
     19! -----------
     20implicit none
     21
    622contains
    7 !=======================================================================
     23!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    824
    925!=======================================================================
    1026SUBROUTINE compute_tend(ngrid,nslope,min_ice,d_ice)
    11 ! Compute the initial tendencies of the ice evolution based on the PCM data
    12 
    13 implicit none
    14 
    15 ! Arguments
     27!-----------------------------------------------------------------------
     28! NAME
     29!     compute_tend
     30!
     31! DESCRIPTION
     32!     Compute initial ice evolution tendencies from PCM data.
     33!
     34! AUTHORS & DATE
     35!     R. Vandemeulebrouck
     36!     L. Lange
     37!     JB Clement, 2023-2025
     38!
     39! NOTES
     40!     Based on minima of ice at each point for the PCM years.
     41!-----------------------------------------------------------------------
     42
     43! DECLARATION
     44! -----------
     45implicit none
     46
     47! ARGUMENTS
    1648! ---------
    17 integer,                         intent(in) :: ngrid   ! # of grid points
    18 integer,                         intent(in) :: nslope  ! # of subslopes
    19 real, dimension(ngrid,nslope,2), intent(in) :: min_ice ! Minima of ice at each point for the PCM years
    20 real, dimension(ngrid,nslope), intent(out) :: d_ice ! Difference between the minima = evolution of perennial ice
    21 
    22 ! Code
     49integer,                         intent(in)  :: ngrid
     50integer,                         intent(in)  :: nslope
     51real, dimension(ngrid,nslope,2), intent(in)  :: min_ice ! Minima of ice at each point for the PCM years
     52real, dimension(ngrid,nslope),   intent(out) :: d_ice   ! Evolution of perennial ice
     53
     54! CODE
    2355! ----
    2456! We compute the difference
    … …  
    3769SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, &
    3870                           vmr_co2_PCM,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global)
    39 ! To compute the evolution of the tendency for co2 ice
    40 
     71!-----------------------------------------------------------------------
     72! NAME
     73!     recomp_tend_co2
     74!
     75! DESCRIPTION
     76!     Recompute CO2 ice tendency based on pressure and atmospheric changes.
     77!
     78! AUTHORS & DATE
     79!     L. Lange
     80!     JB Clement, 2023-2025
     81!
     82! NOTES
     83!     Adjusts CO2 ice evolution based on Clausius-Clapeyron changes.
     84!-----------------------------------------------------------------------
     85
     86! DEPENDENCIES
     87! ------------
    4188use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol
    4289
    43 implicit none
    44 
    45 ! Arguments
     90! DECLARATION
     91! -----------
     92implicit none
     93
     94! ARGUMENTS
    4695! ---------
    47 integer,                        intent(in) :: timelen, ngrid, nslope
    48 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM       ! physical point field: Volume mixing ratio of co2 in the first layer
    49 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM       ! physical point field: Volume mixing ratio of co2 in the first layer
    50 real, dimension(ngrid,timelen), intent(in) :: ps_PCM            ! physical point field: Surface pressure in the PCM
    51 real,                           intent(in) :: ps_avg_global_ini ! global averaged pressure at previous timestep
    52 real,                           intent(in) :: ps_avg_global     ! global averaged pressure at current timestep
    53 real, dimension(ngrid,nslope),  intent(in) :: d_co2ice_ini      ! physical point field: Evolution of perennial ice over one year
    54 real, dimension(ngrid,nslope),  intent(in) :: co2ice            ! CO2 ice per mesh and sub-grid slope (kg/m^2)
    55 real, dimension(ngrid,nslope),  intent(in) :: emissivity        ! Emissivity per mesh and sub-grid slope(1)
    56 real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! physical point field: Evolution of perennial ice over one year
    57 
    58 ! Local variables
     96integer,                        intent(in)    :: timelen, ngrid, nslope ! Time length, # of grid points and slopes
     97real, dimension(ngrid,timelen), intent(in)    :: vmr_co2_PCM            ! CO2 VMR in PCM first layer
     98real, dimension(ngrid,timelen), intent(in)    :: vmr_co2_PEM            ! CO2 VMR in PEM first layer
     99real, dimension(ngrid,timelen), intent(in)    :: ps_PCM                 ! Surface pressure in PCM
     100real,                           intent(in)    :: ps_avg_global_ini      ! Global average pressure (initial)
     101real,                           intent(in)    :: ps_avg_global          ! Global average pressure (current)
     102real, dimension(ngrid,nslope),  intent(in)    :: d_co2ice_ini           ! Initial CO2 ice evolution
     103real, dimension(ngrid,nslope),  intent(in)    :: co2ice                 ! CO2 ice surface [kg/m^2]
     104real, dimension(ngrid,nslope),  intent(in)    :: emissivity             ! Emissivity
     105real, dimension(ngrid,nslope),  intent(inout) :: d_co2ice_phys          ! Updated CO2 ice evolution
     106
     107! LOCAL VARIABLES
    59108! ---------------
    60109integer :: i, t, islope
    61110real    :: coef, avg
    62111
    63 ! Code
     112! CODE
    64113! ----
    65114write(*,*) "> Updating the CO2 ice tendency for the new pressure"
    … …  
    86135!=======================================================================
    87136SUBROUTINE recomp_tend_h2o(h2oice_depth_old,h2oice_depth_new,tsurf,tsoil_PEM_timeseries_old,tsoil_PEM_timeseries_new,d_h2oice)
    88 ! To compute the evolution of the tendency for h2o ice
    89 
     137!-----------------------------------------------------------------------
     138! NAME
     139!     recomp_tend_h2o
     140!
     141! DESCRIPTION
     142!     Recompute H2O ice tendency based on soil depth and temperature changes.
     143!
     144! AUTHORS & DATE
     145!     JB Clement, 2025 (following E. Vos's work)
     146!
     147! NOTES
     148!
     149!-----------------------------------------------------------------------
     150
     151! DEPENDENCIES
     152! ------------
    90153use soil_temp,      only: itp_tsoil
    91154use subsurface_ice, only: psv
    92155
    93 implicit none
    94 
    95 ! Arguments
     156! DECLARATION
     157! -----------
     158implicit none
     159
     160! ARGUMENTS
    96161! ---------
    97 real,                 intent(in) :: h2oice_depth_old, h2oice_depth_new, tsurf
    98 real, dimension(:,:), intent(in) :: tsoil_PEM_timeseries_old, tsoil_PEM_timeseries_new
    99 real, intent(inout) :: d_h2oice ! Evolution of perennial ice over one year
    100 
    101 ! Local variables
     162real,                 intent(in)    :: h2oice_depth_old ! Old H2O ice depth
     163real,                 intent(in)    :: h2oice_depth_new ! New H2O ice depth
     164real,                 intent(in)    :: tsurf            ! Surface temperature
     165real, dimension(:,:), intent(in)    :: tsoil_PEM_timeseries_old ! Old soil temperature time series
     166real, dimension(:,:), intent(in)    :: tsoil_PEM_timeseries_new ! New soil temperature time series
     167real,                 intent(inout) :: d_h2oice         ! Evolution of perennial ice
     168
     169! LOCAL VARIABLES
    102170! ---------------
    103171real            :: Rz_old, Rz_new, R_dec, hum_dec, psv_max_old, psv_max_new
    … …  
    106174real, parameter :: zcdv = 0.0325     ! Drag coefficient
    107175
    108 ! Code
     176! CODE
    109177! ----
    110178! Higher resistance due to growing lag layer (higher depth)

  • trunk/LMDZ.COMMON/libf/evolution/tracers.F90

    r3989 r3991  
    11MODULE tracers
     2!-----------------------------------------------------------------------
     3! NAME
     4!     tracers
     5!
     6! DESCRIPTION
     7!     Tracer species management for tracking.
     8!
     9! AUTHORS & DATE
     10!     JB Clement, 12/2025
     11!
     12! NOTES
     13!
     14!-----------------------------------------------------------------------
    215
     16! DECLARATION
     17! -----------
    318implicit none
    419
    5 ! Indices for tracers taken from the PCM
    6 integer :: iPCM_qh2o
     20! MODULE VARIABLES
     21! ----------------
     22integer                         :: iPCM_qh2o ! Index for H2O vapor tracer from PCM
     23real, dimension(:), allocatable :: mmol      ! Molar masses of tracers [g/mol]
    724
    8 ! Molar masses of tracers
    9 real, dimension(:), allocatable :: mmol
    10 
    11 !=======================================================================
    1225contains
    13 !=======================================================================
     26!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    1427
    1528!=======================================================================
    1629SUBROUTINE ini_tracers_id(nqtot,noms)
     30!-----------------------------------------------------------------------
     31! NAME
     32!     ini_tracers_id
     33!
     34! DESCRIPTION
     35!     Initialize tracer indices from PCM tracer names.
     36!
     37! AUTHORS & DATE
     38!     JB Clement, 12/2025
     39!
     40! NOTES
     41!
     42!-----------------------------------------------------------------------
    1743
     44! DECLARATION
     45! -----------
    1846implicit none
    1947
    20 ! Arguments
    21 !----------
    22 integer,                        intent(in) :: nqtot
    23 character(*), dimension(nqtot), intent(in) :: noms
     48! ARGUMENTS
     49! ---------
     50integer,                        intent(in) :: nqtot ! Total number of tracers
     51character(*), dimension(nqtot), intent(in) :: noms  ! Names of tracers
    2452
    25 ! Local variables
    26 !----------------
     53! LOCAL VARIABLES
     54! ---------------
    2755integer :: i
    2856
    29 ! Code
    30 !-----
     57! CODE
     58! ----
    3159! Allocation
    3260call ini_tracers(nqtot)
    … …  
    5179!=======================================================================
    5280SUBROUTINE ini_tracers(nqtot)
     81!-----------------------------------------------------------------------
     82! NAME
     83!     ini_tracers
     84!
     85! DESCRIPTION
     86!     Allocate tracer molar mass array.
     87!
     88! AUTHORS & DATE
     89!     JB Clement, 12/2025
     90!
     91! NOTES
     92!
     93!-----------------------------------------------------------------------
    5394
     95! DECLARATION
     96! -----------
    5497implicit none
    5598
    56 ! Arguments
    57 !----------
    58 integer, intent(in) :: nqtot
     99! ARGUMENTS
     100! ---------
     101integer, intent(in) :: nqtot ! Total number of tracers
    59102
    60 ! Local variables
    61 !----------------
    62 
    63 ! Code
    64 !-----
     103! CODE
     104! ----
    65105if (.not. allocated(mmol)) allocate(mmol(nqtot))
    66106
    … …  
    70110!=======================================================================
    71111SUBROUTINE end_tracers()
     112!-----------------------------------------------------------------------
     113! NAME
     114!     end_tracers
     115!
     116! DESCRIPTION
     117!     Deallocate tracer molar mass array.
     118!
     119! AUTHORS & DATE
     120!     JB Clement, 12/2025
     121!
     122! NOTES
     123!
     124!-----------------------------------------------------------------------
    72125
     126! DECLARATION
     127! -----------
    73128implicit none
    74129
    75 ! Arguments
    76 !----------
    77 
    78 ! Local variables
    79 !----------------
    80 
    81 ! Code
    82 !-----
     130! CODE
     131! ----
    83132if (allocated(mmol)) deallocate(mmol)
    84133

  • trunk/LMDZ.COMMON/libf/evolution/xios_data.F90

    r3989 r3991  
    11MODULE xios_data
    2 
     2!-----------------------------------------------------------------------
     3! NAME
     4!     xios_data
     5!
     6! DESCRIPTION
     7!     Read XIOS output data and process it for PEM initialization.
     8!
     9! AUTHORS & DATE
     10!     JB Clement, 2025
     11!
     12! NOTES
     13!
     14!-----------------------------------------------------------------------
     15
     16! DEPENDENCIES
     17! ------------
    318use netcdf, only: nf90_open, nf90_close, nf90_inquire_dimension, nf90_inq_dimid, nf90_noerr, nf90_nowrite, nf90_get_var, nf90_inq_varid
    419
    5 implicit none
    6 
    7 character(19), parameter :: file1_daily  = "Xoutdaily4pem_Y1.nc"
    8 character(19), parameter :: file2_daily  = "Xoutdaily4pem_Y2.nc"
    9 character(20), parameter :: file1_yearly = "Xoutyearly4pem_Y1.nc"
    10 character(20), parameter :: file2_yearly = "Xoutyearly4pem_Y2.nc"
    11 character(256)           :: msg      ! For reading
    12 integer                  :: fID, vID ! For reading
    13 
     20! DECLARATION
     21! -----------
     22implicit none
     23
     24! MODULE VARIABLES
     25! ----------------
     26character(19), parameter :: file1_daily  = "Xoutdaily4pem_Y1.nc"   ! XIOS daily output file, year 1
     27character(19), parameter :: file2_daily  = "Xoutdaily4pem_Y2.nc"   ! XIOS daily output file, year 2
     28character(20), parameter :: file1_yearly = "Xoutyearly4pem_Y1.nc"  ! XIOS yearly output file, year 1
     29character(20), parameter :: file2_yearly = "Xoutyearly4pem_Y2.nc"  ! XIOS yearly output file, year 2
     30character(256)           :: msg                                    ! Message for reading errors
     31integer                  :: fID, vID                               ! File and variable IDs for reading
     32
     33! INTERFACES
     34! ----------
    1435interface get_var
    1536    module procedure get_var_1d, get_var_2d, get_var_3d, get_var_4d
    1637end interface get_var
    1738
    18 !=======================================================================
    1939contains
    20 !=======================================================================
     40!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    2141
    2242!=======================================================================
    2343SUBROUTINE load_xios_data(ngrid,nslope,nsoil_PCM,nsol,h2ofrost_PCM,co2frost_PCM,ps_avg,tsurf_avg,tsurf_avg_y1,tsoil_avg,tsoil_ts,watersurf_density_avg,d_h2oice,d_co2ice, &
    2444                         ps_ts,q_h2o_ts,q_co2_ts,watersoil_density_ts,min_h2oice,min_co2ice)
    25 
     45!-----------------------------------------------------------------------
     46! NAME
     47!     load_xios_data
     48!
     49! DESCRIPTION
     50!     Reads yearly and daily XIOS data, computes frost and ice tendencies.
     51!
     52! AUTHORS & DATE
     53!     JB Clement, 2025
     54!
     55! NOTES
     56!
     57!-----------------------------------------------------------------------
     58
     59! DEPENDENCIES
     60! ------------
    2661use grid_conversion, only: lonlat2vect
    2762use soil,            only: do_soil
    … …  
    2964use metamorphism,    only: compute_frost
    3065
    31 implicit none
    32 
    33 ! Arguments
    34 !----------
    35 integer,                       intent(in) :: ngrid, nslope, nsoil_PCM, nsol
    36 real, dimension(ngrid,nslope), intent(in) :: h2ofrost_PCM, co2frost_PCM
    37 real, dimension(ngrid),                       intent(out) :: ps_avg
    38 real, dimension(ngrid,nslope),                intent(out) :: tsurf_avg, tsurf_avg_y1, watersurf_density_avg, d_h2oice, d_co2ice, min_h2oice, min_co2ice
    39 real, dimension(ngrid,nsoil_PCM,nslope),      intent(out) :: tsoil_avg
    40 real, dimension(ngrid,nsol),                  intent(out) :: ps_ts, q_h2o_ts, q_co2_ts
    41 real, dimension(ngrid,nsoil_PCM,nslope,nsol), intent(out) :: tsoil_ts, watersoil_density_ts
    42 
    43 ! Local variables
    44 !----------------
     66! DECLARATION
     67! -----------
     68implicit none
     69
     70! ARGUMENTS
     71! ---------
     72integer,                                      intent(in)  :: ngrid, nslope, nsoil_PCM, nsol ! Grid dimensions
     73real, dimension(ngrid,nslope),                intent(in)  :: h2ofrost_PCM, co2frost_PCM     ! PCM frost fields
     74real, dimension(ngrid),                       intent(out) :: ps_avg                         ! Average surface pressure
     75real, dimension(ngrid,nslope),                intent(out) :: tsurf_avg, tsurf_avg_y1        ! Surface temperature
     76real, dimension(ngrid,nslope),                intent(out) :: watersurf_density_avg          ! Water density
     77real, dimension(ngrid,nslope),                intent(out) :: d_h2oice, d_co2ice             ! Ice tendencies
     78real, dimension(ngrid,nslope),                intent(out) :: min_h2oice, min_co2ice         ! Ice minima
     79real, dimension(ngrid,nsoil_PCM,nslope),      intent(out) :: tsoil_avg                      ! Soil temperature
     80real, dimension(ngrid,nsol),                  intent(out) :: ps_ts, q_h2o_ts, q_co2_ts      ! Time series
     81real, dimension(ngrid,nsoil_PCM,nslope,nsol), intent(out) :: tsoil_ts, watersoil_density_ts ! Soil time series
     82
     83! LOCAL VARIABLES
     84! ---------------
    4585integer                               :: islope, isoil, isol, nlon, nlat
    4686real, dimension(:,:),     allocatable :: var_read_2d
    … …  
    5090real, dimension(ngrid,nslope,2)       :: min_h2operice, min_co2perice, min_h2ofrost, min_co2frost
    5191
    52 ! Code
    53 !-----
     92! CODE
     93! ----
    5494! Initialization
    5595min_h2operice = 0.
    … …  
    187227!=======================================================================
    188228SUBROUTINE get_timelen(filename,timelen)
    189 
     229!-----------------------------------------------------------------------
     230! NAME
     231!     get_timelen
     232!
     233! DESCRIPTION
     234!     Get time dimension length from a NetCDF file.
     235!
     236! AUTHORS & DATE
     237!     JB Clement, 2025
     238!
     239! NOTES
     240!
     241!-----------------------------------------------------------------------
     242
     243! DEPENDENCIES
     244! ------------
    190245use netcdf
    191246
    192 implicit none
    193 
    194 ! Arguments
    195 ! ---------
    196 character(*), intent(in)  :: filename
    197 integer,      intent(out) :: timelen
    198 
    199 ! Local variables
     247! DECLARATION
     248! -----------
     249implicit none
     250
     251! ARGUMENTS
     252! ---------
     253character(*), intent(in)  :: filename ! NetCDF filename
     254integer,      intent(out) :: timelen  ! Length of time dimension
     255
     256! LOCAL VARIABLES
    200257! ---------------
    201258integer :: ncid  ! File ID
    202259integer :: dimid ! Dimension ID
    203 integer :: ierr  ! Return codes
    204 
    205 ! Code
     260integer :: ierr  ! Return code
     261
     262! CODE
    206263! ----
    207264! Open the NetCDF file
    … …  
    237294
    238295!=======================================================================
    239 SUBROUTINE error_msg(ierr,typ,nam)
    240 
    241 implicit none
    242 
    243 integer,      intent(in) :: ierr !--- NetCDF error code
    244 character(*), intent(in) :: typ  !--- Type of operation
    245 character(*), intent(in) :: nam  !--- Field/File name
     296SUBROUTINE error_msg(ierr,typ,nam)
     297!-----------------------------------------------------------------------
     298! NAME
     299!     error_msg
     300!
     301! DESCRIPTION
     302!     Handle and report NetCDF errors.
     303!
     304! AUTHORS & DATE
     305!     JB Clement, 2025
     306!
     307! NOTES
     308!
     309!-----------------------------------------------------------------------
     310
     311! DECLARATION
     312! -----------
     313implicit none
     314
     315! ARGUMENTS
     316! ---------
     317integer,      intent(in) :: ierr ! NetCDF error code
     318character(*), intent(in) :: typ  ! Type of operation (inq, get, put, open, close)
     319character(*), intent(in) :: nam  ! Field/file name
     320
     321! CODE
     322! ----
    246323
    247324if (ierr == nf90_noerr) return
    … …  
    263340!=======================================================================
    264341SUBROUTINE get_var_1d(var,v)
    265 
    266 implicit none
    267 
    268 character(*), intent(in) :: var
    269 real, dimension(:), intent(out) :: v
    270 
     342!-----------------------------------------------------------------------
     343! NAME
     344!     get_var_1d
     345!
     346! DESCRIPTION
     347!     Read a 1D variable from open NetCDF file.
     348!
     349! AUTHORS & DATE
     350!     JB Clement, 2025
     351!
     352! NOTES
     353!
     354!-----------------------------------------------------------------------
     355
     356! DECLARATION
     357! -----------
     358implicit none
     359
     360! ARGUMENTS
     361! ---------
     362character(*),       intent(in)  :: var ! Variable name
     363real, dimension(:), intent(out) :: v   ! Output array
     364
     365! CODE
     366! ----
    271367call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
    272368call error_msg(nf90_get_var(fID,vID,v),"get",var)
    … …  
    277373!=======================================================================
    278374SUBROUTINE get_var_2d(var,v)
    279 
    280 implicit none
    281 
    282 character(*), intent(in) :: var
    283 real, dimension(:,:), intent(out) :: v
     375!-----------------------------------------------------------------------
     376! NAME
     377!     get_var_2d
     378!
     379! DESCRIPTION
     380!     Read a 2D variable from open NetCDF file.
     381!
     382! AUTHORS & DATE
     383!     JB Clement, 2025
     384!
     385! NOTES
     386!
     387!-----------------------------------------------------------------------
     388
     389! DECLARATION
     390! -----------
     391implicit none
     392
     393! ARGUMENTS
     394! ---------
     395character(*),         intent(in)  :: var ! Variable name
     396real, dimension(:,:), intent(out) :: v   ! Output array
     397
     398! CODE
     399! ----
    284400
    285401call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
    … …  
    291407!=======================================================================
    292408SUBROUTINE get_var_3d(var,v)
    293 
    294 implicit none
    295 
    296 character(*), intent(in) :: var
    297 real, dimension(:,:,:), intent(out) :: v
     409!-----------------------------------------------------------------------
     410! NAME
     411!     get_var_3d
     412!
     413! DESCRIPTION
     414!     Read a 3D variable from open NetCDF file.
     415!
     416! AUTHORS & DATE
     417!     JB Clement, 2025
     418!
     419! NOTES
     420!
     421!-----------------------------------------------------------------------
     422
     423! DECLARATION
     424! -----------
     425implicit none
     426
     427! ARGUMENTS
     428! ---------
     429character(*),           intent(in)  ::  var ! Variable name
     430real, dimension(:,:,:), intent(out) ::  v   ! Output array
     431
     432! CODE
     433! ----
    298434
    299435call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
    … …  
    305441!=======================================================================
    306442SUBROUTINE get_var_4d(var,v)
    307 
    308 implicit none
    309 
    310 character(*), intent(in) :: var
    311 real, dimension(:,:,:,:), intent(out) :: v
    312 
     443!-----------------------------------------------------------------------
     444! NAME
     445!     get_var_4d
     446!
     447! DESCRIPTION
     448!     Read a 4D variable from open NetCDF file.
     449!
     450! AUTHORS & DATE
     451!     JB Clement, 2025
     452!
     453! NOTES
     454!
     455!-----------------------------------------------------------------------
     456
     457! DECLARATION
     458! -----------
     459implicit none
     460
     461! ARGUMENTS
     462! ---------
     463character(*),             intent(in)  ::  var ! Variable name
     464real, dimension(:,:,:,:), intent(out) ::  v   ! Output array
     465
     466! CODE
     467! ----
    313468call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
    314469call error_msg(nf90_get_var(fID,vID,v),"get",var)
Note: See TracChangeset for help on using the changeset viewer.