Changeset 3571 for trunk/LMDZ.COMMON

Timestamp:

Jan 10, 2025, 5:45:03 PM (10 days ago)

Author:

jbclement

Message:

PEM:

• New way to manage the pressure: now the PEM manages only the average pressure and keeps the pressure deviation with the instantaneous pressure from the start to reconstruct the pressure at the end ('ps_avg = ps_start + ps_dev'). As a consequence, everything related to pressure in the PEM is modified accordingly.
• Surface temperatures management is now simpler. It follows the strategy for the pressure (and soil temperature) described above.
• Soil temperatures are now adapted to match the surface temperature changes occured during the PEM by modifying the soil temperature deviation at the end.
• Few simplifications/optimizations: notably, the two PCM years are now read in one go in 'read_data_PCM_mod.F90' and only the needed variables are extracted.
• Deletion of unused variables and unnecessary intermediate variables (memory saving and loop deletion in some cases).
• Renaming of variables and subroutines to make everything clearer. In particular, the suffixes: '_avg' = average, '_start' = PCM start file, '_dev' = deviation, '_ini' or '0' = initial, '_dyn' = dynamical grid, '_timeseries' = daily average of last PCM year.
• Cosmetic cleanings for readability.

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• adsorption_mod.F90 (modified) (1 diff)
• changelog.txt (modified) (1 diff)
• evol_ice_mod.F90 (modified) (1 diff)
• get_timelen_PCM_mod.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/evolution/nb_time_step_PCM_mod.F90) (6 diffs)
• glaciers_mod.F90 (modified) (10 diffs)
• ice_table_mod.F90 (modified) (1 diff)
• lask_param_mod.F90 (modified) (1 diff)
• pem.F90 (modified) (49 diffs)
• pemetat0.F90 (modified) (14 diffs)
• pemredem.F90 (modified) (2 diffs)
• read_data_PCM_mod.F90 (modified) (10 diffs)
• recomp_tend_co2_mod.F90 (moved) (moved from trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90) (5 diffs)
• stopping_crit_mod.F90 (modified) (6 diffs)

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

@@ -3 +3 @@
 implicit none
 
-logical                                   :: adsorption_pem     ! True by default, to compute adsorption/desorption. Read in pem.def
-real, save, allocatable, dimension(:,:,:) :: co2_adsorbed_phys ! co2 that is in the regolith [kg/m^2]
-real, save, allocatable, dimension(:,:,:) :: h2o_adsorbed_phys ! h2o that is in the regolith [kg/m^2]
+logical                             :: adsorption_pem     ! True by default, to compute adsorption/desorption. Read in pem.def
+real, allocatable, dimension(:,:,:) :: co2_adsorbed_phys ! co2 that is in the regolith [kg/m^2]
+real, allocatable, dimension(:,:,:) :: h2o_adsorbed_phys ! h2o that is in the regolith [kg/m^2]
 
 !=======================================================================

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

@@ -528 +528 @@
 == 17/12/2024 == JBC
 Correction of Norbert Schorghofer's code due to missing initialization and bad shape array as subroutine argument + some cleanings.
+
+== 10/01/2025 == JBC
+- New way to manage the pressure: now the PEM manages only the average pressure and keeps the pressure deviation with the instantaneous pressure from the start to reconstruct the pressure at the end ('ps_avg = ps_start + ps_dev'). As a consequence, everything related to pressure in the PEM is modified accordingly.
+- Surface temperatures management is now simpler. It follows the strategy for the pressure (and soil temperature) described above.
+- Soil temperatures are now adapted to match the surface temperature changes occured during the PEM by modifying the soil temperature deviation at the end.
+- Few simplifications/optimizations: notably, the two PCM years are now read in one go in 'read_data_PCM_mod.F90' and only the needed variables are extracted.
+- Deletion of unused variables and unnecessary intermediate variables (memory saving and loop deletion in some cases).
+- Renaming of variables and subroutines to make everything clearer. In particular, the suffixes: '_avg' = average, '_start' = PCM start file, '_dev' = deviation, '_ini' or '0' = initial, '_dyn' = dynamical grid, '_timeseries' = daily average of last PCM year.
+- Cosmetic cleanings for readability.

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

@@ -77 +77 @@
 
 ! OUTPUT
-real, dimension(ngrid,nslope), intent(inout) :: h2o_ice     ! Previous and actual density of h2o ice (kg/m^2)
-real, dimension(ngrid,nslope), intent(inout) :: d_h2oice    ! Evolution of perennial ice over one year (kg/m^2/year)
-integer,                       intent(inout) :: stopPEM     ! Stopping criterion code
+real, dimension(ngrid,nslope), intent(inout) :: h2o_ice  ! Previous and actual density of h2o ice (kg/m^2)
+real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Evolution of perennial ice over one year (kg/m^2/year)
+integer,                       intent(inout) :: stopPEM  ! Stopping criterion code
 real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
 

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

@@ -1 @@
-MODULE nb_time_step_PCM_mod
+MODULE get_timelen_PCM_mod
 
 use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
@@ -13 +13 @@
 !=======================================================================
 
-SUBROUTINE nb_time_step_PCM(fichnom,timelen)
+SUBROUTINE get_timelen_PCM(fichnom,timelen)
 
 implicit none
@@ -30 +30 @@
 character(len = *), intent(in) :: fichnom !--- FILE NAME
 !=======================================================================
-!   Local Variables
+!   Local variables
 integer        :: fID, vID, ierr
 integer        :: timelen ! number of times stored in the file
 !-----------------------------------------------------------------------
-modname = "nb_time_step_PCM"
+modname = "get_timelen_PCM"
 
 !  Open initial state NetCDF file
@@ -52 +52 @@
             write(*,*)"read_data_PCM: Le champ <Time> est absent"
             write(*,*)trim(nf90_strerror(ierr))
-            call abort_gcm("nb_time_step_PCM", "", 1)
+            call abort_gcm("get_timelen_PCM", "", 1)
         endif
         ! Get the length of the "Time" dimension
@@ -72 +72 @@
 write(*,*) "The number of timestep of the PCM run data=", timelen
 
-END SUBROUTINE nb_time_step_PCM
+END SUBROUTINE get_timelen_PCM
 
 !=======================================================================
@@ -98 +98 @@
 END SUBROUTINE error_msg
 
-END MODULE nb_time_step_PCM_mod
+END MODULE get_timelen_PCM_mod

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

@@ -21 +21 @@
 !=======================================================================
 
-SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,co2ice,flag_co2flow,flag_co2flow_mesh)
+SUBROUTINE flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,co2ice,flag_co2flow,flag_co2flow_mesh)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -38 +38 @@
 !--------
 integer,                           intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
-real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Physical points x Slopes: Distribution of the subgrid slopes
-real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Physical points: values of the sub grid slope angles
-real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
-real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Physical x Time field: surface pressure given by the PCM [Pa]
-real,                              intent(in) :: global_avg_ps_PCM             ! Global averaged surface pressure from the PCM [Pa]
-real,                              intent(in) :: global_avg_ps_PEM             ! global averaged surface pressure during the PEM iteration [Pa]
+real,    dimension(ngrid,nslope),  intent(in) :: subslope_dist                 ! Grid points x Slopes: Distribution of the subgrid slopes
+real,    dimension(ngrid),         intent(in) :: def_slope_mean                ! Grid points: values of the sub grid slope angles
+real,    dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM                   ! Grid points x Time field : VMR of co2 in the first layer [mol/mol]
+real,    dimension(ngrid,timelen), intent(in) :: ps_PCM                        ! Grid points x Time field: surface pressure given by the PCM [Pa]
+real,                              intent(in) :: ps_avg_global_ini             ! Global averaged surface pressure at the beginning [Pa]
+real,                              intent(in) :: ps_avg_global                 ! Global averaged surface pressure during the PEM iteration [Pa]
 
 ! Ouputs:
 !--------
-real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
+real, dimension(ngrid,nslope), intent(inout) :: co2ice            ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
 real, dimension(ngrid,nslope), intent(inout) :: flag_co2flow      ! flag to see if there is flow on the subgrid slopes
 real, dimension(ngrid),        intent(inout) :: flag_co2flow_mesh ! same but within the mesh
@@ -53 +53 @@
 !------
 real, dimension(ngrid,nslope) :: Tcond ! Physical field: CO2 condensation temperature [K]
-real, dimension(ngrid,nslope) :: hmax  ! Physical x Slope field: maximum thickness for co2  glacier before flow
+real, dimension(ngrid,nslope) :: hmax  ! Grid points x Slope field: maximum thickness for co2  glacier before flow
 
 write(*,*) "Flow of CO2 glaciers"
-call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
+call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
 call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tcond,"co2",hmax)
-call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,"co2",co2ice,flag_co2flow,flag_co2flow_mesh)
+call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tcond,co2ice,flag_co2flow,flag_co2flow_mesh)
 
 END SUBROUTINE flow_co2glaciers
@@ -64 +64 @@
 !=======================================================================
 
-SUBROUTINE flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
+SUBROUTINE flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -82 +82 @@
 
 ! Inputs:
-integer,                       intent(in) :: timelen, ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
-real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Physical points x Slopes : Distribution of the subgrid slopes
+integer,                       intent(in) :: ngrid, nslope, iflat ! number of time sample, physical points, subslopes, index of the flat subslope
+real, dimension(ngrid,nslope), intent(in) :: subslope_dist  ! Grid points x Slopes : Distribution of the subgrid slopes
 real, dimension(ngrid),        intent(in) :: def_slope_mean ! Slopes: values of the sub grid slope angles
 real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice Temperature [K]
 ! Ouputs:
-real, dimension(ngrid,nslope), intent(inout) :: h2oice            ! Physical x Slope field: co2 ice on the subgrid slopes [kg/m^2]
+real, dimension(ngrid,nslope), intent(inout) :: h2oice            ! Grid points x Slope field: co2 ice on the subgrid slopes [kg/m^2]
 real, dimension(ngrid,nslope), intent(inout) :: flag_h2oflow      ! flag to see if there is flow on the subgrid slopes
 real, dimension(ngrid),        intent(inout) :: flag_h2oflow_mesh ! same but within the mesh
 ! Local
-real, dimension(ngrid,nslope) :: hmax ! Physical x Slope field: maximum thickness for co2  glacier before flow
+real, dimension(ngrid,nslope) :: hmax ! Grid points x Slope field: maximum thickness for co2  glacier before flow
 
 write(*,*) "Flow of H2O glaciers"
 call compute_hmaxglaciers(ngrid,nslope,iflat,def_slope_mean,Tice,"h2o",hmax)
-call transfer_ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tice,"h2o",h2oice,flag_h2oflow,flag_h2oflow_mesh)
+call transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,h2oice,flag_h2oflow,flag_h2oflow_mesh)
 
 END SUBROUTINE flow_h2oglaciers
