Changeset 3028


Ignore:
Timestamp:
Aug 14, 2023, 3:22:41 PM (17 months ago)
Author:
jbclement
Message:

Mars PEM:
Big cleaning of main program pem.F90 (indentation, declarations, comments, simplification of conditions/loops, etc).
JBC

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3022 r3028  
    1 
    2 !------------------------
    3 
    4 ! I   Initialisation
    5 !    I_a READ run.def
     1!------------------------
     2! I   Initialization
     3!    I_a READ run.def
    64!    I_b READ of start_evol.nc and starfi_evol.nc
    75!    I_c Subslope parametrisation
    8 !    I_d READ GCM data and convert to the physical grid 
    9 !    I_e Initialisation of the PEM variable and soil
     6!    I_d READ GCM data and convert to the physical grid
     7!    I_e Initialization of the PEM variable and soil
    108!    I_f Compute tendencies & Save initial situation
    119!    I_g Save initial PCM situation
     
    1412
    1513! II  Run
    16 !    II_a update pressure, ice and tracers
     14!    II_a Update pressure, ice and tracers
    1715!    II_b Evolution of the ice
    18 !    II_c CO2 glaciers flows
     16!    II_c CO2 & H2O glaciers flows
    1917!    II_d Update surface and soil temperatures
    20 !    II_e Update the tendencies 
     18!    II_e Update the tendencies
    2119!    II_f Checking the stopping criterion
    2220
    2321! III Output
    2422!    III_a Update surface value for the PCM start files
    25 !    III_b Write start and starfi.nc
    26 !    III_c Write start_pem
    27 
     23!    III_b Write restart_evol.nc and restartfi_evol.nc
     24!    III_c Write restartfi_PEM.nc
    2825!------------------------
    2926
    3027PROGRAM pem
    3128
    32 ! 1: Modules needed for reading and writting startfi:
    33       use phyetat0_mod, only: phyetat0
    34       use phyredem,     only: physdem0, physdem1
    35       use netcdf,       only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
    36                          nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
    37                          nf90_inquire_dimension,nf90_close
    38       ! netcdf is needed to read info like lat and lon in the physiq file
    39       use turb_mod,     only: q2, wstar
    40 ! 1a: Modules specific from the Marsian physiq
     29use phyetat0_mod,               only: phyetat0
     30use phyredem,                   only: physdem0, physdem1
     31use netcdf,                     only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
     32use turb_mod,                   only: q2, wstar
     33use comslope_mod,               only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
     34use logic_mod,                  only: iflag_phys
     35use mod_const_mpi,              only: COMM_LMDZ
     36use comconst_mod,               only: rad, g, cpp, pi
     37use infotrac
     38use geometry_mod,               only: latitude_deg
     39use conf_pem_mod,               only: conf_pem
     40use pemredem,                   only: pemdem0, pemdem1
     41use glaciers_mod,               only: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
     42use criterion_pem_stop_mod,     only: criterion_waterice_stop, criterion_co2_stop
     43use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
     44                                      m_noco2, threshold_water_frost2perenial, threshold_co2_frost2perenial
     45use evol_co2_ice_s_mod,         only: evol_co2_ice_s
     46use evol_h2o_ice_s_mod,         only: evol_h2o_ice_s
     47use 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 GCM
     51                                      water_reservoir                     ! Water ressources
     52use 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
     55use temps_mod_evol,             only: dt_pem, evol_orbit_pem, Max_iter_pem
     56use orbit_param_criterion_mod,  only: orbit_param_criterion
     57use recomp_orb_param_mod,       only: recomp_orb_param
     58use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
     59                                      ini_ice_table_porefilling, computeice_table_equilibrium
     60use soil_thermalproperties_mod, only: update_soil_thermalproperties
     61
    4162#ifndef CPP_STD
    42       use comsoil_h,    only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa,inertiesoil
    43       use surfdat_h,    only: tsurf, emis,&
    44                            qsurf,watercap, ini_surfdat_h, &
    45                            albedodat, zmea, zstd, zsig, zgam, zthe, &
    46                            hmons, summit, base,albedo_h2o_frost, &
    47                            frost_albedo_threshold,emissiv,watercaptag,perenial_co2ice
    48       use dimradmars_mod, only: totcloudfrac, albedo
    49       use dust_param_mod, only: tauscaling
    50       use tracer_mod,     only: noms,igcm_h2o_ice,igcm_co2,mmol,igcm_h2o_vap ! tracer names and molar masses
    51 #else
    52 ! 1b: Modules specific from the Generic physiq
    53       use comsoil_h,    only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
    54       use surfdat_h,    only: albedodat, zmea, zstd, zsig, zgam, zthe, &
    55                            emissiv
    56       use tracer_h,     only: noms,igcm_h2o_ice,igcm_co2 ! tracer names
    57       use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, &
    58                                     PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE,          &
    59                                     ALBEDO_BAREGROUND,ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
    60       use aerosol_mod,  only : iniaerosol
     63    use comsoil_h,          only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil
     64    use surfdat_h,          only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
     65                                  albedodat, zmea, zstd, zsig, zgam, zthe,    &
     66                                  hmons, summit, base,albedo_h2o_frost,        &
     67                                  frost_albedo_threshold, emissiv, watercaptag, perenial_co2ice
     68    use dimradmars_mod,     only: totcloudfrac, albedo
     69    use dust_param_mod,     only: tauscaling
     70    use tracer_mod,         only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
     71    use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
     72    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit
     73    use comcstfi_h,         only: r, mugaz
     74#else
     75    use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
     76    use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT,   &
     77                                  PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
     78                                  ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
     79    use aerosol_mod,  only: iniaerosol
     80    use planete_mod,  only: apoastr, periastr, year_day, peri_day, obliquit
     81    use comcstfi_mod, only: r, mugaz
    6182#endif
    62 ! 1c: Modules specific from the 1d marsian physiq
     83
    6384#ifndef CPP_1D
    64       USE iniphysiq_mod, ONLY: iniphysiq
    65       USE control_mod,   ONLY: iphysiq, day_step,nsplit_phys
     85    use iniphysiq_mod,            only: iniphysiq
     86    use control_mod,              only: iphysiq, day_step, nsplit_phys
    6687#else
    67       use time_phylmdz_mod, only: daysec, iphysiq,day_step
     88    use time_phylmdz_mod,         only: daysec, iphysiq, day_step, dtphys
    6889#endif
    6990
     91#ifdef CPP_1D
     92    use regular_lonlat_mod,       only: init_regular_lonlat
     93    use physics_distribution_mod, only: init_physics_distribution
     94    use mod_grid_phy_lmdz,        only: regular_lonlat
     95    use init_phys_1d_mod,         only: init_phys_1d
     96#endif
     97
     98IMPLICIT NONE
     99
     100include "dimensions.h"
     101include "paramet.h"
     102include "comgeom.h"
     103include "iniprint.h"
     104
     105integer ngridmx
     106parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
     107
     108! Same variable names as in the GCM
     109integer :: ngrid     ! Number of physical grid points
     110integer :: nlayer    ! Number of vertical layer
     111integer :: nq        ! Number of tracer
     112integer :: day_ini   ! First day of the simulation
     113real    :: pday      ! Physical day
     114real    :: time_phys ! Same as GCM
     115real    :: ptimestep ! Same as GCM
     116real    :: ztime_fin ! Same as GCM
     117
     118! Variables to read start.nc
     119character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM
     120
     121! Dynamic variables
     122real                                :: vcov(ip1jm,llm), ucov(ip1jmp1,llm) ! vents covariants
     123real                                :: teta(ip1jmp1,llm)                  ! temperature potentielle
     124real, dimension(:,:,:), allocatable :: q                                  ! champs advectes
     125real                                :: ps(ip1jmp1)                        ! pression  au sol
     126real, dimension(:),     allocatable :: ps_start_GCM                       !(ngrid) pression  au sol
     127real, dimension(:,:),   allocatable :: ps_timeseries                      !(ngrid x timelen) ! pression  au sol instantannées
     128real                                :: masse(ip1jmp1,llm)                 ! masse d'air
     129real                                :: phis(ip1jmp1)                      ! geopotentiel au sol
     130real                                :: time_0
     131
     132! Variables to read starfi.nc
     133character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
     134character*2 str2
     135integer :: ncid, varid,status                       ! Variable for handling opening of files
     136integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files
     137integer :: lonvarid, latvarid, areavarid, sdvarid   ! Variable ID for Netcdf files
     138integer :: apvarid, bpvarid                         ! Variable ID for Netcdf files
     139
     140! Variables to read starfi.nc and write restartfi.nc
     141real, dimension(:), allocatable :: longitude     ! Longitude read in FILE_NAME and written in restartfi
     142real, dimension(:), allocatable :: latitude      ! Latitude read in FILE_NAME and written in restartfi
     143real, dimension(:), allocatable :: ap            ! Coefficient ap read in FILE_NAME_start and written in restart
     144real, dimension(:), allocatable :: bp            ! Coefficient bp read in FILE_NAME_start and written in restart
     145real, dimension(:), allocatable :: cell_area     ! Cell_area read in FILE_NAME and written in restartfi
     146real                            :: Total_surface ! Total surface of the planet
     147
     148! Variables for h2o_ice evolution
     149real                                :: ini_surf_h2o          ! Initial surface of sublimating h2o ice
     150real                                :: ini_surf_co2          ! Initial surface of sublimating co2 ice
     151real                                :: global_ave_press_GCM  ! constant: global average pressure retrieved in the GCM [Pa]
     152real                                :: global_ave_press_old  ! constant: Global average pressure of initial/previous time step
     153real                                :: global_ave_press_new  ! constant: Global average pressure of current time step
     154real, dimension(:,:),   allocatable ::  zplev_new            ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
     155real, dimension(:,:),   allocatable ::  zplev_gcm            ! same but retrieved from the gcm [kg/m^2]
     156real, dimension(:,:,:), allocatable ::  zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
     157real, dimension(:,:,:), allocatable :: zplev_old_timeseries  ! same but with the time series, for oldest time step
     158logical                             :: STOPPING_water        ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
     159logical                             :: STOPPING_1_water      ! Logical : is there still water ice to sublimate?
     160logical                             :: STOPPING_co2          ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
     161logical                             :: STOPPING_pressure     ! Logical : is the criterion (% of change in the surface pressure) reached?
     162integer                             :: criterion_stop        ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
     163real, save                          :: A , B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
     164real, allocatable                   :: vmr_co2_gcm(:,:)      ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
     165real, allocatable                   :: vmr_co2_pem_phys(:,:) ! Physics x Times  co2 volume mixing ratio used in the PEM
     166real, allocatable                   :: q_co2_PEM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]
     167real, allocatable                   :: q_h2o_PEM_phys(:,:)   ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
     168integer                             :: timelen               ! # time samples
     169real                                :: ave                   ! intermediate varibale to compute average
     170real, allocatable                   :: p(:,:)                ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
     171real                                :: extra_mass            ! Intermediate variables Extra mass of a tracer if it is greater than 1
     172
     173! Variables for slopes
     174real, dimension(:,:),   allocatable :: min_co2_ice_1          ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
     175real, dimension(:,:),   allocatable :: min_co2_ice_2          ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
     176real, dimension(:,:),   allocatable :: min_h2o_ice_1          ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
     177real, dimension(:,:),   allocatable :: min_h2o_ice_2          ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
     178real, dimension(:,:,:), allocatable :: co2_ice_GCM            ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
     179real, dimension(:,:),   allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice
     180real, dimension(:,:),   allocatable :: initial_h2o_ice        ! physical point field : Logical array indicating if there is water ice at initial state
     181real, dimension(:,:),   allocatable :: initial_co2_ice        ! physical point field : Logical array indicating if there is co2 ice at initial state
     182real, dimension(:,:),   allocatable :: tendencies_co2_ice     ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
     183real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
     184real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
     185real, dimension(:,:),   allocatable :: flag_co2flow(:,:)      !(ngrid,nslope): Flag where there is a CO2 glacier flow
     186real, dimension(:),     allocatable :: flag_co2flow_mesh(:)   !(ngrid)       : Flag where there is a CO2 glacier flow
     187real, dimension(:,:),   allocatable :: flag_h2oflow(:,:)      !(ngrid,nslope): Flag where there is a H2O glacier flow
     188real, dimension(:),     allocatable :: flag_h2oflow_mesh(:)   !(ngrid)       : Flag where there is a H2O glacier flow
     189
     190! Variables for surface and soil
     191real, allocatable :: tsurf_ave(:,:)              ! Physic x SLOPE field : Averaged Surface Temperature [K]
     192real, allocatable :: tsoil_ave(:,:,:)            ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
     193real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
     194real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
     195real, allocatable :: tsoil_GCM_timeseries(:,:,:,:)      !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
     196real, allocatable :: tsurf_ave_yr1(:,:)                 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
     197real, allocatable :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
     198real, allocatable :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
     199real, allocatable :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
     200real, allocatable :: watersoil_density_timeseries(:,:,:,:)     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
     201real, allocatable :: watersurf_density_ave(:,:)                ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
     202real, allocatable :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
     203real, allocatable :: watersoil_density_PEM_ave(:,:,:)          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
     204real, allocatable :: Tsurfave_before_saved(:,:)                ! Surface temperature saved from previous time step [K]
     205real, allocatable :: delta_co2_adsorbded(:)                    ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
     206real, allocatable :: delta_h2o_adsorbded(:)                    ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
     207real              :: totmassco2_adsorbded                      ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
     208real              :: totmassh2o_adsorbded                      ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
     209logical           :: bool_sublim                               ! logical to check if there is sublimation or not
     210
     211! Some variables for the PEM run
     212real, parameter :: year_step = 1 ! timestep for the pem
     213integer         :: year_iter     ! number of iteration
     214integer         :: year_iter_max ! maximum number of iterations before stopping
     215real            :: timestep      ! timestep [s]
     216real            :: watercap_sum  ! total mass of water cap [kg/m^2]
     217real            :: water_sum     ! total mass of water in the mesh [kg/m^2]
     218
     219#ifdef CPP_STD
     220    real                 :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
     221    real                 :: albedo_h2o_frost              ! albedo of h2o frost
     222    real,    allocatable :: tsurf_read_generic(:)         ! Temporary variable to do the subslope transfert dimensiion when reading form generic
     223    real,    allocatable :: qsurf_read_generic(:,:)       ! Temporary variable to do the subslope transfert dimensiion when reading form generic
     224    real,    allocatable :: tsoil_read_generic(:,:)       ! Temporary variable to do the subslope transfert dimensiion when reading form generic
     225    real,    allocatable :: emis_read_generic(:)          ! Temporary variable to do the subslope transfert dimensiion when reading form generic
     226    real,    allocatable :: albedo_read_generic(:,:)      ! Temporary variable to do the subslope transfert dimensiion when reading form generic
     227    real,    allocatable :: tsurf(:,:)                    ! Subslope variable, only needed in the GENERIC case
     228    real,    allocatable :: qsurf(:,:,:)                  ! Subslope variable, only needed in the GENERIC case
     229    real,    allocatable :: tsoil(:,:,:)                  ! Subslope variable, only needed in the GENERIC case
     230    real,    allocatable :: emis(:,:)                     ! Subslope variable, only needed in the GENERIC case
     231    real,    allocatable :: watercap(:,:)                 ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
     232    logical, allocatable :: WATERCAPTAG(:)                ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
     233    real,    allocatable :: albedo(:,:,:)                 ! Subslope variable, only needed in the GENERIC case
     234    real,    allocatable :: inertiesoil(:,:,:)            ! Subslope variable, only needed in the GENERIC case
     235#endif
     236
     237#ifdef CPP_1D
     238    integer            :: nsplit_phys
     239    integer, parameter :: jjm_value = jjm - 1
     240    integer            :: ierr
     241#else
     242    integer, parameter :: jjm_value = jjm
     243#endif
     244
     245! Loop variables
     246integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil
     247
     248! Parallel variables
    70249#ifndef CPP_STD
    71       use mod_phys_lmdz_para, only: is_parallel, is_sequential, &
    72                                    is_mpi_root, is_omp_root,    &
    73                                    is_master
    74       use planete_h,     only: aphelie, periheli, year_day, peri_day, &
    75                           obliquit
     250    is_sequential = .true.
     251    is_parallel = .false.
     252    is_mpi_root = .true.
     253    is_omp_root = .true.
     254    is_master = .true.
     255#endif
     256
     257day_ini = 0    ! test
     258time_phys = 0. ! test
     259
     260! Some constants
     261ngrid = ngridmx
     262nlayer = llm
     263
     264A = (1/m_co2 - 1/m_noco2)
     265B = 1/m_noco2
     266
     267year_day = 669
     268daysec = 88775.
     269timestep = year_day*daysec/year_step
     270
     271!----------------------------- INITIALIZATION --------------------------
     272!------------------------
     273! I   Initialization
     274!    I_a READ run.def
     275!------------------------
     276#ifndef CPP_1D
     277    dtphys = 0
     278    call conf_gcm(99,.true.)
     279    call infotrac_init
     280    nq = nqtot
     281    allocate(q(ip1jmp1,llm,nqtot))
    76282#else
    77       use planete_mod,   only: apoastr, periastr, year_day, peri_day, &
    78                           obliquit
    79 #endif
    80 
    81       USE comslope_mod,  ONLY: nslope,def_slope,def_slope_mean, &
    82                            subslope_dist,iflat,                &
    83                            major_slope,ini_comslope_h
    84 
    85 #ifndef CPP_STD
    86       USE comcstfi_h, only: r, mugaz
    87 #else
    88       USE comcstfi_mod, only: r, mugaz
    89 #endif
    90 
    91       USE logic_mod,        ONLY: iflag_phys
    92       USE mod_const_mpi,    ONLY: COMM_LMDZ
    93       use time_phylmdz_mod, only: daysec,dtphys
    94       USE comconst_mod,     ONLY: rad,g,cpp,pi
    95       USE infotrac
    96       USE geometry_mod,     only: latitude_deg
    97 
    98       use conf_pem_mod,     only: conf_pem
    99       use pemredem,         only:  pemdem0,pemdem1
    100       use glaciers_mod,        only: co2glaciers_evol,h2oglaciers_evol,co2glaciersflow,h2oglaciersflow
    101       use criterion_pem_stop_mod, only: criterion_waterice_stop,criterion_co2_stop
    102       use constants_marspem_mod,  only: alpha_clap_co2,beta_clap_co2, alpha_clap_h2o,beta_clap_h2o, &
    103                                       m_co2,m_noco2,threshold_water_frost2perenial,threshold_co2_frost2perenial
    104       use evol_co2_ice_s_mod, only: evol_co2_ice_s
    105       use evol_h2o_ice_s_mod, only: evol_h2o_ice_s
    106       use comsoil_h_PEM,      only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
    107                                TI_PEM,inertiedat_PEM,         &                        ! soil thermal inertia         
    108                                tsoil_PEM, mlayer_PEM,layer_PEM,                  &     ! Soil temp, number of subsurface layers, soil mid layer depths
    109                                fluxgeo, &                                              ! geothermal flux for the PEM and GCM
    110                                water_reservoir                                         ! Water ressources
    111       use adsorption_mod,     only : regolith_adsorption,adsorption_pem, &                ! bool to check if adsorption, main subroutine
    112                                ini_adsorption_h_PEM, end_adsorption_h_PEM, &        ! allocate arrays                               
    113                                co2_adsorbded_phys, h2o_adsorbded_phys               ! mass of co2 and h2O adsorbded
    114       USE temps_mod_evol,            ONLY: dt_pem, evol_orbit_pem, Max_iter_pem
    115       use orbit_param_criterion_mod, only : orbit_param_criterion
    116       use recomp_orb_param_mod,      only: recomp_orb_param
    117       use ice_table_mod, only: porefillingice_depth,porefillingice_thickness,&
    118                                end_ice_table_porefilling,ini_ice_table_porefilling, &
    119                                computeice_table_equilibrium
    120       use soil_thermalproperties_mod, only: update_soil_thermalproperties
    121 
    122 #ifdef CPP_1D
    123       use regular_lonlat_mod,       only: init_regular_lonlat
    124       use physics_distribution_mod, only: init_physics_distribution
    125       use mod_grid_phy_lmdz,        only : regular_lonlat
    126       use init_phys_1d_mod,         only : init_phys_1d
    127 #endif
    128 
    129 
    130   IMPLICIT NONE
    131 
    132   include "dimensions.h"
    133   include "paramet.h"
    134 
    135   INTEGER ngridmx
    136   PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
    137 
    138   include "comgeom.h"
    139   include "iniprint.h"
    140 ! Same variable's name as in the GCM
    141 
    142   INTEGER :: ngrid      !Number of physical grid points
    143   INTEGER :: nlayer     !Number of vertical layer
    144   INTEGER :: nq         !Number of tracer
    145   INTEGER :: day_ini    !First day of the simulation
    146   REAL :: pday          !Physical day
    147   REAL :: time_phys     !Same as GCM
    148   REAL :: ptimestep     !Same as GCM
    149   REAL :: ztime_fin     !Same as GCM
    150 
    151 ! Variable for reading start.nc
    152       character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
    153   !   variables dynamiques
    154   REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    155   REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    156   REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
    157   REAL ps(ip1jmp1)                       ! pression  au sol
    158 
    159   REAL, dimension(:),allocatable :: ps_start_GCM    !(ngrid) pression  au sol
    160   REAL, dimension(:,:),allocatable :: ps_timeseries !(ngrid x timelen) ! pression  au sol instantannées
    161 
    162   REAL masse(ip1jmp1,llm)                ! masse d'air
    163   REAL phis(ip1jmp1)                     ! geopotentiel au sol
    164   REAL time_0
    165 
    166 ! Variable for reading starfi.nc
    167 
    168       character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
    169       character*2 str2
    170       integer :: ncid, varid,status                                !Variable for handling opening of files
    171       integer :: phydimid, subdimid, nlayerdimid, nqdimid          !Variable ID for Netcdf files
    172       integer :: lonvarid, latvarid, areavarid,sdvarid             !Variable ID for Netcdf files
    173       integer :: apvarid,bpvarid                                   !Variable ID for Netcdf files
    174 
    175 ! Variable for reading starfi.nc and writting restartfi.nc
    176 
    177       REAL, dimension(:),allocatable :: longitude                  !Longitude read in FILE_NAME and written in restartfi
    178       REAL, dimension(:),allocatable  :: latitude                  !Latitude read in FILE_NAME and written in restartfi
    179       REAL, dimension(:),allocatable :: ap                         !Coefficient ap read in FILE_NAME_start and written in restart
    180       REAL, dimension(:),allocatable :: bp                         !Coefficient bp read in FILE_NAME_start and written in restart
    181       REAL, dimension(:),allocatable  :: cell_area                 !Cell_area read in FILE_NAME and written in restartfi
    182       REAL :: Total_surface                                        !Total surface of the planet
    183 
    184 ! Variable for h2o_ice evolution
    185       REAL :: ini_surf_h2o                                         ! Initial surface of sublimating h2o ice
    186       REAL :: ini_surf_co2                                         ! Initial surface of sublimating co2 ice
    187 
    188       REAL :: global_ave_press_GCM           ! constant: global average pressure retrieved in the GCM [Pa]
    189       REAL :: global_ave_press_old           ! constant: Global average pressure of initial/previous time step
    190       REAL :: global_ave_press_new           ! constant: Global average pressure of current time step
    191 
    192       REAL , dimension(:,:), allocatable ::  zplev_new ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
    193       REAL , dimension(:,:), allocatable ::  zplev_gcm ! same but retrieved from the gcm [kg/m^2]
    194       REAL , dimension(:,:,:), allocatable ::  zplev_new_timeseries  ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
    195       REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries   ! same but with the time series, for oldest time step
    196 
    197       LOGICAL :: STOPPING_water   ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
    198       LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
    199       LOGICAL :: STOPPING_co2     ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
    200       LOGICAL :: STOPPING_pressure! Logical : is the criterion (% of change in the surface pressure) reached?
    201       INTEGER :: criterion_stop   ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
    202 
    203       real,save ::     A , B, mmean              ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
    204       real ,allocatable :: vmr_co2_gcm(:,:)      ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
    205       real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times  co2 volume mixing ratio used in the PEM
    206       real ,allocatable :: q_co2_PEM_phys(:,:)   ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]
    207       REAL ,allocatable :: q_h2o_PEM_phys(:,:)   ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
    208       integer :: timelen                         ! # time samples
    209       REAL :: ave                                ! intermediate varibale to compute average
    210 
    211       REAL, ALLOCATABLE :: p(:,:)  ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
    212       REAL :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
    213 
    214 
    215 !!!!!!!!!!!!!!!!!!!!!!!! SLOPE
    216       REAL, dimension(:,:),  allocatable :: min_co2_ice_1           ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
    217       REAL, dimension(:,:),  allocatable :: min_co2_ice_2           ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
    218       REAL, dimension(:,:),  allocatable :: min_h2o_ice_1           ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
    219       REAL, dimension(:,:),  allocatable :: min_h2o_ice_2           ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
    220       REAL, dimension(:,:,:),allocatable :: co2_ice_GCM             ! Physics x NSLOPE x Times field : co2 ice given by the GCM  [kg/m^2]
    221       REAL, dimension(:,:),  allocatable :: initial_co2_ice_sublim  ! physical point field : Logical array indicating sublimating point of co2 ice
    222       REAL, dimension(:,:),  allocatable :: initial_h2o_ice         ! physical point field : Logical array indicating if there is water ice at initial state
    223       REAL, dimension(:,:),  allocatable :: initial_co2_ice         ! physical point field : Logical array indicating if there is co2 ice at initial state
    224       REAL, dimension(:,:),  allocatable :: tendencies_co2_ice      ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
    225       REAL, dimension(:,:),  allocatable :: tendencies_co2_ice_ini  ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
    226       REAL, dimension(:,:),  allocatable :: tendencies_h2o_ice      ! physical pointx slope  field : Tendency of evolution of perenial h2o ice
    227       REAL, dimension(:,:),  allocatable :: flag_co2flow(:,:)       !(ngrid,nslope): Flag where there is a CO2 glacier flow
    228       REAL, dimension(:),    allocatable :: flag_co2flow_mesh(:)    !(ngrid)       : Flag where there is a CO2 glacier flow
    229       REAL, dimension(:,:),  allocatable :: flag_h2oflow(:,:)       !(ngrid,nslope): Flag where there is a H2O glacier flow
    230       REAL, dimension(:),    allocatable :: flag_h2oflow_mesh(:)    !(ngrid)       : Flag where there is a H2O glacier flow
    231 
    232 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL
    233 
    234      REAL, ALLOCATABLE :: tsurf_ave(:,:)       ! Physic x SLOPE field : Averaged Surface Temperature [K]
    235      REAL, ALLOCATABLE :: tsoil_ave(:,:,:)     ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
    236      REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:)    ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
    237 
    238      REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    239      REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:)      !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
    240      REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:)                 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
    241 
    242      REAL,ALLOCATABLE  :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
    243      REAL,ALLOCATABLE  :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
    244      REAL,ALLOCATABLE  :: Tsurf_locslope(:)   ! Physic x Soil: Intermediate surface temperature  to compute Tsoil [K]
    245      REAL,ALLOCATABLE  :: watersoil_density_timeseries(:,:,:,:)     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
    246      REAL,ALLOCATABLE  :: watersurf_density_ave(:,:)                ! Physic  x Slope, water surface density, yearly averaged [kg/m^3]
    247      REAL,ALLOCATABLE  :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
    248      REAL,ALLOCATABLE  :: watersoil_density_PEM_ave(:,:,:)          ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
    249      REAL,ALLOCATABLE  :: Tsurfave_before_saved(:,:)                ! Surface temperature saved from previous time step [K]
    250      REAL, ALLOCATABLE :: delta_co2_adsorbded(:)                    ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    251      REAL, ALLOCATABLE :: delta_h2o_adsorbded(:)                    ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
    252      REAL :: totmassco2_adsorbded                                   ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
    253      REAL :: totmassh2o_adsorbded                                   ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
    254      LOGICAL :: bool_sublim                                         ! logical to check if there is sublimation or not
    255      
    256 
    257 !! Some parameters for the PEM run
    258     REAL, PARAMETER :: year_step = 1     ! timestep for the pem
    259     INTEGER ::  year_iter                ! number of iteration
    260     INTEGER :: year_iter_max             ! maximum number of iterations before stopping
    261     REAL :: timestep                     ! timestep [s]
    262     REAL :: watercap_sum                 ! total mass of water cap [kg/m^2]
    263     REAL :: water_sum                    ! total mass of water in the mesh [kg/m^2]
    264 
    265 #ifdef CPP_STD
    266      REAL :: frost_albedo_threshold=0.05             ! frost albedo threeshold to convert fresh frost to old ice
    267      REAL :: albedo_h2o_frost                        ! albedo of h2o frost
    268      REAL,allocatable :: tsurf_read_generic(:)       ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    269      REAL,allocatable :: qsurf_read_generic(:,:)     ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    270      REAL,allocatable :: tsoil_read_generic(:,:)     ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    271      REAL,allocatable :: emis_read_generic(:)        ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    272      REAL,allocatable :: albedo_read_generic(:,:)    ! Temporary variable to do the subslope transfert dimensiion when reading form generic
    273      REAL,allocatable :: tsurf(:,:)                  ! Subslope variable, only needed in the GENERIC case
    274      REAL,allocatable :: qsurf(:,:,:)                ! Subslope variable, only needed in the GENERIC case
    275      REAL,allocatable :: tsoil(:,:,:)                ! Subslope variable, only needed in the GENERIC case
    276      REAL,allocatable :: emis(:,:)                   ! Subslope variable, only needed in the GENERIC case
    277      REAL,allocatable :: watercap(:,:)               ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
    278      LOGICAL, allocatable :: WATERCAPTAG(:)          ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
    279      REAL,allocatable :: albedo(:,:,:)               ! Subslope variable, only needed in the GENERIC case
    280      REAL,allocatable :: inertiesoil(:,:,:)          ! Subslope variable, only needed in the GENERIC case
    281 #endif
    282 
    283 #ifdef CPP_1D
    284   INTEGER :: nsplit_phys
    285   integer,parameter:: jjm_value=jjm-1
    286   integer :: ierr
    287 #else
    288   integer,parameter:: jjm_value=jjm
    289 #endif
    290 
    291 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
    292 
    293 ! Loop variable
    294      INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil
    295 #ifndef CPP_STD
    296 ! Parallel variables
    297       is_sequential=.true.
    298       is_parallel=.false.
    299       is_mpi_root=.true.
    300       is_omp_root=.true.
    301       is_master=.true.
    302 #endif
    303 
    304       day_ini=0    !test
    305       time_phys=0. !test
    306 
    307 ! Some constants
    308 
    309       ngrid=ngridmx
    310       nlayer=llm
    311    
    312       A =(1/m_co2 - 1/m_noco2)
    313       B=1/m_noco2
    314 
    315       year_day=669
    316       daysec=88775.
    317       timestep=year_day*daysec/year_step
    318 
    319 !------------------------
    320 
    321 ! I   Initialisation
    322 !    I_a READ run.def
    323 
    324 !------------------------
    325 
    326 !----------------------------READ run.def ---------------------
    327 
    328 #ifndef CPP_1D
    329       dtphys=0
    330       CALL conf_gcm( 99, .TRUE. )
    331       call infotrac_init
    332       nq=nqtot
    333      allocate(q(ip1jmp1,llm,nqtot))
    334 #else
    335       ! load tracer names from file 'traceur.def'
    336         open(90,file='traceur.def',status='old',form='formatted',&
    337             iostat=ierr)
    338         if (ierr.ne.0) then
    339           write(*,*) 'Cannot find required file "traceur.def"'
    340           write(*,*) ' If you want to run with tracers, I need it'
    341           write(*,*) ' ... might as well stop here ...'
    342           stop
    343         else
    344           write(*,*) "pem1d: Reading file traceur.def"
    345           ! read number of tracers:
    346           read(90,*,iostat=ierr) nq
    347           print *, "nq",nq
    348           nqtot=nq ! set value of nqtot (in infotrac module) as nq
    349           if (ierr.ne.0) then
     283    ! load tracer names from file 'traceur.def'
     284    open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
     285    if (ierr /= 0) then
     286        write(*,*) 'Cannot find required file "traceur.def"'
     287        write(*,*) ' If you want to run with tracers, I need it'
     288        write(*,*) ' ... might as well stop here ...'
     289        stop
     290    else
     291        write(*,*) "pem1d: Reading file traceur.def"
     292        ! read number of tracers:
     293        read(90,*,iostat = ierr) nq
     294        write(*,*) "nq",nq
     295        nqtot = nq ! set value of nqtot (in infotrac module) as nq
     296        if (ierr /= 0) then
    350297            write(*,*) "pem1d: error reading number of tracers"
    351298            write(*,*) "   (first line of traceur.def) "
    352299            stop
    353           endif
    354           if (nq<1) then
     300        endif
     301        if (nq < 1) then
    355302            write(*,*) "pem1d: error number of tracers"
    356303            write(*,*) "is nq=",nq," but must be >=1!"
    357304            stop
    358           endif
    359305        endif
    360      nq=nqtot
    361      allocate(q(ip1jmp1,llm,nqtot))
    362      allocate(ap(nlayer+1))
    363      allocate(bp(nlayer+1))
    364      call init_phys_1d(llm,nqtot,vcov,ucov, &
    365                     teta,q,ps, time_0,ap,bp)   
    366      pi=2.E+0*asin(1.E+0)
    367      g=3.72
    368      nsplit_phys=1
     306    endif
     307    nq = nqtot
     308    allocate(q(ip1jmp1,llm,nqtot))
     309    allocate(ap(nlayer + 1))
     310    allocate(bp(nlayer + 1))
     311    call init_phys_1d(llm,nqtot,vcov,ucov,teta,q,ps,time_0,ap,bp)
     312    pi = 2.*asin(1.)
     313    g = 3.72
     314    nsplit_phys = 1
    369315#endif
    370       CALL conf_pem     
    371  
    372 !------------------------
    373 
    374 ! I   Initialisation
    375 !    I_a READ run.def
     316
     317call conf_pem
     318
     319!------------------------
     320! I   Initialization
    376321!    I_b READ of start_evol.nc and starfi_evol.nc
    377 
    378 !------------------------
    379 
    380 !----------------------------Initialisation : READ some constant of startfi_evol.nc ---------------------
    381 
    382 !----------------------------READ start.nc ---------------------
    383 
    384 
     322!------------------------
     323! I_b.1 READ start_evol.nc
     324allocate(ps_start_GCM(ngrid))
    385325#ifndef CPP_1D
    386      CALL dynetat0(FILE_NAME_start,vcov,ucov, &
    387                     teta,q,masse,ps,phis, time_0)
    388 #endif
    389 
    390      allocate(ps_start_GCM(ngrid))
    391 #ifndef CPP_1D
    392      call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
    393 
    394      CALL iniconst !new
    395      CALL inigeom
    396 
    397      allocate(ap(nlayer+1))
    398      allocate(bp(nlayer+1))
    399      status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
    400      status = nf90_inq_varid(ncid, "ap", apvarid)
    401      status = nf90_get_var(ncid, apvarid, ap)
    402      status = nf90_inq_varid(ncid, "bp", bpvarid)
    403      status = nf90_get_var(ncid, bpvarid, bp)
    404      status =nf90_close(ncid)
    405 
     326    call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0)
     327
     328    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
     329
     330    call iniconst !new
     331    call inigeom
     332
     333    allocate(ap(nlayer + 1))
     334    allocate(bp(nlayer + 1))
     335    status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid)
     336    status = nf90_inq_varid(ncid,"ap",apvarid)
     337    status = nf90_get_var(ncid,apvarid,ap)
     338    status = nf90_inq_varid(ncid,"bp",bpvarid)
     339    status = nf90_get_var(ncid,bpvarid,bp)
     340    status = nf90_close(ncid)
     341
     342    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, &
     343                   rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
    406344#else
    407 
    408      ps_start_GCM(1)=ps(1)
    409 
    410 #endif
    411 
    412 
    413 #ifndef CPP_1D
    414      CALL iniphysiq(iim,jjm,llm, &
    415           (jjm-1)*iim+2,comm_lmdz, &
    416           daysec,day_ini,dtphys/nsplit_phys, &
    417           rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
    418           iflag_phys)
     345    ps_start_GCM(1) = ps(1)
    419346#endif
    420347
    421348! In the gcm, these values are given to the physic by the dynamic.
    422349! Here we simply read them in the startfi_evol.nc file
    423       status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
    424 
    425       allocate(longitude(ngrid))
    426       allocate(latitude(ngrid))
    427       allocate(cell_area(ngrid))
    428 
    429       status = nf90_inq_varid(ncid, "longitude", lonvarid)
    430       status = nf90_get_var(ncid, lonvarid, longitude)
    431 
    432       status = nf90_inq_varid(ncid, "latitude", latvarid)
    433       status = nf90_get_var(ncid, latvarid, latitude)
    434 
    435       status = nf90_inq_varid(ncid, "area", areavarid)
    436       status = nf90_get_var(ncid, areavarid, cell_area)
    437 
    438       status = nf90_inq_varid(ncid, "soildepth", sdvarid)
    439       status = nf90_get_var(ncid, sdvarid, mlayer)
    440 
    441       status =nf90_close(ncid)
    442 
    443 !----------------------------READ startfi.nc ---------------------
    444 
     350status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
     351
     352allocate(longitude(ngrid))
     353allocate(latitude(ngrid))
     354allocate(cell_area(ngrid))
     355
     356status = nf90_inq_varid(ncid,"longitude",lonvarid)
     357status = nf90_get_var(ncid,lonvarid,longitude)
     358
     359status = nf90_inq_varid(ncid,"latitude",latvarid)
     360status = nf90_get_var(ncid,latvarid,latitude)
     361
     362status = nf90_inq_varid(ncid,"area",areavarid)
     363status = nf90_get_var(ncid,areavarid,cell_area)
     364
     365status = nf90_inq_varid(ncid,"soildepth",sdvarid)
     366status = nf90_get_var(ncid,sdvarid,mlayer)
     367
     368status = nf90_close(ncid)
     369
     370! I_b.2 READ startfi_evol.nc
    445371! First we read the initial state (starfi.nc)
    446 
    447372#ifndef CPP_STD
    448       CALL phyetat0 (FILE_NAME,0,0, &
    449               nsoilmx,ngrid,nlayer,nq,   &
    450               day_ini,time_phys,         &
    451               tsurf,tsoil,albedo,emis,   &
    452               q2,qsurf,tauscaling,totcloudfrac,wstar,     &
    453               watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
     373    call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,day_ini,time_phys,tsurf, &
     374                  tsoil,albedo,emis,q2,qsurf,tauscaling,totcloudfrac,wstar,      &
     375                  watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
    454376
    455377! Remove unphysical values of surface tracer
    456      DO i=1,ngrid
    457        DO nnq=1,nqtot
    458          DO islope=1,nslope
    459            if(qsurf(i,nnq,islope).LT.0) then
    460              qsurf(i,nnq,islope)=0.
    461            endif
    462          enddo
    463        enddo
    464      enddo
    465    
    466      call surfini(ngrid,qsurf)
    467 
     378    do i = 1,ngrid
     379        do nnq = 1,nqtot
     380            do islope = 1,nslope
     381                if (qsurf(i,nnq,islope) < 0) qsurf(i,nnq,islope) = 0.
     382            enddo
     383        enddo
     384    enddo
     385
     386    call surfini(ngrid,qsurf)
    468387#else
    469         call phys_state_var_init(nq)
    470          IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
    471          call initracer(ngrid,nq)
    472          call iniaerosol()
    473          allocate(tsurf_read_generic(ngrid))
    474          allocate(qsurf_read_generic(ngrid,nq))
    475          allocate(tsoil_read_generic(ngrid,nsoilmx))
    476          allocate(emis_read_generic(ngrid))
    477          allocate(tsurf(ngrid,1))
    478          allocate(qsurf(ngrid,nq,1))
    479          allocate(tsoil(ngrid,nsoilmx,1))
    480          allocate(emis(ngrid,1))
    481          allocate(watercap(ngrid,1))
    482          allocate(watercaptag(ngrid))
    483          allocate(albedo_read_generic(ngrid,2))
    484          allocate(albedo(ngrid,2,1))
    485          allocate(inertiesoil(ngrid,nsoilmx,1))
    486          call phyetat0(.true.,                                 &
    487                        ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,      &
    488                        day_ini,time_phys,tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,qsurf_read_generic,   &
    489                        cloudfrac,totcloudfrac,hice,                   &
    490                        rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
    491          call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
    492 
    493          nslope=1
    494          call ini_comslope_h(ngrid,1)
    495 
    496          qsurf(:,:,1)=qsurf_read_generic(:,:)
    497          tsurf(:,1)=tsurf_read_generic(:)
    498          tsoil(:,:,1)=tsoil_read_generic(:,:)
    499          emis(:,1)=emis_read_generic(:)
    500          watercap(:,1)=0.
    501          watercaptag(:)=.false.
    502          albedo(:,1,1)=albedo_read_generic(:,1)
    503          albedo(:,2,1)=albedo_read_generic(:,2)
    504          inertiesoil(:,:,1)=inertiedat(:,:)
    505 
    506 if (nslope.eq.1) then
    507  def_slope(1) = 0
    508  def_slope(2) = 0
    509  def_slope_mean=0
    510  subslope_dist(:,1) = 1.
    511 endif
     388    call phys_state_var_init(nq)
     389    if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
     390    call initracer(ngrid,nq)
     391    call iniaerosol()
     392    allocate(tsurf_read_generic(ngrid))
     393    allocate(qsurf_read_generic(ngrid,nq))
     394    allocate(tsoil_read_generic(ngrid,nsoilmx))
     395    allocate(emis_read_generic(ngrid))
     396    allocate(tsurf(ngrid,1))
     397    allocate(qsurf(ngrid,nq,1))
     398    allocate(tsoil(ngrid,nsoilmx,1))
     399    allocate(emis(ngrid,1))
     400    allocate(watercap(ngrid,1))
     401    allocate(watercaptag(ngrid))
     402    allocate(albedo_read_generic(ngrid,2))
     403    allocate(albedo(ngrid,2,1))
     404    allocate(inertiesoil(ngrid,nsoilmx,1))
     405    call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,day_ini,time_phys, &
     406                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,     &
     407                  qsurf_read_generic,cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, &
     408                  tslab, tsea_ice,sea_ice)
     409    call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
     410
     411    nslope = 1
     412    call ini_comslope_h(ngrid,1)
     413
     414    qsurf(:,:,1) = qsurf_read_generic(:,:)
     415    tsurf(:,1) = tsurf_read_generic(:)
     416    tsoil(:,:,1) = tsoil_read_generic(:,:)
     417    emis(:,1) = emis_read_generic(:)
     418    watercap(:,1) = 0.
     419    watercaptag(:) = .false.
     420    albedo(:,1,1) = albedo_read_generic(:,1)
     421    albedo(:,2,1) = albedo_read_generic(:,2)
     422    inertiesoil(:,:,1) = inertiedat(:,:)
     423
     424    if (nslope == 1) then
     425        def_slope(1) = 0
     426        def_slope(2) = 0
     427        def_slope_mean = 0
     428        subslope_dist(:,1) = 1.
     429    endif
    512430
    513431! Remove unphysical values of surface tracer
    514      DO i=1,ngrid
    515        DO nnq=1,nqtot
    516          qsurf(i,nnq,1)=qsurf_read_generic(i,nnq)
    517          if(qsurf(i,nnq,1).LT.0) then
    518            qsurf(i,nnq,1)=0.
    519          endif
    520        enddo
    521      enddo
     432    do i = 1,ngrid
     433        do nnq = 1,nqtot
     434            qsurf(i,nnq,1) = qsurf_read_generic(i,nnq)
     435            if (qsurf(i,nnq,1) < 0) qsurf(i,nnq,1) = 0.
     436        enddo
     437    enddo
    522438#endif
    523439
    524      DO nnq=1,nqtot  ! Why not using ini_tracer ?
    525        if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
    526        if(noms(nnq).eq."h2o_vap") then
    527           igcm_h2o_vap = nnq
    528           mmol(igcm_h2o_vap)=18.
    529        endif
    530        if(noms(nnq).eq."co2") igcm_co2 = nnq
    531      ENDDO
    532      r= 8.314511E+0 *1000.E+0/mugaz
    533 !------------------------
    534 
    535 ! I   Initialisation
    536 !    I_a READ run.def
    537 !    I_b READ of start_evol.nc and starfi_evol.nc
     440do nnq = 1,nqtot  ! Why not using ini_tracer ?
     441    if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
     442    if (noms(nnq) == "h2o_vap") then
     443        igcm_h2o_vap = nnq
     444        mmol(igcm_h2o_vap)=18.
     445    endif
     446    if (noms(nnq) == "co2") igcm_co2 = nnq
     447enddo
     448r = 8.314511E+0*1000.E+0/mugaz
     449
     450!------------------------
     451! I   Initialization
    538452!    I_c Subslope parametrisation
    539 
    540 !------------------------
    541 
    542 !----------------------------Subslope parametrisation definition ---------------------
    543 
    544 !     Define some slope statistics
    545      iflat=1
    546      DO islope=2,nslope
    547        IF(abs(def_slope_mean(islope)).lt. &
    548          abs(def_slope_mean(iflat))) THEN
    549          iflat = islope
    550        ENDIF     
    551      ENDDO
    552 
    553      PRINT*,'Flat slope for islope = ',iflat
    554      PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
    555 
    556      allocate(flag_co2flow(ngrid,nslope))
    557      allocate(flag_co2flow_mesh(ngrid))
    558      allocate(flag_h2oflow(ngrid,nslope))
    559      allocate(flag_h2oflow_mesh(ngrid))
    560 
    561      flag_co2flow(:,:) = 0     
    562      flag_co2flow_mesh(:) = 0
    563      flag_h2oflow(:,:) = 0     
    564      flag_h2oflow_mesh(:) = 0
    565 !---------------------------- READ GCM data ---------------------
    566 
    567 ! I   Initialisation
    568 !    I_a READ run.def
    569 !    I_b READ of start_evol.nc and starfi_evol.nc
    570 !    I_c Subslope parametrisation
    571 !    I_d READ GCM data
    572 
    573 !------------------------
    574 
     453!------------------------
     454! Define some slope statistics
     455iflat = 1
     456do islope = 2,nslope
     457    if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
     458enddo
     459
     460write(*,*) 'Flat slope for islope = ',iflat
     461write(*,*) 'corresponding criterium = ',def_slope_mean(iflat)
     462
     463allocate(flag_co2flow(ngrid,nslope))
     464allocate(flag_co2flow_mesh(ngrid))
     465allocate(flag_h2oflow(ngrid,nslope))
     466allocate(flag_h2oflow_mesh(ngrid))
     467
     468flag_co2flow(:,:) = 0     
     469flag_co2flow_mesh(:) = 0
     470flag_h2oflow(:,:) = 0     
     471flag_h2oflow_mesh(:) = 0
     472
     473!------------------------
     474! I   Initialization
     475!    I_d READ GCM data and convert to the physical grid
     476!------------------------
    575477! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value
    576 
    577      call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
    578 
    579      allocate(tsoil_ave(ngrid,nsoilmx,nslope))
    580      allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
    581 
    582      allocate(vmr_co2_gcm(ngrid,timelen))
    583      allocate(ps_timeseries(ngrid,timelen))
    584      allocate(min_co2_ice_1(ngrid,nslope))
    585      allocate(min_h2o_ice_1(ngrid,nslope))
    586      allocate(min_co2_ice_2(ngrid,nslope))
    587      allocate(min_h2o_ice_2(ngrid,nslope))
    588      allocate(tsurf_ave_yr1(ngrid,nslope))
    589      allocate(tsurf_ave(ngrid,nslope))
    590 
    591      allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
    592      allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
    593      allocate(q_co2_PEM_phys(ngrid,timelen))
    594      allocate(q_h2o_PEM_phys(ngrid,timelen))
    595      allocate(co2_ice_GCM(ngrid,nslope,timelen))
    596      allocate(watersurf_density_ave(ngrid,nslope))
    597      allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
    598 
    599      allocate(Tsurfave_before_saved(ngrid,nslope))
    600      allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    601      allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    602      allocate(delta_co2_adsorbded(ngrid))
    603      allocate(delta_h2o_adsorbded(ngrid))
    604      allocate(vmr_co2_pem_phys(ngrid,timelen))
    605 
    606      print *, "Downloading data Y1..."
    607 
    608      call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1,&   
    609                        tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &     
    610                        watersurf_density_ave,watersoil_density_timeseries)
     478call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
     479
     480allocate(tsoil_ave(ngrid,nsoilmx,nslope))
     481allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
     482allocate(vmr_co2_gcm(ngrid,timelen))
     483allocate(ps_timeseries(ngrid,timelen))
     484allocate(min_co2_ice_1(ngrid,nslope))
     485allocate(min_h2o_ice_1(ngrid,nslope))
     486allocate(min_co2_ice_2(ngrid,nslope))
     487allocate(min_h2o_ice_2(ngrid,nslope))
     488allocate(tsurf_ave_yr1(ngrid,nslope))
     489allocate(tsurf_ave(ngrid,nslope))
     490allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
     491allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
     492allocate(q_co2_PEM_phys(ngrid,timelen))
     493allocate(q_h2o_PEM_phys(ngrid,timelen))
     494allocate(co2_ice_GCM(ngrid,nslope,timelen))
     495allocate(watersurf_density_ave(ngrid,nslope))
     496allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
     497allocate(Tsurfave_before_saved(ngrid,nslope))
     498allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
     499allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
     500allocate(delta_co2_adsorbded(ngrid))
     501allocate(delta_h2o_adsorbded(ngrid))
     502allocate(vmr_co2_pem_phys(ngrid,timelen))
     503
     504write(*,*) "Downloading data Y1..."
     505call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
     506                   tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,           &
     507                   co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
     508write(*,*) "Downloading data Y1 done"
    611509
    612510! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
    613 
    614      print *, "Downloading data Y1 done"
    615      print *, "Downloading data Y2"
    616 
    617      call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
    618                   tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM, &     
    619                   watersurf_density_ave,watersoil_density_timeseries)
    620 
    621      print *, "Downloading data Y2 done"
    622 
    623 !------------------------
    624 
    625 ! I   Initialisation
    626 !    I_a READ run.def
    627 !    I_b READ of start_evol.nc and starfi_evol.nc
    628 !    I_c Subslope parametrisation
    629 !    I_d READ GCM data and convert to the physical grid
    630 !    I_e Initialisation of the PEM variable and soil
    631 
    632 !------------------------
    633 
    634 !---------------------------- Initialisation of the PEM soil and values ---------------------
    635 
    636    call end_comsoil_h_PEM
    637    call ini_comsoil_h_PEM(ngrid,nslope)
    638    call end_adsorption_h_PEM
    639    call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
    640    call end_ice_table_porefilling
    641    call ini_ice_table_porefilling(ngrid,nslope)
    642 
    643    if(soil_pem) then
    644       call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
    645       DO l=1,nsoilmx
     511write(*,*) "Downloading data Y2"
     512call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
     513                   tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,              &
     514                   co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
     515write(*,*) "Downloading data Y2 done"
     516
     517!------------------------
     518! I   Initialization
     519!    I_e Initialization of the PEM variables and soil
     520!------------------------
     521call end_comsoil_h_PEM
     522call ini_comsoil_h_PEM(ngrid,nslope)
     523call end_adsorption_h_PEM
     524call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
     525call end_ice_table_porefilling
     526call ini_ice_table_porefilling(ngrid,nslope)
     527
     528if (soil_pem) then
     529    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
     530    do l=1,nsoilmx
    646531        tsoil_PEM(:,l,:)=tsoil_ave(:,l,:)
    647532        tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:)
    648533        watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:)
    649       ENDDO
    650       DO l=nsoilmx+1,nsoilmx_PEM
     534    enddo
     535    do l=nsoilmx+1,nsoilmx_PEM
    651536        tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:)
    652537        watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:)
    653       ENDDO
    654       watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
    655   endif !soil_pem
    656       deallocate(tsoil_ave)
    657       deallocate(tsoil_GCM_timeseries)
    658 
    659 !------------------------
    660 
    661 ! I   Initialisation
    662 !    I_a READ run.def
    663 !    I_b READ of start_evol.nc and starfi_evol.nc
    664 !    I_c Subslope parametrisation
    665 !    I_d READ GCM data and convert to the physical grid
    666 !    I_e Initialisation of the PEM variable and soil
     538    enddo
     539    watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
     540endif !soil_pem
     541deallocate(tsoil_ave,tsoil_GCM_timeseries)
     542
     543!------------------------
     544! I   Initialization
    667545!    I_f Compute tendencies & Save initial situation
    668 
    669 !----- Compute tendencies from the PCM run
    670 
    671      allocate(tendencies_co2_ice(ngrid,nslope))
    672      allocate(tendencies_co2_ice_ini(ngrid,nslope))
    673      allocate(tendencies_h2o_ice(ngrid,nslope))
    674 
    675 !  Compute the tendencies of the evolution of ice over the years
    676 
    677       call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,&
    678              min_co2_ice_2,tendencies_co2_ice)
    679 
    680       tendencies_co2_ice_ini(:,:)=tendencies_co2_ice(:,:)
    681 
    682       call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,&
    683              min_h2o_ice_2,tendencies_h2o_ice)
    684 
    685      deallocate(min_co2_ice_1)
    686      deallocate(min_co2_ice_2)
    687      deallocate(min_h2o_ice_1)
    688      deallocate(min_h2o_ice_2)
    689 
    690 !------------------------
    691 
    692 ! I   Initialisation
    693 !    I_a READ run.def
    694 !    I_b READ of start_evol.nc and starfi_evol.nc
    695 !    I_c Subslope parametrisation
    696 !    I_d READ GCM data and convert to the physical grid
    697 !    I_e Initialisation of the PEM variable and soil
    698 !    I_f Compute tendencies & Save initial situation
     546!------------------------
     547allocate(tendencies_co2_ice(ngrid,nslope))
     548allocate(tendencies_co2_ice_ini(ngrid,nslope))
     549allocate(tendencies_h2o_ice(ngrid,nslope))
     550
     551! Compute the tendencies of the evolution of ice over the years
     552call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice)
     553tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:)
     554call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice)
     555
     556deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2)
     557
     558!------------------------
     559! I   Initialization
    699560!    I_g Save initial PCM situation
    700 
    701 !---------------------------- Save initial PCM situation ---------------------
    702 
    703      allocate(initial_co2_ice_sublim(ngrid,nslope))
    704      allocate(initial_co2_ice(ngrid,nslope))
    705      allocate(initial_h2o_ice(ngrid,nslope))
     561!------------------------
     562allocate(initial_co2_ice_sublim(ngrid,nslope))
     563allocate(initial_co2_ice(ngrid,nslope))
     564allocate(initial_h2o_ice(ngrid,nslope))
    706565
    707566! We save the places where water ice is sublimating
    708567! We compute the surface of water ice sublimating
    709   ini_surf_co2=0.
    710   ini_surf_h2o=0.
    711   Total_surface=0.
    712   do i=1,ngrid
    713     Total_surface=Total_surface+cell_area(i)
    714     do islope=1,nslope
    715       if (tendencies_co2_ice(i,islope).LT.0) then
    716          initial_co2_ice_sublim(i,islope)=1.
    717          ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
    718       else
    719          initial_co2_ice_sublim(i,islope)=0.         
    720       endif
    721       if (qsurf(i,igcm_co2,islope).GT.0) then
    722          initial_co2_ice(i,islope)=1.
    723       else
    724          initial_co2_ice(i,islope)=0.         
    725       endif
    726       if (tendencies_h2o_ice(i,islope).LT.0) then
    727          initial_h2o_ice(i,islope)=1.
    728          ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
    729       else
    730          initial_h2o_ice(i,islope)=0.         
    731       endif
    732     enddo
    733   enddo
    734 
    735      print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
    736      print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
    737      print *, "Total surface of the planet=", Total_surface
    738    
    739      allocate(zplev_gcm(ngrid,nlayer+1))
    740 
    741      DO l=1,nlayer+1
    742        DO ig=1,ngrid
    743          zplev_gcm(ig,l) = ap(l)  + bp(l)*ps_start_GCM(ig)
    744        ENDDO
    745      ENDDO
    746 
    747      global_ave_press_old=0.
    748      do i=1,ngrid
    749        global_ave_press_old=global_ave_press_old+cell_area(i)*ps_start_GCM(i)/Total_surface
    750      enddo
    751 
    752      global_ave_press_GCM=global_ave_press_old
    753      global_ave_press_new=global_ave_press_old
    754      print *, "Initial global average pressure=", global_ave_press_GCM
    755 
    756 !------------------------
    757 
    758 ! I   Initialisation
    759 !    I_a READ run.def
    760 !    I_b READ of start_evol.nc and starfi_evol.nc
    761 !    I_c Subslope parametrisation
    762 !    I_d READ GCM data and convert to the physical grid
    763 !    I_e Initialisation of the PEM variable and soil
    764 !    I_f Compute tendencies & Save initial situation
    765 !    I_g Save initial PCM situation
     568ini_surf_co2 = 0.
     569ini_surf_h2o = 0.
     570Total_surface = 0.
     571do i = 1,ngrid
     572    Total_surface = Total_surface+cell_area(i)
     573    do islope = 1,nslope
     574        if (tendencies_co2_ice(i,islope) < 0) then
     575            initial_co2_ice_sublim(i,islope) = 1.
     576            ini_surf_co2 = ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
     577        else
     578            initial_co2_ice_sublim(i,islope) = 0.
     579        endif
     580        if (qsurf(i,igcm_co2,islope) > 0) then
     581            initial_co2_ice(i,islope) = 1.
     582        else
     583            initial_co2_ice(i,islope) = 0.
     584        endif
     585        if (tendencies_h2o_ice(i,islope) < 0) then
     586            initial_h2o_ice(i,islope) = 1.
     587            ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
     588        else
     589            initial_h2o_ice(i,islope) = 0.
     590        endif
     591    enddo
     592enddo
     593
     594write(*,*) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
     595write(*,*) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
     596write(*,*) "Total surface of the planet=", Total_surface
     597allocate(zplev_gcm(ngrid,nlayer + 1))
     598
     599do l = 1,nlayer + 1
     600    do ig = 1,ngrid
     601        zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
     602    enddo
     603enddo
     604
     605global_ave_press_old = 0.
     606do i = 1,ngrid
     607       global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
     608enddo
     609
     610global_ave_press_GCM = global_ave_press_old
     611global_ave_press_new = global_ave_press_old
     612write(*,*) "Initial global average pressure=", global_ave_press_GCM
     613
     614!------------------------
     615! I   Initialization
    766616!    I_h Read the PEMstart
    767 
    768 !---------------------------- Read the PEMstart ---------------------
    769 
    770       call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth,porefillingice_thickness, &
    771       tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,&
    772       tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,&
    773       watersurf_density_ave,watersoil_density_PEM_ave, &
    774       co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
     617!------------------------
     618call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
     619              porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,          &
     620              tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:),                      &
     621              qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave,             &
     622              co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
     623
     624do ig = 1,ngrid
     625    do islope = 1,nslope
     626        qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
     627    enddo
     628enddo
     629
     630if (adsorption_pem) then
     631    totmassco2_adsorbded = 0.
     632    totmassh2o_adsorbded = 0.
     633    do ig = 1,ngrid
     634        do islope = 1, nslope
     635            do l = 1,nsoilmx_PEM - 1
     636                totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     637                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
     638                totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     639                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
     640            enddo
     641        enddo
     642    enddo
     643
     644    write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
     645    write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
     646endif ! adsorption
     647deallocate(tsurf_ave_yr1)
     648
     649!------------------------
     650! I   Initialization
     651!    I_i Compute orbit criterion
     652!------------------------
     653#ifndef CPP_STD
     654     call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
     655#else
     656     call iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
     657#endif
     658
     659if (evol_orbit_pem) then
     660    call orbit_param_criterion(year_iter_max)
     661else
     662    year_iter_max = Max_iter_pem
     663endif
     664!-------------------------- END INITIALIZATION -------------------------
     665
     666!-------------------------------- RUN ----------------------------------
     667!------------------------
     668! II  Run
     669!    II_a Update pressure, ice and tracers
     670!------------------------
     671year_iter = 0
     672do while (year_iter < year_iter_max)
     673
     674! II.a.1. Compute updated global pressure
     675    write(*,*) "Recomputing the new pressure..."
     676
     677    do i = 1,ngrid
     678        do islope = 1,nslope
     679            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
     680        enddo
     681    enddo
     682    write(*,*) 'Global average pressure old time step',global_ave_press_old
     683    call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
     684     
     685    if (adsorption_pem) then
     686        do i = 1,ngrid
     687            global_ave_press_new = global_ave_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
     688        enddo
     689        write(*,*) 'Global average pressure old time step',global_ave_press_old
     690        write(*,*) 'Global average pressure new time step',global_ave_press_new
     691    endif
     692
     693! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
     694    allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
     695    write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
     696    do l = 1,nlayer + 1
     697        do ig = 1,ngrid
     698            zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
     699        enddo
     700    enddo
     701
     702! II.a.3. Surface pressure timeseries
     703    write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..."
     704    do ig = 1,ngrid
     705        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
     706    enddo
     707
     708! II.a.4. New pressure levels timeseries
     709    allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
     710    write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
     711    do l = 1,nlayer + 1
     712        do ig = 1,ngrid
     713            zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
     714        enddo
     715    enddo
     716
     717! II.a.5. Tracers timeseries
     718    write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..."
     719
     720    l = 1
     721    do ig = 1,ngrid
     722        do t = 1,timelen
     723            q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
     724                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
     725            if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30
     726            if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1.
     727        enddo
     728    enddo
     729
     730    do ig = 1,ngrid
     731        do t = 1, timelen
     732            q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
     733                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))  &
     734                                + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))  &
     735                                -  (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
     736                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
     737            if (q_co2_PEM_phys(ig,t) < 0) then
     738                q_co2_PEM_phys(ig,t) = 1.e-30
     739            elseif (q_co2_PEM_phys(ig,t) > 1) then
     740                q_co2_PEM_phys(ig,t) = 1.
     741            endif
     742            mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
     743            vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
     744        enddo
     745    enddo
     746
     747    deallocate(zplev_new_timeseries,zplev_old_timeseries)
     748
     749!------------------------
     750! II  Run
     751!    II_b Evolution of the ice
     752!------------------------
     753    write(*,*) "Evolution of h2o ice"
     754    call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)
     755
     756    write(*,*) "Evolution of co2 ice"
     757    call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
     758
     759    do islope=1, nslope
     760        write(str2(1:2),'(i2.2)') islope
     761        call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
     762        call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
     763        call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
     764        call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
     765    enddo
     766
     767!------------------------
     768! II  Run
     769!    II_c CO2 & H2O glaciers flows
     770!------------------------
     771    write(*,*) "CO2 glacier flows"
     772
     773    if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
     774                                               global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
    775775 
    776       do ig = 1,ngrid
     776    write(*,*) "H2O glacier flows"
     777
     778    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)
     779
     780    do islope = 1,nslope
     781        write(str2(1:2),'(i2.2)') islope
     782        call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
     783        call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
     784        call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
     785    enddo
     786
     787!------------------------
     788! II  Run
     789!    II_d Update surface and soil temperatures
     790!------------------------
     791! II_d.1 Update Tsurf
     792    write(*,*) "Updating the new Tsurf"
     793    bool_sublim = .false.
     794    Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
     795    do ig = 1,ngrid
    777796        do islope = 1,nslope
    778           qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    779         enddo
    780       enddo
    781 
    782     if(adsorption_pem) then
    783      totmassco2_adsorbded = 0.
    784      totmassh2o_adsorbded = 0.
    785      do ig = 1,ngrid
    786       do islope =1, nslope
    787        do l = 1,nsoilmx_PEM - 1
    788           totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
    789           subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
    790           cell_area(ig)
    791           totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
    792           subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
    793           cell_area(ig)
    794        enddo
    795       enddo
    796      enddo
    797 
    798      write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
    799      write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
    800      endif ! adsorption
    801      deallocate(tsurf_ave_yr1)
    802 
    803 !------------------------
    804 
    805 ! I   Initialisation
    806 !    I_a READ run.def
    807 !    I_b READ of start_evol.nc and starfi_evol.nc
    808 !    I_c Subslope parametrisation
    809 !    I_d READ GCM data and convert to the physical grid
    810 !    I_e Initialisation of the PEM variable and soil
    811 !    I_f Compute tendencies & Save initial situation
    812 !    I_g Save initial PCM situation
    813 !    I_h Read the PEMstar
    814 !    I_i Compute orbit criterion
    815 
    816 #ifndef CPP_STD
    817      CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
    818 #else
    819      CALL iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
    820 #endif
    821 
    822      if (evol_orbit_pem) then
    823        call orbit_param_criterion(year_iter_max)
    824      else
    825        year_iter_max = Max_iter_pem
    826      endif
    827 
    828 !--------------------------- END INITIALISATION ---------------------
    829 
    830 !---------------------------- RUN ---------------------
    831 
    832 !------------------------
    833 
     797            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
     798                if (latitude_deg(ig) > 0) then
     799                    do ig_loop = ig,ngrid
     800                        do islope_loop = islope,iflat,-1
     801                            if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
     802                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
     803                                bool_sublim = .true.
     804                                exit
     805                            endif
     806                        enddo
     807                        if (bool_sublim) exit
     808                    enddo
     809                else
     810                    do ig_loop = ig,1,-1
     811                        do islope_loop = islope,iflat
     812                            if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
     813                                tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
     814                                bool_sublim = .true.
     815                                exit
     816                            endif
     817                        enddo
     818                        if (bool_sublim) exit
     819                    enddo
     820                endif
     821                initial_co2_ice(ig,islope) = 0
     822                if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
     823                    albedo(ig,1,islope) = albedo_h2o_frost
     824                    albedo(ig,2,islope) = albedo_h2o_frost
     825                else
     826                    albedo(ig,1,islope) = albedodat(ig)
     827                    albedo(ig,2,islope) = albedodat(ig)
     828                    emis(ig,islope) = emissiv
     829                endif
     830            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
     831                ave = 0.
     832                do t = 1,timelen
     833                    if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
     834                        ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
     835                    else
     836                        ave = ave + tsurf_GCM_timeseries(ig,islope,t)
     837                    endif
     838                enddo
     839                tsurf_ave(ig,islope) = ave/timelen
     840            endif
     841        enddo
     842    enddo
     843
     844    do t = 1,timelen
     845        tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:))
     846    enddo
     847    ! for the start
     848    do ig = 1,ngrid
     849        do islope = 1,nslope
     850            tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
     851        enddo
     852    enddo
     853
     854    do islope = 1,nslope
     855        write(str2(1:2),'(i2.2)') islope
     856        call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
     857    enddo
     858
     859    if (soil_pem) then
     860
     861! II_d.2 Update soil temperature
     862        allocate(TI_locslope(ngrid,nsoilmx_PEM))
     863        allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
     864        allocate(Tsurf_locslope(ngrid))
     865        write(*,*)"Updating soil temperature"
     866
     867        ! Soil averaged
     868        do islope = 1,nslope
     869            TI_locslope(:,:) = TI_PEM(:,:,islope)
     870            do t = 1,timelen
     871                Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
     872                Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
     873                call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
     874                call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
     875                tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
     876                do ig = 1,ngrid
     877                    do isoil = 1,nsoilmx_PEM
     878                        watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r)
     879                        if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1)
     880                    enddo
     881                enddo
     882            enddo
     883        enddo
     884        tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
     885        watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
     886
     887        write(*,*) "Update of soil temperature done"
     888
     889        deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
     890        write(*,*) "Compute ice table"
     891
     892! II_d.3 Update the ice table
     893        call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
     894        write(*,*) "Update soil propreties"
     895
     896! II_d.4 Update the soil thermal properties
     897        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, &
     898                                           porefillingice_depth,porefillingice_thickness,TI_PEM)
     899
     900! II_d.5 Update the mass of the regolith adsorbded
     901        if (adsorption_pem) then
     902            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
     903                                     qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, &
     904                                     q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
     905
     906            totmassco2_adsorbded = 0.
     907            totmassh2o_adsorbded = 0.
     908            do ig = 1,ngrid
     909                do islope =1, nslope
     910                    do l = 1,nsoilmx_PEM - 1
     911                        totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     912                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
     913                        cell_area(ig)
     914                        totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
     915                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
     916                        cell_area(ig)
     917                    enddo
     918                enddo
     919            enddo
     920            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
     921            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
     922        endif
     923    endif !soil_pem
     924
     925!------------------------
    834926! II  Run
    835 !    II_a update pressure,ice and tracers
    836 
    837 !------------------------
    838      year_iter=0
    839      do while (year_iter.LT.year_iter_max)
    840 
    841 ! II.a.1. Compute updated global pressure
    842      print *, "Recomputing the new pressure..."
    843 
    844      do i=1,ngrid
    845        do islope=1,nslope
    846            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
    847       enddo
    848      enddo
    849      print *, 'Global average pressure old time step',global_ave_press_old
    850      call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
    851      
    852      if(adsorption_pem) then
    853       do i=1,ngrid
    854         global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
    855        enddo
    856        print *, 'Global average pressure old time step',global_ave_press_old
    857        print *, 'Global average pressure new time step',global_ave_press_new
    858      endif
    859 
    860 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
    861      allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
    862      print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
    863      DO l=1,nlayer+1
    864        DO ig=1,ngrid
    865          zplev_old_timeseries(ig,l,:) = ap(l)  + bp(l)*ps_timeseries(ig,:)
    866        ENDDO
    867      ENDDO
    868 
    869 ! II.a.3. Surface pressure timeseries
    870      print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
    871      do ig = 1,ngrid
    872        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
    873      enddo
    874 
    875 ! II.a.4. New pressure levels timeseries
    876      allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
    877      print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
    878      do l=1,nlayer+1
    879        do ig=1,ngrid
    880          zplev_new_timeseries(ig,l,:)  = ap(l)  + bp(l)*ps_timeseries(ig,:)
    881        enddo
    882      enddo
    883 
    884 ! II.a.5. Tracers timeseries
    885       print *, "Recomputing of tracer VMR timeseries for the new pressure..."
    886 
    887      l=1
    888      DO ig=1,ngrid
    889        DO t = 1, timelen
    890          q_h2o_PEM_phys(ig,t)=q_h2o_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
    891          if(q_h2o_PEM_phys(ig,t).LT.0) then
    892            q_h2o_PEM_phys(ig,t)=1E-30
    893          endif
    894          if(q_h2o_PEM_phys(ig,t).GT.1) then
    895            q_h2o_PEM_phys(ig,t)=1.
    896          endif
    897        enddo
    898      enddo
    899 
    900      DO ig=1,ngrid
    901        DO t = 1, timelen
    902          q_co2_PEM_phys(ig,t)=q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  &
    903                                 + (   (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))  -       &
    904                                       (zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t))     )  / &
    905                                       (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
    906          if (q_co2_PEM_phys(ig,t).LT.0) then
    907               q_co2_PEM_phys(ig,t)=1E-30
    908          elseif (q_co2_PEM_phys(ig,t).GT.1) then
    909               q_co2_PEM_phys(ig,t)=1.
    910          endif
    911          mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
    912          vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
    913        ENDDO
    914      ENDDO
    915    
    916      deallocate(zplev_new_timeseries)
    917      deallocate(zplev_old_timeseries)
    918 
     927!    II_e Update the tendencies
     928!------------------------
     929    write(*,*) "Adaptation of the new co2 tendencies to the current pressure"
     930    call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
     931                               global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
     932
     933!------------------------
    919934! II  Run
    920 !    II_a update pressure, ice and tracers
    921 !    II_b Evolution of the ice
    922 
    923 ! II.b. Evolution of the ice
    924      print *, "Evolution of h2o ice"
    925      call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)
    926      
    927      print *, "Evolution of co2 ice"
    928      call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
    929 
    930      DO islope=1, nslope
    931        write(str2(1:2),'(i2.2)') islope
    932        call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
    933        call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
    934        call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
    935        call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
    936      ENDDO
    937 
    938 !------------------------
    939 
    940 ! II  Run
    941 !    II_a update pressure, ice and tracers
    942 !    II_b Evolution of the ice
    943 !    II_c CO2  & H2O glaciers flows
    944 
    945 !------------------------
    946 
    947       print *, "CO2 glacier flows"
    948 
    949       if (co2glaciersflow) then
    950        call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,&
    951                          global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
    952       endif
    953 
    954       print *, "H2O glacier flows"
    955 
    956       if (h2oglaciersflow) then
    957         call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
    958       endif
    959      DO islope=1, nslope
    960        write(str2(1:2),'(i2.2)') islope
    961        call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
    962        call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
    963        call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
    964      ENDDO
    965 
    966 !------------------------
    967 
    968 ! II  Run
    969 !    II_a update pressure, ice and tracers
    970 !    II_b Evolution of the ice
    971 !    II_c CO2 glaciers flows
    972 !    II_d Update surface and soil temperatures
    973 
    974 !------------------------
    975      
    976 ! II_d.1 Update Tsurf
    977 
    978       print *, "Updating the new Tsurf"
    979       bool_sublim=.false.
    980       Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
    981       DO ig = 1,ngrid
    982         DO islope = 1,nslope
    983           if(initial_co2_ice(ig,islope).gt.0.5 .and. qsurf(ig,igcm_co2,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice
    984             if(latitude_deg(ig).gt.0) then
    985               do ig_loop=ig,ngrid
    986                 DO islope_loop=islope,iflat,-1
    987                   if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then
    988                     tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
    989                     bool_sublim=.true.
    990                     exit
    991                   endif
    992                 enddo
    993                 if (bool_sublim.eqv. .true.) then
    994                   exit
    995                 endif
    996               enddo
    997             else
    998               do ig_loop=ig,1,-1
    999                 DO islope_loop=islope,iflat
    1000                   if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then
    1001                     tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
    1002                     bool_sublim=.true.
    1003                     exit
    1004                   endif
    1005                 enddo
    1006                 if (bool_sublim.eqv. .true.) then
    1007                   exit
    1008                 endif
    1009               enddo
    1010             endif
    1011             initial_co2_ice(ig,islope)=0
    1012             if ((qsurf(ig,igcm_co2,islope).lt.1e-10).and. (qsurf(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then   
    1013               albedo(ig,1,islope) = albedo_h2o_frost
    1014               albedo(ig,2,islope) = albedo_h2o_frost
    1015             else
    1016               albedo(ig,1,islope) = albedodat(ig)
    1017               albedo(ig,2,islope) = albedodat(ig)
    1018               emis(ig,islope) = emissiv         
    1019             endif
    1020           elseif( (qsurf(ig,igcm_co2,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
    1021             ave=0.
    1022             do t=1,timelen
    1023               if(co2_ice_GCM(ig,islope,t).gt.1e-3) then
    1024                 ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
    1025               else
    1026                 ave = ave + tsurf_GCM_timeseries(ig,islope,t)
    1027               endif
    1028             enddo
    1029             tsurf_ave(ig,islope)=ave/timelen
    1030           endif
    1031         enddo
    1032       enddo
    1033 
    1034       do t = 1,timelen
    1035         tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) +( tsurf_ave(:,:) -Tsurfave_before_saved(:,:))
    1036       enddo
    1037       ! for the start
    1038       do ig = 1,ngrid
    1039         do islope = 1,nslope
    1040           tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope))
    1041         enddo
    1042       enddo
    1043 
    1044      DO islope=1, nslope
    1045        write(str2(1:2),'(i2.2)') islope
    1046        call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
    1047      ENDDO
    1048 
    1049     if(soil_pem) then
    1050 
    1051 ! II_d.2 Update soil temperature
    1052 
    1053   allocate(TI_locslope(ngrid,nsoilmx_PEM))
    1054   allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
    1055   allocate(Tsurf_locslope(ngrid))
    1056   print *,"Updating soil temperature"
    1057 
    1058 ! Soil averaged
    1059   do islope = 1,nslope
    1060      TI_locslope(:,:) = TI_PEM(:,:,islope)
    1061      do t = 1,timelen
    1062        Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
    1063        Tsurf_locslope(:) =  tsurf_GCM_timeseries(:,islope,t)
    1064        call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    1065        call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
    1066        tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
    1067        do ig = 1,ngrid
    1068          do isoil = 1,nsoilmx_PEM
    1069           watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r)
    1070           if(isnan(Tsoil_locslope(ig,isoil))) then
    1071             call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1)
    1072           endif
    1073          enddo
    1074        enddo
    1075      enddo
    1076   enddo
    1077      tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
    1078      watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
    1079 
    1080   print *, "Update of soil temperature done"
    1081 
    1082   deallocate(TI_locslope)
    1083   deallocate(Tsoil_locslope)
    1084   deallocate(Tsurf_locslope)
    1085   write(*,*) "Compute ice table"
    1086 
    1087 ! II_d.3 Update the ice table
    1088      call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
    1089      print *, "Update soil propreties"
    1090 
    1091 ! II_d.4 Update the soil thermal properties
    1092       call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, &
    1093         porefillingice_depth,porefillingice_thickness,TI_PEM)
    1094 
    1095 ! II_d.5 Update the mass of the regolith adsorbded
    1096      if(adsorption_pem) then
    1097        call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:), &
    1098                              tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
    1099                              h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
    1100  
    1101          
    1102       totmassco2_adsorbded = 0.
    1103       totmassh2o_adsorbded = 0.
    1104       do ig = 1,ngrid
    1105        do islope =1, nslope
    1106         do l = 1,nsoilmx_PEM - 1
    1107            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
    1108            subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
    1109            cell_area(ig)
    1110            totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
    1111            subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
    1112            cell_area(ig)
    1113         enddo
    1114        enddo
    1115       enddo
    1116       write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
    1117       write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded               
    1118       endif
    1119   endif !soil_pem
    1120 
    1121 !------------------------
    1122 
    1123 ! II  Run
    1124 !    II_a update pressure, ice and tracers
    1125 !    II_b Evolution of the ice
    1126 !    II_c CO2 glaciers flows
    1127 !    II_d Update surface and soil temperatures
    1128 !    II_e Update the tendencies
    1129 
    1130 !------------------------
    1131 
    1132       print *, "Adaptation of the new co2 tendencies to the current pressure"
    1133       call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,&     
    1134                                global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
    1135 
    1136 !------------------------
    1137 
    1138 ! II  Run
    1139 !    II_a update pressure, ice and tracers
    1140 !    II_b Evolution of the ice
    1141 !    II_c CO2 glaciers flows
    1142 !    II_d Update surface and soil temperatures
    1143 !    II_e Update the tendencies
    1144935!    II_f Checking the stopping criterion
    1145 
    1146 !------------------------
    1147       call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
    1148 
    1149       call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
    1150                                    initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
    1151 
    1152       year_iter=year_iter+dt_pem     
    1153 
    1154       print *, "Checking all the stopping criterion."
    1155       if (STOPPING_water) then
    1156         print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
    1157         criterion_stop=1
    1158       endif
    1159       if (STOPPING_1_water) then
    1160         print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
    1161         criterion_stop=1
    1162       endif
    1163       if (STOPPING_co2) then
    1164         print *, "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
    1165         criterion_stop=2
    1166       endif
    1167       if (STOPPING_pressure) then
    1168         print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
    1169         criterion_stop=3
    1170       endif
    1171       if (year_iter.ge.year_iter_max) then
    1172         print *, "STOPPING because maximum number of iterations reached"
    1173         criterion_stop=4
    1174       endif
    1175 
    1176      if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
     936!------------------------
     937    call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
     938
     939    call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
     940                            initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
     941
     942    year_iter = year_iter + dt_pem
     943
     944    write(*,*) "Checking all the stopping criterion."
     945    if (STOPPING_water) then
     946        write(*,*) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
     947        criterion_stop = 1
     948    endif
     949    if (STOPPING_1_water) then
     950        write(*,*) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
     951        criterion_stop = 1
     952    endif
     953    if (STOPPING_co2) then
     954        write(*,*) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
     955        criterion_stop = 2
     956    endif
     957    if (STOPPING_pressure) then
     958        write(*,*) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
     959        criterion_stop = 3
     960    endif
     961    if (year_iter >= year_iter_max) then
     962        write(*,*) "STOPPING because maximum number of iterations reached"
     963        criterion_stop = 4
     964    endif
     965
     966    if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure)  then
    1177967        exit
    1178       else
    1179         print *, "We continue!"
    1180         print *, "Number of iteration of the PEM : year_iter=", year_iter
    1181       endif
    1182 
    1183       global_ave_press_old=global_ave_press_new
    1184 
    1185       enddo  ! big time iteration loop of the pem
    1186 
    1187 
    1188 !---------------------------- END RUN PEM ---------------------
    1189 
    1190 !---------------------------- OUTPUT ---------------------
    1191 
    1192 !------------------------
    1193 
     968    else
     969        write(*,*) "We continue!"
     970        write(*,*) "Number of iteration of the PEM : year_iter=", year_iter
     971    endif
     972
     973    global_ave_press_old=global_ave_press_new
     974
     975enddo  ! big time iteration loop of the pem
     976!------------------------------ END RUN --------------------------------
     977
     978!------------------------------- OUTPUT --------------------------------
     979!------------------------
    1194980! III Output
    1195981!    III_a Update surface value for the PCM start files
    1196 
    1197 !------------------------
    1198 
     982!------------------------
    1199983! III_a.1 Ice update (for startfi)
    1200984
    1201985! H2O ice
    1202      DO ig=1,ngrid
    1203        if(watercaptag(ig)) then
    1204          watercap_sum=0.
    1205          DO islope=1,nslope
    1206            if(qsurf(ig,igcm_h2o_ice,islope).GT. (watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy
    1207              qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost
    1208            else
    1209 !             2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
    1210              watercap(ig,islope)=watercap(ig,islope)+qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
    1211              qsurf(ig,igcm_h2o_ice,islope)=0.
     986do ig = 1,ngrid
     987    if (watercaptag(ig)) then
     988        watercap_sum = 0.
     989        do islope = 1,nslope
     990            if (qsurf(ig,igcm_h2o_ice,islope) > (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy
     991                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)) ! put back ancien frost
     992            else
     993!               2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
     994                watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.))
     995                qsurf(ig,igcm_h2o_ice,islope)=0.
     996            endif
     997            watercap_sum = watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
     998            watercap(ig,islope) = 0.
     999        enddo
     1000        water_reservoir(ig) = water_reservoir(ig) + watercap_sum
     1001    endif
     1002enddo
     1003
     1004do ig = 1,ngrid
     1005    water_sum = 0.
     1006    do islope = 1,nslope
     1007        water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
     1008    enddo
     1009    if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
     1010        if (water_sum > threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
     1011            watercaptag(ig) = .true.
     1012            water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perenial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
     1013            do islope = 1,nslope
     1014                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.)
     1015            enddo
     1016        endif
     1017    else ! let's check that the infinite source of water disapear
     1018        if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then
     1019            watercaptag(ig) = .false.
     1020            do islope = 1,nslope
     1021                qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
     1022            enddo
     1023            water_reservoir(ig) = 0.
     1024        endif
     1025    endif
     1026enddo
     1027
     1028! CO2 ice
     1029do ig = 1,ngrid
     1030    do islope = 1,nslope
     1031        if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then
     1032            perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
     1033            qsurf(ig,igcm_co2,islope)  = 0.5*qsurf(ig,igcm_co2,islope)
     1034        endif
     1035    enddo
     1036enddo
     1037
     1038! III_a.2  Tsoil update (for startfi)
     1039if (soil_pem) then
     1040    call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
     1041    tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)     
     1042endif
     1043
     1044! III_a.4 Pressure (for start)
     1045do i = 1,ip1jmp1
     1046    ps(i) = ps(i)*global_ave_press_new/global_ave_press_GCM
     1047enddo
     1048
     1049do i = 1,ngrid
     1050    ps_start_GCM(i) = ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM
     1051enddo
     1052
     1053! III_a.5 Tracer (for start)
     1054allocate(zplev_new(ngrid,nlayer + 1))
     1055
     1056do l = 1,nlayer + 1
     1057    do ig = 1,ngrid
     1058        zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
     1059    enddo
     1060enddo
     1061
     1062do nnq = 1,nqtot
     1063    if (noms(nnq) /= "co2") then
     1064        do l = 1,llm - 1
     1065            do ig = 1,ngrid
     1066                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))
     1067            enddo
     1068            q(:,llm,nnq) = q(:,llm - 1,nnq)
     1069        enddo
     1070    else
     1071        do l = 1,llm - 1
     1072            do ig = 1,ngrid
     1073                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)) &
     1074                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ &
     1075                              (zplev_new(ig,l) - zplev_new(ig,l + 1))
     1076            enddo
     1077            q(:,llm,nnq) = q(:,llm - 1,nnq)
     1078        enddo
     1079    endif
     1080enddo
     1081
     1082! Conserving the tracers mass for GCM start files
     1083do nnq = 1,nqtot
     1084    do ig = 1,ngrid
     1085        do l = 1,llm - 1
     1086            if (q(ig,l,nnq) > 1 .and. (noms(nnq) /= "dust_number") .and. (noms(nnq) /= "ccn_number") .and. (noms(nnq) /= "stormdust_number") .and. (noms(nnq) /= "topdust_number")) then
     1087                extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
     1088                q(ig,l,nnq)=1.
     1089                q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
     1090                write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
    12121091           endif
    1213            watercap_sum=watercap_sum+watercap(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    1214            watercap(ig,islope)=0.
    1215          enddo
    1216            water_reservoir(ig)=water_reservoir(ig)+watercap_sum
    1217        endif
    1218      enddo
    1219 
    1220      DO ig=1,ngrid
    1221        water_sum = 0.
    1222        DO islope=1,nslope
    1223          water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
    1224        ENDDO
    1225        if(watercaptag(ig).eqv..false.) then ! let's check if we have an 'infinite' source of water that has been forming.
    1226          if(water_sum.gt.threshold_water_frost2perenial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
    1227             watercaptag(ig)=.true.
    1228             water_reservoir(ig)=water_reservoir(ig)+threshold_water_frost2perenial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
    1229             DO islope = 1,nslope
    1230               qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.*cos(pi*def_slope_mean(islope)/180.)
    1231             ENDDO   
    1232          endif
    1233        else ! let's check that the infinite source of water disapear
    1234          if((water_sum + water_reservoir(ig)).lt.threshold_water_frost2perenial) then
    1235            watercaptag(ig)=.false.
    1236            DO islope = 1,nslope
    1237               qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
    1238            ENDDO
    1239            water_reservoir(ig) = 0.
    1240          endif
    1241        endif
    1242      enddo
    1243 
    1244 ! CO2 ice
    1245 
    1246      DO ig=1,ngrid
    1247        DO islope=1,nslope
    1248          if(qsurf(ig,igcm_co2,islope).gt.threshold_co2_frost2perenial) then
    1249            perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
    1250            qsurf(ig,igcm_co2,islope)  = 0.5*qsurf(ig,igcm_co2,islope)
    1251          endif
    1252        ENDDO
    1253      ENDDO
    1254        
    1255 
    1256 ! III_a.2  Tsoil update (for startfi)
    1257 
    1258    if(soil_pem) then
    1259      call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
    1260      tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)     
    1261    endif !soil_pem
    1262 
    1263 ! III_a.4 Pressure (for start)
    1264      do i=1,ip1jmp1
    1265        ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
    1266      enddo
    1267 
    1268      do i = 1,ngrid
    1269        ps_start_GCM(i)=ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM
    1270      enddo
    1271 
    1272 ! III_a.5 Tracer (for start)
    1273      allocate(zplev_new(ngrid,nlayer+1))
    1274 
    1275      do l=1,nlayer+1
    1276        do ig=1,ngrid
    1277          zplev_new(ig,l) = ap(l)  + bp(l)*ps_start_GCM(ig)
    1278        enddo
    1279      enddo
    1280 
    1281      DO nnq=1,nqtot
    1282        if (noms(nnq).NE."co2") then
    1283          DO l=1,llm-1
    1284            DO ig=1,ngrid
    1285              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))
    1286            ENDDO
    1287            q(:,llm,nnq)=q(:,llm-1,nnq)
    1288          ENDDO
    1289        else
    1290          DO l=1,llm-1
    1291            DO ig=1,ngrid
    1292              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))  &
    1293                              + (   (zplev_new(ig,l)-zplev_new(ig,l+1))  -       &
    1294                                    (zplev_gcm(ig,l)-zplev_gcm(ig,l+1))     )  / &
    1295                                    (zplev_new(ig,l)-zplev_new(ig,l+1))
    1296            ENDDO
    1297            q(:,llm,nnq)=q(:,llm-1,nnq)
    1298          ENDDO
    1299        endif
    1300      ENDDO
    1301 
    1302 ! Conserving the tracers's mass for GCM start files
    1303      DO nnq=1,nqtot
    1304        DO ig=1,ngrid
    1305          DO l=1,llm-1
    1306            if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") .and. (noms(nnq).NE."stormdust_number") .and. (noms(nnq).NE."topdust_number")) then
    1307              extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
    1308              q(ig,l,nnq)=1.
    1309              q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
    1310              write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number"
    1311            endif
    1312            if(q(ig,l,nnq).LT.0) then
    1313              q(ig,l,nnq)=1E-30
    1314            endif
    1315          ENDDO
    1316        enddo
    1317      enddo
    1318    
    1319 !------------------------
    1320      if(evol_orbit_pem) then
    1321       call recomp_orb_param(year_iter)
    1322      endif
    1323 
     1092            if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
     1093        enddo
     1094    enddo
     1095enddo
     1096
     1097if (evol_orbit_pem) call recomp_orb_param(year_iter)
     1098
     1099!------------------------
    13241100! III Output
    1325 !    III_a Update surface value for the PCM start files
    1326 !    III_b Write start and starfi.nc
    1327 
    1328 !------------------------
    1329 
    1330 ! III_b.1 WRITE restart.nc
    1331 
    1332       ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
    1333       pday=day_ini
    1334       ztime_fin=0.
    1335 
    1336      allocate(p(ip1jmp1,nlayer+1))
     1101!    III_b Write restart_evol.nc and restartfi_evol.nc
     1102!------------------------
     1103! III_b.1 Write restart_evol.nc
     1104ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys
     1105pday = day_ini
     1106ztime_fin = 0.
     1107
     1108allocate(p(ip1jmp1,nlayer + 1))
    13371109#ifndef CPP_1D
    1338      CALL pression (ip1jmp1,ap,bp,ps,p)
    1339      CALL massdair(p,masse)
    1340 
    1341       CALL dynredem0("restart_evol.nc", day_ini, phis)
    1342 
    1343       CALL dynredem1("restart_evol.nc",   &
    1344                 time_0,vcov,ucov,teta,q,masse,ps)
    1345       print *, "restart_evol.nc has been written"
    1346 
     1110    call pression (ip1jmp1,ap,bp,ps,p)
     1111    call massdair(p,masse)
     1112    call dynredem0("restart_evol.nc", day_ini, phis)
     1113    call dynredem1("restart_evol.nc",   &
     1114    time_0,vcov,ucov,teta,q,masse,ps)
     1115    write(*,*) "restart_evol.nc has been written"
    13471116#else
    1348      DO nnq = 1, nqtot
    1349         call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf)
    1350      ENDDO
     1117    do nnq = 1, nqtot
     1118        call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf)
     1119    enddo
    13511120#endif
    13521121
    1353 ! III_b.2 WRITE restartfi.nc
     1122! III_b.2 Write restartfi_evol.nc
    13541123#ifndef CPP_STD
    1355       call physdem0("restartfi_evol.nc",longitude,latitude, &
    1356                         nsoilmx,ngrid,nlayer,nq,              &
    1357                         ptimestep,pday,0.,cell_area,          &
    1358                         albedodat,inertiedat,def_slope,   &
    1359                         subslope_dist)
    1360 
    1361      call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
    1362                      ptimestep,ztime_fin,                        &
    1363                      tsurf,tsoil,inertiesoil,albedo,             &
    1364                      emis,q2,qsurf,tauscaling,totcloudfrac,wstar,&
    1365                      watercap,perenial_co2ice)
     1124    call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
     1125                  nlayer,nq,ptimestep,pday,0.,cell_area,albedodat,      &
     1126                  inertiedat,def_slope,subslope_dist)
     1127
     1128    call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,  &
     1129                  ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,  &
     1130                  albedo,emis,q2,qsurf,tauscaling,totcloudfrac, &
     1131                  wstar,watercap,perenial_co2ice)
    13661132#else
    1367      call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
    1368                          ptimestep,pday,time_phys,cell_area,          &
    1369                          albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
    1370 
    1371      call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
    1372                           ptimestep,ztime_fin,                    &
    1373                           tsurf,tsoil,emis,q2,qsurf,         &
    1374                           cloudfrac,totcloudfrac,hice,            &
    1375                           rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
     1133    call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
     1134                  nlayer,nq,ptimestep,pday,time_phys,cell_area,         &
     1135                  albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
     1136
     1137    call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,   &
     1138                  ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf, &
     1139                  cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,   &
     1140                  tslab,tsea_ice,sea_ice)
    13761141#endif
    1377 
    1378       print *, "restartfi_evol.nc has been written"
    1379 !------------------------
    1380 
     1142write(*,*) "restartfi_evol.nc has been written"
     1143
     1144!------------------------
    13811145! III Output
    1382 !    III_a Update surface value for the PCM start files
    1383 !    III_b Write start and starfi.nc
    1384 !    III_c Write start_pem
    1385 
    1386 !------------------------
    1387         call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
    1388                          float(day_ini),0.,nslope,def_slope,subslope_dist)
    1389 
    1390 
    1391         call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
    1392                       tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, &
    1393                       co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
    1394 
    1395         call info_run_PEM(year_iter, criterion_stop)
    1396 
    1397       print *, "restartfi_PEM.nc has been written"
    1398       print *, "The PEM had run for ", year_iter, " years."
    1399       print *, "LL & RV : So far so good"
    1400 
    1401      deallocate(vmr_co2_gcm)
    1402      deallocate(ps_timeseries)
    1403      deallocate(tsurf_GCM_timeseries)
    1404      deallocate(q_co2_PEM_phys)
    1405      deallocate(q_h2o_PEM_phys)
    1406      deallocate(co2_ice_GCM)
    1407      deallocate(watersurf_density_ave)
    1408      deallocate(watersoil_density_timeseries)
    1409      deallocate(Tsurfave_before_saved)
    1410      deallocate(tsoil_phys_PEM_timeseries)
    1411      deallocate(watersoil_density_PEM_timeseries)
    1412      deallocate(watersoil_density_PEM_ave)
    1413      deallocate(delta_co2_adsorbded)
    1414      deallocate(delta_h2o_adsorbded)
    1415      deallocate(vmr_co2_pem_phys)
     1146!    III_c Write restartfi_PEM.nc
     1147!------------------------
     1148call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
     1149             float(day_ini),0.,nslope,def_slope,subslope_dist)
     1150
     1151
     1152call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope ,           &
     1153             tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, &
     1154             co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
     1155
     1156call info_run_PEM(year_iter,criterion_stop)
     1157
     1158write(*,*) "restartfi_PEM.nc has been written"
     1159write(*,*) "The PEM had run for ", year_iter, " years."
     1160write(*,*) "LL & RV : So far so good"
     1161
     1162deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
     1163deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
     1164deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
     1165deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys)
     1166!----------------------------- END OUTPUT ------------------------------
    14161167
    14171168END PROGRAM pem
     1169
  • trunk/LMDZ.MARS/changelog.txt

    r3026 r3028  
    41694169== 09/08/2023 == JBC
    41704170Small fixes to be able to run the Mars PCM 1D without "water" + Improvements/addition of scripts in deftank/pem to run the PEM 1D model according to Laskar orbital parameters.
     4171
     4172== 14/08/2023 == JBC
     4173Related to commit r3026: improvement of error message in initracer.F (now it gives correctly the only identified tracers) + one small correction to run PCM 1D without water.
     4174Big cleaning of main program pem.F90 (indentation, declarations, comments, simplification of conditions/loops, etc).
Note: See TracChangeset for help on using the changeset viewer.