Ignore:
Timestamp:
Dec 6, 2023, 4:02:06 PM (12 months ago)
Author:
jbclement
Message:

PEM:

  • Simplification of the algorithm managing the stopping criteria;
  • Complete rework of the ice management in the PEM (H2O & CO2);

    Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
    Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
    Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
    Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
    Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
    Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
    New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').

  • Renaming of variables/subroutines for clarity;
  • Some cleanings throughout the code;
  • Small updates in files of the deftank.

JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
1 added
2 deleted
18 edited
3 moved

Legend:

Unmodified
Added
Removed

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

    r3145 r3149  
    154154- Addition of a script in LMDZ.MARS/deftank/pem/ to launch a chained simulation of 1D PCM runs which follow, year by year, the orbital parameters (obliquity, eccentricity, Ls perihelion) given in a specified file.
    155155- Small changes to other files of the deftank directory (check and cosmetic).
     156
     157== 06/12/2023 == JBC
     158- Simplification of the algorithm managing the stopping criteria;
     159- Complete rework of the ice management in the PEM (H2O & CO2);
     160      > Subroutines to evolve the H2O and CO2 ice are now in the same module "evol_ice_mod.F90".
     161      > Tendencies are computed from the variation of "ice + frost" between the 2 PCM runs.
     162      > Evolving ice in the PEM is now called 'h2o_ice' or 'co2_ice' (not anymore in 'qsurf' and free of 'water_reservoir').
     163      > Default value 'ini_h2o_bigreservoir' (= 10 m) initializes the H2O ice of the first PEM run where there is 'watercap'. For the next PEM runs, initialization is done with the value kept in "startpem.nc". CO2 ice is taken from 'perennial_co2ice' of the PCM (paleoclimate flag must be true).
     164      > Simplification of the condition to compute the surface ice cover needed for the stopping criteria.
     165      > Frost ('qsurf') is not evolved by the PEM and given back to the PCM.
     166      > New default threshold value 'inf_h2oice_threshold' (= 2 m) to decide at the end of the PEM run if the H2O ice should be 'watercap' or not for the next PCM runs. If H2O ice cannot be 'watercap', then the remaining H2O ice is transferred to the frost ('qsurf').
     167- Renaming of variables/subroutines for clarity;
     168- Some cleanings throughout the code;
     169- Small updates in files of the deftank.

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

    r3148 r3149  
    1 MODULE compute_tendencies_slope_mod
     1MODULE compute_tend_mod
    22
    33implicit none
    … …  
    77!=======================================================================
    88
    9 SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,min_ice_Y2,tendencies_ice)
     9SUBROUTINE compute_tend(ngrid,nslope,min_ice,tendencies_ice)
    1010
    1111implicit none
    … …  
    1313!=======================================================================
    1414!
    15 !  Compute the tendencies of the evolution of water ice over the years
     15! Compute the tendencies of the evolution of water ice over the years
    1616!
    1717!=======================================================================
    … …  
    2020!   ----------
    2121!   INPUT
    22 integer,                       intent(in) :: ngrid, nslope ! # of grid points along longitude/latitude/ total
    23 real, dimension(ngrid,nslope), intent(in) :: min_ice_Y1    ! LON x LAT field : minimum of water ice at each point for the first year
    24 real, dimension(ngrid,nslope), intent(in) :: min_ice_Y2    ! LON x LAT field : minimum of water ice at each point for the second year
     22integer,                         intent(in) :: ngrid   ! # of grid points
     23integer,                         intent(in) :: nslope  ! # of subslopes
     24real, dimension(ngrid,nslope,2), intent(in) :: min_ice ! Minima of ice at each point for the PCM years
    2525
    2626!   OUTPUT
    27 real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perennial ice
     27real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! Difference between the minima = evolution of perennial ice
    2828!=======================================================================
     29! We compute the difference
     30tendencies_ice = min_ice(:,:,2) - min_ice(:,:,1)
    2931
    30 !  We compute the difference
    31 tendencies_ice = min_ice_Y2 - min_ice_Y1
    32 
    33 !  If the difference is too small; there is no evolution
     32! If the difference is too small, then there is no evolution
    3433where (abs(tendencies_ice) < 1.e-10) tendencies_ice = 0.
    3534
    36 END SUBROUTINE compute_tendencies_slope
     35END SUBROUTINE compute_tend
    3736
    38 END MODULE compute_tendencies_slope_mod
    39 
     37END MODULE compute_tend_mod

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

    r3143 r3149  
    33implicit none
    44
    5 integer, parameter                        :: nsoilmx_PEM = 68    ! number of layers in the PEM
    6 real, allocatable, dimension(:),     save :: layer_PEM           ! soil layer depths [m]
    7 real, allocatable, dimension(:),     save :: mlayer_PEM          ! soil mid-layer depths [m]
    8 real, allocatable, dimension(:,:,:), save :: TI_PEM              ! soil thermal inertia [SI]
    9 real, allocatable, dimension(:,:),   save :: inertiedat_PEM      ! soil thermal inertia saved as reference for current climate [SI]
     5integer, parameter                        :: nsoilmx_PEM = 68     ! number of layers in the PEM
     6real, allocatable, dimension(:),     save :: layer_PEM            ! soil layer depths [m]
     7real, allocatable, dimension(:),     save :: mlayer_PEM           ! soil mid-layer depths [m]
     8real, allocatable, dimension(:,:,:), save :: TI_PEM               ! soil thermal inertia [SI]
     9real, allocatable, dimension(:,:),   save :: inertiedat_PEM       ! soil thermal inertia saved as reference for current climate [SI]
    1010! Variables (FC: built in firstcall in soil.F)
    11 real, allocatable, dimension(:,:,:), save :: tsoil_PEM           ! sub-surface temperatures [K]
    12 real, allocatable, dimension(:,:),   save :: mthermdiff_PEM      ! (FC) mid-layer thermal diffusivity [SI]
    13 real, allocatable, dimension(:,:),   save :: thermdiff_PEM       ! (FC) inter-layer thermal diffusivity [SI]
    14 real, allocatable, dimension(:),     save :: coefq_PEM           ! (FC) q_{k+1/2} coefficients [SI]
    15 real, allocatable, dimension(:,:),   save :: coefd_PEM           ! (FC) d_k coefficients [SI]
    16 real, allocatable, dimension(:,:),   save :: alph_PEM            ! (FC) alpha_k coefficients [SI]
    17 real, allocatable, dimension(:,:),   save :: beta_PEM            ! beta_k coefficients [SI]
    18 real,                                save :: mu_PEM              ! mu coefficient [SI]
    19 real,                                save :: fluxgeo             ! Geothermal flux [W/m^2]
    20 real,                                save :: depth_breccia       ! Depth at which we have breccia [m]
    21 real,                                save :: depth_bedrock       ! Depth at which we have bedrock [m]
    22 integer,                             save :: index_breccia       ! last index of the depth grid before having breccia
    23 integer,                             save :: index_bedrock       ! last index of the depth grid before having bedrock
    24 logical                                   :: soil_pem            ! True by default, to run with the subsurface physic. Read in pem.def
    25 real, allocatable, dimension(:),     save :: water_reservoir     ! Large reserve of water [kg/m^2]
    26 real,                                save :: water_reservoir_nom ! Nominal value for the large reserve of water [kg/m^2]
    27 logical,                             save :: reg_thprop_dependp  ! thermal properites of the regolith vary with the surface pressure
     11real, allocatable, dimension(:,:,:), save :: tsoil_PEM            ! sub-surface temperatures [K]
     12real, allocatable, dimension(:,:),   save :: mthermdiff_PEM       ! (FC) mid-layer thermal diffusivity [SI]
     13real, allocatable, dimension(:,:),   save :: thermdiff_PEM        ! (FC) inter-layer thermal diffusivity [SI]
     14real, allocatable, dimension(:),     save :: coefq_PEM            ! (FC) q_{k+1/2} coefficients [SI]
     15real, allocatable, dimension(:,:),   save :: coefd_PEM            ! (FC) d_k coefficients [SI]
     16real, allocatable, dimension(:,:),   save :: alph_PEM             ! (FC) alpha_k coefficients [SI]
     17real, allocatable, dimension(:,:),   save :: beta_PEM             ! beta_k coefficients [SI]
     18real,                                save :: mu_PEM               ! mu coefficient [SI]
     19real,                                save :: fluxgeo              ! Geothermal flux [W/m^2]
     20real,                                save :: depth_breccia        ! Depth at which we have breccia [m]
     21real,                                save :: depth_bedrock        ! Depth at which we have bedrock [m]
     22integer,                             save :: index_breccia        ! last index of the depth grid before having breccia
     23integer,                             save :: index_bedrock        ! last index of the depth grid before having bedrock
     24logical                                   :: soil_pem             ! True by default, to run with the subsurface physic. Read in pem.def
     25real,                                save :: ini_h2o_bigreservoir ! Initial value for the large reservoir of water [kg/m^2]
     26logical,                             save :: reg_thprop_dependp   ! thermal properites of the regolith vary with the surface pressure
    2827
    2928!=======================================================================
    … …  
    3938
    4039allocate(layer_PEM(nsoilmx_PEM))
    41 allocate(mlayer_PEM(0:nsoilmx_PEM-1))
     40allocate(mlayer_PEM(0:nsoilmx_PEM - 1))
    4241allocate(TI_PEM(ngrid,nsoilmx_PEM,nslope))
    4342allocate(tsoil_PEM(ngrid,nsoilmx_PEM,nslope))
    44 allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM-1))
    45 allocate(thermdiff_PEM(ngrid,nsoilmx_PEM-1))
    46 allocate(coefq_PEM(0:nsoilmx_PEM-1))
    47 allocate(coefd_PEM(ngrid,nsoilmx_PEM-1))
    48 allocate(alph_PEM(ngrid,nsoilmx_PEM-1))
    49 allocate(beta_PEM(ngrid,nsoilmx_PEM-1))
     43allocate(mthermdiff_PEM(ngrid,0:nsoilmx_PEM - 1))
     44allocate(thermdiff_PEM(ngrid,nsoilmx_PEM - 1))
     45allocate(coefq_PEM(0:nsoilmx_PEM - 1))
     46allocate(coefd_PEM(ngrid,nsoilmx_PEM - 1))
     47allocate(alph_PEM(ngrid,nsoilmx_PEM - 1))
     48allocate(beta_PEM(ngrid,nsoilmx_PEM - 1))
    5049allocate(inertiedat_PEM(ngrid,nsoilmx_PEM))
    51 allocate(water_reservoir(ngrid))
    5250
    5351END SUBROUTINE ini_comsoil_h_PEM
    … …  
    7068if (allocated(beta_PEM)) deallocate(beta_PEM)
    7169if (allocated(inertiedat_PEM)) deallocate(inertiedat_PEM)
    72 if (allocated(water_reservoir)) deallocate(water_reservoir)
    7370
    7471END SUBROUTINE end_comsoil_h_PEM

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

    r3096 r3149  
    2525use time_evol_mod,    only: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, Max_iter_pem, &
    2626                          evol_orbit_pem, var_obl, var_ecc, var_lsp, convert_years
    27 use comsoil_h_pem,    only: soil_pem, fluxgeo, water_reservoir_nom, depth_breccia, depth_bedrock, reg_thprop_dependp
     27use comsoil_h_pem,    only: soil_pem, fluxgeo, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, reg_thprop_dependp
    2828use adsorption_mod,   only: adsorption_pem
    29 use glaciers_mod,     only: co2glaciersflow, h2oglaciersflow
     29use glaciers_mod,     only: co2_ice_flow, h2o_ice_flow
    3030use ice_table_mod,    only: icetable_equilibrium, icetable_dynamic
    3131use time_phylmdz_mod, only: ecritphy
    … …  
    7171var_ecc = .true.
    7272call getin('var_ecc',var_ecc)
    73 write(*,*) 'Does excentricity vary ?',var_ecc
     73write(*,*) 'Does eccentricity vary ?',var_ecc
    7474
    7575var_lsp = .true.
    … …  
    100100call getin('adsorption_pem',adsorption_pem)
    101101
    102 co2glaciersflow = .true.
    103 call getin('co2glaciersflow',co2glaciersflow)
     102co2_ice_flow = .true.
     103call getin('co2_ice_flow',co2_ice_flow)
    104104
    105 h2oglaciersflow = .true.
    106 call getin('h2oglaciersflow',h2oglaciersflow)
     105h2o_ice_flow = .true.
     106call getin('h2o_ice_flow',h2o_ice_flow)
    107107
    108108reg_thprop_dependp = .false.
    … …  
    156156endif
    157157
    158 water_reservoir_nom = 1.e4
    159 call getin('water_reservoir_nom',water_reservoir_nom)
     158ini_h2o_bigreservoir = 1.e4
     159call getin('ini_h2o_bigreservoir',ini_h2o_bigreservoir)
    160160
    161161END SUBROUTINE conf_pem

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

    r3143 r3149  
    1212real, parameter :: m_h2o = 18.01528e-3 ! Molecular weight of h2o (kg/mol)
    1313
    14 !     Coefficient for Clapeyron law for CO2 condensation temperature (Tco2 = beta/(alpha-log(vmr)),following James et al. 1992
     14! Coefficient for Clapeyron law for CO2 condensation temperature (Tco2 = beta/(alpha-log(vmr)),following James et al. 1992
    1515real, parameter :: alpha_clap_co2 = 23.3494 ! Uniteless, James et al. 1992
    1616real, parameter :: beta_clap_co2 = 3182.48  ! Kelvin, James et al. 1992
    1717
    18 !     Coefficient for Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
     18! Coefficient for Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005
    1919real, parameter :: alpha_clap_h2o = 28.9074 ! Uniteless, Murphy and Koop 2005
    2020real, parameter :: beta_clap_h2o = -6143.7  ! Kelvin, Murphy and Koop 2005
    2121
    22 !     Density of the regolith (Zent et al., 1995, Buhler and Piqueux 2021)     
     22! Density of the regolith (Zent et al., 1995, Buhler and Piqueux 2021)     
    2323real, parameter :: rho_regolith = 2000. ! kg/m^3
    2424
    25 !     Average  Thermal inertia of the surface, breccia, bedrock, following Mellon et al., 2000., Wood et al., 2008
     25! Average thermal inertia of the surface, breccia, bedrock, following Mellon et al., 2000., Wood et al., 2008
    2626real, parameter :: TI_regolith_avg = 250. ! Averaged of the observed thermal inertia for regolith following Mellon et al., 2000[SI]
    2727real, parameter :: TI_breccia = 750.      ! Thermal inertia of Breccia following Wood 2009 [SI]
    2828real, parameter :: TI_bedrock = 2300.     ! Thermal inertia of Bedrock following Wood 2009 [SI]
    2929
    30 !     Porosity of the soil
     30! Porosity of the soil
    3131real, parameter :: porosity = 0.4 ! porosity of the martian soil, correspond to the value for a random loose packing of monodiperse sphere (Scott, 1960)
    3232
    33 !     Stefan Boltzmann constant
     33! Stefan Boltzmann constant
    3434real, parameter :: sigmaB = 5.678e-8
    3535
    36 !     Latent heat of CO2
     36! Latent heat of CO2
    3737real, parameter :: Lco2 =  5.71e5 ! Pilorget and Forget 2016
    3838
    39 ! Conversion H2O/CO2 frost to perennial frost and vice versa
    40 real, parameter :: threshold_h2o_frost2perennial = 1000.  !~ 1 m
    41 real, parameter :: threshold_co2_frost2perennial = 16000. !~ 10 m
     39! Threshold to consider the amount of H2O ice as an infinite reservoir
     40real, parameter :: inf_h2oice_threshold = 2000.  !~ 1 m
    4241
    4342END MODULE constants_marspem_mod

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

    r3082 r3149  
    1 module glaciers_mod
    2        
    3  implicit none
    4  LOGICAL  co2glaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
    5  LOGICAL  h2oglaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
    6 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7 !!!
    8 !!! Purpose: Compute CO2 glacier flows
    9 !!!
    10 !!! Author: LL
    11 !!!
    12 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    13 
     1MODULE glaciers_mod
     2
     3implicit none
     4
     5logical :: co2_ice_flow ! True by default, to compute co2 ice flow. Read in run_PEM.def
     6logical :: h2o_ice_flow ! True by default, to compute h2o ice flow. Read in run_PEM.def
     7
     8!=======================================================================
    149contains
    15 
    16 subroutine co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
     10!=======================================================================
     11
     12SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
    1713
    1814!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    2622!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2723
    28 IMPLICIT NONE
     24implicit none
    2925
    3026! arguments
    … …  
    3228
    3329! Inputs:
    34       INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
    35       REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
    36       REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
    37       REAL,INTENT(IN) :: ps_GCM(ngrid,timelen)      ! Physical x Time field: surface pressure given by the GCM [Pa]
    38       REAL,INTENT(IN) :: global_ave_ps_GCM          ! Global averaged surface pressure from the GCM [Pa]
    39       REAL,INTENT(IN) :: global_ave_ps_PEM          ! global averaged surface pressure during the PEM iteration [Pa]
     30      integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
     31      real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
     32      real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
     33      real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
     34      real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
     35      real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
     36      real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
    4037     
    4138! Ouputs:
    42       REAL,INTENT(INOUT) :: co2ice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    43       REAL,INTENT(INOUT) :: flag_co2flow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
    44       REAL,INTENT(INOUT) :: flag_co2flow_mesh(ngrid)  ! same but within the mesh
    45 
     39      real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     40      real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
     41      real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
    4642
    4743! Local
    48       REAL :: Tcond(ngrid,nslope) ! Physical field: CO2 condensation temperature [K]
    49       REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
     44      real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
     45      real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    5046
    5147!-----------------------------
    52       call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
    53 
    54       call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
    55 
    56       call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
    57    RETURN   
    58 end subroutine
    59 
    60 
    61 
    62 
    63 
    64 subroutine h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
     48write(*,*) "Flow of CO2 glacier"
     49   
     50call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
     51call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
     52call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
     53
     54END SUBROUTINE flow_co2glaciers
     55
     56!=======================================================================
     57
     58SUBROUTINE flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
    6559
    6660!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    7468!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7569
    76 
    77 IMPLICIT NONE
     70implicit none
    7871
    7972! arguments
    … …  
    8174
    8275! Inputs:
    83       INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
    84       REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
    85       REAL,INTENT(IN) :: Tice(ngrid,nslope) ! Ice Temperature [K]
     76      integer,intent(in) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
     77      real,intent(in) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
     78      real,intent(in) :: Tice(ngrid,nslope) ! Ice Temperature [K]
    8679! Ouputs:
    87       REAL,INTENT(INOUT) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
    88       REAL,INTENT(INOUT) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
    89       REAL,INTENT(INOUT) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
     80      real,intent(inout) :: h2oice(ngrid,nslope) ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
     81      real,intent(inout) :: flag_h2oflow(ngrid,nslope) ! flag to see if there is flow on the subgrid slopes
     82      real,intent(inout) :: flag_h2oflow_mesh(ngrid)  ! same but within the mesh
    9083! Local
    91       REAL :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
     84      real :: hmax(ngrid,nslope)  ! Physical x Slope field: maximum thickness for co2  glacier before flow
    9285
    9386!-----------------------------
    94 
    95       call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
    96       call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
    97 
    98    RETURN
    99 end subroutine
    100 
    101 
    102 subroutine compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
     87write(*,*) "Flow of H2O glaciers"
     88
     89call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
     90call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
     91
     92END SUBROUTINE flow_h2oglaciers
     93
     94!=======================================================================
     95
     96SUBROUTINE compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,name_ice,hmax)
    10397
    10498use abort_pem_mod, only: abort_pem
    … …  
    118112!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    119113
    120    IMPLICIT NONE
     114implicit none
    121115
    122116! arguments
    … …  
    124118
    125119! Inputs
    126       INTEGER,INTENT(IN) :: ngrid,nslope ! # of grid points and subslopes
    127       INTEGER,INTENT(IN) :: iflat        ! index of the flat subslope
    128       REAL,INTENT(IN) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
    129       REAL,INTENT(IN) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
    130       character(len=3), INTENT(IN) :: name_ice ! Nature of the ice
     120      integer,intent(in) :: ngrid,nslope ! # of grid points and subslopes
     121      integer,intent(in) :: iflat        ! index of the flat subslope
     122      real,intent(in) :: def_slope_mean(nslope) ! Slope field: Values of the subgrid slope angles [deg]
     123      real,intent(in) :: Tice(ngrid,nslope)     ! Physical field:  ice temperature [K]
     124      character(len=3), intent(in) :: name_ice ! Nature of the ice
    131125! Outputs
    132       REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
     126      real,intent(out) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum  thickness before flaw [m]
    133127! Local
    134128      DOUBLE PRECISION :: 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
    135       REAL :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
    136       INTEGER :: ig,islope ! loop variables
    137       REAL :: slo_angle
     129      real :: rho(ngrid,nslope) ! co2 ice density [kg/m^3]
     130      integer :: ig,islope ! loop variables
     131      real :: slo_angle
    138132
    139133! 1. Compute rho
    … …  
    167161       ENDDO
    168162    ENDDO
    169 RETURN
    170 
    171 end subroutine
    172 
    173 
    174 
    175 
    176 subroutine transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
     163END SUBROUTINE compute_hmaxglaciers
     164
     165!=======================================================================
     166
     167SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
    177168!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    178169!!!
    … …  
    196187
    197188! Inputs
    198       INTEGER, INTENT(IN) :: ngrid,nslope               !# of physical points and subslope
    199       INTEGER, INTENT(IN) :: iflat                      ! index of the flat subslope
    200       REAL, INTENT(IN) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
    201       REAL, INTENT(IN) :: def_slope_mean(nslope)        ! values of the subgrid slopes
    202       REAL, INTENT(IN) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
    203       REAL, INTENT(IN) :: Tice(ngrid,nslope)            ! Ice temperature[K]
    204       character(len=3), INTENT(IN) :: name_ice              ! Nature of the ice
     189      integer, intent(in) :: ngrid,nslope               !# of physical points and subslope
     190      integer, intent(in) :: iflat                      ! index of the flat subslope
     191      real, intent(in) :: subslope_dist(ngrid,nslope)   ! Distribution of the subgrid slopes within the mesh
     192      real, intent(in) :: def_slope_mean(nslope)        ! values of the subgrid slopes
     193      real, intent(in) :: hmax(ngrid,nslope)            ! maximum height of the  glaciers before initiating flow [m]
     194      real, intent(in) :: Tice(ngrid,nslope)            ! Ice temperature[K]
     195      character(len=3), intent(in) :: name_ice              ! Nature of the ice
    205196
    206197! Outputs
    207       REAL, INTENT(INOUT) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
    208       REAL, INTENT(INOUT) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
    209       REAL, INTENT(INOUT) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
     198      real, intent(inout) :: qice(ngrid,nslope)      ! CO2 in the subslope [kg/m^2]
     199      real, intent(inout) :: flag_flow(ngrid,nslope) ! boolean to check if there is flow on a subgrid slope
     200      real, intent(inout) :: flag_flowmesh(ngrid)    ! boolean to check if there is flow in the mesh
    210201! Local
    211       INTEGER ig,islope ! loop
    212       REAL rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
    213       INTEGER iaval ! ice will be transfered here
     202      integer ig,islope ! loop
     203      real rho(ngrid,nslope) ! density of ice, temperature dependant [kg/m^3]
     204      integer iaval ! ice will be transfered here
    214205
    215206! 0. Compute rho
    … …  
    268259        ENDDO !islope
    269260       ENDDO !ig
    270 RETURN
    271 end subroutine
    272 
    273      subroutine computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_GCM,global_ave_ps_GCM,global_ave_ps_PEM,Tcond)
     261END SUBROUTINE
     262
     263!=======================================================================
     264
     265SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
    274266!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    275267!!!
    … …  
    288280
    289281! INPUT
    290       INTEGER,INTENT(IN) :: timelen, ngrid,nslope ! # of timesample,physical points, subslopes
    291       REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
    292       REAL,INTENT(IN) :: ps_GCM(ngrid,timelen) ! Physical points x times field: surface pressure in the GCM [Pa]
    293       REAL,INTENT(IN) :: global_ave_ps_GCM ! Global averaged surfacepressure in the GCM [Pa]
    294       REAL, INTENT(IN) :: global_ave_ps_PEM ! Global averaged surface pressure computed during the PEM iteration
     282integer,                        intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes
     283real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
     284real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Physical points x times field: surface pressure in the PCM [Pa]
     285real,                           intent(in) :: global_avg_ps_PCM      ! Global averaged surfacepressure in the PCM [Pa]
     286real,                           intent(in) :: global_avg_ps_PEM      ! Global averaged surface pressure computed during the PEM iteration
     287
    295288! OUTPUT
    296       REAL,INTENT(OUT) :: Tcond(ngrid,nslope) ! Physical points : condensation temperature of CO2, yearly averaged
     289real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged
    297290
    298291! LOCAL
    299 
    300       INTEGER :: ig,it,islope ! for loop
    301       REAL :: ave ! intermediate to compute average
     292integer :: ig, it ! For loop
     293real    :: ave    ! Intermediate to compute average
    302294
    303295!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    304 
    305 
    306       DO ig = 1,ngrid
    307         ave = 0
    308         DO it = 1,timelen
    309            ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(ig,it)*ps_GCM(ig,it)*global_ave_ps_GCM/global_ave_ps_PEM/100))
    310         ENDDO
    311         DO islope = 1,nslope
    312           Tcond(ig,islope) = ave/timelen
    313         ENDDO
    314       ENDDO
    315 RETURN
    316 
    317 
    318 end subroutine
    319 end module
     296do ig = 1,ngrid
     297    ave = 0
     298    do it = 1,timelen
     299        ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*global_avg_ps_PCM/global_avg_ps_PEM/100))
     300    enddo
     301    Tcond(ig,:) = ave/timelen
     302enddo
     303
     304END SUBROUTINE computeTcondCO2
     305
     306END MODULE glaciers_mod

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

    r3097 r3149  
    77!=======================================================================
    88
    9 SUBROUTINE info_PEM(year_iter,criterion_stop,i_myear,n_myear)
     9SUBROUTINE info_PEM(year_iter,stopPEM,i_myear,n_myear)
    1010
    1111!=======================================================================
    … …  
    2222
    2323!----- Arguments
    24 integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
     24integer, intent(in) :: year_iter, stopPEM ! # of year and reason to stop
    2525integer, intent(in) :: i_myear, n_myear          ! Current simulated Martian year and maximum number of Martian years to be simulated
    2626
    … …  
    4343    endif
    4444    open(1,file = 'info_PEM.txt',status = "old",position = "append",action = "write")
    45     write(1,*) year_iter, i_myear, criterion_stop
     45    write(1,*) year_iter, i_myear, stopPEM
    4646    close(1)
    4747else

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

    r3148 r3149  
    1 MODULE interpolate_TIPEM_TIGCM_mod
     1MODULE interpol_TI_PEM2PCM_mod
    22
    33implicit none
    … …  
    77!=======================================================================
    88
    9 SUBROUTINE interpolate_TIPEM_TIGCM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_PEM,TI_GCM)
     9SUBROUTINE interpol_TI_PEM2PCM(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PEM,TI_PCM)
    1010
    1111!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1212!!!
    13 !!! Purpose: Transfer the thermal inertia from the PEM vertical  grid to the GCM vertical grid
     13!!! Purpose: Transfer the thermal inertia from the PEM vertical grid to the PCM vertical grid
    1414!!!
    1515!!!
    … …  
    2727integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
    2828integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
    29 integer,                                 intent(in) :: nsoil_GCM ! # of soil layers in the GCM
     29integer,                                 intent(in) :: nsoil_PCM ! # of soil layers in the GCM
    3030real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: TI_PEM    ! Thermal inertia in the PEM vertical grid [J/m^2/K/s^{1/2}]
    3131
    32 real, dimension(ngrid,nsoil_GCM,nslope), intent(inout) :: TI_GCM ! Thermal inertia in the GCM vertical grid [J/m^2/K/s^{1/2}]
     32real, dimension(ngrid,nsoil_PCM,nslope), intent(inout) :: TI_PCM ! Thermal inertia in the PCM vertical grid [J/m^2/K/s^{1/2}]
    3333
    3434!----- Code
    35 TI_GCM(:,:,:) = TI_PEM(:,:nsoil_GCM,:)
     35TI_PCM = TI_PEM(:,:nsoil_PCM,:)
    3636
    37 END SUBROUTINE interpolate_TIPEM_TIGCM
     37END SUBROUTINE interpol_TI_PEM2PCM
    3838
    39 END MODULE interpolate_TIPEM_TIGCM_mod
     39END MODULE interpol_TI_PEM2PCM_mod

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

    r3050 r3149  
    55!!!
    66!!!
    7 !!! Author: LL, inspired by iostart from the GCM
     7!!! Author: LL, inspired by iostart from the PCM
    88!!!
    99!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

    r3076 r3149  
    1313real, dimension(:), allocatable, save :: yearlask   ! year before present from Laskar et al. Tab
    1414real, dimension(:), allocatable, save :: obllask    ! obliquity    [deg]
    15 real, dimension(:), allocatable, save :: ecclask    ! excentricity [deg]
     15real, dimension(:), allocatable, save :: ecclask    ! eccentricity [deg]
    1616real, dimension(:), allocatable, save :: lsplask    ! ls perihelie [deg]
    1717integer,                         save :: last_ilask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year

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

    r3096 r3149  
    1919!=======================================================================
    2020!
    21 ! Purpose: Read in the data_PCM_Yr*.nc the number of time step
     21! Purpose: Read in the data_PCM_Y*.nc the number of time steps
    2222!
    2323! Author: RV

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

    r3082 r3149  
    9999
    100100call getin('max_change_ecc',max_change_ecc)
    101 max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ex < 1.
    102 min_ecc = max(e_elips - max_change_ecc,0.) ! ex >= 0.
     101max_ecc = min(e_elips + max_change_ecc,1. - 1.e-10) ! ecc < 1.
     102min_ecc = max(e_elips - max_change_ecc,0.) ! ecc >= 0.
    103103write(*,*) 'Maximum eccentricity accepted =', max_ecc
    104104write(*,*) 'Minimum eccentricity accepted =', min_ecc

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

    r3143 r3149  
    1313! II  Run
    1414!    II_a Update pressure, ice and tracers
    15 !    II_b Evolution of the ice
    16 !    II_c CO2 & H2O glaciers flows
     15!    II_b Evolution of ice
     16!    II_c Flow of glaciers
    1717!    II_d Update surface and soil temperatures
    1818!    II_e Outputs
    … …  
    2828PROGRAM pem
    2929
    30 use phyetat0_mod,                 only: phyetat0
    31 use phyredem,                     only: physdem0, physdem1
    32 use netcdf,                       only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
    33 use turb_mod,                     only: q2, wstar
    34 use comslope_mod,                 only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
    35 use logic_mod,                    only: iflag_phys
    36 use mod_const_mpi,                only: COMM_LMDZ
     30use phyetat0_mod,                only: phyetat0
     31use phyredem,                    only: physdem0, physdem1
     32use netcdf,                      only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
     33use turb_mod,                    only: q2, wstar
     34use comslope_mod,                only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
     35use logic_mod,                   only: iflag_phys
     36use mod_const_mpi,               only: COMM_LMDZ
    3737use infotrac
    38 use geometry_mod,                 only: latitude_deg
    39 use conf_pem_mod,                 only: conf_pem
    40 use pemredem,                     only: pemdem0, pemdem1
    41 use glaciers_mod,                 only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
    42 use criterion_pem_stop_mod,       only: criterion_waterice_stop, criterion_co2_stop
    43 use constants_marspem_mod,        only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
    44                                         m_noco2, threshold_h2o_frost2perennial, threshold_co2_frost2perennial
    45 use evol_co2_ice_s_mod,           only: evol_co2_ice_s
    46 use evol_h2o_ice_s_mod,           only: evol_h2o_ice_s
    47 use comsoil_h_PEM,                only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
    48                                         TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
    49                                         tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
    50                                         fluxgeo,                          & ! Geothermal flux for the PEM and PCM
    51                                         water_reservoir                     ! Water ressources
    52 use adsorption_mod,               only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
    53                                         ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
    54                                         co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
    55 use time_evol_mod,                only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
    56 use orbit_param_criterion_mod,    only: orbit_param_criterion
    57 use recomp_orb_param_mod,         only: recomp_orb_param
    58 use ice_table_mod,                only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
    59                                         ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
    60 use soil_thermalproperties_mod,   only: update_soil_thermalproperties
    61 use time_phylmdz_mod,             only: daysec, dtphys, day_end
    62 use abort_pem_mod,                only: abort_pem
    63 use soil_settings_PEM_mod,        only: soil_settings_PEM
    64 use compute_tendencies_slope_mod, only: compute_tendencies_slope
    65 use info_PEM_mod,                 only: info_PEM
    66 use interpolate_TIPEM_TIGCM_mod,  only: interpolate_TIPEM_TIGCM
    67 use nb_time_step_PCM_mod,         only: nb_time_step_PCM
    68 use pemetat0_mod,                 only: pemetat0
    69 use read_data_PCM_mod,            only: read_data_PCM
    70 use recomp_tend_co2_slope_mod,    only: recomp_tend_co2_slope
    71 use soil_pem_compute_mod,         only: soil_pem_compute
    72 use writediagpem_mod,             only: writediagpem
     38use geometry_mod,                only: latitude_deg
     39use conf_pem_mod,                only: conf_pem
     40use pemredem,                    only: pemdem0, pemdem1
     41use glaciers_mod,                only: flow_co2glaciers, flow_h2oglaciers, co2_ice_flow, h2o_ice_flow
     42use stopping_crit_mod,           only: stopping_crit_h2o_ice, stopping_crit_co2
     43use constants_marspem_mod,       only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
     44                                       m_noco2, inf_h2oice_threshold
     45use evol_ice_mod,                only: evol_co2_ice, evol_h2o_ice
     46use comsoil_h_PEM,               only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
     47                                       TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
     48                                       tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
     49                                       fluxgeo                             ! Geothermal flux for the PEM and PCM
     50use adsorption_mod,              only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
     51                                       ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
     52                                       co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
     53use time_evol_mod,               only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
     54use orbit_param_criterion_mod,   only: orbit_param_criterion
     55use recomp_orb_param_mod,        only: recomp_orb_param
     56use ice_table_mod,               only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
     57                                       ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
     58use soil_thermalproperties_mod,  only: update_soil_thermalproperties
     59use time_phylmdz_mod,            only: daysec, dtphys, day_end
     60use abort_pem_mod,               only: abort_pem
     61use soil_settings_PEM_mod,       only: soil_settings_PEM
     62use compute_tend_mod,            only: compute_tend
     63use info_PEM_mod,                only: info_PEM
     64use interpol_TI_PEM2PCM_mod, only: interpol_TI_PEM2PCM
     65use nb_time_step_PCM_mod,        only: nb_time_step_PCM
     66use pemetat0_mod,                only: pemetat0
     67use read_data_PCM_mod,           only: read_data_PCM
     68use recomp_tend_co2_slope_mod,   only: recomp_tend_co2_slope
     69use soil_pem_compute_mod,        only: soil_pem_compute
     70use writediagpem_mod,            only: writediagpem
    7371
    7472#ifndef CPP_STD
    … …  
    7775                                  albedodat, zmea, zstd, zsig, zgam, zthe,     &
    7876                                  hmons, summit, base,albedo_h2o_frost,        &
    79                                   frost_albedo_threshold, emissiv, watercaptag, perennial_co2ice, &
    80                                   emisice, albedice
     77                                  frost_albedo_threshold, emissiv, watercaptag, perennial_co2ice, emisice, albedice
    8178    use dimradmars_mod,     only: totcloudfrac, albedo
    8279    use dust_param_mod,     only: tauscaling
    83     use tracer_mod,         only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
     80    use tracer_mod,         only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
    8481    use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
    8582    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
    86     use paleoclimate_mod,   only: albedo_perennialco2
     83    use paleoclimate_mod,   only: albedo_perennialco2, paleoclimate
    8784    use comcstfi_h,         only: pi, rad, g, cpp, mugaz, r
    8885    use surfini_mod,        only: surfini
    … …  
    127124integer            :: day_ini      ! First day of the simulation
    128125real               :: pday         ! Physical day
    129 real               :: time_phys    ! Same as PCM
    130 real               :: ptimestep    ! Same as PCM
    131 real               :: ztime_fin    ! Same as PCM
     126real               :: time_phys    ! Same as in PCM
     127real               :: ptimestep    ! Same as in PCM
     128real               :: ztime_fin    ! Same as in PCM
    132129
    133130! Variables to read start.nc
    … …  
    140137real, dimension(:,:,:), allocatable :: q             ! champs advectes
    141138real, dimension(ip1jmp1)            :: ps            ! pression au sol
    142 real, dimension(:),     allocatable :: ps_start_GCM  ! (ngrid) pression  au sol
    143 real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées
     139real, dimension(:),     allocatable :: ps_start_PCM  ! (ngrid) surface pressure
     140real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure
    144141real, dimension(ip1jmp1,llm)        :: masse         ! masse d'air
    145142real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
    … …  
    161158
    162159! Variables for h2o_ice evolution
    163 real                                :: ini_surf_h2o         ! Initial surface of sublimating h2o ice
    164 real                                :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa]
    165 real                                :: global_ave_press_old ! constant: Global average pressure of initial/previous time step
    166 real                                :: global_ave_press_new ! constant: Global average pressure of current time step
    167 real, dimension(:,:),   allocatable :: zplev_new            ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
    168 real, dimension(:,:),   allocatable :: zplev_gcm            ! same but retrieved from the PCM [kg/m^2]
     160real, dimension(:,:),   allocatable :: h2o_ice              ! h2o ice in the PEM
     161real, dimension(:,:),   allocatable :: h2o_ice_ini          ! Initial amount of h2o ice in the PEM
     162real, dimension(:,:,:), allocatable :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
     163real                                :: h2o_surf_ini         ! Initial surface of sublimating h2o ice
     164real                                :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
     165real                                :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
     166real                                :: global_avg_press_new ! constant: Global average pressure of current time step
     167real, dimension(:,:),   allocatable :: zplev_new            ! Physical x Atmospheric field: mass of the atmospheric  layers in the pem at current time step [kg/m^2]
     168real, dimension(:,:),   allocatable :: zplev_PCM            ! same but retrieved from the PCM [kg/m^2]
    169169real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
    170170real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
    171 logical                             :: STOPPING_water       ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?
    172 logical                             :: STOPPING_1_water     ! Logical: is there still water ice to sublimate?
    173 logical                             :: STOPPING_pressure    ! Logical: is the criterion (% of change in the surface pressure) reached?
    174 integer                             :: criterion_stop       ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
     171integer                             :: stopPEM              ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
    175172real, save                          :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
    176173real, dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
    … …  
    181178
    182179! Variables for co2_ice evolution
    183 real                              :: ini_surf_co2         ! Initial surface of sublimating co2 ice
    184 logical                           :: STOPPING_co2         ! Logical: is the criterion (% of change in the surface of sublimating co2 ice) reached?
    185 real, dimension(:,:), allocatable :: vmr_co2_gcm          ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
    186 real, dimension(:,:), allocatable :: vmr_co2_pem_phys     ! Physics x Times  co2 volume mixing ratio used in the PEM
    187 real, dimension(:,:), allocatable :: q_co2_PEM_phys       ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
    188 real, dimension(:,:), allocatable :: co2frost_ini         ! Initial amount of co2 frost (at the start of the PEM run)
    189 real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run)
     180real, dimension(:,:),   allocatable :: co2_ice          ! co2 ice in the PEM
     181real, dimension(:,:),   allocatable :: co2_ice_ini      ! Initial amount of co2 ice in the PEM
     182real, dimension(:,:,:), allocatable :: min_co2_ice      ! Minimum of co2 ice at each point for the first year [kg/m^2]
     183real                                :: co2_surf_ini     ! Initial surface of sublimating co2 ice
     184real, dimension(:,:),   allocatable :: vmr_co2_PCM      ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
     185real, dimension(:,:),   allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM
     186real, dimension(:,:),   allocatable :: q_co2_PEM_phys   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
    190187
    191188! Variables for slopes
    192 real, dimension(:,:),   allocatable :: min_co2_ice_1          ! ngrid field: minimum of co2 ice at each point for the first year [kg/m^2]
    193 real, dimension(:,:),   allocatable :: min_co2_ice_2          ! ngrid field: minimum of co2 ice at each point for the second year [kg/m^2]
    194 real, dimension(:,:),   allocatable :: min_h2o_ice_1          ! ngrid field: minimum of water ice at each point for the first year [kg/m^2]
    195 real, dimension(:,:),   allocatable :: min_h2o_ice_2          ! ngrid field: minimum of water ice at each point for the second year [kg/m^2]
    196 real, dimension(:,:,:), allocatable :: co2_ice_GCM            ! Physics x NSLOPE x Times field: co2 ice given by the PCM  [kg/m^2]
    197 real, dimension(:,:),   allocatable :: initial_co2_ice_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
    198 real, dimension(:,:),   allocatable :: initial_h2o_ice        ! physical point field: Logical array indicating if there is water ice at initial state
    199 real, dimension(:,:),   allocatable :: initial_co2_ice        ! physical point field: Logical array indicating if there is co2 ice at initial state
    200 real, dimension(:,:),   allocatable :: tendencies_co2_ice     ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
    201 real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM
    202 real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical point x slope  field: Tendency of evolution of perennial h2o ice
    203 real, dimension(:,:),   allocatable :: flag_co2flow           ! (ngrid,nslope): Flag where there is a CO2 glacier flow
    204 real, dimension(:),     allocatable :: flag_co2flow_mesh      ! (ngrid)       : Flag where there is a CO2 glacier flow
    205 real, dimension(:,:),   allocatable :: flag_h2oflow           ! (ngrid,nslope): Flag where there is a H2O glacier flow
    206 real, dimension(:),     allocatable :: flag_h2oflow_mesh      ! (ngrid)       : Flag where there is a H2O glacier flow
     189real, dimension(:,:,:), allocatable :: co2_ice_PCM        ! Physics x NSLOPE x Times field: co2 ice given by the PCM [kg/m^2]
     190real, dimension(:,:),   allocatable :: co2_ice_ini_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
     191real, dimension(:,:),   allocatable :: initial_h2o_ice    ! physical point field: Logical array indicating if there is water ice at initial state
     192real, dimension(:,:),   allocatable :: tend_co2_ice       ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
     193real, dimension(:,:),   allocatable :: tend_co2_ice_ini   ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM
     194real, dimension(:,:),   allocatable :: tend_h2o_ice       ! physical point x slope  field: Tendency of evolution of perennial h2o ice
     195real, dimension(:,:),   allocatable :: flag_co2flow       ! (ngrid,nslope): Flag where there is a CO2 glacier flow
     196real, dimension(:),     allocatable :: flag_co2flow_mesh  ! (ngrid)       : Flag where there is a CO2 glacier flow
     197real, dimension(:,:),   allocatable :: flag_h2oflow       ! (ngrid,nslope): Flag where there is a H2O glacier flow
     198real, dimension(:),     allocatable :: flag_h2oflow_mesh  ! (ngrid)       : Flag where there is a H2O glacier flow
    207199
    208200! Variables for surface and soil
    209 real, dimension(:,:),     allocatable :: tsurf_ave                          ! Physic x SLOPE field : Averaged Surface Temperature [K]
    210 real, dimension(:,:,:),   allocatable :: tsoil_ave                          ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
    211 real, dimension(:,:,:),   allocatable :: tsurf_GCM_timeseries               ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
    212 real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    213 real, dimension(:,:,:,:), allocatable :: tsoil_GCM_timeseries               ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    214 real, dimension(:,:),     allocatable :: tsurf_ave_yr1                      ! Physic x SLOPE field : Averaged Surface Temperature of first call of the PCM [K]
     201real, dimension(:,:),     allocatable :: tsurf_ave                          ! Physic x SLOPE field: Averaged Surface Temperature [K]
     202real, dimension(:,:,:),   allocatable :: tsoil_ave                          ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
     203real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries               ! ngrid x SLOPE XTULES field: Surface Temperature in timeseries [K]
     204real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SLOPE XTULES field: Non averaged Soil Temperature [K]
     205real, dimension(:,:,:,:), allocatable :: tsoil_PCM_timeseries               ! IG x SLOPE XTULES field: Non averaged Soil Temperature [K]
     206real, dimension(:,:),     allocatable :: tsurf_avg_yr1                      ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
    215207real, dimension(:,:),     allocatable :: TI_locslope                        ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
    216208real, dimension(:,:),     allocatable :: Tsoil_locslope                     ! Physic x Soil: intermediate when computing Tsoil [K]
    217 real, dimension(:),       allocatable :: Tsurf_locslope                     ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
     209real, dimension(:),       allocatable :: Tsurf_locslope                     ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
    218210real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries       ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
    219211real, dimension(:,:),     allocatable :: watersurf_density_ave              ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
    220212real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries   ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
    221213real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_ave          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
    222 real, dimension(:,:),     allocatable :: Tsurfave_before_saved              ! Surface temperature saved from previous time step [K]
     214real, dimension(:,:),     allocatable :: Tsurfavg_before_saved              ! Surface temperature saved from previous time step [K]
    223215real, dimension(:),       allocatable :: delta_co2_adsorbded                ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    224216real, dimension(:),       allocatable :: delta_h2o_adsorbded                ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
    … …  
    236228integer         :: n_myear       ! Maximum number of Martian years of the chained simulations
    237229real            :: timestep      ! timestep [s]
    238 real            :: watercap_sum  ! total mass of water cap [kg/m^2]
    239 real            :: water_sum     ! total mass of water in the mesh [kg/m^2]
    240230
    241231#ifdef CPP_STD
    … …  
    341331!------------------------
    342332! I_b.1 Read start_evol.nc
    343 allocate(ps_start_GCM(ngrid))
     333allocate(ps_start_PCM(ngrid))
    344334#ifndef CPP_1D
    345335    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
    346336
    347     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
     337    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
    348338
    349339    call iniconst !new
    … …  
    359349    status = nf90_close(ncid)
    360350
    361     call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, &
    362                    rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
     351    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
    363352#else
    364     ps_start_GCM(1) = ps(1)
     353    ps_start_PCM(1) = ps(1)
    365354#endif
    366355
    … …  
    387376#ifndef CPP_STD
    388377    call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
    389                   tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,       &
     378                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,          &
    390379                  watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
    391380
    … …  
    414403    allocate(inertiesoil(ngrid,nsoilmx,1))
    415404    call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
    416                   tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,            &
    417                   qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,     &
     405                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,               &
     406                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,        &
    418407                  rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
    419408    call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
    … …  
    422411    call ini_comslope_h(ngrid,1)
    423412
    424     qsurf(:,:,1) = qsurf_read_generic(:,:)
    425     tsurf(:,1) = tsurf_read_generic(:)
    426     tsoil(:,:,1) = tsoil_read_generic(:,:)
    427     emis(:,1) = emis_read_generic(:)
     413    qsurf(:,:,1) = qsurf_read_generic
     414    tsurf(:,1) = tsurf_read_generic
     415    tsoil(:,:,1) = tsoil_read_generic
     416    emis(:,1) = emis_read_generic
    428417    watercap(:,1) = 0.
    429418    watercaptag(:) = .false.
    430419    albedo(:,1,1) = albedo_read_generic(:,1)
    431420    albedo(:,2,1) = albedo_read_generic(:,2)
    432     inertiesoil(:,:,1) = inertiedat(:,:)
     421    inertiesoil(:,:,1) = inertiedat
    433422
    434423    if (nslope == 1) then
    … …  
    440429
    441430    ! Remove unphysical values of surface tracer
    442     qsurf(:,:,1) = qsurf_read_generic(:,:)
     431    qsurf(:,:,1) = qsurf_read_generic
    443432    where (qsurf < 0.) qsurf = 0.
    444433#endif
    … …  
    472461allocate(flag_h2oflow_mesh(ngrid))
    473462
    474 flag_co2flow(:,:) = 0
    475 flag_co2flow_mesh(:) = 0
    476 flag_h2oflow(:,:) = 0
    477 flag_h2oflow_mesh(:) = 0
     463flag_co2flow = 0
     464flag_co2flow_mesh = 0
     465flag_h2oflow = 0
     466flag_h2oflow_mesh = 0
    478467
    479468!------------------------
    … …  
    486475allocate(tsoil_ave(ngrid,nsoilmx,nslope))
    487476allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
    488 allocate(vmr_co2_gcm(ngrid,timelen))
     477allocate(vmr_co2_PCM(ngrid,timelen))
    489478allocate(ps_timeseries(ngrid,timelen))
    490 allocate(min_co2_ice_1(ngrid,nslope))
    491 allocate(min_h2o_ice_1(ngrid,nslope))
    492 allocate(min_co2_ice_2(ngrid,nslope))
    493 allocate(min_h2o_ice_2(ngrid,nslope))
    494 allocate(tsurf_ave_yr1(ngrid,nslope))
     479allocate(min_co2_ice(ngrid,nslope,2))
     480allocate(min_h2o_ice(ngrid,nslope,2))
     481allocate(tsurf_avg_yr1(ngrid,nslope))
    495482allocate(tsurf_ave(ngrid,nslope))
    496 allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
    497 allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
     483allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
     484allocate(tsoil_PCM_timeseries(ngrid,nsoilmx,nslope,timelen))
    498485allocate(q_co2_PEM_phys(ngrid,timelen))
    499486allocate(q_h2o_PEM_phys(ngrid,timelen))
    500 allocate(co2_ice_GCM(ngrid,nslope,timelen))
     487allocate(co2_ice_PCM(ngrid,nslope,timelen))
    501488allocate(watersurf_density_ave(ngrid,nslope))
    502489allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
    503 allocate(Tsurfave_before_saved(ngrid,nslope))
     490allocate(Tsurfavg_before_saved(ngrid,nslope))
    504491allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    505492allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    … …  
    508495allocate(delta_h2o_icetablesublim(ngrid))
    509496allocate(delta_h2o_adsorbded(ngrid))
    510 allocate(vmr_co2_pem_phys(ngrid,timelen))
     497allocate(vmr_co2_PEM_phys(ngrid,timelen))
    511498
    512499write(*,*) "Downloading data Y1..."
    513 call read_data_PCM("data_PCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
    514                    tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,           &
    515                    co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
     500call read_data_PCM("data_PCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), &
     501                   tsurf_avg_yr1,tsoil_ave,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,                      &
     502                   co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries)
    516503write(*,*) "Downloading data Y1 done"
    517504
    518505! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value
    519506write(*,*) "Downloading data Y2"
    520 call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
    521                    tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,              &
    522                    co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
     507call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,2),min_h2o_ice(:,:,2), &
     508                   tsurf_ave,tsoil_ave,tsurf_PCM_timeseries,tsoil_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,                         &
     509                   co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries)
    523510write(*,*) "Downloading data Y2 done"
    524511
    … …  
    536523if (soil_pem) then
    537524    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
    538     do l = 1,nsoilmx
    539         tsoil_PEM(:,l,:) = tsoil_ave(:,l,:)
    540         tsoil_phys_PEM_timeseries(:,l,:,:) = tsoil_GCM_timeseries(:,l,:,:)
    541         watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,l,:,:)
    542     enddo
     525    tsoil_PEM(:,1:nsoilmx,:) = tsoil_ave
     526    tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_PCM_timeseries
     527    watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
    543528    do l = nsoilmx + 1,nsoilmx_PEM
    544529        tsoil_PEM(:,l,:) = tsoil_ave(:,nsoilmx,:)
    545530        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
    546531    enddo
    547     watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
     532    watersoil_density_PEM_ave = sum(watersoil_density_PEM_timeseries,4)/timelen
    548533endif !soil_pem
    549 deallocate(tsoil_ave,tsoil_GCM_timeseries)
     534deallocate(tsoil_ave,tsoil_PCM_timeseries)
    550535
    551536!------------------------
    … …  
    553538!    I_f Compute tendencies
    554539!------------------------
    555 allocate(tendencies_co2_ice(ngrid,nslope))
    556 allocate(tendencies_co2_ice_ini(ngrid,nslope))
    557 allocate(tendencies_h2o_ice(ngrid,nslope))
     540allocate(tend_co2_ice(ngrid,nslope),tend_h2o_ice(ngrid,nslope))
     541allocate(tend_co2_ice_ini(ngrid,nslope))
    558542
    559543! Compute the tendencies of the evolution of ice over the years
    560 call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice)
    561 tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:)
    562 call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice)
    563 
    564 deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2)
     544call compute_tend(ngrid,nslope,min_co2_ice,tend_co2_ice)
     545call compute_tend(ngrid,nslope,min_h2o_ice,tend_h2o_ice)
     546tend_co2_ice_ini = tend_co2_ice
    565547
    566548!------------------------
    … …  
    568550!    I_g Save initial PCM situation
    569551!------------------------
    570 allocate(initial_co2_ice_sublim(ngrid,nslope))
    571 allocate(initial_co2_ice(ngrid,nslope))
    572 allocate(initial_h2o_ice(ngrid,nslope))
    573 
    574552! We save the places where water ice is sublimating
    575553! We compute the surface of water ice sublimating
    576 ini_surf_co2 = 0.
    577 ini_surf_h2o = 0.
     554co2_surf_ini = 0.
     555h2o_surf_ini = 0.
    578556Total_surface = 0.
    579557do i = 1,ngrid
    580558    Total_surface = Total_surface + cell_area(i)
    581559    do islope = 1,nslope
    582         if (tendencies_co2_ice(i,islope) < 0) then
    583             initial_co2_ice_sublim(i,islope) = 1.
    584             ini_surf_co2 = ini_surf_co2 + cell_area(i)*subslope_dist(i,islope)
    585         else
    586             initial_co2_ice_sublim(i,islope) = 0.
    587         endif
    588         if (qsurf(i,igcm_co2,islope) > 0) then
    589             initial_co2_ice(i,islope) = 1.
    590         else
    591             initial_co2_ice(i,islope) = 0.
    592         endif
    593         if (tendencies_h2o_ice(i,islope) < 0) then
    594             initial_h2o_ice(i,islope) = 1.
    595             ini_surf_h2o = ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
    596         else
    597             initial_h2o_ice(i,islope) = 0.
    598         endif
     560        if (tend_co2_ice(i,islope) < 0.) co2_surf_ini = co2_surf_ini + cell_area(i)*subslope_dist(i,islope)
     561        if (tend_h2o_ice(i,islope) < 0.) h2o_surf_ini = h2o_surf_ini + cell_area(i)*subslope_dist(i,islope)
    599562    enddo
    600563enddo
    601564
    602 write(*,*) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2
    603 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o
     565write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2_surf_ini
     566write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2o_surf_ini
    604567write(*,*) "Total surface of the planet =", Total_surface
    605 allocate(zplev_gcm(ngrid,nlayer + 1))
     568allocate(zplev_PCM(ngrid,nlayer + 1))
    606569
    607570do ig = 1,ngrid
    608     zplev_gcm(ig,:) = ap(:) + bp(:)*ps_start_GCM(ig)
     571    zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
    609572enddo
    610573
    611 global_ave_press_old = 0.
    612 do i = 1,ngrid
    613     global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
    614 enddo
    615 
    616 global_ave_press_GCM = global_ave_press_old
    617 global_ave_press_new = global_ave_press_old
    618 write(*,*) "Initial global average pressure =", global_ave_press_GCM
     574global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
     575global_avg_press_PCM = global_avg_press_old
     576global_avg_press_new = global_avg_press_old
     577write(*,*) "Initial global average pressure =", global_avg_press_PCM
    619578
    620579!------------------------
    … …  
    622581!    I_h Read the startpem.nc
    623582!------------------------
     583allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
     584deallocate(min_co2_ice,min_h2o_ice)
     585
    624586call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
    625               porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,       &
    626               tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),                   &
    627               qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave,          &
    628               co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
    629 
    630 delta_h2o_icetablesublim(:) = 0.
    631 
    632 ! We save the initial values for the co2 frost and perennial ice
    633 allocate(co2frost_ini(ngrid,nslope),perennial_co2ice_ini(ngrid,nslope))
    634 co2frost_ini = qsurf(:,igcm_co2,:)
    635 perennial_co2ice_ini = perennial_co2ice
    636 
    637 do ig = 1,ngrid
    638     do islope = 1,nslope
    639         qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    640         qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope)
    641     enddo
    642 enddo
     587              porefillingice_thickness,tsurf_avg_yr1,tsurf_ave,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,          &
     588              tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM,              &
     589              watersurf_density_ave,watersoil_density_PEM_ave,co2_adsorbded_phys,delta_co2_adsorbded,                &
     590              h2o_adsorbded_phys,delta_h2o_adsorbded)
     591
     592allocate(co2_ice_ini(ngrid,nslope))
     593co2_ice_ini = co2_ice
     594
     595delta_h2o_icetablesublim = 0.
    643596
    644597if (adsorption_pem) then
    … …  
    648601        do islope = 1,nslope
    649602            do l = 1,nsoilmx_PEM - 1
    650                 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     603                totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
    651604                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    652                 totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     605                totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
    653606                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    654607            enddo
    655608        enddo
    656609    enddo
    657 
    658610    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
    659611    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
    660612endif ! adsorption
    661 deallocate(tsurf_ave_yr1)
     613deallocate(tsurf_avg_yr1)
    662614
    663615!------------------------
    … …  
    684636!------------------------
    685637year_iter = 0
     638stopPEM = 0
    686639
    687640do while (year_iter < year_iter_max .and. i_myear < n_myear)
    … …  
    690643    do i = 1,ngrid
    691644        do islope = 1,nslope
    692             global_ave_press_new = global_ave_press_new - g*cell_area(i)*tendencies_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
     645            global_avg_press_new = global_avg_press_new - g*cell_area(i)*tend_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
    693646        enddo
    694647    enddo
    … …  
    696649    if (adsorption_pem) then
    697650        do i = 1,ngrid
    698             global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
     651            global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
    699652        enddo
    700653    endif
    701     write(*,*) 'Global average pressure old time step',global_ave_press_old
    702     write(*,*) 'Global average pressure new time step',global_ave_press_new
     654    write(*,*) 'Global average pressure old time step',global_avg_press_old
     655    write(*,*) 'Global average pressure new time step',global_avg_press_new
    703656
    704657! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
    … …  
    714667    write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..."
    715668    do ig = 1,ngrid
    716         ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
     669        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
    717670    enddo
    718671
    719672! II.a.4. New pressure levels timeseries
    720     allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
     673    allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
    721674    write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
    722675    do l = 1,nlayer + 1
    … …  
    755708            endif
    756709            mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
    757             vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
     710            vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
    758711        enddo
    759712    enddo
    … …  
    763716!------------------------
    764717! II  Run
    765 !    II_b Evolution of the ice
    766 !------------------------
    767     write(*,*) "Evolution of h2o ice"
    768     call evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,STOPPING_1_water)
    769 
    770     write(*,*) "Evolution of co2 ice"
    771     call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
     718!    II_b Evolution of ice
     719!------------------------
     720    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
     721    call evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
    772722
    773723!------------------------
    774724! II  Run
    775 !    II_c CO2 & H2O glaciers flows
    776 !------------------------
    777     write(*,*) "CO2 glacier flows"
    778     if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
    779                                                global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
    780 
    781     write(*,*) "H2O glacier flows"
    782     if (h2oglaciersflow) call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
     725!    II_c Flow of glaciers
     726!------------------------
     727    if (co2_ice_flow) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
     728                                            global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
     729    if (h2o_ice_flow) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
    783730
    784731!------------------------
    … …  
    789736    write(*,*) "Updating the new Tsurf"
    790737    bool_sublim = .false.
    791     Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
     738    Tsurfavg_before_saved = tsurf_ave
    792739    do ig = 1,ngrid
    793740        do islope = 1,nslope
    794             if (initial_co2_ice(ig,islope) > 0.5 .and. qsurf(ig,igcm_co2,islope) < 1.e-10) then !co2ice disappeared, look for closest point without co2ice
     741            if (co2_ice_ini(ig,islope) > 0.5 .and. co2_ice(ig,islope) < 1.e-10) then ! co2 ice disappeared, look for closest point without co2ice
    795742                if (latitude_deg(ig) > 0) then
    796743                    do ig_loop = ig,ngrid
    797744                        do islope_loop = islope,iflat,-1
    798                             if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
     745                            if (co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
    799746                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
    800747                                bool_sublim = .true.
    … …  
    807754                    do ig_loop = ig,1,-1
    808755                        do islope_loop = islope,iflat
    809                             if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
     756                            if(co2_ice_ini(ig_loop,islope_loop) < 0.5 .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
    810757                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
    811758                                bool_sublim = .true.
    … …  
    816763                    enddo
    817764                endif
    818                 initial_co2_ice(ig,islope) = 0
    819                 if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
     765                co2_ice_ini(ig,islope) = 0
     766                if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
    820767                    albedo(ig,1,islope) = albedo_h2o_frost
    821768                    albedo(ig,2,islope) = albedo_h2o_frost
    … …  
    825772                    emis(ig,islope) = emissiv
    826773                endif
    827             else if ((qsurf(ig,igcm_co2,islope) > 1.e-3) .and. (tendencies_co2_ice(ig,islope) > 1.e-10)) then !Put tsurf as tcond co2
     774            else if ((co2_ice(ig,islope) > 1.e-3) .and. (tend_co2_ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
    828775                ave = 0.
    829776                do t = 1,timelen
    830                     if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
    831                         ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
     777                    if (co2_ice_PCM(ig,islope,t) > 1.e-3) then
     778                        ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
    832779                    else
    833                         ave = ave + tsurf_GCM_timeseries(ig,islope,t)
     780                        ave = ave + tsurf_PCM_timeseries(ig,islope,t)
    834781                    endif
    835782                enddo
    … …  
    849796
    850797    do t = 1,timelen
    851         tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:))
     798        tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_ave - Tsurfavg_before_saved
    852799    enddo
    853800    ! for the start
    854801    do ig = 1,ngrid
    855802        do islope = 1,nslope
    856             tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
     803            tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_ave(ig,islope))
    857804        enddo
    858805    enddo
    … …  
    868815        ! Soil averaged
    869816        do islope = 1,nslope
    870             TI_locslope(:,:) = TI_PEM(:,:,islope)
     817            TI_locslope = TI_PEM(:,:,islope)
    871818            do t = 1,timelen
    872                 Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
    873                 Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
     819                Tsoil_locslope = tsoil_phys_PEM_timeseries(:,:,islope,t)
     820                Tsurf_locslope = tsurf_PCM_timeseries(:,islope,t)
    874821                call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    875822                call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    876                 tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
     823                tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope
    877824                do ig = 1,ngrid
    878825                    do isoil = 1,nsoilmx_PEM
    … …  
    883830            enddo
    884831        enddo
    885         tsoil_PEM(:,:,:) = sum(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
    886         watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
     832        tsoil_PEM = sum(tsoil_phys_PEM_timeseries,4)/timelen
     833        watersoil_density_PEM_ave = sum(watersoil_density_PEM_timeseries,4)/timelen
    887834
    888835        write(*,*) "Update of soil temperature done"
    … …  
    892839! II_d.3 Update the ice table
    893840        write(*,*) "Compute ice table"
    894         porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:)
     841        porefillingice_thickness_prev_iter = porefillingice_thickness
    895842        call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
    896843        call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
    897844
    898845! II_d.4 Update the soil thermal properties
    899         write(*,*) "Update soil propreties"
    900         call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
     846        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
    901847
    902848! II_d.5 Update the mass of the regolith adsorbed
    903849        if (adsorption_pem) then
    904             call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
    905                                      qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, &
    906                                      q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
     850            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,      &
     851                                     h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
     852                                     h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
    907853
    908854            totmassco2_adsorbded = 0.
    … …  
    911857                do islope =1, nslope
    912858                    do l = 1,nsoilmx_PEM - 1
    913                         totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     859                        totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
    914860                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    915                         totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     861                        totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
    916862                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    917863                    enddo
    … …  
    927873!    II_e Outputs
    928874!------------------------
    929     call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_ave_press_new/))
     875    call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_avg_press_new/))
    930876    do islope = 1,nslope
    931877        write(str2(1:2),'(i2.2)') islope
    932         call writediagpem(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
    933         call writediagpem(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
    934         call writediagpem(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
    935         call writediagpem(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
     878        call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
     879        call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
     880        call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
     881        call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
    936882        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
    937883        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
    … …  
    942888!    II_f Update the tendencies
    943889!------------------------
    944     write(*,*) "Adaptation of the new co2 tendencies to the current pressure"
    945     call recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),emis,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
    946                                global_ave_press_GCM,global_ave_press_new)
     890    call recomp_tend_co2_slope(ngrid,nslope,timelen,tend_co2_ice,tend_co2_ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries, &
     891                               global_avg_press_PCM,global_avg_press_new)
    947892
    948893!------------------------
    … …  
    950895!    II_g Checking the stopping criterion
    951896!------------------------
    952     call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
    953 
    954     call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
    955                             initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
     897    call stopping_crit_h2o_ice(cell_area,h2o_surf_ini,h2o_ice,stopPEM,ngrid)
     898    call stopping_crit_co2(cell_area,co2_surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
    956899
    957900    year_iter = year_iter + dt_pem
    958901    i_myear = i_myear + dt_pem
    959902
    960     write(*,*) "Checking all the stopping criteria..."
    961     if (STOPPING_water) then
    962         write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
    963         criterion_stop = 1
    964     endif
    965     if (STOPPING_1_water) then
    966         write(*,*) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
    967         criterion_stop = 1
    968     endif
    969     if (STOPPING_co2) then
    970         write(*,*) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
    971         criterion_stop = 2
    972     endif
    973     if (STOPPING_pressure) then
    974         write(*,*) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
    975         criterion_stop = 3
    976     endif
    977     if (year_iter >= year_iter_max) then
    978         write(*,*) "STOPPING because maximum number of iterations reached"
    979         criterion_stop = 4
    980     endif
    981     if (i_myear >= n_myear) then
    982         write(*,*) "STOPPING because maximum number of Martian years to be simulated reached"
    983         criterion_stop = 5
    984     endif
    985 
    986     if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
     903    write(*,*) "Checking the stopping criteria..."
     904    if (year_iter >= year_iter_max) stopPEM = 5
     905    if (i_myear >= n_myear) stopPEM = 6
     906    if (stopPEM > 0) then
     907        select case (stopPEM)
     908            case(1)
     909                write(*,*) "STOPPING because surface of water ice sublimating is too low:", stopPEM, "See message above."
     910            case(2)
     911                write(*,*) "STOPPING because tendencies on water ice = 0:", stopPEM, "See message above."
     912            case(3)
     913                write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
     914            case(4)
     915                write(*,*) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
     916            case(5)
     917                write(*,*) "STOPPING because maximum number of iterations due to orbital parameters is reached:", stopPEM
     918            case(6)
     919                write(*,*) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
     920            case default
     921                write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM
     922        end select
    987923        exit
    988924    else
    … …  
    992928    endif
    993929
    994     global_ave_press_old = global_ave_press_new
    995 
    996 enddo  ! big time iteration loop of the pem
     930    global_avg_press_old = global_avg_press_new
     931
     932enddo ! big time iteration loop of the pem
    997933!------------------------------ END RUN --------------------------------
    998934
    … …  
    1005941
    1006942! H2O ice
     943watercap = 0.
    1007944do ig = 1,ngrid
    1008     if (watercaptag(ig)) then
    1009         watercap_sum = 0.
    1010         do islope = 1,nslope
    1011             ! water_reservoir and watercap have not changed since PCM call: here we check if we have accumulate frost or not.
    1012             if (qsurf(ig,igcm_h2o_ice,islope) > water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) then
    1013                 ! 1st case: more ice than initialy. Put back ancien frost
    1014                 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    1015             else
    1016                 ! 2nd case: ice sublimated. Then let's put qsurf = 0 and add the difference in watercap
    1017                 watercap(ig,islope) = qsurf(ig,igcm_h2o_ice,islope)
    1018                 qsurf(ig,igcm_h2o_ice,islope) = 0.
    1019             endif
    1020             watercap_sum = watercap_sum + watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    1021             watercap(ig,islope) = 0.
    1022         enddo
    1023         water_reservoir(ig) = water_reservoir(ig) + watercap_sum
     945    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
     946    ! If there is enough ice to be considered as an infinite reservoir
     947        watercaptag(ig) = .true.
     948    else
     949    ! If there too little ice to be considered as an infinite reservoir
     950    ! Ice is transferred to the frost
     951        watercaptag(ig) = .false.
     952        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
     953        h2o_ice(ig,:) = 0.
    1024954    endif
    1025955enddo
    1026956
    1027 do ig = 1,ngrid
    1028     water_sum = 0.
    1029     do islope = 1,nslope
    1030         water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    1031     enddo
    1032     if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
    1033         if (water_sum > threshold_h2o_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
    1034             watercaptag(ig) = .true.
    1035             water_reservoir(ig) = water_reservoir(ig) + threshold_h2o_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
    1036             do islope = 1,nslope
    1037                 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_h2o_frost2perennial/2.*cos(pi*def_slope_mean(islope)/180.)
    1038             enddo
    1039         endif
    1040     else ! let's check that the infinite source of water disapear
    1041         if ((water_sum + water_reservoir(ig)) < threshold_h2o_frost2perennial) then
    1042             watercaptag(ig) = .false.
    1043             do islope = 1,nslope
    1044                 qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    1045             enddo
    1046             water_reservoir(ig) = 0.
    1047         endif
    1048     endif
    1049 enddo
    1050 
    1051 ! CO2 ice
    1052 do ig = 1,ngrid
    1053     do islope = 1,nslope
    1054         if (qsurf(ig,igcm_co2,islope) == threshold_co2_frost2perennial) then
    1055         ! If co2 ice is equal to the threshold, then everything is transformed into perennial ice
    1056             perennial_co2ice(ig,islope) = threshold_co2_frost2perennial
    1057             qsurf(ig,igcm_co2,islope) = 0.
    1058         else if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perennial) then
    1059         ! If co2 ice is superior to the threshold, then co2 frost is equal to the threshold (max possible value)
    1060         ! and the leftover is transformed into perennial ice
    1061             perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope) - threshold_co2_frost2perennial
    1062             qsurf(ig,igcm_co2,islope) = threshold_co2_frost2perennial
    1063         else
    1064         ! If co2 ice is inferior to the threshold, then we compare with the initial state
    1065             if (qsurf(ig,igcm_co2,islope) > perennial_co2ice_ini(ig,islope)) then
    1066             ! If co2 ice higher than the initial perennial ice, then the change is affected only to the frost
    1067                 perennial_co2ice(ig,islope) = perennial_co2ice_ini(ig,islope)
    1068                 qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) - perennial_co2ice_ini(ig,islope)
    1069             else
    1070             ! If co2 ice is lower than the initial perennial ice, then there is no frost anymore
    1071                 perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope)
    1072                 qsurf(ig,igcm_co2,islope) = 0.
    1073             endif
    1074         endif
    1075     enddo
    1076 enddo
    1077 
    1078957! III_a.2  Tsoil update (for startfi)
    1079958if (soil_pem) then
    1080     call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
    1081     tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
     959    call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
     960    tsoil = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
    1082961endif
    1083962
    1084963! III_a.4 Pressure (for start)
    1085 ps(:) = ps(:)*global_ave_press_new/global_ave_press_GCM
    1086 ps_start_GCM(:) = ps_start_GCM(:)*global_ave_press_new/global_ave_press_GCM
     964ps = ps*global_avg_press_new/global_avg_press_PCM
     965ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
    1087966
    1088967! III_a.5 Tracer (for start)
    … …  
    1090969
    1091970do l = 1,nlayer + 1
    1092     zplev_new(:,l) = ap(l) + bp(l)*ps_start_GCM(:)
     971    zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
    1093972enddo
    1094973
    … …  
    1097976        do l = 1,llm - 1
    1098977            do ig = 1,ngrid
    1099                 q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
     978                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
    1100979            enddo
    1101980            q(:,llm,nnq) = q(:,llm - 1,nnq)
    … …  
    1104983        do l = 1,llm - 1
    1105984            do ig = 1,ngrid
    1106                 q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
    1107                               + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
     985                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
     986                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_PCM(ig,l) - zplev_PCM(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
    1108987            enddo
    1109988            q(:,llm,nnq) = q(:,llm - 1,nnq)
    … …  
    11781057call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
    11791058             TI_PEM, porefillingice_depth,porefillingice_thickness,      &
    1180              co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
     1059             co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice)
    11811060write(*,*) "restartpem.nc has been written"
    1182 call info_PEM(year_iter,criterion_stop,i_myear,n_myear)
     1061
     1062call info_PEM(year_iter,stopPEM,i_myear,n_myear)
    11831063
    11841064write(*,*) "The PEM has run for", year_iter, "Martian years."
    … …  
    11871067write(*,*) "LL & RV & JBC: so far, so good!"
    11881068
    1189 deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
    1190 deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
     1069deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
     1070deallocate(co2_ice_PCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfavg_before_saved)
    11911071deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
    1192 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
    1193 deallocate(co2frost_ini,perennial_co2ice_ini)
     1072deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
     1073deallocate(co2_ice_ini)
    11941074!----------------------------- END OUTPUT ------------------------------
    11951075

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

    r3143 r3149  
    77!=======================================================================
    88
    9 SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness,  &
    10                     tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
    11                     global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,                &
    12                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
     9SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table,ice_table_thickness, &
     10                    tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,      &
     11                    global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,              &
     12                    m_h2o_regolith_phys,deltam_h2o_regolith_phys)
    1313
    1414use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
    15 use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
    16 use comsoil_h,                  only: volcapa,inertiedat
     15use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_h2o_bigreservoir, depth_breccia, depth_bedrock, index_breccia, index_bedrock
     16use comsoil_h,                  only: volcapa, inertiedat
    1717use adsorption_mod,             only: regolith_adsorption, adsorption_pem
    1818use ice_table_mod,              only: computeice_table_equilibrium
    … …  
    2626#ifndef CPP_STD
    2727    use comcstfi_h,   only: r, mugaz
    28     use surfdat_h,    only: watercaptag
     28    use surfdat_h,    only: watercaptag, perennial_co2ice
    2929#else
    3030    use comcstfi_mod, only: r, mugaz
    … …  
    3737character(len=*),               intent(in) :: filename            ! name of the startfi_PEM.nc
    3838integer,                        intent(in) :: ngrid               ! # of physical grid points
    39 integer,                        intent(in) :: nsoil_GCM           ! # of vertical grid points in the GCM
     39integer,                        intent(in) :: nsoil_PCM           ! # of vertical grid points in the PCM
    4040integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
    4141integer,                        intent(in) :: nslope              ! # of sub-grid slopes
    4242integer,                        intent(in) :: timelen             ! # time samples
    4343real,                           intent(in) :: timestep            ! time step [s]
    44 real,                           intent(in) :: global_ave_pressure ! mean average pressure on the planet [Pa]
    45 real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr1       ! surface temperature at the first year of GCM call [K]
    46 real, dimension(ngrid,nslope),  intent(in) :: tsurf_ave_yr2       ! surface temperature at the second  year of GCM call [K]
     44real,                           intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
     45real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1       ! surface temperature at the first year of PCM call [K]
     46real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second  year of PCM call [K]
    4747real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
    4848real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
    4949real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
    50 real, dimension(ngrid,nslope),  intent(in) :: tend_h2oglaciers    ! tendencies on h2o glaciers
    51 real, dimension(ngrid,nslope),  intent(in) :: tend_co2glaciers    ! tendencies on co2 glaciers
    52 real, dimension(ngrid,nslope),  intent(in) :: co2ice              ! co2 ice amount [kg/m^2]
    53 real, dimension(ngrid,nslope),  intent(in) :: waterice            ! water ice amount [kg/m^2]
    54 real, dimension(ngrid,nslope),  intent(in) :: watersurf_ave       ! surface water ice density, yearly averaged  (kg/m^3)
     50real, dimension(ngrid,nslope),  intent(in) :: tend_h2o_ice        ! tendencies on h2o ice
     51real, dimension(ngrid,nslope),  intent(in) :: tend_co2_ice        ! tendencies on co2 ice
     52real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg       ! surface water ice density, yearly averaged (kg/m^3)
    5553! outputs
    5654real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
    5755real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
     56real, dimension(ngrid,nslope),                   intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
     57real, dimension(ngrid,nslope),                   intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
    5858real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
    5959real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
    … …  
    6363real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
    6464real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
    65 real, dimension(ngrid),                          intent(inout) :: water_reservoir     ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
    66 real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_ave       ! surface water ice density, yearly averaged (kg/m^3)
     65real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
    6766! local
    6867real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM               ! soil temperature saved in the start [K]
    … …  
    7574character(2)                            :: num                          ! intermediate string to read PEM start sloped variables
    7675real, dimension(ngrid,nsoil_PEM)        :: tsoil_saved                  ! saved soil temperature [K]
    77 real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr1[K]
     76real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                ! intermediate soil temperature during yr 1 [K]
    7877real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                ! intermediate soil temperature during yr 2 [K]
    7978real, dimension(ngrid,nsoil_PEM - 1)    :: alph_tmp                     ! Intermediate for tsoil computation []
    … …  
    8281
    8382#ifdef CPP_STD
    84    logical, dimension(ngrid) :: watercaptag
    85    watercaptag(:) = .false.
     83    logical, dimension(ngrid) :: watercaptag
     84    watercaptag(:) = .false.
    8685#endif
    8786
    … …  
    9190!!!
    9291!!! Order: 0. Previous year of the PEM run
     92!!!           Ice initialization
    9393!!!        1. Thermal Inertia
    9494!!!        2. Soil Temperature
    9595!!!        3. Ice table
    9696!!!        4. Mass of CO2 & H2O adsorbed
    97 !!!        5. Water reservoir
    9897!!!
    9998!!! /!\ This order must be respected !
    … …  
    113112    call open_startphy(filename)
    114113
     114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     115#ifndef CPP_STD
     116    ! h2o ice
     117    h2o_ice = 0.
     118    call get_field("h2o_ice",h2o_ice,found)
     119    if (.not. found) then
     120        write(*,*)'Pemetat0: failed loading <h2o_ice>'
     121        write(*,*)'will reconstruct the values from watercaptag'
     122        write(*,*)'with default value ''ini_h2o_bigreservoir'''
     123        do ig = 1,ngrid
     124            if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
     125        enddo
     126    endif
     127
     128    ! co2 ice
     129    co2_ice = perennial_co2ice
     130#endif
     131
    115132    if (soil_pem) then
    116133
     134!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    117135!1. Thermal Inertia
    118136! a. General case
    119     do islope = 1,nslope
    120         write(num,'(i2.2)') islope
    121         call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
     137        do islope = 1,nslope
     138            write(num,'(i2.2)') islope
     139            call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
     140            if (.not. found) then
     141                write(*,*)'PEM settings: failed loading < TI_PEM_slope'//num//'>'
     142                write(*,*)'will reconstruct the values of TI_PEM'
     143
     144                do ig = 1,ngrid
     145                    if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
     146                        !!! transition
     147                        delta = depth_breccia
     148                        TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
     149                                                                 (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
     150                                                                 ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
     151                        do iloop=index_breccia+2,index_bedrock
     152                            TI_PEM(ig,iloop,islope) = TI_breccia
     153                        enddo
     154                    else ! we keep the high ti values
     155                        do iloop = index_breccia + 1,index_bedrock
     156                            TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
     157                        enddo
     158                    endif ! TI PEM and breccia comparison
     159                    !! transition
     160                    delta = depth_bedrock
     161                    TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
     162                                                               (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
     163                                                               ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
     164                    do iloop = index_bedrock + 2,nsoil_PEM
     165                        TI_PEM(ig,iloop,islope) = TI_bedrock
     166                    enddo
     167                enddo
     168            else ! found
     169                do iloop = nsoil_PCM + 1,nsoil_PEM
     170                    TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
     171                enddo
     172            endif ! not found
     173        enddo ! islope
     174
     175        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
     176
     177! b. Special case for inertiedat, inertiedat_PEM
     178        call get_field("inertiedat_PEM",inertiedat_PEM,found)
     179        if (.not.found) then
     180            do iloop = 1,nsoil_PCM
     181                inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
     182            enddo
     183        !!! zone de transition
     184            delta = depth_breccia
     185            do ig = 1,ngrid
     186                if (inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
     187                    inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
     188                                                              (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
     189                                                              ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
     190
     191                    do iloop = index_breccia+2,index_bedrock
     192                        inertiedat_PEM(ig,iloop) = TI_breccia
     193                    enddo
     194                else
     195                    do iloop=index_breccia+1,index_bedrock
     196                        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_PCM)
     197                    enddo
     198                endif ! comparison ti breccia
     199            enddo ! ig
     200
     201            !!! zone de transition
     202            delta = depth_bedrock
     203            do ig = 1,ngrid
     204                inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
     205                                                          (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
     206                                                          ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
     207            enddo
     208
     209            do iloop = index_bedrock + 2, nsoil_PEM
     210                do ig = 1,ngrid
     211                    inertiedat_PEM(ig,iloop) = TI_bedrock
     212                enddo
     213            enddo
     214        endif ! not found
     215
     216!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     217!2. Soil Temperature
     218        do islope=1,nslope
     219            write(num,fmt='(i2.2)') islope
     220            call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
     221            if (.not.found) then
     222                write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
     223                write(*,*)'will reconstruct the values of Tsoil'
     224!                do ig = 1,ngrid
     225!                    kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
     226!                    tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
     227!                    do iloop=index_breccia+2,index_bedrock
     228!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
     229!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
     230!                    enddo
     231!                    kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
     232!                    tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
     233!
     234!                    do iloop=index_bedrock+2,nsoil_PEM
     235!                        kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
     236!                        tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
     237!                    enddo
     238!                enddo
     239                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
     240                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
     241            else
     242! predictor corrector: restart from year 1 of the PCM and build the evolution of
     243! tsoil at depth
     244                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
     245                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
     246                call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
     247                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
     248                call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
     249
     250                do iloop = nsoil_PCM+1,nsoil_PEM
     251                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
     252                enddo
     253            endif !found
     254
     255            do it = 1,timelen
     256                do isoil = nsoil_PCM+1,nsoil_PEM
     257                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
     258                enddo
     259            enddo
     260            do isoil = nsoil_PCM+1,nsoil_PEM
     261                do ig = 1,ngrid
     262                    watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
     263                enddo
     264            enddo
     265        enddo ! islope
     266        write(*,*) 'PEMETAT0: SOIL TEMP done'
     267
     268!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     269!3. Ice Table
     270        call get_field("ice_table",ice_table,found)
    122271        if (.not. found) then
    123             write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
    124             write(*,*)'will reconstruct the values of TI_PEM'
    125 
     272            write(*,*)'PEM settings: failed loading <ice_table>'
     273            write(*,*)'will reconstruct the values of the ice table given the current state'
     274            ice_table = -1  ! by default, no ice table
     275            ice_table_thickness = -1  ! by default, no ice table
     276            call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg, TI_PEM(:,1,:),ice_table,ice_table_thickness)
     277            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
     278            do islope = 1,nslope
     279                call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
     280            enddo
     281        endif
     282
     283        write(*,*) 'PEMETAT0: ICE TABLE done'
     284
     285!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     286!4. CO2 & H2O Adsorption
     287        if (adsorption_pem) then
     288            do islope = 1,nslope
     289                write(num,fmt='(i2.2)') islope
     290                call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
     291                if (.not. found) then
     292                    m_co2_regolith_phys = 0.
     293                    exit
     294                endif
     295            enddo
     296            do islope=1,nslope
     297                write(num,fmt='(i2.2)') islope
     298                call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
     299                if (.not. found2) then
     300                    m_h2o_regolith_phys = 0.
     301                    exit
     302                endif
     303            enddo
     304
     305            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
     306                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
     307
     308            if (.not. found) deltam_co2_regolith_phys = 0.
     309            if (.not.found2) deltam_h2o_regolith_phys = 0.
     310
     311            write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
     312        endif ! adsorption_pem
     313    endif ! soil_pem
     314
     315    call close_startphy
     316
     317else !No startfi, let's build all by hand
     318
     319    ! h2o ice
     320    h2o_ice = 0.
     321    write(*,*)'There is no "startpem.nc" so ''h2o_ice'' is reconstructed from watercaptag with default value ''ini_h2o_bigreservoir''.'
     322        do ig = 1,ngrid
     323            if (watercaptag(ig)) h2o_ice(ig,:) = ini_h2o_bigreservoir
     324        enddo
     325
     326    ! co2 ice
     327    co2_ice = perennial_co2ice
     328
     329    if (soil_pem) then
     330
     331!a) Thermal inertia
     332        do islope = 1,nslope
    126333            do ig = 1,ngrid
    127334                if (TI_PEM(ig,index_breccia,islope) < TI_breccia) then
    128335                    !!! transition
    129336                    delta = depth_breccia
    130                     TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/ &
    131                                                         (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
    132                                                         ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
    133                     do iloop=index_breccia+2,index_bedrock
     337                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
     338                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
     339                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
     340                    do iloop=index_breccia + 2,index_bedrock
    134341                        TI_PEM(ig,iloop,islope) = TI_breccia
    135342                    enddo
    136343                else ! we keep the high ti values
    137                     do iloop = index_breccia + 1,index_bedrock
     344                    do iloop=index_breccia + 1,index_bedrock
    138345                        TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
    139346                    enddo
    140                 endif ! TI PEM and breccia comparison
     347                endif
    141348                !! transition
    142349                delta = depth_bedrock
    143                 TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/ &
    144                                                         (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
    145                                                         ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
     350                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
     351                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
     352                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
    146353                do iloop = index_bedrock + 2,nsoil_PEM
    147354                    TI_PEM(ig,iloop,islope) = TI_bedrock
    148355                enddo
    149356            enddo
    150         else ! found
    151             do iloop = nsoil_GCM + 1,nsoil_PEM
    152                 TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
    153             enddo
    154         endif ! not found
    155     enddo ! islope
    156 
    157     write(*,*) 'PEMETAT0: THERMAL INERTIA done'
    158 
    159 ! b. Special case for inertiedat, inertiedat_PEM
    160 call get_field("inertiedat_PEM",inertiedat_PEM,found)
    161 if(.not.found) then
    162  do iloop = 1,nsoil_GCM
    163    inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
    164  enddo
    165 !!! zone de transition
    166 delta = depth_breccia
    167 do ig = 1,ngrid
    168 if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
    169 inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    170             (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
    171                       ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
    172 
    173  do iloop = index_breccia+2,index_bedrock
    174        inertiedat_PEM(ig,iloop) = TI_breccia
    175   enddo
    176 
    177 else
    178    do iloop=index_breccia+1,index_bedrock
    179       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
    180    enddo
    181 endif ! comparison ti breccia
    182 enddo!ig
    183 
    184 !!! zone de transition
    185 delta = depth_bedrock
    186 do ig = 1,ngrid
    187 inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
    188             (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
    189                       ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
    190 enddo
    191 
    192  do iloop = index_bedrock+2, nsoil_PEM
    193    do ig = 1,ngrid
    194       inertiedat_PEM(ig,iloop) = TI_bedrock
    195    enddo
    196  enddo
    197 endif ! not found
    198 
    199 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    200 !2. Soil Temperature
    201 do islope=1,nslope
    202   write(num,fmt='(i2.2)') islope
    203    call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
    204   if(.not.found) then
    205       write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
    206       write(*,*)'will reconstruct the values of Tsoil'
    207 !      do ig = 1,ngrid
    208 !        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
    209 !        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
    210 !       do iloop=index_breccia+2,index_bedrock
    211 !            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
    212 !            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
    213 !      enddo
    214 !        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
    215 !        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
     357        enddo
     358
     359        do iloop = 1,nsoil_PCM
     360            inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
     361        enddo
     362        !!! zone de transition
     363        delta = depth_breccia
     364        do ig = 1,ngrid
     365            if (inertiedat_PEM(ig,index_breccia) < TI_breccia) then
     366                inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
     367                                                            (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
     368                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
     369                do iloop = index_breccia + 2,index_bedrock
     370                    inertiedat_PEM(ig,iloop) = TI_breccia
     371                enddo
     372            else
     373                do iloop = index_breccia + 1,index_bedrock
     374                    inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
     375                enddo
     376            endif
     377        enddo
     378
     379        !!! zone de transition
     380        delta = depth_bedrock
     381        do ig = 1,ngrid
     382            inertiedat_PEM(ig,index_bedrock + 1) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                    &
     383                                                        (((delta - layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2)) + &
     384                                                        ((layer_PEM(index_bedrock + 1) - delta)/(TI_bedrock**2))))
     385        enddo
     386
     387        do iloop = index_bedrock + 2,nsoil_PEM
     388            do ig = 1,ngrid
     389                inertiedat_PEM(ig,iloop) = TI_bedrock
     390            enddo
     391        enddo
     392
     393        write(*,*) 'PEMETAT0: TI done'
     394
     395!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     396!b) Soil temperature
     397        do islope = 1,nslope
     398!            do ig = 1,ngrid
     399!                kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
     400!                tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
    216401!
    217 !      do iloop=index_bedrock+2,nsoil_PEM
    218 !            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
    219 !            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
    220 !      enddo
    221 !      enddo
    222 
    223 
    224      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    225      call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    226 
    227 
    228    else
    229 ! predictor corrector: restart from year 1 of the GCM and build the evolution of
    230 ! tsoil at depth
    231 
    232     tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
    233     call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
    234     call soil_pem_compute(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
    235     tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
    236     call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
    237 
    238 
    239      do iloop = nsoil_GCM+1,nsoil_PEM
    240        tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
    241      enddo
    242    endif !found
    243 
    244     do it = 1,timelen
    245         do isoil = nsoil_GCM+1,nsoil_PEM
    246         tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
    247         enddo
    248      enddo
    249       do isoil = nsoil_GCM+1,nsoil_PEM
    250         do ig = 1,ngrid
    251         watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    252         enddo
    253       enddo
    254 enddo ! islope
    255 
    256 write(*,*) 'PEMETAT0: SOIL TEMP done'
    257 
    258 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    259 !3. Ice Table
    260   call get_field("ice_table",ice_table,found)
    261    if(.not.found) then
    262       write(*,*)'PEM settings: failed loading <ice_table>'
    263       write(*,*)'will reconstruct the values of the ice table given the current state'
    264       ice_table(:,:) = -1  ! by default, no ice table 
    265       ice_table_thickness(:,:) = -1  ! by default, no ice table
    266      call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, TI_PEM(:,1,:),ice_table,ice_table_thickness)
    267      call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
    268      do islope = 1,nslope
    269      call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    270      enddo
    271    endif
    272 
    273 write(*,*) 'PEMETAT0: ICE TABLE done'
    274 
    275 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    276 !4. CO2 & H2O Adsorption
    277  if(adsorption_pem) then
    278   do islope=1,nslope
    279    write(num,fmt='(i2.2)') islope
    280    call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
    281     if((.not.found)) then
    282        m_co2_regolith_phys(:,:,:) = 0.
    283        exit
    284     endif
    285 
    286   enddo
    287 
    288   do islope=1,nslope
    289    write(num,fmt='(i2.2)') islope
    290    call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
    291     if((.not.found2)) then
    292        m_h2o_regolith_phys(:,:,:) = 0.
    293       exit
    294     endif
    295 
    296   enddo
    297 
    298     call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
    299                                 m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
    300 
    301     if((.not.found)) then
    302       deltam_co2_regolith_phys(:) = 0.
    303     endif
    304     if((.not.found2)) then
    305        deltam_h2o_regolith_phys(:) = 0.
    306     endif
    307 write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
    308     endif
    309 endif ! soil_pem
    310 
    311 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    312 !. 5 water reservoir
    313 #ifndef CPP_STD
    314    water_reservoir = 0.
    315    call get_field("water_reservoir",water_reservoir,found)
    316    if (.not. found) then
    317        write(*,*)'Pemetat0: failed loading <water_reservoir>'
    318        write(*,*)'will reconstruct the values from watercaptag'
    319        where (watercaptag) water_reservoir = water_reservoir_nom
    320    endif
    321 #endif
    322 
    323 call close_startphy
    324 
    325 else !No startfi, let's build all by hand
    326 
    327     if (soil_pem) then
    328 
    329 !a) Thermal inertia
    330    do islope = 1,nslope
    331       do ig = 1,ngrid
    332 
    333           if(TI_PEM(ig,index_breccia,islope).lt.TI_breccia) then
    334 !!! transition
    335              delta = depth_breccia
    336              TI_PEM(ig,index_breccia+1,islope) =sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    337             (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
    338                       ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
    339 
    340           do iloop=index_breccia+2,index_bedrock
    341             TI_PEM(ig,iloop,islope) = TI_breccia
    342          enddo
    343          else ! we keep the high ti values
    344            do iloop=index_breccia+1,index_bedrock
    345                   TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
    346            enddo
    347          endif
    348 
    349 !! transition
    350              delta = depth_bedrock
    351              TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
    352             (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
    353                       ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
    354           do iloop=index_bedrock+2,nsoil_PEM
    355             TI_PEM(ig,iloop,islope) = TI_bedrock
    356          enddo
    357       enddo
    358 enddo
    359 
    360  do iloop = 1,nsoil_GCM
    361    inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
    362  enddo
    363 !!! zone de transition
    364 delta = depth_breccia
    365 do ig = 1,ngrid
    366       if(inertiedat_PEM(ig,index_breccia).lt.TI_breccia) then
    367 inertiedat_PEM(ig,index_breccia+1) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
    368             (((delta-layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2))+ &
    369                       ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
    370 
    371 
    372  do iloop = index_breccia+2,index_bedrock
    373 
    374        inertiedat_PEM(ig,iloop) = TI_breccia
    375 
    376  enddo
    377 else
    378    do iloop = index_breccia+1,index_bedrock
    379        inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,index_breccia)
    380     enddo
    381 
    382 endif
    383 enddo
    384 
    385 !!! zone de transition
    386 delta = depth_bedrock
    387 do ig = 1,ngrid
    388 inertiedat_PEM(ig,index_bedrock+1) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
    389             (((delta-layer_PEM(index_bedrock))/(inertiedat_PEM(ig,index_bedrock)**2))+ &
    390                       ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
    391 enddo
    392 
    393 
    394 
    395  do iloop = index_bedrock+2, nsoil_PEM
    396    do ig = 1,ngrid
    397       inertiedat_PEM(ig,iloop) = TI_bedrock
    398    enddo
    399  enddo
    400 
    401 write(*,*) 'PEMETAT0: TI done'
    402 
    403 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    404 !b) Soil temperature
    405     do islope = 1,nslope
    406 !     do ig = 1,ngrid
    407 !        kcond = (TI_PEM(ig,index_breccia+1,islope)*TI_PEM(ig,index_breccia+1,islope))/volcapa
    408 !        tsoil_PEM(ig,index_breccia+1,islope) = tsoil_PEM(ig,index_breccia,islope) + fluxgeo/kcond*(mlayer_PEM(index_breccia)-mlayer_PEM(index_breccia-1))
     402!                do iloop=index_breccia+2,index_bedrock
     403!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
     404!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
     405!                enddo
     406!                kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
     407!                tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
    409408!
    410 !       do iloop=index_breccia+2,index_bedrock
    411 !            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
    412 !            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_breccia+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_breccia))
    413 !      enddo
    414 !        kcond = (TI_PEM(ig,index_bedrock+1,islope)*TI_PEM(ig,index_bedrock+1,islope))/volcapa
    415 !        tsoil_PEM(ig,index_bedrock+1,islope) = tsoil_PEM(ig,index_bedrock,islope) + fluxgeo/kcond*(mlayer_PEM(index_bedrock)-mlayer_PEM(index_bedrock-1))
    416 
    417 !       do iloop=index_bedrock+2,nsoil_PEM
    418 !            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
    419 !            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
    420 !      enddo
    421 
    422 !       enddo
    423 
    424       call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    425       call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    426  
    427 ! First raw initialization
    428       do it = 1,timelen
    429         do isoil = nsoil_GCM+1,nsoil_PEM
    430           tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
    431         enddo
    432       enddo
    433 
    434       do it = 1,timelen
    435         do isoil = nsoil_GCM+1,nsoil_PEM
    436           call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_inst(:,:,islope,it))
    437         enddo
    438        enddo
    439 
    440        do isoil = nsoil_GCM+1,nsoil_PEM
    441          do ig = 1,ngrid
    442            watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
    443          enddo
    444        enddo
    445 enddo !islope
    446 write(*,*) 'PEMETAT0: TSOIL done'
     409!                do iloop=index_bedrock+2,nsoil_PEM
     410!                    kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
     411!                    tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,index_bedrock+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(index_bedrock))
     412!                enddo
     413!            enddo
     414            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
     415            call soil_pem_compute(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
     416
     417! First raw initialization
     418            do it = 1,timelen
     419                do isoil = nsoil_PCM + 1,nsoil_PEM
     420                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
     421                enddo
     422            enddo
     423
     424            do it = 1,timelen
     425                do isoil = nsoil_PCM + 1,nsoil_PEM
     426                    call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
     427                enddo
     428            enddo
     429
     430            do isoil = nsoil_PCM + 1,nsoil_PEM
     431                do ig = 1,ngrid
     432                    watersoil_avg(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
     433                enddo
     434            enddo
     435        enddo !islope
     436        write(*,*) 'PEMETAT0: TSOIL done'
    447437
    448438!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    449439!c) Ice table
    450        ice_table(:,:) = -1.  ! by default, no ice table
    451        ice_table_thickness(:,:) = -1.
    452        call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
    453        call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
    454        do islope = 1,nslope
    455          call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    456        enddo
    457        write(*,*) 'PEMETAT0: Ice table done'
     440        ice_table = -1.  ! by default, no ice table
     441        ice_table_thickness = -1.
     442        call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),ice_table,ice_table_thickness)
     443        call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2o_ice,h2o_ice,global_avg_pressure,ice_table,ice_table_thickness,TI_PEM)
     444        do islope = 1,nslope
     445            call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
     446        enddo
     447        write(*,*) 'PEMETAT0: Ice table done'
    458448
    459449!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    460450!d) Regolith adsorbed
    461 if (adsorption_pem) then
    462     m_co2_regolith_phys(:,:,:) = 0.
    463     m_h2o_regolith_phys(:,:,:) = 0.
    464 
    465     call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2oglaciers,tend_co2glaciers,waterice,co2ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
    466                              m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
    467 
    468     deltam_co2_regolith_phys(:) = 0.
    469     deltam_h2o_regolith_phys(:) = 0.
    470 endif
    471 
    472 write(*,*) 'PEMETAT0: CO2 adsorption done'
    473 endif !soil_pem
    474 
    475 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    476 !. e) water reservoir
    477 #ifndef CPP_STD
    478     water_reservoir = 0.
    479     where (watercaptag) water_reservoir = water_reservoir_nom
    480 #endif
    481 
     451        if (adsorption_pem) then
     452            m_co2_regolith_phys = 0.
     453            m_h2o_regolith_phys = 0.
     454            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,tend_h2o_ice,tend_co2_ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
     455                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
     456            deltam_co2_regolith_phys = 0.
     457            deltam_h2o_regolith_phys = 0.
     458        endif
     459
     460        write(*,*) 'PEMETAT0: CO2 adsorption done'
     461    endif !soil_pem
    482462endif ! of if (startphy_file)
    483463
    … …  
    495475
    496476END MODULE pemetat0_mod
    497 

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

    r3082 r3149  
    1 module pemredem
     1MODULE pemredem
     2
    23!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    34!!!
    4 !!! Purpose: Write specific netcdf restart for the PEM 
    5 !!! 
    6 !!! 
    7 !!! Author: LL, inspired by phyredem  from the GCM
     5!!! Purpose: Write specific netcdf restart for the PEM
     6!!!
     7!!!
     8!!! Author: LL, inspired by phyredem from the PCM
    89!!!
    910!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    1112implicit none
    1213
     14!=======================================================================
    1315contains
     16!=======================================================================
    1417
    15 subroutine pemdem0(filename,lonfi,latfi,cell_area,nsoil_PEM,ngrid,day_ini,time,nslope,def_slope,subslope_dist)
     18SUBROUTINE pemdem0(filename,lonfi,latfi,cell_area,nsoil_PEM,ngrid,day_ini,time,nslope,def_slope,subslope_dist)
     19
    1620! create physics restart file and write time-independent variables
    17   use comsoil_h_PEM,     only: soil_pem,mlayer_PEM,fluxgeo,inertiedat_PEM
    18   use iostart_PEM,       only: open_restartphy, close_restartphy, put_var, put_field, length
    19   use mod_grid_phy_lmdz, only: klon_glo
    20   use time_phylmdz_mod,  only: daysec
     21use comsoil_h_PEM,     only: soil_pem,mlayer_PEM,fluxgeo,inertiedat_PEM
     22use iostart_PEM,       only: open_restartphy, close_restartphy, put_var, put_field, length
     23use mod_grid_phy_lmdz, only: klon_glo
     24use time_phylmdz_mod,  only: daysec
     25
    2126#ifndef CPP_STD
    22   use planete_h,    only: aphelie, emin_turb, lmixmin, obliquit, peri_day, periheli, year_day
    23   use comcstfi_h,   only: g, mugaz, omeg, rad, rcp
     27    use planete_h,    only: aphelie, emin_turb, lmixmin, obliquit, peri_day, periheli, year_day
     28    use comcstfi_h,   only: g, mugaz, omeg, rad, rcp
    2429#else
    25   use planete_mod,  only: apoastr, emin_turb, lmixmin, obliquit, peri_day, periastr, year_day
    26   use comcstfi_mod, only: g, mugaz, omeg, rad, rcp
     30    use planete_mod,  only: apoastr, emin_turb, lmixmin, obliquit, peri_day, periastr, year_day
     31    use comcstfi_mod, only: g, mugaz, omeg, rad, rcp
    2732#endif
    2833
    29   implicit none
    30  
    31   character(len=*), intent(in) :: filename
    32   real, intent(in)             :: lonfi(ngrid)
    33   real, intent(in)             :: latfi(ngrid)
    34   integer, intent(in)          :: nsoil_PEM
    35   integer, intent(in)          :: ngrid
    36   real, intent(in)             :: day_ini
    37   real, intent(in)             :: time
    38   real, intent(in)             :: cell_area(ngrid) !boundaries for bining of the slopes
    39   integer, intent(in)          :: nslope
    40   real, intent(in)             :: def_slope(nslope+1) !boundaries for bining of the slopes
    41   real, intent(in)             :: subslope_dist(ngrid,nslope) !undermesh statistics
    42  
    43   ! Create physics start file
    44   call open_restartphy(filename)
    45  
    46   ! Write the mid-layer depths
    47   call put_var("soildepth","Soil mid-layer depth",mlayer_PEM)
    48  
    49   ! Write longitudes
    50   call put_field("longitude","Longitudes of physics grid",lonfi)
    51  
    52   ! Write latitudes
    53   call put_field("latitude","Latitudes of physics grid",latfi)
    54  
    55   ! Write mesh areas
    56   call put_field("area","Mesh area",cell_area)
    57      
    58   ! Multidimensionnal variables (nogcm undermesh slope statistics)
    59   call put_var("def_slope","slope criterium stages",def_slope)
    60   call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
     34implicit none
    6135
     36character(*),                  intent(in) :: filename
     37integer,                       intent(in) :: nsoil_PEM, ngrid, nslope
     38real, dimension(ngrid),        intent(in) :: lonfi, latfi
     39real,                          intent(in) :: day_ini, time
     40real, dimension(ngrid),        intent(in) :: cell_area     ! boundaries for bining of the slopes
     41real, dimension(ngrid + 1),    intent(in) :: def_slope     ! boundaries for bining of the slopes
     42real, dimension(ngrid,nslope), intent(in) :: subslope_dist ! undermesh statistics
    6243
    63   ! Close file
    64   call close_restartphy
    65  
    66 end subroutine pemdem0
     44! Create physics start file
     45call open_restartphy(filename)
    6746
    68 subroutine pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope, &
    69                     tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table_depth,ice_table_thickness, &
    70                     m_co2_regolith,m_h2o_regolith,water_reservoir)
     47! Write the mid-layer depths
     48call put_var("soildepth","Soil mid-layer depth",mlayer_PEM)
    7149
     50! Write longitudes
     51call put_field("longitude","Longitudes of physics grid",lonfi)
    7252
    73   ! write time-dependent variable to restart file
    74   use iostart_PEM,   only: open_restartphy, close_restartphy, put_var, put_field
    75   use comsoil_h_PEM, only: mlayer_PEM,fluxgeo, inertiedat_PEM,soil_pem
    76   use time_evol_mod, only: year_bp_ini, convert_years
    77                
    78   implicit none
    79  
     53! Write latitudes
     54call put_field("latitude","Latitudes of physics grid",latfi)
     55
     56! Write mesh areas
     57call put_field("area","Mesh area",cell_area)
     58
     59! Multidimensionnal variables (nopcm undermesh slope statistics)
     60call put_var("def_slope","slope criterium stages",def_slope)
     61call put_field("subslope_dist","under mesh slope distribution",subslope_dist)
     62
     63! Close file
     64call close_restartphy
     65
     66END SUBROUTINE pemdem0
     67
     68!=======================================================================
     69
     70SUBROUTINE pemdem1(filename,i_myear,nsoil_PEM,ngrid,nslope,tsoil_slope_PEM,inertiesoil_slope_PEM, &
     71                   ice_table_depth,ice_table_thickness,m_co2_regolith,m_h2o_regolith,h2o_ice)
     72
     73! write time-dependent variable to restart file
     74use iostart_PEM,   only: open_restartphy, close_restartphy, put_var, put_field
     75use comsoil_h_PEM, only: mlayer_PEM,fluxgeo, inertiedat_PEM, soil_pem
     76use time_evol_mod, only: year_bp_ini, convert_years
     77
     78implicit none
     79
    8080#ifndef CPP_STD
    81   include "callkeys.h"
     81    include "callkeys.h"
    8282#endif
    83  
    84   character(len=*), intent(in) :: filename
    85   integer, intent(in)          :: nsoil_PEM, ngrid, nslope, i_myear
    86   real, intent(in)             :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
    87   real, intent(in)             :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
    88   real, intent(in)             :: ice_table_depth(ngrid,nslope) !under mesh bining according to slope
    89   real, intent(in)             :: ice_table_thickness(ngrid,nslope) !under mesh bining according to slope
    90   real, intent(in)             :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
    91   real, intent(in)             :: m_h2o_regolith(ngrid,nsoil_PEM,nslope)
    92   real, intent(in)             :: water_reservoir(ngrid)
    9383
    94   integer           :: islope
    95   character*2       :: num
    96   integer           :: iq
    97   character(len=30) :: txt  ! To store some text
    98   real              :: Year ! Year of the simulation
     84character(*),                            intent(in) :: filename
     85integer,                                 intent(in) :: nsoil_PEM, ngrid, nslope, i_myear
     86real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_slope_PEM       ! under mesh bining according to slope
     87real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: inertiesoil_slope_PEM ! under mesh bining according to slope
     88real, dimension(ngrid,nslope),           intent(in) :: ice_table_depth       ! under mesh bining according to slope
     89real, dimension(ngrid,nslope),           intent(in) :: ice_table_thickness   ! under mesh bining according to slope
     90real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith
     91real, dimension(ngrid,nslope),           intent(in) :: h2o_ice
    9992
    100   ! Open file
    101   call open_restartphy(filename)
    102  
    103   ! First variable to write must be "Time", in order to correctly
    104   ! set time counter in file
    105   Year = (year_bp_ini + i_myear)*convert_years
    106   call put_var("Time","Year of simulation",Year)
    107   call put_field('water_reservoir','water_reservoir',water_reservoir,Year)
     93integer       :: islope
     94character(2)  :: num
     95integer       :: iq
     96character(30) :: txt  ! To store some text
     97real          :: Year ! Year of the simulation
    10898
    109     if(soil_pem) then
    110  
     99! Open file
     100call open_restartphy(filename)
     101
     102! First variable to write must be "Time", in order to correctly
     103! set time counter in file
     104Year = (year_bp_ini + i_myear)*convert_years
     105call put_var("Time","Year of simulation",Year)
     106call put_field('h2o_ice','h2o_ice',h2o_ice,Year)
     107
     108if (soil_pem) then
    111109  ! Multidimensionnal variables (undermesh slope statistics)
    112 do islope = 1,nslope
    113   write(num,fmt='(i2.2)') islope
    114   call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type", &
    115                  tsoil_slope_PEM(:,:,islope),Year)
    116   call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type", &
    117                  inertiesoil_slope_PEM(:,:,islope),Year)
     110    do islope = 1,nslope
     111        write(num,fmt='(i2.2)') islope
     112        call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type",tsoil_slope_PEM(:,:,islope),Year)
     113        call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type",inertiesoil_slope_PEM(:,:,islope),Year)
     114        call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith",m_co2_regolith(:,:,islope),Year)
     115        call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith",m_h2o_regolith(:,:,islope),Year)
     116    enddo
     117    call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
     118    call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
     119    call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
     120endif ! soil_pem
    118121
    119   call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbded in the regolith", &
    120                 m_co2_regolith(:,:,islope),Year)
    121   call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbded in the regolith", &
    122                  m_h2o_regolith(:,:,islope),Year)
    123 enddo
     122! Close file
     123call close_restartphy
    124124
    125   call put_field("ice_table_depth","Depth of ice table",ice_table_depth,Year)
     125END SUBROUTINE pemdem1
    126126
    127   call put_field("ice_table_thickness","Depth of ice table",ice_table_thickness,Year)
    128 
    129   call put_field("inertiedat_PEM","Thermal inertie of PEM ",inertiedat_PEM,Year)
    130 
    131    endif !soil_pem
    132 
    133   ! Close file
    134   call close_restartphy
    135  
    136 end subroutine pemdem1
    137 
    138 end module pemredem
     127END MODULE pemredem

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

    r3143 r3149  
    1515!=======================================================================
    1616
    17 SUBROUTINE read_data_PCM(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,          &
    18                          min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
    19                          watersurf_density_ave,watersoil_density)
     17SUBROUTINE read_data_PCM(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_PCM_phys,ps_timeseries,          &
     18                         min_co2_ice,min_h2o_ice,tsurf_avg,tsoil_avg,tsurf_PCM,tsoil_PCM,q_co2,q_h2o,co2_ice_slope, &
     19                         watersurf_density_avg,watersoil_density)
    2020use comsoil_h,             only: nsoilmx
    2121use comsoil_h_PEM,         only: soil_pem
    … …  
    3636! Arguments:
    3737
    38 character(len = *), intent(in) :: filename                            ! File name
    39 integer,            intent(in) :: timelen                             ! number of times stored in the file
    40 integer                        :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
     38character(*), intent(in) :: filename                            ! File name
     39integer,      intent(in) :: timelen                             ! number of times stored in the file
     40integer,      intent(in) :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
    4141! Ouputs
    4242real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice  per slope of the year [kg/m^2]
    4343real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
    44 real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
     44real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_PCM_phys ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
    4545real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
    4646real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
    4747real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
    48 real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per  slope of the year [kg/m^2]
     48real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per slope of the year [kg/m^2]
    4949!SOIL
    50 real, dimension(ngrid,nslope),                 intent(out) :: tsurf_ave             ! Average surface temperature of the concatenated file [K]
    51 real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
    52 real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
    53 real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
    54 real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
     50real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg             ! Average surface temperature of the concatenated file [K]
     51real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_avg             ! Average soil temperature of the concatenated file [K]
     52real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_PCM             ! Surface temperature of the concatenated file, time series [K]
     53real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_PCM             ! Soil temperature of the concatenated file, time series [K]
     54real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_avg ! Water density at the surface [kg/m^3]
    5555real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
    5656!=======================================================================
    … …  
    6767real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
    6868real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: perennial_co2ice
    69 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm               ! CO2 volume mixing ratio in the first layer [mol/m^3]
    70 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                    ! Surface Pressure [Pa]
     69real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_PCM               ! CO2 volume mixing ratio in the first layer [mol/m^3]
     70real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_PCM                    ! Surface Pressure [Pa]
    7171real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
    7272real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
    73 real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn             ! Average surface temperature of the concatenated file [K]
    74 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn             ! Average soil temperature of the concatenated file [K]
    75 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn             ! Surface temperature of the concatenated file, time series [K]
    76 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn             ! Soil temperature of the concatenated file, time series [K]
     73real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_avg_dyn             ! Average surface temperature of the concatenated file [K]
     74real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_avg_dyn             ! Average soil temperature of the concatenated file [K]
     75real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_PCM_dyn             ! Surface temperature of the concatenated file, time series [K]
     76real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_PCM_dyn             ! Soil temperature of the concatenated file, time series [K]
    7777real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn                 ! CO2 mass mixing ratio in the first layer [kg/m^3]
    7878real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn                 ! H2O mass mixing ratio in the first layer [kg/m^3]
    7979real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn         ! co2 ice amount per  slope of the year [kg/m^2]
    8080real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn     ! Water density at the surface, time series [kg/m^3]
    81 real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_ave ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
     81real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: watersurf_density_dyn_avg ! Water density at the surface, dynamic grid, yearly averaged [kg/m^3]
    8282real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn     ! Water density in the soil layer, time series [kg/m^3]
    8383real, dimension(ngrid,nslope,timelen)                               :: watersurf_density         ! Water density at the surface, physical grid, time series [kg/m^3]
    … …  
    8888
    8989! Dowload the data from the file
    90 call get_var3("ps",ps_GCM)
     90call get_var3("ps",ps_PCM)
    9191write(*,*) "Data for surface pressure downloaded!"
    9292
    … …  
    104104    write(*,*) "Data for h2o_ice_s downloaded!"
    105105
    106     call get_var3("tsurf",tsurf_gcm_dyn(:,:,1,:))
     106    call get_var3("tsurf",tsurf_PCM_dyn(:,:,1,:))
    107107    write(*,*) "Data for tsurf downloaded!"
    108108
    … …  
    115115
    116116    if (soil_pem) then
    117         call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
     117        call get_var4("soiltemp",tsoil_PCM_dyn(:,:,:,1,:))
    118118        write(*,*) "Data for soiltemp downloaded!"
    119119
    … …  
    140140    do islope = 1,nslope
    141141        write(num,'(i2.2)') islope
    142         call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
     142        call get_var3("tsurf_slope"//num,tsurf_PCM_dyn(:,:,islope,:))
    143143    enddo
    144144    write(*,*) "Data for tsurf downloaded!"
    … …  
    160160        do islope = 1,nslope
    161161            write(num,'(i2.2)') islope
    162             call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
     162            call get_var4("soiltemp_slope"//num,tsoil_PCM_dyn(:,:,:,islope,:))
    163163        enddo
    164164        write(*,*) "Data for soiltemp downloaded!"
    … …  
    182182write(*,*) "Computing the min of h2o_ice..."
    183183where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
    184 min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap,4)
     184min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
    185185write(*,*) "Computing the min of co2_ice..."
    186186where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
    187 min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn + perennial_co2ice,4)
     187min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
    188188
    189189! Compute averages over the year for each point
    190190write(*,*) "Computing the average of tsurf..."
    191 tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
     191tsurf_avg_dyn = sum(tsurf_PCM_dyn,4)/timelen
    192192
    193193if (soil_pem) then
    194194    write(*,*) "Computing average of tsoil..."
    195     tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
     195    tsoil_avg_dyn = sum(tsoil_PCM_dyn,5)/timelen
    196196    write(*,*) "Computing average of waterdensity_surface..."
    197     watersurf_density_dyn_ave(:,:,:) = sum(watersurf_density_dyn(:,:,:,:),4)/timelen
     197    watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen
    198198endif
    199199
    … …  
    215215            endif
    216216            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
    217             vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
     217            vmr_co2_PCM(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
    218218        enddo
    219219    enddo
    … …  
    221221
    222222#ifndef CPP_1D
    223     call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
    224     call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,ps_timeseries)
     223    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_PCM,vmr_co2_PCM_phys)
     224    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_PCM,ps_timeseries)
    225225    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
    226226    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
    … …  
    231231        if (soil_pem) then
    232232            do l = 1,nsoilmx
    233                 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
     233                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope))
    234234                do t = 1,timelen
    235                     call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
     235                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_PCM_dyn(:,:,l,islope,t),tsoil_PCM(:,l,islope,t))
    236236                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
    237237                enddo
    238238            enddo
    239             call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_ave(:,:,islope),watersurf_density_ave(:,islope))
     239            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope))
    240240        endif ! soil_pem
    241241        do t = 1,timelen
    242             call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
     242            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_PCM_dyn(:,:,islope,t),tsurf_PCM(:,islope,t))
    243243            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
    244244        enddo
    245245    enddo
    246246
    247     call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_ave_dyn,tsurf_ave)
     247    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
    248248#else
    249     vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
    250     ps_timeseries(1,:) = ps_GCM(1,1,:)
     249    vmr_co2_PCM_phys(1,:) = vmr_co2_PCM(1,1,:)
     250    ps_timeseries(1,:) = ps_PCM(1,1,:)
    251251    q_co2(1,:) = q_co2_dyn(1,1,:)
    252252    q_h2o(1,:) = q_h2o_dyn(1,1,:)
    … …  
    254254    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
    255255    if (soil_pem) then
    256         tsoil_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
    257         tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
     256        tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:)
     257        tsoil_PCM(1,:,:,:) = tsoil_PCM_dyn(1,1,:,:,:)
    258258        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
    259         watersurf_density_ave(1,:) = watersurf_density_dyn_ave(1,1,:)
     259        watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:)
    260260    endif ! soil_pem
    261     tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
     261    tsurf_PCM(1,:,:) = tsurf_PCM_dyn(1,1,:,:)
    262262    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
    263     tsurf_ave(1,:) = tsurf_ave_dyn(1,1,:)
     263    tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:)
    264264#endif
    265265

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

    r3130 r3149  
    88
    99SUBROUTINE recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope,emissivity_slope, &
    10                                  vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new)
     10                                 vmr_co2_PCM,vmr_co2_pem,ps_PCM_2,global_avg_press_PCM,global_avg_press_new)
    1111
    12 use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB,Lco2, sols_per_my, sec_per_sol
     12use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol
    1313
    1414implicit none
    … …  
    2424!   INPUT
    2525integer,                        intent(in) :: timelen, ngrid, nslope
    26 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_gcm                 ! physical point field : Volume mixing ratio of co2 in the first layer
     26real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM                 ! physical point field : Volume mixing ratio of co2 in the first layer
    2727real, dimension(ngrid,timelen), intent(in) :: vmr_co2_pem                 ! physical point field : Volume mixing ratio of co2 in the first layer
    28 real, dimension(ngrid,timelen), intent(in) :: ps_GCM_2                    ! physical point field : Surface pressure in the GCM
    29 real,                           intent(in) :: global_ave_press_GCM        ! global averaged pressure at previous timestep
    30 real,                           intent(in) :: global_ave_press_new        ! global averaged pressure at current timestep
     28real, dimension(ngrid,timelen), intent(in) :: ps_PCM_2                    ! physical point field : Surface pressure in the PCM
     29real,                           intent(in) :: global_avg_press_PCM        ! global averaged pressure at previous timestep
     30real,                           intent(in) :: global_avg_press_new        ! global averaged pressure at current timestep
    3131real, dimension(ngrid,nslope),  intent(in) :: tendencies_co2_ice_phys_ini ! physical point field : Evolution of perennial ice over one year
    32 real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope(kg/m^2)
     32real, dimension(ngrid,nslope),  intent(in) :: co2ice_slope                ! CO2 ice per mesh and sub-grid slope (kg/m^2)
    3333real, dimension(ngrid,nslope),  intent(in) :: emissivity_slope            ! Emissivity per mesh and sub-grid slope(1)
    3434!   OUTPUT
    … …  
    4040real    :: coef, ave
    4141
     42write(*,*) "Update of the CO2 tendency from the current pressure"
     43   
    4244! Evolution of the water ice for each physical point
    4345do i = 1,ngrid
    … …  
    4749        if (co2ice_slope(i,islope) > 1.e-4 .and. abs(tendencies_co2_ice_phys(i,islope)) > 1.e-5) then
    4850            do t=1,timelen
    49                 ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4 &
    50                       - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4
     51                ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM_2(i,t)/100.)))**4 &
     52                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem(i,t)*ps_PCM_2(i,t)*(global_avg_press_new/global_avg_press_PCM)/100.)))**4
    5153            enddo
    5254        endif

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

    r3123 r3149  
    77!=======================================================================
    88
    9 SUBROUTINE soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_GCM,TI_GCM,TI_PEM)
     9SUBROUTINE soil_settings_PEM(ngrid,nslope,nsoil_PEM,nsoil_PCM,TI_PCM,TI_PEM)
    1010
    1111!      use netcdf
    … …  
    3939integer,                                 intent(in) :: nslope    ! # of subslope wihtin the mesh
    4040integer,                                 intent(in) :: nsoil_PEM ! # of soil layers in the PEM
    41 integer,                                 intent(in) :: nsoil_GCM ! # of soil layers in the GCM
    42 real, dimension(ngrid,nsoil_GCM,nslope), intent(in) :: TI_GCM    ! Thermal inertia  in the GCM [SI]
     41integer,                                 intent(in) :: nsoil_PCM ! # of soil layers in the PCM
     42real, dimension(ngrid,nsoil_PCM,nslope), intent(in) :: TI_PCM    ! Thermal inertia  in the PCM [SI]
    4343
    4444real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM ! Thermal inertia   in the PEM [SI]
    … …  
    5757lay1 = 2.e-4
    5858alpha = 2
    59 do iloop = 0,nsoil_GCM - 1
     59do iloop = 0,nsoil_PCM - 1
    6060    mlayer_PEM(iloop) = lay1*(1+iloop**2.9*(1-exp(-real(iloop)/20.)))
    6161enddo
    6262
    6363do iloop = 0, nsoil_PEM-1
    64     if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_GCM-1)) then
     64    if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM-1)) then
    6565        index_powerlaw = iloop
    6666        exit
    … …  
    6868enddo
    6969
    70 do iloop = nsoil_GCM, nsoil_PEM-1
    71     mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop-nsoil_GCM)-0.5))
     70do iloop = nsoil_PCM, nsoil_PEM-1
     71    mlayer_PEM(iloop) = lay1*(alpha**(index_powerlaw + (iloop-nsoil_PCM)-0.5))
    7272enddo
    7373
    … …  
    8383do ig = 1,ngrid
    8484    do islope = 1,nslope
    85         do iloop = 1,nsoil_GCM
    86             TI_PEM(ig,iloop,islope) = TI_GCM(ig,iloop,islope)
     85        do iloop = 1,nsoil_PCM
     86            TI_PEM(ig,iloop,islope) = TI_PCM(ig,iloop,islope)
    8787        enddo
    88         if (nsoil_PEM > nsoil_GCM) then
    89             do iloop = nsoil_GCM + 1,nsoil_PEM
    90                 TI_PEM(ig,iloop,islope) = TI_GCM(ig,nsoil_GCM,islope)
     88        if (nsoil_PEM > nsoil_PCM) then
     89            do iloop = nsoil_PCM + 1,nsoil_PEM
     90                TI_PEM(ig,iloop,islope) = TI_PCM(ig,nsoil_PCM,islope)
    9191            enddo
    9292        endif

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

    r3130 r3149  
    8989 REAL :: ice_bottom_depth
    9090
     91write(*,*) "Update soil propreties"
     92       
    9193! 1. Modification of the regolith thermal inertia.
    92 
    93 
    94 
    9594 do islope = 1,nslope
    9695   regolith_inertia(:,islope) = inertiedat_PEM(:,1)
    … …  
    112111
    113112! 2. Build new Thermal Inertia
    114 
    115113do  islope=1,nslope
    116114   do ig = 1,ngrid

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

    r3148 r3149  
    1 MODULE criterion_pem_stop_mod
     1MODULE stopping_crit_mod
    22
    33implicit none
    44
    5 ! ----------------------------------------------------------------------
     5!=======================================================================
    66contains
    7 ! ----------------------------------------------------------------------
     7!=======================================================================
    88
    99!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    1414!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1515
    16 SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)
     16SUBROUTINE stopping_crit_h2o_ice(cell_area,surf_ini,h2o_ice,stopPEM,ngrid)
    1717
    1818use time_evol_mod, only: water_ice_criterion
    … …  
    2323!=======================================================================
    2424!
    25 !  Routine that checks if the water ice criterion to stop the PEM is reached
     25! Routine to check if the h2o ice criterion to stop the PEM is reached
    2626!
    2727!=======================================================================
    … …  
    3030!   ----------
    3131!   INPUT
    32 integer,                       intent(in) :: ngrid           ! # of physical grid points
    33 real, dimension(ngrid),        intent(in) :: cell_area       ! Area of the cells
    34 real, dimension(ngrid,nslope), intent(in) :: qsurf           ! Actual density of water ice
    35 real,                          intent(in) :: ini_surf        ! Initial surface of h2o ice that was sublimating
    36 real, dimension(ngrid,nslope), intent(in) :: initial_h2o_ice ! Grid point that initialy were covered by h2o_ice
     32integer,                       intent(in) :: ngrid     ! # of physical grid points
     33real, dimension(ngrid),        intent(in) :: cell_area ! Area of the cells
     34real, dimension(ngrid,nslope), intent(in) :: h2o_ice   ! Actual density of water ice
     35real,                          intent(in) :: surf_ini  ! Initial surface of h2o ice that was sublimating
    3736!   OUTPUT
    38 logical, intent(out) :: STOPPING ! Is the criterion reached?
     37integer, intent(inout) :: stopPEM ! Stopping criterion code
    3938!   local:
    4039!   ------
    41 integer :: i, islope    ! Loop
    42 real    :: present_surf ! Initial/Actual surface of water ice
     40integer :: i, islope ! Loop
     41real    :: surf_now  ! Current surface of h2o ice
    4342
    4443!=======================================================================
    45 ! Initialisation to false
    46 STOPPING = .false.
    47 
    4844! Computation of the present surface of h2o ice
    49 present_surf = 0.
     45surf_now = 0.
    5046do i = 1,ngrid
    5147    do islope = 1,nslope
    52         if (initial_h2o_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
     48        if (h2o_ice(i,islope) > 0.) surf_now = surf_now + cell_area(i)*subslope_dist(i,islope)
    5349    enddo
    5450enddo
    5551
    5652! Check of the criterion
    57 if (present_surf < ini_surf*(1. - water_ice_criterion)) then
    58     STOPPING = .true.
    59     write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
    60     write(*,*) "present_surf < ini_surf*(1. - water_ice_criterion)", present_surf < ini_surf*(1. - water_ice_criterion)
    61     write(*,*) "Current surface of water ice sublimating =", present_surf
    62     write(*,*) "Initial surface of water ice sublimating =", ini_surf
     53if (surf_now < surf_ini*(1. - water_ice_criterion)) then
     54    stopPEM = 1
     55    write(*,*) "Reason of stopping: the surface of water ice sublimating reaches the threshold"
     56    write(*,*) "surf_now < surf_ini*(1. - water_ice_criterion)", surf_now < surf_ini*(1. - water_ice_criterion)
     57    write(*,*) "Current surface of water ice sublimating =", surf_now
     58    write(*,*) "Initial surface of water ice sublimating =", surf_ini
    6359    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
    64 endif
    65 if (present_surf > ini_surf*(1. + water_ice_criterion)) then
    66     STOPPING = .true.
    67     write(*,*) "Reason of stopping: the surface of water ice sublimating reach the threshold"
    68     write(*,*) "present_surf > ini_surf*(1. + water_ice_criterion)", present_surf > ini_surf*(1. + water_ice_criterion)
    69     write(*,*) "Current surface of water ice sublimating =", present_surf
    70     write(*,*) "Initial surface of water ice sublimating =", ini_surf
     60else if (surf_now > surf_ini*(1. + water_ice_criterion)) then
     61    stopPEM = 1
     62    write(*,*) "Reason of stopping: the surface of water ice sublimating reaches the threshold"
     63    write(*,*) "surf_now > surf_ini*(1. + water_ice_criterion)", surf_now > surf_ini*(1. + water_ice_criterion)
     64    write(*,*) "Current surface of water ice sublimating =", surf_now
     65    write(*,*) "Initial surface of water ice sublimating =", surf_ini
    7166    write(*,*) "Percentage of change accepted =", water_ice_criterion*100
    7267endif
    7368
    74 if (ini_surf < 1.e-5 .and. ini_surf > -1.e-5) STOPPING = .false.
     69if (abs(surf_ini) < 1.e-5) stopPEM = 0
    7570
    76 END SUBROUTINE criterion_waterice_stop
     71END SUBROUTINE stopping_crit_h2o_ice
    7772
    78 ! ----------------------------------------------------------------------
     73!=======================================================================
    7974
    80 SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_co2_ice,global_ave_press_GCM,global_ave_press_new,nslope)
     75SUBROUTINE stopping_crit_co2(cell_area,surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
    8176
    8277use time_evol_mod, only: co2_ice_criterion, ps_criterion
    … …  
    8782!=======================================================================
    8883!
    89 !  Routine that checks if the co2 and pressure criteria to stop the PEM are reached
     84! Routine to check if the co2 and pressure criteria to stop the PEM are reached
    9085!
    9186!=======================================================================
    … …  
    9691integer,                       intent(in) :: ngrid, nslope        ! # of grid physical grid points
    9792real, dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
    98 real, dimension(ngrid,nslope), intent(in) :: qsurf                ! Actual density of water ice
    99 real,                          intent(in) :: ini_surf             ! Initial surface of co2 ice that was sublimating
    100 real, dimension(ngrid,nslope), intent(in) :: initial_co2_ice      ! Grid point that initialy were covered by co2_ice
    101 real,                          intent(in) :: global_ave_press_GCM ! Planet average pressure from the PCM start files
    102 real,                          intent(in) :: global_ave_press_new ! Planet average pressure from the PEM computations
     93real, dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of water ice
     94real,                          intent(in) :: surf_ini             ! Initial surface of co2 ice that was sublimating
     95real,                          intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
     96real,                          intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations
    10397!   OUTPUT
    104 logical, intent(out) :: STOPPING_ice ! Is the criterion for co2 ice reached?
    105 logical, intent(out) :: STOPPING_ps  ! Is the criterion for pressure reached?
     98integer, intent(inout) :: stopPEM ! Stopping criterion code
     99
    106100!   local:
    107101!   ------
    108 integer :: i, islope    ! Loop
    109 real    :: present_surf ! Initial/Actual surface of water ice
     102integer :: i, islope ! Loop
     103real    :: surf_now  ! Current surface of co2 ice
    110104
    111105!=======================================================================
    112 ! Initialisation to false
    113 STOPPING_ice = .false.
    114 STOPPING_ps = .false.
    115 
    116106! Computation of the present surface of co2 ice
    117 present_surf = 0.
     107surf_now = 0.
    118108do i = 1,ngrid
    119109    do islope = 1,nslope
    120         if (initial_co2_ice(i,islope) > 0.5 .and. qsurf(i,islope) > 0.) present_surf = present_surf + cell_area(i)*subslope_dist(i,islope)
     110        if (co2_ice(i,islope) > 0.) surf_now = surf_now + cell_area(i)*subslope_dist(i,islope)
    121111    enddo
    122112enddo
    123113
    124114! Check of the criterion
    125 if (present_surf < ini_surf*(1. - co2_ice_criterion)) then
    126     STOPPING_ice = .true.
    127     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reached the threshold"
    128     write(*,*) "present_surf < ini_surf*(1. - co2_ice_criterion)", present_surf < ini_surf*(1. - co2_ice_criterion)
    129     write(*,*) "Current surface of co2 ice sublimating =", present_surf
    130     write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
     115if (surf_now < surf_ini*(1. - co2_ice_criterion)) then
     116    stopPEM = 3
     117    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
     118    write(*,*) "surf_now < surf_ini*(1. - co2_ice_criterion)", surf_now < surf_ini*(1. - co2_ice_criterion)
     119    write(*,*) "Current surface of co2 ice sublimating =", surf_now
     120    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
    131121    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
    132 endif
    133 if (present_surf > ini_surf*(1. + co2_ice_criterion)) then
    134     STOPPING_ice = .true.
    135     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reach the threshold"
    136     write(*,*) "present_surf > ini_surf*(1. + co2_ice_criterion)", present_surf > ini_surf*(1. + co2_ice_criterion)
    137     write(*,*) "Current surface of co2 ice sublimating =", present_surf
    138     write(*,*) "Initial surface of co2 ice sublimating =", ini_surf
     122else if (surf_now > surf_ini*(1. + co2_ice_criterion)) then
     123    stopPEM = 3
     124    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
     125    write(*,*) "surf_now > surf_ini*(1. + co2_ice_criterion)", surf_now > surf_ini*(1. + co2_ice_criterion)
     126    write(*,*) "Current surface of co2 ice sublimating =", surf_now
     127    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
    139128    write(*,*) "Percentage of change accepted =", co2_ice_criterion*100.
    140129endif
    141130
    142 if (abs(ini_surf) < 1.e-5) STOPPING_ice = .false.
     131if (abs(surf_ini) < 1.e-5) stopPEM = .false.
    143132
    144 if (global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)) then
    145     STOPPING_ps = .true.
    146     write(*,*) "Reason of stopping: the global pressure reach the threshold"
    147     write(*,*) "global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)", global_ave_press_new < global_ave_press_GCM*(1. - ps_criterion)
    148     write(*,*) "Current global pressure =", global_ave_press_new
    149     write(*,*) "PCM global pressure =", global_ave_press_GCM
     133if (global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)) then
     134    stopPEM = 4
     135    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
     136    write(*,*) "global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)", global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)
     137    write(*,*) "Current global pressure =", global_avg_press_new
     138    write(*,*) "PCM global pressure =", global_avg_press_PCM
    150139    write(*,*) "Percentage of change accepted =", ps_criterion*100.
    151 endif
    152 if (global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)) then
    153     STOPPING_ps = .true.
    154     write(*,*) "Reason of stopping: the global pressure reach the threshold"
    155     write(*,*) "global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)", global_ave_press_new > global_ave_press_GCM*(1. + ps_criterion)
    156     write(*,*) "Current global pressure =", global_ave_press_new
    157     write(*,*) "PCM global pressure =", global_ave_press_GCM
     140else if (global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)) then
     141    stopPEM = 4
     142    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
     143    write(*,*) "global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)", global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)
     144    write(*,*) "Current global pressure =", global_avg_press_new
     145    write(*,*) "PCM global pressure =", global_avg_press_PCM
    158146    write(*,*) "Percentage of change accepted =", ps_criterion*100.
    159147endif
    160148
    161 END SUBROUTINE criterion_co2_stop
     149END SUBROUTINE stopping_crit_co2
    162150
    163151END MODULE

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

    r3076 r3149  
    1212logical :: evol_orbit_pem      ! True if we want to follow the orbital parameters of obl_ecc_lsp.asc, read in evol.def
    1313logical :: var_obl             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for obliquity
    14 logical :: var_ecc             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for excenticity
     14logical :: var_ecc             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for eccenticity
    1515logical :: var_lsp             ! True if we want the PEM to follow obl_ecc_lsp.asc parameters for ls perihelie
    1616
Note: See TracChangeset for help on using the changeset viewer.