@@ -161 +161 @@
 !=======================================================================
 
-SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,name_ice,qice,flag_flow,flag_flowmesh)
+SUBROUTINE transfer_ice_duringflow(ngrid,nslope,iflat,subslope_dist,def_slope_mean,hmax,Tice,qice,flag_flow,flag_flowmesh)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -188 +188 @@
 real, dimension(ngrid,nslope), intent(in) :: hmax           ! maximum height of the  glaciers before initiating flow [m]
 real, dimension(ngrid,nslope), intent(in) :: Tice           ! Ice temperature[K]
-character(3),                  intent(in) :: name_ice       ! Nature of the ice
 ! Outputs
 !--------
@@ -232 +231 @@
 !=======================================================================
 
-SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,global_avg_ps_PCM,global_avg_ps_PEM,Tcond)
+SUBROUTINE computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global,Tcond)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -241 +240 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-use constants_marspem_mod, only: alpha_clap_co2,beta_clap_co2
+use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2
 
 implicit none
@@ -250 +249 @@
 ! INPUT
 integer,                        intent(in) :: timelen, ngrid, nslope ! # of timesample, physical points, subslopes
-real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Physical points x times field: VMR of CO2 in the first layer [mol/mol]
-real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Physical points x times field: surface pressure in the PCM [Pa]
-real,                           intent(in) :: global_avg_ps_PCM      ! Global averaged surfacepressure in the PCM [Pa]
-real,                           intent(in) :: global_avg_ps_PEM      ! Global averaged surface pressure computed during the PEM iteration
+real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM            ! Grid points x times field: VMR of CO2 in the first layer [mol/mol]
+real, dimension(ngrid,timelen), intent(in) :: ps_PCM                 ! Grid points x times field: surface pressure in the PCM [Pa]
+real,                           intent(in) :: ps_avg_global_ini      ! Global averaged surfacepressure in the PCM [Pa]
+real,                           intent(in) :: ps_avg_global          ! Global averaged surface pressure computed during the PEM iteration
 
 ! OUTPUT
-real, dimension(ngrid,nslope), intent(out) :: Tcond ! Physical points: condensation temperature of CO2, yearly averaged
+real, dimension(ngrid,nslope), intent(out) :: Tcond ! Grid points: condensation temperature of CO2, yearly averaged
 
 ! LOCAL
 integer :: ig, it ! For loop
-real    :: ave    ! Intermediate to compute average
 
 do ig = 1,ngrid
-    ave = 0
+    Tcond(ig,:) = 0
     do it = 1,timelen
-        ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*global_avg_ps_PCM/global_avg_ps_PEM/100))
+        Tcond(ig,:) = Tcond(ig,:) + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM(ig,it)*ps_PCM(ig,it)*ps_avg_global_ini/ps_avg_global/100))
     enddo
-    Tcond(ig,:) = ave/timelen
+    Tcond(ig,:) = Tcond(ig,:)/timelen
 enddo
 

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

@@ -383 +383 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-use comsoil_h_PEM, only: layer_PEM, mlayer_PEM
+use comsoil_h_PEM, only: mlayer_PEM
 use abort_pem_mod, only: abort_pem
 

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

@@ -11 +11 @@
 implicit none
 
-real, dimension(:), allocatable, save :: yearlask   ! year before present from Laskar et al. Tab
-real, dimension(:), allocatable, save :: obllask    ! obliquity    [deg]
-real, dimension(:), allocatable, save :: ecclask    ! eccentricity [deg]
-real, dimension(:), allocatable, save :: lsplask    ! ls perihelie [deg]
-integer,                         save :: last_ilask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
+real, dimension(:), allocatable :: yearlask   ! year before present from Laskar et al. Tab
+real, dimension(:), allocatable :: obllask    ! obliquity    [deg]
+real, dimension(:), allocatable :: ecclask    ! eccentricity [deg]
+real, dimension(:), allocatable :: lsplask    ! ls perihelie [deg]
+integer                         :: last_ilask ! Index of the line in the file year_obl_lask.asc corresponding to the closest lower year to the current year
 
 !=======================================================================

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

@@ -7 +7 @@
 !    I_e Initialization of the PEM variable and soil
 !    I_f Compute tendencies
-!    I_g Save the initial situation
+!    I_g Compute global surface pressure
 !    I_h Read the "startpem.nc"
 !    I_i Compute orbit criterion
@@ -40 +40 @@
 use pemredem,                   only: pemdem0, pemdem1
 use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
-                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice
+                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
 use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2
 use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
 use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
-use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, mlayer_PEM, &
-                                      TI_PEM,               & ! soil thermal inertia
+use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
+                                      TI_PEM,               & ! Soil thermal inertia
                                       tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
                                       fluxgeo                 ! Geothermal flux for the PEM and PCM
 use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
                                       ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
-                                      co2_adsorbed_phys, h2o_adsorbed_phys        ! Mass of co2 and h2O adsorbed
+                                      co2_adsorbed_phys, h2o_adsorbed_phys          ! Mass of co2 and h2O adsorbed
 use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
 use orbit_param_criterion_mod,  only: orbit_param_criterion
@@ -62 +62 @@
 use compute_tend_mod,           only: compute_tend
 use info_PEM_mod,               only: info_PEM
-use nb_time_step_PCM_mod,       only: nb_time_step_PCM
+use get_timelen_PCM_mod,        only: get_timelen_PCM
 use pemetat0_mod,               only: pemetat0
 use read_data_PCM_mod,          only: read_data_PCM
-use recomp_tend_co2_slope_mod,  only: recomp_tend_co2
+use recomp_tend_co2_mod,        only: recomp_tend_co2
 use compute_soiltemp_mod,       only: compute_tsoil_pem, shift_tsoil2surf
 use writediagpem_mod,           only: writediagpem, writediagsoilpem
@@ -77 +77 @@
                                   albedodat, zmea, zstd, zsig, zgam, zthe,     &
                                   albedo_h2o_frost,frost_albedo_threshold,     &
-                                  emissiv, watercaptag, perennial_co2ice, emisice, albedice
+                                  emissiv, watercaptag, perennial_co2ice
     use dimradmars_mod,     only: totcloudfrac, albedo
     use dust_param_mod,     only: tauscaling
@@ -128 +128 @@
 real               :: time_phys    ! Same as in PCM
 real               :: ptimestep    ! Same as in PCM
-real               :: ztime_fin     ! Same as in PCM
-
-! Variables to read start.nc
+real               :: ztime_fin    ! Same as in PCM
+
+! Variables to read "start.nc"
 character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM
 
@@ -136 +136 @@
 real, dimension(ip1jm,llm)          :: vcov          ! vents covariants
 real, dimension(ip1jmp1,llm)        :: ucov          ! vents covariants
-real, dimension(ip1jmp1,llm)        :: teta          ! temperature potentielle
+real, dimension(ip1jmp1,llm)        :: teta          ! Potential temperature
 real, dimension(:,:,:), allocatable :: q             ! champs advectes
-real, dimension(ip1jmp1)            :: ps            ! pression au sol
-real, dimension(ip1jmp1)            :: ps0           ! pression au sol initiale
-real, dimension(:),     allocatable :: ps_start_PCM  ! (ngrid) surface pressure
-real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure
-real, dimension(ip1jmp1,llm)        :: masse         ! masse d'air
+real, dimension(ip1jmp1)            :: ps_start_dyn  ! surface pressure in the start file (dynamic grid)
+real, dimension(:),     allocatable :: ps_start      ! surface pressure in the start file
+real, dimension(:),     allocatable :: ps_start0     ! surface pressure in the start file at the beginning
+real, dimension(:),     allocatable :: ps_avg        ! (ngrid) Averaged surface pressure
+real, dimension(:),     allocatable :: ps_dev        ! (ngrid x timelen) Surface pressure deviation
+real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure
+real, dimension(ip1jmp1,llm)        :: masse         ! Air mass
 real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
 real                                :: time_0
@@ -157 +159 @@
 real, dimension(:), allocatable :: latitude      ! Latitude read in startfi_name and written in restartfi
 real, dimension(:), allocatable :: cell_area     ! Cell_area read in startfi_name and written in restartfi
-real                            :: Total_surface ! Total surface of the planet
+real                            :: total_surface ! Total surface of the planet
 
 ! Variables for h2o_ice evolution
 real, dimension(:,:),    allocatable  :: h2o_ice              ! h2o ice in the PEM
+real, dimension(:,:),    allocatable  :: d_h2oice             ! physical point x slope field: Tendency of evolution of perennial h2o ice
 real, dimension(:,:,:),  allocatable  :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
 real                                  :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
-logical, dimension(:,:), allocatable  :: ini_h2oice_sublim    ! Logical array to know if h2o ice is sublimating
-real                                  :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
-real                                  :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
-real                                  :: global_avg_press_new ! constant: Global average pressure of current time step
-real,   dimension(:,:),   allocatable :: zplev_new            ! Physical x Atmospheric field: mass of the atmospheric  layers in the pem at current time step [kg/m^2]
-real,   dimension(:,:),   allocatable :: zplev_PCM            ! same but retrieved from the PCM [kg/m^2]
-real,   dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
-real,   dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
+logical, dimension(:,:), allocatable  :: is_h2oice_sublim_ini ! Logical array to know if h2o ice is sublimating
+real                                  :: ps_avg_global_ini    ! constant: Global average pressure at initialization [Pa]
+real                                  :: ps_avg_global_old    ! constant: Global average pressure of previous time step
+real                                  :: ps_avg_global_new    ! constant: Global average pressure of current time step
+real,   dimension(:,:),   allocatable :: zplev_new            ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
+real,   dimension(:,:),   allocatable :: zplev_start0         ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
+real,   dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
+real,   dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
 integer                               :: stopPEM              ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu
-real, save                            :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
-real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
+real                                  :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
+real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Grid points x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
 integer                               :: timelen              ! # time samples
-real                                  :: ave                  ! intermediate varibale to compute average
-real,   dimension(:,:),   allocatable :: p                    ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
+real,   dimension(:,:),   allocatable :: p                    ! Grid points x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
 real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
 
 ! Variables for co2_ice evolution
-real,    dimension(:,:), allocatable :: co2_ice           ! co2 ice in the PEM
-logical, dimension(:,:), allocatable :: is_co2ice_ini     ! Was there co2 ice initially in the PEM?
-real,  dimension(:,:,:), allocatable :: min_co2_ice       ! Minimum of co2 ice at each point for the first year [kg/m^2]
-real                                 :: co2ice_ini_surf   ! Initial surface of sublimating co2 ice
-logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating
-real,    dimension(:,:), allocatable :: vmr_co2_PCM       ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
-real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys  ! Physics x Times co2 volume mixing ratio used in the PEM
-real,    dimension(:,:), allocatable :: q_co2_PEM_phys    ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
+real,    dimension(:,:), allocatable :: co2_ice                ! co2 ice in the PEM
+real,    dimension(:,:), allocatable :: d_co2ice               ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
+real,    dimension(:,:), allocatable :: d_co2ice_ini           ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
+logical, dimension(:,:), allocatable :: is_co2ice_ini          ! Was there co2 ice initially in the PEM?
+real,  dimension(:,:,:), allocatable :: min_co2_ice            ! Minimum of co2 ice at each point for the first year [kg/m^2]
+real                                 :: co2ice_sublim_surf_ini ! Initial surface of sublimating co2 ice
+logical, dimension(:,:), allocatable :: is_co2ice_sublim_ini   ! Logical array to know if co2 ice is sublimating
+real,    dimension(:,:), allocatable :: vmr_co2_PCM            ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
+real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys       ! Grid points x Times co2 volume mixing ratio used in the PEM
+real,    dimension(:,:), allocatable :: q_co2_PEM_phys         ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
 
 ! Variables for stratification (layering) evolution
@@ -195 +199 @@
 
 ! Variables for slopes
-real, dimension(:,:),   allocatable :: d_co2ice          ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
-real, dimension(:,:),   allocatable :: d_co2ice_ini      ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
-real, dimension(:,:),   allocatable :: d_h2oice          ! physical point x slope field: Tendency of evolution of perennial h2o ice
 real, dimension(:,:),   allocatable :: flag_co2flow      ! (ngrid,nslope): Flag where there is a CO2 glacier flow
 real, dimension(:),     allocatable :: flag_co2flow_mesh ! (ngrid)       : Flag where there is a CO2 glacier flow
@@ -204 +205 @@
 
 ! Variables for surface and soil
-real, dimension(:,:),     allocatable :: tsurf_avg                        ! Physic x SLOPE field: Averaged Surface Temperature [K]
-real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
-real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries             ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries        ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_anom                       ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K]
-real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
-real, dimension(:,:),     allocatable :: Tsoil_locslope                   ! Physic x Soil: Intermediate when computing Tsoil [K]
-real, dimension(:),       allocatable :: Tsurf_locslope                   ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
-real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
-real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
-real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
-real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
-real, dimension(:,:),     allocatable :: Tsurfavg_before_saved            ! Surface temperature saved from previous time step [K]
+real, dimension(:,:),     allocatable :: tsurf_avg                        ! Grid points x Slope field: Averaged surface temperature [K]
+real, dimension(:,:),     allocatable :: tsurf_dev                        ! ngrid x Slope x Times field: Surface temperature deviation [K]
+real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Grid points x Slope field: Averaged surface temperature of first call of the PCM [K]
+real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Grid points x Soil x Slope field: Averaged Soil Temperature [K]
+real, dimension(:,:),     allocatable :: tsoil_avg_old                    ! Grid points x Soil field: Averaged Soil Temperature at the previous time step [K]
+real, dimension(:,:,:),   allocatable :: tsoil_dev                        ! Grid points x Soil x Slope field: Soil temperature deviation [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_timeseries                 ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries             ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K]
+real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
+real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Grid points x Slope: Averaged  water surface density [kg/m^3]
+real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
+real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Averaged water soil density [kg/m^3]
+real, dimension(:,:),     allocatable :: tsurf_avg_old                    ! Surface temperature saved from previous time step [K]
 real, dimension(:),       allocatable :: delta_co2_adsorbed               ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
 real, dimension(:),       allocatable :: delta_h2o_adsorbed               ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
 real                                  :: totmassco2_adsorbed              ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
 real                                  :: totmassh2o_adsorbed              ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
-logical                               :: bool_sublim                      ! logical to check if there is sublimation or not
 logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
 real, dimension(:,:),     allocatable :: icetable_thickness_old           ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
@@ -233 +233 @@
 ! Some variables for the PEM run
 real, parameter :: year_step = 1   ! Timestep for the pem
-real            :: i_myear_leg       ! Number of iteration
+real            :: i_myear_leg     ! Number of iteration
 real            :: n_myear_leg     ! Maximum number of iterations before stopping
 real            :: i_myear         ! Global number of Martian years of the chained simulations
@@ -249 +249 @@
 
 #ifdef CPP_STD
-    real                                :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
-    real                                :: albedo_h2o_frost              ! albedo of h2o frost
+    real                                :: frost_albedo_threshold = 0.05 ! Frost albedo threeshold to convert fresh frost to old ice
+    real                                :: albedo_h2o_frost              ! Albedo of h2o frost
     real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
     real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
@@ -286 +286 @@
 
 ! Loop variables
-integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
+integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil
 
 ! Elapsed time with system clock
@@ -379 +379 @@
     endif
 
-    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,        &
-                         time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
+    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,            &
+                         time_0,ps_start_dyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
                          play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
-    ps(2) = ps(1)
+    ps_start_dyn(2) = ps_start_dyn(1)
     nsplit_phys = 1
     day_step = steps_per_sol
@@ -391 +391 @@
 !------------------------
 ! I   Initialization
-!    I_b Read of the "start.nc" and starfi_evol.nc
+!    I_b Read of the "start.nc" and "starfi.nc"
 !------------------------
 ! I_b.1 Read "start.nc"
-allocate(ps_start_PCM(ngrid))
+allocate(ps_start(ngrid),ps_start0(ngrid))
 #ifndef CPP_1D
-    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
-
-    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
+    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps_start_dyn,phis,time_0)
+
+    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps_start_dyn,ps_start)
 
     call iniconst !new
@@ -412 +412 @@
     status = nf90_close(ncid)
 
-    call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
+    call iniphysiq(iim,jjm,llm,(jjm - 1)*iim + 2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
 #else
-    ps_start_PCM(1) = ps(1)
+    ps_start(1) = ps_start_dyn(1)
 #endif
 
 ! Save initial surface pressure
-ps0 = ps
+ps_start0 = ps_start
 
 ! In the PCM, these values are given to the physic by the dynamic.
 ! Here we simply read them in the "startfi.nc" file
-status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
+status = nf90_open(startfi_name,NF90_NOWRITE,ncid)
 
 status = nf90_inq_varid(ncid,"longitude",lonvarid)
@@ -537 +537 @@
 !------------------------
 ! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value
-call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
-
-allocate(tsoil_avg(ngrid,nsoilmx,nslope))
+call get_timelen_PCM("data_PCM_Y1.nc",timelen)
+
 allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope))
 allocate(vmr_co2_PCM(ngrid,timelen))
-allocate(ps_timeseries(ngrid,timelen))
+allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid),ps_dev(ngrid))
 allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
-allocate(tsurf_avg_yr1(ngrid,nslope))
-allocate(tsurf_avg(ngrid,nslope))
-allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
-allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen))
+allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope),tsurf_avg_old(ngrid,nslope),tsurf_dev(ngrid,nslope))
+allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
+allocate(tsoil_timeseries(ngrid,nsoilmx,nslope,timelen),tsoil_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
 allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
 allocate(q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
 allocate(watersurf_density_avg(ngrid,nslope))
-allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
-allocate(Tsurfavg_before_saved(ngrid,nslope))
-allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
-allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
+allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen),watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
 allocate(delta_co2_adsorbed(ngrid))
 allocate(co2ice_disappeared(ngrid,nslope))
@@ -561 +556 @@
 allocate(vmr_co2_PEM_phys(ngrid,timelen))
 
-write(*,*) "Downloading data Y1..."
-call read_data_PCM("data_PCM_Y1.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), &
-                   tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys,                     &
-                   watersurf_density_avg,watersoil_density_timeseries)
-write(*,*) "Downloading data Y1 done!"
-
-! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value
-write(*,*) "Downloading data Y2..."
-call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,2),min_h2o_ice(:,:,2), &
-                   tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys,                         &
-                   watersurf_density_avg,watersoil_density_timeseries)
-write(*,*) "Downloading data Y2 done!"
+call read_data_PCM("data_PCM_Y1.nc","data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,ps_avg,tsurf_avg_yr1,tsurf_avg, &
+                   tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2_PEM_phys,q_h2o_PEM_phys,watersurf_density_avg,watersoil_density_timeseries)
+
+! Compute the deviation from the average
+ps_dev = ps_start - ps_avg
+tsurf_dev = tsurf - tsurf_avg
+tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:)
 
 !------------------------
@@ -586 +576 @@
 
 if (soil_pem) then
-    do t = 1,timelen
-        tsoil_anom(:,:,:,t) = tsoil_anom(:,:,:,t) - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - <Tsoil> to recompute properly tsoil in the restart
-    enddo
     call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
     tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
-    tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_anom
     watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
+    tsoil_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_timeseries
     do l = nsoilmx + 1,nsoilmx_PEM
         tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
         watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
+        tsoil_PEM_timeseries(:,l,:,:) = tsoil_timeseries(:,nsoilmx,:,:)
     enddo
     watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
 endif !soil_pem
-deallocate(tsoil_avg)
+deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries)
 
 !------------------------
@@ -616 +604 @@
 !------------------------
 ! I   Initialization
-!    I_g Save the initial situation
-!------------------------
-allocate(zplev_PCM(ngrid,nlayer + 1))
-Total_surface = 0.
-do ig = 1,ngrid
-    Total_surface = Total_surface + cell_area(ig)
-    zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
-enddo
-global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
-global_avg_press_PCM = global_avg_press_old
-global_avg_press_new = global_avg_press_old
-write(*,*) "Total surface of the planet =", Total_surface
-write(*,*) "Initial global average pressure =", global_avg_press_PCM
+!    I_g Compute global surface pressure
+!------------------------
+total_surface = sum(cell_area)
+ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
+ps_avg_global_new = ps_avg_global_ini
+write(*,*) "Total surface of the planet =", total_surface
+write(*,*) "Initial global averaged pressure =", ps_avg_global_ini
 
 !------------------------
@@ -635 +617 @@
 !------------------------
 allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
-
 allocate(stratif(ngrid,nslope))
 if (layering_algo) then
@@ -647 +628 @@
 call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
               icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
-              ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice,                       &
-              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,         &
-              delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
+              ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice,watersurf_density_avg,         &
+              watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
+deallocate(tsurf_avg_yr1)
 
 ! We save the places where h2o ice is sublimating
 ! We compute the surface of h2o ice sublimating
-allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
-co2ice_ini_surf = 0.
+allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
+co2ice_sublim_surf_ini = 0.
 h2oice_ini_surf = 0.
-ini_co2ice_sublim = .false.
-ini_h2oice_sublim = .false.
+is_co2ice_sublim_ini = .false.
+is_h2oice_sublim_ini = .false.
 is_co2ice_ini = .false.
 co2ice_disappeared = .false.
@@ -664 +645 @@
         if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
         if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
-            ini_co2ice_sublim(i,islope) = .true.
-            co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
+            is_co2ice_sublim_ini(i,islope) = .true.
+            co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
         endif
         if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
-            ini_h2oice_sublim(i,islope) = .true.
+            is_h2oice_sublim_ini(i,islope) = .true.
             h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
         endif
     enddo
 enddo
-write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf
+write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_sublim_surf_ini
 write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
 
@@ -684 +665 @@
         do islope = 1,nslope
             do l = 1,nsoilmx_PEM - 1
-                if (l==1) then
+                if (l == 1) then
                    totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
@@ -701 +682 @@
     write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbed
 endif ! adsorption
-deallocate(tsurf_avg_yr1)
 
 !------------------------
@@ -741 +721 @@
 ! II.a.1. Compute updated global pressure
     write(*,*) "Recomputing the new pressure..."
+    ps_avg_global_old = ps_avg_global_new
     do i = 1,ngrid
         do islope = 1,nslope
-            global_avg_press_new = global_avg_press_new - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
+            ps_avg_global_new = ps_avg_global_old - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/total_surface
         enddo
     enddo
@@ -749 +730 @@
     if (adsorption_pem) then
         do i = 1,ngrid
-            global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbed(i)/Total_surface
-        enddo
-    endif
-    write(*,*) 'Global average pressure old time step',global_avg_press_old
-    write(*,*) 'Global average pressure new time step',global_avg_press_new
-
-! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
-    allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
-    write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
+            ps_avg_global_new = ps_avg_global_old - g*cell_area(i)*delta_co2_adsorbed(i)/total_surface
+        enddo
+    endif
+    write(*,*) 'Global average pressure old time step',ps_avg_global_old
+    write(*,*) 'Global average pressure new time step',ps_avg_global_new
+
+! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
+    allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
+    write(*,*) "Recomputing the pressure levels timeseries adapted to the old pressure..."
     do l = 1,nlayer + 1
         do ig = 1,ngrid
-            zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
-        enddo
-    enddo
-
-! II.a.3. Surface pressure timeseries
-    write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..."
+            zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
+        enddo
+    enddo
+    write(*,*) "Recomputing the surface pressure timeseries adapted to the new pressure..."
     do ig = 1,ngrid
-        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
-    enddo
-
-! II.a.4. New pressure levels timeseries
-    allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
-    write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
+        ps_timeseries(ig,:) = ps_timeseries(ig,:)*ps_avg_global_new/ps_avg_global_old
+    enddo
+    write(*,*) "Recomputing the pressure levels timeseries adapted to the new pressure..."
+    allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
     do l = 1,nlayer + 1
         do ig = 1,ngrid
-            zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
-        enddo
-    enddo
-
-! II.a.5. Tracers timeseries
+            zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
+        enddo
+    enddo
+
+! II.a.3. Tracers timeseries
     write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..."
-
     l = 1
     do ig = 1,ngrid
         do t = 1,timelen
-            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))
+            ! H2O
+            q_h2o_PEM_phys(ig,t) = q_h2o_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ &
+                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
             if (q_h2o_PEM_phys(ig,t) < 0) then
                 q_h2o_PEM_phys(ig,t) = 1.e-30
@@ -792 +769 @@
                 q_h2o_PEM_phys(ig,t) = 1.
             endif
-        enddo
-    enddo
-
-    do ig = 1,ngrid
-        do t = 1,timelen
-            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))                       &
-                                + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,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))
+            ! CO2
+            q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t))/ &
+                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))                       &
+                                + ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))                       &
+                                -  (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/                     &
+                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
             if (q_co2_PEM_phys(ig,t) < 0) then
                 q_co2_PEM_phys(ig,t) = 1.e-30
@@ -807 +780 @@
                 q_co2_PEM_phys(ig,t) = 1.
             endif
-            mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
+            mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
             vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
         enddo
     enddo
-
-    deallocate(zplev_new_timeseries,zplev_old_timeseries)
+    deallocate(zplev_timeseries_new,zplev_timeseries_old)
 
 !------------------------
@@ -834 +806 @@
 !    II_c Flow of glaciers
 !------------------------
-    if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
-                                                            global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
-    if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
+    if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
+                                                            ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
+    if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow,flag_h2oflow_mesh)
 
 !------------------------
@@ -844 +816 @@
 ! II_d.1 Update Tsurf
     write(*,*) "Updating the new Tsurf"
-    bool_sublim = .false.
-    Tsurfavg_before_saved = tsurf_avg
     do ig = 1,ngrid
         do islope = 1,nslope
-            if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then ! co2 ice disappeared, look for closest point without co2ice
+            ! CO2 ice disappeared so we look for the closest point without CO2 ice
+            if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then
                 co2ice_disappeared(ig,islope) = .true.
                 if (latitude_deg(ig) > 0) then
-                    do ig_loop = ig,ngrid
+                    outer1: do ig_loop = ig,ngrid
                         do islope_loop = islope,iflat,-1
                             if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
                                 tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
-                                bool_sublim = .true.
-                                exit
+                                exit outer1
                             endif
                         enddo
-                        if (bool_sublim) exit
-                    enddo
+                    enddo outer1
                 else
-                    do ig_loop = ig,1,-1
+                    outer2: do ig_loop = ig,1,-1
                         do islope_loop = islope,iflat
                             if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
                                 tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
-                                bool_sublim = .true.
-                                exit
+                                exit outer2
                             endif
                         enddo
-                        if (bool_sublim) exit
-                    enddo
+                    enddo outer2
                 endif
-                if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
+                if (h2o_ice(ig,islope) > frost_albedo_threshold) then
                     albedo(ig,1,islope) = albedo_h2o_frost
                     albedo(ig,2,islope) = albedo_h2o_frost
@@ -881 +848 @@
                     emis(ig,islope) = emissiv
                 endif
-            else if ((co2_ice(ig,islope) > 1.e-3) .and. (d_co2ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
-                ave = 0.
-                do t = 1,timelen
-                    ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
-                enddo
-                tsurf_avg(ig,islope) = ave/timelen
+            else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as tcond CO2
+                call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_ini,ps_avg_global_new,tsurf_avg)
             endif
-            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBLIMATION??? !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        enddo
-    enddo
-
-    do t = 1,timelen
-        tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_avg - Tsurfavg_before_saved
-    enddo
-    ! for the start
-    do ig = 1,ngrid
-        do islope = 1,nslope
-            tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_avg(ig,islope))
         enddo
     enddo
@@ -908 +860 @@
 ! II_d.3 Update soil temperature
         write(*,*)"Updating soil temperature"
-        allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
+        allocate(tsoil_avg_old(ngrid,nsoilmx_PEM))
         do islope = 1,nslope
+            tsoil_avg_old = tsoil_PEM(:,:,islope)
             call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
             call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
 
             do t = 1,timelen
-                Tsoil_locslope(:,1:nsoilmx) = tsoil_PEM(:,1:nsoilmx,islope) + tsoil_anom(:,:,islope,t)
-                Tsoil_locslope(:,nsoilmx + 1:) = tsoil_PEM(:,nsoilmx + 1:,islope)
                 do ig = 1,ngrid
                     do isoil = 1,nsoilmx_PEM
-                        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)
+            ! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
+            tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
+            ! Update of watersoil density
+                        watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/tsoil_PEM_timeseries(ig,isoil,islope,t) + alpha_clap_h2o)/tsoil_PEM_timeseries(ig,isoil,islope,t)*mmol(igcm_h2o_vap)/(mugaz*r)
                         if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
                     enddo
@@ -925 +879 @@
         enddo
         watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
-
+        deallocate(tsoil_avg_old)
         write(*,*) "Update of soil temperature done"
-
-        deallocate(Tsoil_locslope)
 
 ! II_d.4 Update the ice table
@@ -942 +894 @@
             do ig = 1,ngrid
                 do islope = 1,nslope
-                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth)
+                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps_avg(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth)
                     icetable_depth(ig,islope) = ssi_depth
                     ice_porefilling(ig,:,islope) = porefill
@@ -952 +904 @@
 
 ! II_d.5 Update the soil thermal properties
-        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,ps_avg_global_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
 
 ! II_d.6 Update the mass of the regolith adsorbed
         if (adsorption_pem) then
-            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                           &
-                                     h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
+            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
+                                     tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,       &
                                      h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
 
@@ -988 +940 @@
 !    II_e Outputs
 !------------------------
-    call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/global_avg_press_new/))
+    call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
     do islope = 1,nslope
         write(str2(1:2),'(i2.2)') islope
@@ -996 +948 @@
         call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
         call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
-        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
+        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
         if (icetable_equilibrium) then
             call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
@@ -1005 +957 @@
 
         if (soil_pem) then
-            call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
-            call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
+            call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
+            call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
             if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
             if (adsorption_pem) then
@@ -1019 +971 @@
 !    II_f Update the tendencies
 !------------------------
-    call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,global_avg_press_PCM,global_avg_press_new)
+    call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new)
 
 !------------------------
@@ -1026 +978 @@
 !------------------------
     write(*,*) "Checking the stopping criteria..."
-    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
-    call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
+    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
+    call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope)
     i_myear_leg = i_myear_leg + dt
     i_myear = i_myear + dt
@@ -1061 +1013 @@
         write(*,*) "Number of simulated Martian years: i_myear =", i_myear
     endif
-
-    global_avg_press_old = global_avg_press_new
-
 enddo ! big time iteration loop of the pem
 !------------------------------ END RUN --------------------------------
@@ -1072 +1021 @@
 !    III_a Update surface value for the PCM start files
 !------------------------
-! III_a.1 Ice update (for startfi)
+! III_a.1 Ice update for start file
 watercap = 0.
 perennial_co2ice = co2_ice
@@ -1100 +1049 @@
 enddo
 
-! III_a.2 Tsoil update (for startfi)
+! III.a.2. Tsurf update for start file
+tsurf = tsurf_avg + tsurf_dev
+
+! III_a.3 Tsoil update for start file
 if (soil_pem) then
     inertiesoil = TI_PEM(:,:nsoilmx,:)
-    tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen)
+    ! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value
+    do isoil = 1,nsoilmx
+        tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:)
+    enddo
+    tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev
 #ifndef CPP_STD
     flux_geo = fluxgeo
 #endif
 endif
-deallocate(tsoil_anom)
-
-! III_a.4 Pressure (for start)
-ps = ps*global_avg_press_new/global_avg_press_PCM
-ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
-
-! III_a.5 Tracer (for start)
-allocate(zplev_new(ngrid,nlayer + 1))
-
+
+! III_a.4 Pressure update for start file
+ps_start = ps_avg + ps_dev
+
+! III_a.5 Tracers update for start file
+allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
 do l = 1,nlayer + 1
-    zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
+    zplev_start0(:,l) = ap(l) + bp(l)*ps_start0
+    zplev_new(:,l) = ap(l) + bp(l)*ps_start
 enddo
 
@@ -1125 +1079 @@
         do l = 1,llm - 1
             do ig = 1,ngrid
-                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
+                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
             enddo
             q(:,llm,nnq) = q(:,llm - 1,nnq)
@@ -1132 +1086 @@
         do l = 1,llm - 1
             do ig = 1,ngrid
-                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
-                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_PCM(ig,l) - zplev_PCM(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
+                q(ig,l,nnq) = q(ig,l,nnq)*(zplev_start0(ig,l) - zplev_start0(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
+                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_start0(ig,l) - zplev_start0(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
             enddo
             q(:,llm,nnq) = q(:,llm - 1,nnq)
@@ -1140 +1094 @@
 enddo
 
-! Conserving the tracers mass for PCM start files
+! Conserving the tracers mass for start file
 do nnq = 1,nqtot
     do ig = 1,ngrid
@@ -1154 +1108 @@
     enddo
 enddo
+deallocate(zplev_start0,zplev_new)
 
 if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
@@ -1171 +1126 @@
     do l = 1,nlayer
         do i = 1,ip1jmp1
-            teta(i,l)= teta(i,l)*(ps0(i)/ps(i))**kappa
+            teta(i,l)= teta(i,l)*(ps_start0(i)/ps_start(i))**kappa
         enddo
     enddo
     ! Correction on atmospheric pressure
-    call pression(ip1jmp1,ap,bp,ps,p)
+    call pression(ip1jmp1,ap,bp,ps_start,p)
     ! Correction on the mass of atmosphere
     call massdair(p,masse)
     call dynredem0("restart.nc",day_ini,phis)
-    call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps)
+    call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps_start)
     write(*,*) "restart.nc has been written"
 #else
-    call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
+    call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
     write(*,*) "restart1D.txt has been written"
 #endif
@@ -1231 +1186 @@
     deallocate(new_str,new_lag1,new_lag2,current1,current2)
 endif
-deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
-deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
-deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
-deallocate(delta_co2_adsorbed,delta_h2o_adsorbed,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
+deallocate(ps_start,ps_start0,ps_timeseries,ps_avg,ps_dev)
+deallocate(tsurf_avg,tsurf_dev,tsurf_avg_old)
+deallocate(tsoil_PEM,tsoil_dev,tsoil_PEM_timeseries)
+deallocate(vmr_co2_PCM,vmr_co2_PEM_phys,q_co2_PEM_phys,q_h2o_PEM_phys)
+deallocate(watersurf_density_avg,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
+deallocate(delta_co2_adsorbed,delta_h2o_adsorbed,delta_h2o_icetablesublim)
 deallocate(icetable_thickness_old,ice_porefilling_old,zshift_surf,zlag)
-deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
+deallocate(is_co2ice_ini,co2ice_disappeared,is_co2ice_sublim_ini,is_h2oice_sublim_ini,stratif)
 !----------------------------- END OUTPUT ------------------------------
 

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

@@ -8 +8 @@
 
 SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness, &
-                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,d_h2oice,d_co2ice,co2_ice,h2o_ice,    &
-                    global_avg_pressure,watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                    &
-                    m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
+                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,                           &
+                    watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                                      &
+                    d_h2oice,d_co2ice,co2_ice,h2o_ice,m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
 
 use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
-use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
-use comsoil_h,                  only: volcapa, inertiedat
+use comsoil_h_PEM,              only: soil_pem, layer_PEM, fluxgeo, inertiedat_PEM, ini_huge_h2oice, depth_breccia, depth_bedrock, index_breccia, index_bedrock
+use comsoil_h,                  only: inertiedat
 use adsorption_mod,             only: regolith_adsorption, adsorption_pem
 use ice_table_mod,              only: computeice_table_equilibrium, icetable_equilibrium, icetable_dynamic
@@ -22 +22 @@
 use abort_pem_mod,              only: abort_pem
 use compute_soiltemp_mod,       only: ini_tsoil_pem, compute_tsoil_pem
-use comslope_mod,               only: def_slope_mean, subslope_dist
 use layering_mod,               only: layering, array2stratif, nb_str_max, layering_algo
 
@@ -36 +35 @@
 include "callkeys.h"
 
-character(*),                   intent(in) :: filename            ! name of the startfi_PEM.nc
-integer,                        intent(in) :: ngrid               ! # of physical grid points
-integer,                        intent(in) :: nsoil_PCM           ! # of vertical grid points in the PCM
-integer,                        intent(in) :: nsoil_PEM           ! # of vertical grid points in the PEM
-integer,                        intent(in) :: nslope              ! # of sub-grid slopes
-integer,                        intent(in) :: timelen             ! # time samples
-real,                           intent(in) :: timestep            ! time step [s]
-real,                           intent(in) :: global_avg_pressure ! mean average pressure on the planet [Pa]
-real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1       ! surface temperature at the first year of PCM call [K]
-real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2       ! surface temperature at the second year of PCM call [K]
-real, dimension(ngrid,timelen), intent(in) :: q_co2               ! MMR tracer co2 [kg/kg]
-real, dimension(ngrid,timelen), intent(in) :: q_h2o               ! MMR tracer h2o [kg/kg]
-real, dimension(ngrid,timelen), intent(in) :: ps_inst             ! surface pressure [Pa]
-real, dimension(ngrid,nslope),  intent(in) :: d_h2oice        ! tendencies on h2o ice
-real, dimension(ngrid,nslope),  intent(in) :: d_co2ice        ! tendencies on co2 ice
-real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg       ! surface water ice density, yearly averaged (kg/m^3)
+character(*),                   intent(in) :: filename      ! name of the startfi_PEM.nc
+integer,                        intent(in) :: ngrid         ! # of physical grid points
+integer,                        intent(in) :: nsoil_PCM     ! # of vertical grid points in the PCM
+integer,                        intent(in) :: nsoil_PEM     ! # of vertical grid points in the PEM
+integer,                        intent(in) :: nslope        ! # of sub-grid slopes
+integer,                        intent(in) :: timelen       ! # time samples
+real,                           intent(in) :: timestep      ! time step [s]
+real,                           intent(in) :: ps_avg_global ! mean average pressure on the planet [Pa]
+real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr1 ! surface temperature at the first year of PCM call [K]
+real, dimension(ngrid,nslope),  intent(in) :: tsurf_avg_yr2 ! surface temperature at the second year of PCM call [K]
+real, dimension(ngrid,timelen), intent(in) :: q_co2         ! MMR tracer co2 [kg/kg]
+real, dimension(ngrid,timelen), intent(in) :: q_h2o         ! MMR tracer h2o [kg/kg]
+real, dimension(ngrid,timelen), intent(in) :: ps_timeseries ! surface pressure [Pa]
+real, dimension(ngrid,nslope),  intent(in) :: d_h2oice      ! tendencies on h2o ice
+real, dimension(ngrid,nslope),  intent(in) :: d_co2ice      ! tendencies on co2 ice
+real, dimension(ngrid,nslope),  intent(in) :: watersurf_avg ! surface water ice density, yearly averaged (kg/m^3)
 ! outputs
-real, dimension(ngrid),                          intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
-real, dimension(ngrid),                          intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
-real, dimension(ngrid,nslope),                   intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
-real, dimension(ngrid,nslope),                   intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
-type(layering), dimension(ngrid,nslope),         intent(inout) :: stratif             ! stratification (layerings)
-real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
-real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
-real, dimension(ngrid,nslope),                   intent(inout) :: icetable_depth      ! Ice table depth [m]
-real, dimension(ngrid,nslope),                   intent(inout) :: icetable_thickness  ! Ice table thickness [m]
-real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: ice_porefilling     ! Subsurface ice pore filling [m3/m3]
-real, dimension(ngrid,nsoil_PEM,nslope,timelen), intent(inout) :: tsoil_inst          ! instantaneous soil (mid-layer) temperature [K]
-real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
-real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
-real, dimension(ngrid,nsoil_PEM,nslope),         intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
+real, dimension(ngrid),                  intent(out) :: deltam_co2_regolith_phys ! mass of co2 that is exchanged due to adsorption desorption [kg/m^2]
+real, dimension(ngrid),                  intent(out) :: deltam_h2o_regolith_phys ! mass of h2o that is exchanged due to adsorption desorption [kg/m^2]
+real, dimension(ngrid,nslope),           intent(out) :: h2o_ice                  ! h2o ice amount [kg/m^2]
+real, dimension(ngrid,nslope),           intent(out) :: co2_ice                  ! co2 ice amount [kg/m^2]
+type(layering), dimension(ngrid,nslope), intent(inout) :: stratif             ! stratification (layerings)
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: TI_PEM              ! soil (mid-layer) thermal inertia in the PEM grid [SI]
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: tsoil_PEM           ! soil (mid-layer) temperature [K]
+real, dimension(ngrid,nslope),           intent(inout) :: icetable_depth      ! Ice table depth [m]
+real, dimension(ngrid,nslope),           intent(inout) :: icetable_thickness  ! Ice table thickness [m]
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: ice_porefilling     ! Subsurface ice pore filling [m3/m3]
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_co2_regolith_phys ! mass of co2 adsorbed [kg/m^2]
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: m_h2o_regolith_phys ! mass of h2o adsorbed [kg/m^2]
+real, dimension(ngrid,nsoil_PEM,nslope), intent(inout) :: watersoil_avg       ! surface water ice density, yearly averaged (kg/m^3)
 ! local
-real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startPEM                  ! soil temperature saved in the start [K]
-real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM                     ! soil thermal inertia saved in the start [SI]
-logical                                 :: found, found2                   ! check if variables are found in the start
-integer                                 :: iloop, ig, islope, it, isoil, k ! index for loops
-real                                    :: kcond                           ! Thermal conductivity, intermediate variable [SI]
-real                                    :: delta                           ! Depth of the interface regolith-breccia, breccia -bedrock [m]
-character(2)                            :: num                             ! intermediate string to read PEM start sloped variables
-real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr1                   ! intermediate soil temperature during yr 1 [K]
-real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_tmp_yr2                   ! intermediate soil temperature during yr 2 [K]
-logical                                 :: startpem_file                   ! boolean to check if we read the startfile or not
-real, dimension(:,:,:,:), allocatable   :: stratif_array                   ! Array for stratification (layerings)
+real, dimension(ngrid,nsoil_PEM,nslope) :: tsoil_startpem           ! soil temperature saved in the start [K]
+real, dimension(ngrid,nsoil_PEM,nslope) :: TI_startPEM              ! soil thermal inertia saved in the start [SI]
+logical                                 :: found, found2            ! check if variables are found in the start
+integer                                 :: iloop, ig, islope, isoil ! index for loops
+real                                    :: delta                    ! Depth of the interface regolith-breccia, breccia -bedrock [m]
+character(2)                            :: num                      ! intermediate string to read PEM start sloped variables
+logical                                 :: startpem_file            ! boolean to check if we read the startfile or not
+real, dimension(:,:,:,:), allocatable   :: stratif_array            ! Array for stratification (layerings)
 
 #ifdef CPP_STD
@@ -210 +205 @@
             else ! found
                 do iloop = nsoil_PCM + 1,nsoil_PEM
-                    TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
+                    TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope) ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
                 enddo
             endif ! not found
@@ -258 +253 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !2. Soil Temperature
-        do islope=1,nslope
+        do islope = 1,nslope
             write(num,fmt='(i2.2)') islope
-            call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
+            call get_field("tsoil_PEM_slope"//num,tsoil_startpem(:,:,islope),found)
             if (.not. found) then
                 write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
@@ -267 +262 @@
                 call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
             else
-! predictor corrector: restart from year 1 of the PCM and build the evolution of tsoil at depth
-                tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
-                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
-                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_tmp_yr1(:,:,islope))
-                tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)
-                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_tmp_yr2(:,:,islope))
+! predictor-corrector: due to changes of surface temperature in the PCM, the tsoil_startpem is adapted firstly
+!                      for PCM year 1 and then for PCM year 2 in order to build the evolution of the profile at depth
+                call compute_tsoil_pem(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
+                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr1(:,islope),tsoil_startpem(:,:,islope))
+                call compute_tsoil_pem(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg_yr2(:,islope),tsoil_startpem(:,:,islope))
 
                 do iloop = nsoil_PCM + 1,nsoil_PEM
-                    tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
+                    tsoil_PEM(:,iloop,islope) = tsoil_startpem(:,iloop,islope)
                 enddo
             endif !found
 
-            do it = 1,timelen
-                tsoil_inst(:,nsoil_PCM + 1:nsoil_PEM,islope,it) = tsoil_PEM(:,nsoil_PCM + 1:nsoil_PEM,islope)
-            enddo
             do isoil = nsoil_PCM + 1,nsoil_PEM
                 watersoil_avg(:,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(:,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(:,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
@@ -296 +287 @@
                 write(*,*)'will reconstruct the values of the ice table given the current state'
                 call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
-                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
                 do islope = 1,nslope
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -309 +300 @@
                 write(*,*)'PEM settings: failed loading <icetable_depth>'
                 write(*,*)'will reconstruct the values of the ice table given the current state'
-                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+                call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
                 do islope = 1,nslope
                     call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -328 +319 @@
                 endif
             enddo
-            do islope=1,nslope
+            do islope = 1,nslope
                 write(num,fmt='(i2.2)') islope
                 call get_field("mh2o_reg_ads_slope"//num,m_h2o_regolith_phys(:,:,islope),found2)
@@ -337 +328 @@
             enddo
 
-            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
-                                        m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
+            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
+                                     m_h2o_regolith_phys,deltam_h2o_regolith_phys,m_co2_regolith_phys,deltam_co2_regolith_phys)
 
             if (.not. found) deltam_co2_regolith_phys = 0.
@@ -441 +432 @@
 
 ! First raw initialization
-            do it = 1,timelen
-                do isoil = nsoil_PCM + 1,nsoil_PEM
-                    tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
-                enddo
-            enddo
-
-            do it = 1,timelen
-                do isoil = nsoil_PCM + 1,nsoil_PEM
-                    call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_inst(:,:,islope,it))
-                enddo
-            enddo
-
             do isoil = nsoil_PCM + 1,nsoil_PEM
                 do ig = 1,ngrid
@@ -465 +444 @@
         if (icetable_equilibrium) then
             call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_avg,watersoil_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
             do islope = 1,nslope
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -471 +450 @@
             write(*,*) 'PEMETAT0: Ice table done'
         else if (icetable_dynamic) then
-            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,global_avg_pressure,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
+            call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
             do islope = 1,nslope
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
@@ -483 +462 @@
             m_co2_regolith_phys = 0.
             m_h2o_regolith_phys = 0.
-            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_inst,q_co2,q_h2o, &
+            call regolith_adsorption(ngrid,nslope,nsoil_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2,q_h2o, &
                                      m_h2o_regolith_phys,deltam_h2o_regolith_phys, m_co2_regolith_phys,deltam_co2_regolith_phys)
             deltam_co2_regolith_phys = 0.

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

@@ -77 +77 @@
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: tsoil_slope_PEM       ! under mesh bining according to slope
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: inertiesoil_slope_PEM ! under mesh bining according to slope
-real, dimension(ngrid,nslope),           intent(in) :: icetable_depth       ! under mesh bining according to slope
-real, dimension(ngrid,nslope),           intent(in) :: icetable_thickness   ! under mesh bining according to slope
+real, dimension(ngrid,nslope),           intent(in) :: icetable_depth        ! under mesh bining according to slope
+real, dimension(ngrid,nslope),           intent(in) :: icetable_thickness    ! under mesh bining according to slope
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: ice_porefilling       ! under mesh bining according to slope
 real, dimension(ngrid,nsoil_PEM,nslope), intent(in) :: m_co2_regolith, m_h2o_regolith
@@ -119 +119 @@
         call put_field("tsoil_PEM_slope"//num,"Soil temperature by slope type",tsoil_slope_PEM(:,:,islope),Year)
         call put_field("TI_PEM_slope"//num,"Soil Thermal Inertia by slope type",inertiesoil_slope_PEM(:,:,islope),Year)
-        call put_field("mco2_reg_ads_slope"//num, "Mass of co2 adsorbed in the regolith",m_co2_regolith(:,:,islope),Year)
-        call put_field("mh2o_reg_ads_slope"//num, "Mass of h2o adsorbed in the regolith",m_h2o_regolith(:,:,islope),Year)
+        call put_field("mco2_reg_ads_slope"//num, "Mass of CO2 adsorbed in the regolith",m_co2_regolith(:,:,islope),Year)
+        call put_field("mh2o_reg_ads_slope"//num, "Mass of H2O adsorbed in the regolith",m_h2o_regolith(:,:,islope),Year)
     enddo
     call put_field("icetable_depth","Depth of ice table",icetable_depth,Year)

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

@@ -15 +15 @@
 !=======================================================================
 
-SUBROUTINE read_data_PCM(filename,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_PCM_phys,ps_timeseries, &
-                         min_co2_ice,min_h2o_ice,tsurf_avg,tsoil_avg,tsurf_PCM,tsoil_PCM,q_co2,q_h2o,      &
-                         watersurf_density_avg,watersoil_density)
+SUBROUTINE read_data_PCM(filename1,filename2,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_PCM,ps_timeseries,ps_avg,tsurf_avg_yr1,tsurf_avg, &
+                         tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2,q_h2o,watersurf_density_avg,watersoil_density_timeseries)
 use comsoil_h,             only: nsoilmx
 use comsoil_h_PEM,         only: soil_pem
@@ -34 +33 @@
 
 !=======================================================================
-! Arguments:
-character(*), intent(in) :: filename                            ! File name
-integer,      intent(in) :: timelen                             ! number of times stored in the file
-integer,      intent(in) :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
+! Inputs:
+character(*), intent(in) :: filename1, filename2                ! File names
+integer,      intent(in) :: timelen                             ! Number of times stored in the file
+integer,      intent(in) :: iim_input, jjm_input, ngrid, nslope ! Number of points in the lat x lon dynamical grid, number of subgrid slopes
 ! Ouputs
-real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice per slope of the year [kg/m^2]
-real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
-real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_PCM_phys ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
-real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
-real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
-real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
-!SOIL
-real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg             ! Average surface temperature of the concatenated file [K]
-real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_avg             ! Average soil temperature of the concatenated file [K]
-real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_PCM             ! Surface temperature of the concatenated file, time series [K]
-real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_PCM             ! Soil temperature of the concatenated file, time series [K]
-real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_avg ! Water density at the surface [kg/m^3]
-real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
-!=======================================================================
-! Local Variables
-integer                         :: i, j, l, t                                                    ! loop variables
-real                            :: A, B, mmean                                                   ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
-integer                         :: islope                                                        ! loop for variables
-character(2)                    :: num                                                           ! for reading sloped variables
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn             ! h2o ice per slope of the concatenated file [kg/m^2]
+real, dimension(ngrid,nslope,2),               intent(out) :: min_co2_ice                  ! Minimum of co2 ice per slope of the year [kg/m^2]
+real, dimension(ngrid,nslope,2),               intent(out) :: min_h2o_ice                  ! Minimum of h2o ice per slope of the year [kg/m^2]
+real, dimension(ngrid,timelen),                intent(out) :: vmr_co2_PCM                  ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
+real, dimension(ngrid,timelen),                intent(out) :: q_co2                        ! CO2 mass mixing ratio in the first layer [kg/m^3]
+real, dimension(ngrid,timelen),                intent(out) :: q_h2o                        ! H2O mass mixing ratio in the first layer [kg/m^3]
+real, dimension(ngrid,timelen),                intent(out) :: ps_timeseries                ! Surface pressure timeseries [Pa]
+real, dimension(ngrid),                        intent(out) :: ps_avg                       ! Averaged surface pressure [K]
+real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg                    ! Averaged surface temperature [K]
+real, dimension(ngrid,nslope),                 intent(out) :: tsurf_avg_yr1                ! Averaged surface temperature for year 1 [K]
+real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_avg                    ! Averaged soil temperature for year 2 [K]
+real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_timeseries             ! Soil temperature timeseries [K]
+real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_avg        ! Water density at the surface [kg/m^3]
+real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density_timeseries ! Water density timeseries in the soil layer [kg/m^3]
+! Local variables
+integer                                                             :: i, j, l, t, islope        ! Loop variables
+real                                                                :: A, B, mmean               ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
+character(2)                                                        :: num                       ! Hor reading sloped variables
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn             ! H2O ice per slope of the concatenated file [kg/m^2]
 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap
 real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: perennial_co2ice
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_PCM               ! CO2 volume mixing ratio in the first layer [mol/m^3]
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_PCM                    ! Surface Pressure [Pa]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_dyn               ! CO2 volume mixing ratio in the first layer [mol/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_dyn                   ! Surface pressure [Pa]
+real, dimension(iim_input + 1,jjm_input + 1)                        :: ps_avg_dyn                ! Averaged surface pressure of the concatenated file [Pa]
 real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
 real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
-real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_avg_dyn             ! Average surface temperature of the concatenated file [K]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_avg_dyn             ! Average soil temperature of the concatenated file [K]
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_PCM_dyn             ! Surface temperature of the concatenated file, time series [K]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_PCM_dyn             ! Soil temperature of the concatenated file, time series [K]
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_avg_dyn             ! Averaged surface temperature of the concatenated file [K]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_avg_dyn             ! Averaged soil temperature of the concatenated file [K]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_dyn                ! Surface temperature of the concatenated file, time series [K]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_dyn                ! Soil temperature of the concatenated file, time series [K]
 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn                 ! CO2 mass mixing ratio in the first layer [kg/m^3]
 real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn                 ! H2O mass mixing ratio in the first layer [kg/m^3]
@@ -76 +74 @@
 real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn     ! Water density in the soil layer, time series [kg/m^3]
 !-----------------------------------------------------------------------
+
+!------------------ Year 1 ------------------
 ! Open the NetCDF file
-write(*,*) "Opening "//filename//"..."
-call error_msg(NF90_OPEN(filename,NF90_NOWRITE,fID),"open",filename)
-
-! Download the data from the file
-call get_var3("ps",ps_PCM)
-write(*,*) "Data for surface pressure downloaded!"
-
-call get_var3("co2_layer1",q_co2_dyn)
-write(*,*) "Data for vmr co2 downloaded!"
-
-call get_var3("h2o_layer1",q_h2o_dyn)
-write(*,*) "Data for vmr h2o downloaded!"
-
+write(*,*) "Opening "//filename1//"..."
+call error_msg(NF90_OPEN(filename1,NF90_NOWRITE,fID),"open",filename1)
+
+! Compute averages over the year for each point
+write(*,*) "Computing the average of tsurf..."
 if (nslope == 1) then ! There is no slope
     call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
@@ -97 +89 @@
     write(*,*) "Data for h2o_ice_s downloaded!"
 
-    call get_var3("tsurf",tsurf_PCM_dyn(:,:,1,:))
+    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
     write(*,*) "Data for tsurf downloaded!"
 
@@ -106 +98 @@
     call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
     write(*,*) "Data for perennial_co2ice downloaded!"
-
-    if (soil_pem) then
-        call get_var4("soiltemp",tsoil_PCM_dyn(:,:,:,1,:))
-        write(*,*) "Data for soiltemp downloaded!"
-
-        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
-        write(*,*) "Data for waterdensity_soil downloaded!"
-
-        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:))
-        write(*,*) "Data for waterdensity_surface downloaded!"
-    endif !soil_pem
 #endif
 else ! We use slopes
@@ -133 +114 @@
     do islope = 1,nslope
         write(num,'(i2.2)') islope
-        call get_var3("tsurf_slope"//num,tsurf_PCM_dyn(:,:,islope,:))
+        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
     enddo
     write(*,*) "Data for tsurf downloaded!"
@@ -149 +130 @@
     enddo
     write(*,*) "Data for perennial_co2ice downloaded!"
-
-    if (soil_pem) then
-        do islope = 1,nslope
-            write(num,'(i2.2)') islope
-            call get_var4("soiltemp_slope"//num,tsoil_PCM_dyn(:,:,:,islope,:))
-        enddo
-        write(*,*) "Data for soiltemp downloaded!"
-
-        do islope = 1,nslope
-            write(num,'(i2.2)') islope
-            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
-        enddo
-        write(*,*) "Data for waterdensity_soil downloaded!"
-
-        do islope = 1,nslope
-            write(num,'(i2.2)') islope
-            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
-        enddo
-        write(*,*) "Data for waterdensity_surface downloaded!"
-    endif !soil_pem
 #endif
 endif
 
 ! Close the NetCDF file
-write(*,*) "Closing "//filename//"..."
-call error_msg(nf90_close(fID),"close",filename)
+write(*,*) "Closing "//filename1//"..."
+call error_msg(nf90_close(fID),"close",filename1)
 
 ! Compute the minimum over the year for each point
@@ -186 +147 @@
 ! Compute averages over the year for each point
 write(*,*) "Computing the average of tsurf..."
-tsurf_avg_dyn = sum(tsurf_PCM_dyn,4)/timelen
+tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen
+
+#ifndef CPP_1D
+    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
+    do islope = 1,nslope
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope,1))
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,1))
+    enddo
+#else
+    tsurf_avg_yr1(1,:) = tsurf_avg_dyn(1,1,:)
+    min_co2_ice(1,:,1) = min_co2_ice_dyn(1,1,:)
+    min_h2o_ice(1,:,1) = min_h2o_ice_dyn(1,1,:)
+#endif
+write(*,*) "Downloading data Y1 done!"
+
+!------------------ Year 2 ------------------
+
+! Open the NetCDF file
+write(*,*) "Opening "//filename2//"..."
+call error_msg(NF90_OPEN(filename2,NF90_NOWRITE,fID),"open",filename2)
+
+! Download the data from the file
+call get_var3("ps",ps_dyn)
+write(*,*) "Data for surface pressure downloaded!"
+
+call get_var3("co2_layer1",q_co2_dyn)
+write(*,*) "Data for vmr co2 downloaded!"
+
+call get_var3("h2o_layer1",q_h2o_dyn)
+write(*,*) "Data for vmr h2o downloaded!"
+
+if (nslope == 1) then ! There is no slope
+    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
+    write(*,*) "Data for co2_ice downloaded!"
+
+    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
+    write(*,*) "Data for h2o_ice_s downloaded!"
+
+    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
+    write(*,*) "Data for tsurf downloaded!"
+
+#ifndef CPP_STD
+    call get_var3("watercap",watercap(:,:,1,:))
+    write(*,*) "Data for watercap downloaded!"
+
+    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
+    write(*,*) "Data for perennial_co2ice downloaded!"
+
+    if (soil_pem) then
+        call get_var4("soiltemp",tsoil_dyn(:,:,:,1,:))
+        write(*,*) "Data for soiltemp downloaded!"
+
+        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
+        write(*,*) "Data for waterdensity_soil downloaded!"
+
+        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:))
+        write(*,*) "Data for waterdensity_surface downloaded!"
+    endif !soil_pem
+#endif
+else ! We use slopes
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
+    enddo
+    write(*,*) "Data for co2_ice downloaded!"
+
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
+    enddo
+    write(*,*) "Data for h2o_ice_s downloaded!"
+
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
+    enddo
+    write(*,*) "Data for tsurf downloaded!"
+
+#ifndef CPP_STD
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
+    enddo
+    write(*,*) "Data for watercap downloaded!"
+
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
+    enddo
+    write(*,*) "Data for perennial_co2ice downloaded!"
+
+    if (soil_pem) then
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var4("soiltemp_slope"//num,tsoil_dyn(:,:,:,islope,:))
+        enddo
+        write(*,*) "Data for soiltemp downloaded!"
+
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
+        enddo
+        write(*,*) "Data for waterdensity_soil downloaded!"
+
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
+        enddo
+        write(*,*) "Data for waterdensity_surface downloaded!"
+    endif !soil_pem
+#endif
+endif
+
+! Close the NetCDF file
+write(*,*) "Closing "//filename2//"..."
+call error_msg(nf90_close(fID),"close",filename2)
+
+! Compute the minimum over the year for each point
+write(*,*) "Computing the min of h2o_ice..."
+where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
+min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
+write(*,*) "Computing the min of co2_ice..."
+where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
+min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
+
+! Compute averages over the year for each point
+write(*,*) "Computing the average of tsurf..."
+tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen
+
+write(*,*) "Computing the average of Ps..."
+ps_avg_dyn = sum(ps_dyn,3)/timelen
 
 if (soil_pem) then
-    write(*,*) "Computing average of tsoil..."
-    tsoil_avg_dyn = sum(tsoil_PCM_dyn,5)/timelen
-    write(*,*) "Computing average of waterdensity_surface..."
+    write(*,*) "Computing the average of tsoil..."
+    tsoil_avg_dyn = sum(tsoil_dyn,5)/timelen
+    write(*,*) "Computing the average of waterdensity_surface..."
     watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen
 endif
@@ -212 +303 @@
             endif
             mmean = 1/(A*q_co2_dyn(i,j,t) + B)
-            vmr_co2_PCM(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
+            vmr_co2_dyn(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
         enddo
     enddo
@@ -218 +309 @@
 
 #ifndef CPP_1D
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_PCM,vmr_co2_PCM_phys)
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_PCM,ps_timeseries)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_dyn,vmr_co2_PCM)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_dyn,ps_timeseries)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_avg_dyn,ps_avg)
     call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
     call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
-
-    do islope = 1,nslope
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
+    do islope = 1,nslope
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope,2))
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,2))
         if (soil_pem) then
             do l = 1,nsoilmx
                 call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope))
                 do t = 1,timelen
-                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_PCM_dyn(:,:,l,islope,t),tsoil_PCM(:,l,islope,t))
-                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
+                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_dyn(:,:,l,islope,t),tsoil_timeseries(:,l,islope,t))
+                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density_timeseries(:,l,islope,t))
                 enddo
             enddo
             call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope))
         endif ! soil_pem
-        do t = 1,timelen
-            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_PCM_dyn(:,:,islope,t),tsurf_PCM(:,islope,t))
-        enddo
-    enddo
-
+    enddo
     call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
 #else
-    vmr_co2_PCM_phys(1,:) = vmr_co2_PCM(1,1,:)
-    ps_timeseries(1,:) = ps_PCM(1,1,:)
+    vmr_co2_PCM(1,:) = vmr_co2_dyn(1,1,:)
+    ps_timeseries(1,:) = ps_dyn(1,1,:)
+    ps_avg(1) = ps_avg_dyn(1,1)
     q_co2(1,:) = q_co2_dyn(1,1,:)
     q_h2o(1,:) = q_h2o_dyn(1,1,:)
-    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
-    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
+    min_co2_ice(1,:,2) = min_co2_ice_dyn(1,1,:)
+    min_h2o_ice(1,:,2) = min_h2o_ice_dyn(1,1,:)
     if (soil_pem) then
+        tsoil_timeseries(1,:,:,:) = tsoil_dyn(1,1,:,:,:)
         tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:)
-        tsoil_PCM(1,:,:,:) = tsoil_PCM_dyn(1,1,:,:,:)
-        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
+        watersoil_density_timeseries(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
         watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:)
     endif ! soil_pem
-    tsurf_PCM(1,:,:) = tsurf_PCM_dyn(1,1,:,:)
     tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:)
 #endif
+write(*,*) "Downloading data Y2 done!"
 
 END SUBROUTINE read_data_PCM

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

@@ -1 @@
-MODULE recomp_tend_co2_slope_mod
+MODULE recomp_tend_co2_mod
 
 implicit none
@@ -8 +8 @@
 
 SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, &
-                           vmr_co2_PCM,vmr_co2_PEM,ps_PCM_2,global_avg_press_PCM,global_avg_press_new)
+                           vmr_co2_PCM,vmr_co2_PEM,ps_PCM,ps_avg_global_ini,ps_avg_global)
 
 use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol
@@ -20 +20 @@
 !=======================================================================
 
-!   arguments:
-!   ----------
-!   INPUT
+! Inputs
+! ------
 integer,                        intent(in) :: timelen, ngrid, nslope
-real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM          ! physical point field: Volume mixing ratio of co2 in the first layer
-real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM          ! physical point field: Volume mixing ratio of co2 in the first layer
-real, dimension(ngrid,timelen), intent(in) :: ps_PCM_2             ! physical point field: Surface pressure in the PCM
-real,                           intent(in) :: global_avg_press_PCM ! global averaged pressure at previous timestep
-real,                           intent(in) :: global_avg_press_new ! global averaged pressure at current timestep
-real, dimension(ngrid,nslope),  intent(in) :: d_co2ice_ini         ! physical point field: Evolution of perennial ice over one year
-real, dimension(ngrid,nslope),  intent(in) :: co2ice               ! CO2 ice per mesh and sub-grid slope (kg/m^2)
-real, dimension(ngrid,nslope),  intent(in) :: emissivity           ! Emissivity per mesh and sub-grid slope(1)
-!   OUTPUT
+real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM       ! physical point field: Volume mixing ratio of co2 in the first layer
+real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM       ! physical point field: Volume mixing ratio of co2 in the first layer
+real, dimension(ngrid,timelen), intent(in) :: ps_PCM            ! physical point field: Surface pressure in the PCM
+real,                           intent(in) :: ps_avg_global_ini ! global averaged pressure at previous timestep
+real,                           intent(in) :: ps_avg_global     ! global averaged pressure at current timestep
+real, dimension(ngrid,nslope),  intent(in) :: d_co2ice_ini      ! physical point field: Evolution of perennial ice over one year
+real, dimension(ngrid,nslope),  intent(in) :: co2ice            ! CO2 ice per mesh and sub-grid slope (kg/m^2)
+real, dimension(ngrid,nslope),  intent(in) :: emissivity        ! Emissivity per mesh and sub-grid slope(1)
+! Outputs
+! -------
 real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! physical point field: Evolution of perennial ice over one year
 
-!   local:
-!   ------
+! Local:
+! ------
 integer :: i, t, islope
-real    :: coef, ave
+real    :: coef, avg
 
 write(*,*) "Update of the CO2 tendency from the current pressure"
@@ -46 +46 @@
     do islope = 1,nslope
         coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2
-        ave = 0.
+        avg = 0.
         if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
             do t = 1,timelen
-                ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM_2(i,t)/100.)))**4 &
-                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM_2(i,t)*(global_avg_press_new/global_avg_press_PCM)/100.)))**4
+                avg = avg + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM(i,t)/100.)))**4 &
+                      - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM(i,t)*(ps_avg_global/ps_avg_global_ini)/100.)))**4
             enddo
-            if (ave < 1.e-4) ave = 0.
-            d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*ave/timelen
+            if (avg < 1.e-4) avg = 0.
+            d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*avg/timelen
         endif
     enddo
@@ -59 +59 @@
 
 END SUBROUTINE recomp_tend_co2
+!=======================================================================
 
-END MODULE recomp_tend_co2_slope_mod
+SUBROUTINE recomp_tend_h2o(ngrid,nslope,timelen,d_h2oice,PCM_temp,PEM_temp)
+
+implicit none
+
+!=======================================================================
+!
+!  Routine that compute the evolution of the tendencie for h2o ice
+!
+!=======================================================================
+
+! Inputs
+! ------
+integer,            intent(in) :: timelen, ngrid, nslope
+real, dimension(:), intent(in) :: PCM_temp, PEM_temp
+! Outputs
+! -------
+real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! physical point field: Evolution of perennial ice over one year
+
+! Local:
+! ------
+
+write(*,*) "Update of the H2O tendency due to lag layer"
+
+! Flux correction due to lag layer
+!~ Rz_old = h2oice_depth_old*0.0325/4.e-4              ! resistance from PCM
+!~ Rz_new = h2oice_depth_new*0.0325/4.e-4              ! new resistance based on new depth
+!~ R_dec = (1./Rz_old)/(1./Rz_new)                     ! decrease because of resistance
+!~ soil_psv_old = psv(max(PCM_temp(h2oice_depth_old))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the GCM run temperature at the old ice location
+!~ soil_psv_new = psv(max(PEM_temp(h2oice_depth_new))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the PEM run temperature at the new ice location
+!~ hum_dec = soil_psv_old/soil_psv_new                 ! decrease because of lower water vapor pressure at the new depth
+!~ d_h2oice = d_h2oice*R_dec*hum_dec                   ! decrease of flux
+
+END SUBROUTINE recomp_tend_h2o
+
+END MODULE recomp_tend_co2_mod

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

@@ -14 +14 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
+SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
 
 use time_evol_mod, only: h2o_ice_crit
@@ -28 +28 @@
 ! Inputs
 !-------
-integer,                          intent(in) :: ngrid             ! # of physical grid points
-real,    dimension(ngrid),        intent(in) :: cell_area         ! Area of the cells
-real,    dimension(ngrid,nslope), intent(in) :: h2o_ice           ! Actual density of h2o ice
-real,                             intent(in) :: h2oice_ini_surf   ! Initial surface of sublimating h2o ice
-logical, dimension(ngrid,nslope), intent(in) :: ini_h2oice_sublim ! Grid points where h2o ice was initially sublimating
+integer,                          intent(in) :: ngrid                ! # of physical grid points
+real,    dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
+real,    dimension(ngrid,nslope), intent(in) :: h2o_ice              ! Actual density of h2o ice
+real,                             intent(in) :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
+logical, dimension(ngrid,nslope), intent(in) :: is_h2oice_sublim_ini ! Grid points where h2o ice was initially sublimating
 ! Outputs
 !--------
@@ -48 +48 @@
 do i = 1,ngrid
     do islope = 1,nslope
-        if (ini_h2oice_sublim(i,islope) .and. h2o_ice(i,islope) > 0.) h2oice_now_surf = h2oice_now_surf + cell_area(i)*subslope_dist(i,islope)
+        if (is_h2oice_sublim_ini(i,islope) .and. h2o_ice(i,islope) > 0.) h2oice_now_surf = h2oice_now_surf + cell_area(i)*subslope_dist(i,islope)
     enddo
 enddo
@@ -75 +75 @@
 !=======================================================================
 
-SUBROUTINE stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
+SUBROUTINE stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global,nslope)
 
 use time_evol_mod, only: co2_ice_crit, ps_criterion
@@ -89 +89 @@
 ! Inputs
 !-------
-integer,                          intent(in) :: ngrid, nslope        ! # of grid physical grid points
-real,    dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
-real,    dimension(ngrid,nslope), intent(in) :: co2_ice              ! Actual density of co2 ice
-real,                             intent(in) :: co2ice_ini_surf      ! Initial surface of sublimatingco2 ice
-logical, dimension(ngrid,nslope), intent(in) :: ini_co2ice_sublim    ! Grid points where co2 ice was initially sublimating
-real,                             intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
-real,                             intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations
+integer,                          intent(in) :: ngrid, nslope          ! # of grid physical grid points
+real,    dimension(ngrid),        intent(in) :: cell_area              ! Area of the cells
+real,    dimension(ngrid,nslope), intent(in) :: co2_ice                ! Actual density of co2 ice
+real,                             intent(in) :: co2ice_sublim_surf_ini ! Initial surface of sublimatingco2 ice
+logical, dimension(ngrid,nslope), intent(in) :: is_co2ice_sublim_ini   ! Grid points where co2 ice was initially sublimating
+real,                             intent(in) :: ps_avg_global_ini      ! Planet average pressure from the PCM start files
+real,                             intent(in) :: ps_avg_global          ! Planet average pressure from the PEM computations
 ! Outputs
 !--------
@@ -112 +112 @@
 do i = 1,ngrid
     do islope = 1,nslope
-        if (ini_co2ice_sublim(i,islope) .and. co2_ice(i,islope) > 0.) co2ice_now_surf = co2ice_now_surf + cell_area(i)*subslope_dist(i,islope)
+        if (is_co2ice_sublim_ini(i,islope) .and. co2_ice(i,islope) > 0.) co2ice_now_surf = co2ice_now_surf + cell_area(i)*subslope_dist(i,islope)
     enddo
 enddo
 
 ! Check of the criterion
-if (co2ice_now_surf < co2ice_ini_surf*(1. - co2_ice_crit)) then
+if (co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)) then
     stopPEM = 3
     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
-    write(*,*) "co2ice_now_surf < co2ice_ini_surf*(1. - co2_ice_crit)", co2ice_now_surf < co2ice_ini_surf*(1. - co2_ice_crit)
-    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_ini_surf
+    write(*,*) "co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)", co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)
+    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_sublim_surf_ini
     write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf
     write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
-else if (co2ice_now_surf > co2ice_ini_surf*(1. + co2_ice_crit)) then
+else if (co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)) then
     stopPEM = 3
     write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
-    write(*,*) "co2ice_now_surf > co2ice_ini_surf*(1. + co2_ice_crit)", co2ice_now_surf > co2ice_ini_surf*(1. + co2_ice_crit)
+    write(*,*) "co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)", co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)
     write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf
-    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_ini_surf
+    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_sublim_surf_ini
     write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
 endif
 
-if (abs(co2ice_ini_surf) < 1.e-5) stopPEM = 0
+if (abs(co2ice_sublim_surf_ini) < 1.e-5) stopPEM = 0
 
-if (global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)) then
+if (ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)) then
     stopPEM = 4
     write(*,*) "Reason of stopping: the global pressure reaches the threshold"
-    write(*,*) "global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)", global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)
-    write(*,*) "Initial global pressure =", global_avg_press_PCM
-    write(*,*) "Current global pressure =", global_avg_press_new
+    write(*,*) "ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)", ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)
+    write(*,*) "Initial global pressure =", ps_avg_global_ini
+    write(*,*) "Current global pressure =", ps_avg_global
     write(*,*) "Percentage of change accepted =", ps_criterion*100.
-else if (global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)) then
+else if (ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)) then
     stopPEM = 4
     write(*,*) "Reason of stopping: the global pressure reaches the threshold"
-    write(*,*) "global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)", global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)
-    write(*,*) "Initial global pressure =", global_avg_press_PCM
-    write(*,*) "Current global pressure =", global_avg_press_new
+    write(*,*) "ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)", ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)
+    write(*,*) "Initial global pressure =", ps_avg_global_ini
+    write(*,*) "Current global pressure =", ps_avg_global
     write(*,*) "Percentage of change accepted =", ps_criterion*100.
 endif