Ignore:
Timestamp:
Jan 10, 2025, 5:45:03 PM (5 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

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3554 r3571  
    77!    I_e Initialization of the PEM variable and soil
    88!    I_f Compute tendencies
    9 !    I_g Save the initial situation
     9!    I_g Compute global surface pressure
    1010!    I_h Read the "startpem.nc"
    1111!    I_i Compute orbit criterion
    … …  
    4040use pemredem,                   only: pemdem0, pemdem1
    4141use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
    42                                       metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice
     42                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice, computeTcondCO2
    4343use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2
    4444use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
    4545use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
    46 use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, mlayer_PEM, &
    47                                       TI_PEM,               & ! soil thermal inertia
     46use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
     47                                      TI_PEM,               & ! Soil thermal inertia
    4848                                      tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
    4949                                      fluxgeo                 ! Geothermal flux for the PEM and PCM
    5050use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
    5151                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
    52                                       co2_adsorbed_phys, h2o_adsorbed_phys        ! Mass of co2 and h2O adsorbed
     52                                      co2_adsorbed_phys, h2o_adsorbed_phys          ! Mass of co2 and h2O adsorbed
    5353use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
    5454use orbit_param_criterion_mod,  only: orbit_param_criterion
    … …  
    6262use compute_tend_mod,           only: compute_tend
    6363use info_PEM_mod,               only: info_PEM
    64 use nb_time_step_PCM_mod,       only: nb_time_step_PCM
     64use get_timelen_PCM_mod,        only: get_timelen_PCM
    6565use pemetat0_mod,               only: pemetat0
    6666use read_data_PCM_mod,          only: read_data_PCM
    67 use recomp_tend_co2_slope_mod,  only: recomp_tend_co2
     67use recomp_tend_co2_mod,        only: recomp_tend_co2
    6868use compute_soiltemp_mod,       only: compute_tsoil_pem, shift_tsoil2surf
    6969use writediagpem_mod,           only: writediagpem, writediagsoilpem
    … …  
    7777                                  albedodat, zmea, zstd, zsig, zgam, zthe,     &
    7878                                  albedo_h2o_frost,frost_albedo_threshold,     &
    79                                   emissiv, watercaptag, perennial_co2ice, emisice, albedice
     79                                  emissiv, watercaptag, perennial_co2ice
    8080    use dimradmars_mod,     only: totcloudfrac, albedo
    8181    use dust_param_mod,     only: tauscaling
    … …  
    128128real               :: time_phys    ! Same as in PCM
    129129real               :: ptimestep    ! Same as in PCM
    130 real               :: ztime_fin     ! Same as in PCM
    131 
    132 ! Variables to read start.nc
     130real               :: ztime_fin    ! Same as in PCM
     131
     132! Variables to read "start.nc"
    133133character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM
    134134
    … …  
    136136real, dimension(ip1jm,llm)          :: vcov          ! vents covariants
    137137real, dimension(ip1jmp1,llm)        :: ucov          ! vents covariants
    138 real, dimension(ip1jmp1,llm)        :: teta          ! temperature potentielle
     138real, dimension(ip1jmp1,llm)        :: teta          ! Potential temperature
    139139real, dimension(:,:,:), allocatable :: q             ! champs advectes
    140 real, dimension(ip1jmp1)            :: ps            ! pression au sol
    141 real, dimension(ip1jmp1)            :: ps0           ! pression au sol initiale
    142 real, dimension(:),     allocatable :: ps_start_PCM  ! (ngrid) surface pressure
    143 real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure
    144 real, dimension(ip1jmp1,llm)        :: masse         ! masse d'air
     140real, dimension(ip1jmp1)            :: ps_start_dyn  ! surface pressure in the start file (dynamic grid)
     141real, dimension(:),     allocatable :: ps_start      ! surface pressure in the start file
     142real, dimension(:),     allocatable :: ps_start0     ! surface pressure in the start file at the beginning
     143real, dimension(:),     allocatable :: ps_avg        ! (ngrid) Averaged surface pressure
     144real, dimension(:),     allocatable :: ps_dev        ! (ngrid x timelen) Surface pressure deviation
     145real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) Instantaneous surface pressure
     146real, dimension(ip1jmp1,llm)        :: masse         ! Air mass
    145147real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
    146148real                                :: time_0
    … …  
    157159real, dimension(:), allocatable :: latitude      ! Latitude read in startfi_name and written in restartfi
    158160real, dimension(:), allocatable :: cell_area     ! Cell_area read in startfi_name and written in restartfi
    159 real                            :: Total_surface ! Total surface of the planet
     161real                            :: total_surface ! Total surface of the planet
    160162
    161163! Variables for h2o_ice evolution
    162164real, dimension(:,:),    allocatable  :: h2o_ice              ! h2o ice in the PEM
     165real, dimension(:,:),    allocatable  :: d_h2oice             ! physical point x slope field: Tendency of evolution of perennial h2o ice
    163166real, dimension(:,:,:),  allocatable  :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
    164167real                                  :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
    165 logical, dimension(:,:), allocatable  :: ini_h2oice_sublim    ! Logical array to know if h2o ice is sublimating
    166 real                                  :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
    167 real                                  :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
    168 real                                  :: global_avg_press_new ! constant: Global average pressure of current time step
    169 real,   dimension(:,:),   allocatable :: zplev_new            ! Physical x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
    170 real,   dimension(:,:),   allocatable :: zplev_PCM            ! same but retrieved from the PCM [kg/m^2]
    171 real,   dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
    172 real,   dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
     168logical, dimension(:,:), allocatable  :: is_h2oice_sublim_ini ! Logical array to know if h2o ice is sublimating
     169real                                  :: ps_avg_global_ini    ! constant: Global average pressure at initialization [Pa]
     170real                                  :: ps_avg_global_old    ! constant: Global average pressure of previous time step
     171real                                  :: ps_avg_global_new    ! constant: Global average pressure of current time step
     172real,   dimension(:,:),   allocatable :: zplev_new            ! Grid points x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2]
     173real,   dimension(:,:),   allocatable :: zplev_start0         ! Grid points x Atmospheric field: mass of the atmospheric layers in the start [kg/m^2]
     174real,   dimension(:,:,:), allocatable :: zplev_timeseries_new ! Grid points x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
     175real,   dimension(:,:,:), allocatable :: zplev_timeseries_old ! same but with the time series, for previous time step
    173176integer                               :: 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
    174 real, save                            :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
    175 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]
     177real                                  :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
     178real,   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]
    176179integer                               :: timelen              ! # time samples
    177 real                                  :: ave                  ! intermediate varibale to compute average
    178 real,   dimension(:,:),   allocatable :: p                    ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
     180real,   dimension(:,:),   allocatable :: p                    ! Grid points x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
    179181real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
    180182
    181183! Variables for co2_ice evolution
    182 real,    dimension(:,:), allocatable :: co2_ice           ! co2 ice in the PEM
    183 logical, dimension(:,:), allocatable :: is_co2ice_ini     ! Was there co2 ice initially in the PEM?
    184 real,  dimension(:,:,:), allocatable :: min_co2_ice       ! Minimum of co2 ice at each point for the first year [kg/m^2]
    185 real                                 :: co2ice_ini_surf   ! Initial surface of sublimating co2 ice
    186 logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating
    187 real,    dimension(:,:), allocatable :: vmr_co2_PCM       ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
    188 real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys  ! Physics x Times co2 volume mixing ratio used in the PEM
    189 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]
     184real,    dimension(:,:), allocatable :: co2_ice                ! co2 ice in the PEM
     185real,    dimension(:,:), allocatable :: d_co2ice               ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
     186real,    dimension(:,:), allocatable :: d_co2ice_ini           ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
     187logical, dimension(:,:), allocatable :: is_co2ice_ini          ! Was there co2 ice initially in the PEM?
     188real,  dimension(:,:,:), allocatable :: min_co2_ice            ! Minimum of co2 ice at each point for the first year [kg/m^2]
     189real                                 :: co2ice_sublim_surf_ini ! Initial surface of sublimating co2 ice
     190logical, dimension(:,:), allocatable :: is_co2ice_sublim_ini   ! Logical array to know if co2 ice is sublimating
     191real,    dimension(:,:), allocatable :: vmr_co2_PCM            ! Grid points x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
     192real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys       ! Grid points x Times co2 volume mixing ratio used in the PEM
     193real,    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]
    190194
    191195! Variables for stratification (layering) evolution
    … …  
    195199
    196200! Variables for slopes
    197 real, dimension(:,:),   allocatable :: d_co2ice          ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
    198 real, dimension(:,:),   allocatable :: d_co2ice_ini      ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
    199 real, dimension(:,:),   allocatable :: d_h2oice          ! physical point x slope field: Tendency of evolution of perennial h2o ice
    200201real, dimension(:,:),   allocatable :: flag_co2flow      ! (ngrid,nslope): Flag where there is a CO2 glacier flow
    201202real, dimension(:),     allocatable :: flag_co2flow_mesh ! (ngrid)       : Flag where there is a CO2 glacier flow
    … …  
    204205
    205206! Variables for surface and soil
    206 real, dimension(:,:),     allocatable :: tsurf_avg                        ! Physic x SLOPE field: Averaged Surface Temperature [K]
    207 real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
    208 real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries             ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
    209 real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries        ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K]
    210 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]
    211 real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
    212 real, dimension(:,:),     allocatable :: Tsoil_locslope                   ! Physic x Soil: Intermediate when computing Tsoil [K]
    213 real, dimension(:),       allocatable :: Tsurf_locslope                   ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
    214 real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
    215 real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
    216 real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
    217 real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
    218 real, dimension(:,:),     allocatable :: Tsurfavg_before_saved            ! Surface temperature saved from previous time step [K]
     207real, dimension(:,:),     allocatable :: tsurf_avg                        ! Grid points x Slope field: Averaged surface temperature [K]
     208real, dimension(:,:),     allocatable :: tsurf_dev                        ! ngrid x Slope x Times field: Surface temperature deviation [K]
     209real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Grid points x Slope field: Averaged surface temperature of first call of the PCM [K]
     210real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Grid points x Soil x Slope field: Averaged Soil Temperature [K]
     211real, dimension(:,:),     allocatable :: tsoil_avg_old                    ! Grid points x Soil field: Averaged Soil Temperature at the previous time step [K]
     212real, dimension(:,:,:),   allocatable :: tsoil_dev                        ! Grid points x Soil x Slope field: Soil temperature deviation [K]
     213real, dimension(:,:,:,:), allocatable :: tsoil_timeseries                 ! Grid points x Soil x Slope x Times field: Soil temperature timeseries [K]
     214real, dimension(:,:,:,:), allocatable :: tsoil_PEM_timeseries             ! Grid points x Soil x Slope x Times field: Soil temperature timeseries for PEM [K]
     215real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Grid points x Soil x Slope x Times Water soil density timeseries [kg /m^3]
     216real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Grid points x Slope: Averaged  water surface density [kg/m^3]
     217real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Grid points x Soil x Slope x Times: Water soil density timeseries for PEM [kg/m^3]
     218real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Grid points x Soil x Slopes: Averaged water soil density [kg/m^3]
     219real, dimension(:,:),     allocatable :: tsurf_avg_old                    ! Surface temperature saved from previous time step [K]
    219220real, dimension(:),       allocatable :: delta_co2_adsorbed               ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    220221real, dimension(:),       allocatable :: delta_h2o_adsorbed               ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
    221222real                                  :: totmassco2_adsorbed              ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
    222223real                                  :: totmassh2o_adsorbed              ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
    223 logical                               :: bool_sublim                      ! logical to check if there is sublimation or not
    224224logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
    225225real, dimension(:,:),     allocatable :: icetable_thickness_old           ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
    … …  
    233233! Some variables for the PEM run
    234234real, parameter :: year_step = 1   ! Timestep for the pem
    235 real            :: i_myear_leg       ! Number of iteration
     235real            :: i_myear_leg     ! Number of iteration
    236236real            :: n_myear_leg     ! Maximum number of iterations before stopping
    237237real            :: i_myear         ! Global number of Martian years of the chained simulations
    … …  
    249249
    250250#ifdef CPP_STD
    251     real                                :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
    252     real                                :: albedo_h2o_frost              ! albedo of h2o frost
     251    real                                :: frost_albedo_threshold = 0.05 ! Frost albedo threeshold to convert fresh frost to old ice
     252    real                                :: albedo_h2o_frost              ! Albedo of h2o frost
    253253    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
    254254    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
    … …  
    286286
    287287! Loop variables
    288 integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
     288integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil
    289289
    290290! Elapsed time with system clock
    … …  
    379379    endif
    380380
    381     call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,        &
    382                          time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
     381    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,            &
     382                         time_0,ps_start_dyn(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
    383383                         play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
    384     ps(2) = ps(1)
     384    ps_start_dyn(2) = ps_start_dyn(1)
    385385    nsplit_phys = 1
    386386    day_step = steps_per_sol
    … …  
    391391!------------------------
    392392! I   Initialization
    393 !    I_b Read of the "start.nc" and starfi_evol.nc
     393!    I_b Read of the "start.nc" and "starfi.nc"
    394394!------------------------
    395395! I_b.1 Read "start.nc"
    396 allocate(ps_start_PCM(ngrid))
     396allocate(ps_start(ngrid),ps_start0(ngrid))
    397397#ifndef CPP_1D
    398     call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
    399 
    400     call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
     398    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps_start_dyn,phis,time_0)
     399
     400    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps_start_dyn,ps_start)
    401401
    402402    call iniconst !new
    … …  
    412412    status = nf90_close(ncid)
    413413
    414     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)
     414    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)
    415415#else
    416     ps_start_PCM(1) = ps(1)
     416    ps_start(1) = ps_start_dyn(1)
    417417#endif
    418418
    419419! Save initial surface pressure
    420 ps0 = ps
     420ps_start0 = ps_start
    421421
    422422! In the PCM, these values are given to the physic by the dynamic.
    423423! Here we simply read them in the "startfi.nc" file
    424 status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
     424status = nf90_open(startfi_name,NF90_NOWRITE,ncid)
    425425
    426426status = nf90_inq_varid(ncid,"longitude",lonvarid)
    … …  
    537537!------------------------
    538538! 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
    539 call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
    540 
    541 allocate(tsoil_avg(ngrid,nsoilmx,nslope))
     539call get_timelen_PCM("data_PCM_Y1.nc",timelen)
     540
    542541allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope))
    543542allocate(vmr_co2_PCM(ngrid,timelen))
    544 allocate(ps_timeseries(ngrid,timelen))
     543allocate(ps_timeseries(ngrid,timelen),ps_avg(ngrid),ps_dev(ngrid))
    545544allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
    546 allocate(tsurf_avg_yr1(ngrid,nslope))
    547 allocate(tsurf_avg(ngrid,nslope))
    548 allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
    549 allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen))
     545allocate(tsurf_avg_yr1(ngrid,nslope),tsurf_avg(ngrid,nslope),tsurf_avg_old(ngrid,nslope),tsurf_dev(ngrid,nslope))
     546allocate(tsoil_avg(ngrid,nsoilmx,nslope),tsoil_dev(ngrid,nsoilmx,nslope))
     547allocate(tsoil_timeseries(ngrid,nsoilmx,nslope,timelen),tsoil_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    550548allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
    551549allocate(q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
    552550allocate(watersurf_density_avg(ngrid,nslope))
    553 allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
    554 allocate(Tsurfavg_before_saved(ngrid,nslope))
    555 allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    556 allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
     551allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen),watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    557552allocate(delta_co2_adsorbed(ngrid))
    558553allocate(co2ice_disappeared(ngrid,nslope))
    … …  
    561556allocate(vmr_co2_PEM_phys(ngrid,timelen))
    562557
    563 write(*,*) "Downloading data Y1..."
    564 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), &
    565                    tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys,                     &
    566                    watersurf_density_avg,watersoil_density_timeseries)
    567 write(*,*) "Downloading data Y1 done!"
    568 
    569 ! 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
    570 write(*,*) "Downloading data Y2..."
    571 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), &
    572                    tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys,                         &
    573                    watersurf_density_avg,watersoil_density_timeseries)
    574 write(*,*) "Downloading data Y2 done!"
     558call 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, &
     559                   tsoil_avg,tsoil_timeseries,min_co2_ice,min_h2o_ice,q_co2_PEM_phys,q_h2o_PEM_phys,watersurf_density_avg,watersoil_density_timeseries)
     560
     561! Compute the deviation from the average
     562ps_dev = ps_start - ps_avg
     563tsurf_dev = tsurf - tsurf_avg
     564tsoil_dev = tsoil - tsoil_avg(:,1:nsoilmx,:)
    575565
    576566!------------------------
    … …  
    586576
    587577if (soil_pem) then
    588     do t = 1,timelen
    589         tsoil_anom(:,:,:,t) = tsoil_anom(:,:,:,t) - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - <Tsoil> to recompute properly tsoil in the restart
    590     enddo
    591578    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
    592579    tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
    593     tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_anom
    594580    watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
     581    tsoil_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_timeseries
    595582    do l = nsoilmx + 1,nsoilmx_PEM
    596583        tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
    597584        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
     585        tsoil_PEM_timeseries(:,l,:,:) = tsoil_timeseries(:,nsoilmx,:,:)
    598586    enddo
    599587    watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
    600588endif !soil_pem
    601 deallocate(tsoil_avg)
     589deallocate(tsoil_avg,watersoil_density_timeseries,tsoil_timeseries)
    602590
    603591!------------------------
    … …  
    616604!------------------------
    617605! I   Initialization
    618 !    I_g Save the initial situation
    619 !------------------------
    620 allocate(zplev_PCM(ngrid,nlayer + 1))
    621 Total_surface = 0.
    622 do ig = 1,ngrid
    623     Total_surface = Total_surface + cell_area(ig)
    624     zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
    625 enddo
    626 global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
    627 global_avg_press_PCM = global_avg_press_old
    628 global_avg_press_new = global_avg_press_old
    629 write(*,*) "Total surface of the planet =", Total_surface
    630 write(*,*) "Initial global average pressure =", global_avg_press_PCM
     606!    I_g Compute global surface pressure
     607!------------------------
     608total_surface = sum(cell_area)
     609ps_avg_global_ini = sum(cell_area*ps_avg)/total_surface
     610ps_avg_global_new = ps_avg_global_ini
     611write(*,*) "Total surface of the planet =", total_surface
     612write(*,*) "Initial global averaged pressure =", ps_avg_global_ini
    631613
    632614!------------------------
    … …  
    635617!------------------------
    636618allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
    637 
    638619allocate(stratif(ngrid,nslope))
    639620if (layering_algo) then
    … …  
    647628call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
    648629              icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
    649               ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice,                       &
    650               global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,         &
    651               delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
     630              ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice,watersurf_density_avg,         &
     631              watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
     632deallocate(tsurf_avg_yr1)
    652633
    653634! We save the places where h2o ice is sublimating
    654635! We compute the surface of h2o ice sublimating
    655 allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
    656 co2ice_ini_surf = 0.
     636allocate(is_co2ice_sublim_ini(ngrid,nslope),is_h2oice_sublim_ini(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
     637co2ice_sublim_surf_ini = 0.
    657638h2oice_ini_surf = 0.
    658 ini_co2ice_sublim = .false.
    659 ini_h2oice_sublim = .false.
     639is_co2ice_sublim_ini = .false.
     640is_h2oice_sublim_ini = .false.
    660641is_co2ice_ini = .false.
    661642co2ice_disappeared = .false.
    … …  
    664645        if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
    665646        if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
    666             ini_co2ice_sublim(i,islope) = .true.
    667             co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
     647            is_co2ice_sublim_ini(i,islope) = .true.
     648            co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
    668649        endif
    669650        if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
    670             ini_h2oice_sublim(i,islope) = .true.
     651            is_h2oice_sublim_ini(i,islope) = .true.
    671652            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
    672653        endif
    673654    enddo
    674655enddo
    675 write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf
     656write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_sublim_surf_ini
    676657write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
    677658
    … …  
    684665        do islope = 1,nslope
    685666            do l = 1,nsoilmx_PEM - 1
    686                 if (l==1) then
     667                if (l == 1) then
    687668                   totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* &
    688669                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    … …  
    701682    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbed
    702683endif ! adsorption
    703 deallocate(tsurf_avg_yr1)
    704684
    705685!------------------------
    … …  
    741721! II.a.1. Compute updated global pressure
    742722    write(*,*) "Recomputing the new pressure..."
     723    ps_avg_global_old = ps_avg_global_new
    743724    do i = 1,ngrid
    744725        do islope = 1,nslope
    745             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
     726            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
    746727        enddo
    747728    enddo
    … …  
    749730    if (adsorption_pem) then
    750731        do i = 1,ngrid
    751             global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbed(i)/Total_surface
    752         enddo
    753     endif
    754     write(*,*) 'Global average pressure old time step',global_avg_press_old
    755     write(*,*) 'Global average pressure new time step',global_avg_press_new
    756 
    757 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
    758     allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
    759     write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
     732            ps_avg_global_new = ps_avg_global_old - g*cell_area(i)*delta_co2_adsorbed(i)/total_surface
     733        enddo
     734    endif
     735    write(*,*) 'Global average pressure old time step',ps_avg_global_old
     736    write(*,*) 'Global average pressure new time step',ps_avg_global_new
     737
     738! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
     739    allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
     740    write(*,*) "Recomputing the pressure levels timeseries adapted to the old pressure..."
    760741    do l = 1,nlayer + 1
    761742        do ig = 1,ngrid
    762             zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
    763         enddo
    764     enddo
    765 
    766 ! II.a.3. Surface pressure timeseries
    767     write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..."
     743            zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
     744        enddo
     745    enddo
     746    write(*,*) "Recomputing the surface pressure timeseries adapted to the new pressure..."
    768747    do ig = 1,ngrid
    769         ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
    770     enddo
    771 
    772 ! II.a.4. New pressure levels timeseries
    773     allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
    774     write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
     748        ps_timeseries(ig,:) = ps_timeseries(ig,:)*ps_avg_global_new/ps_avg_global_old
     749    enddo
     750    write(*,*) "Recomputing the pressure levels timeseries adapted to the new pressure..."
     751    allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
    775752    do l = 1,nlayer + 1
    776753        do ig = 1,ngrid
    777             zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
    778         enddo
    779     enddo
    780 
    781 ! II.a.5. Tracers timeseries
     754            zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
     755        enddo
     756    enddo
     757
     758! II.a.3. Tracers timeseries
    782759    write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..."
    783 
    784760    l = 1
    785761    do ig = 1,ngrid
    786762        do t = 1,timelen
    787             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))/ &
    788                                    (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
     763            ! H2O
     764            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))/ &
     765                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
    789766            if (q_h2o_PEM_phys(ig,t) < 0) then
    790767                q_h2o_PEM_phys(ig,t) = 1.e-30
    … …  
    792769                q_h2o_PEM_phys(ig,t) = 1.
    793770            endif
    794         enddo
    795     enddo
    796 
    797     do ig = 1,ngrid
    798         do t = 1,timelen
    799             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))/ &
    800                                    (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))                       &
    801                                 + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))                       &
    802                                 -  (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/                     &
    803                                    (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
     771            ! CO2
     772            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))/ &
     773                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))                       &
     774                                + ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))                       &
     775                                -  (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/                     &
     776                                   (zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
    804777            if (q_co2_PEM_phys(ig,t) < 0) then
    805778                q_co2_PEM_phys(ig,t) = 1.e-30
    … …  
    807780                q_co2_PEM_phys(ig,t) = 1.
    808781            endif
    809             mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
     782            mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
    810783            vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
    811784        enddo
    812785    enddo
    813 
    814     deallocate(zplev_new_timeseries,zplev_old_timeseries)
     786    deallocate(zplev_timeseries_new,zplev_timeseries_old)
    815787
    816788!------------------------
    … …  
    834806!    II_c Flow of glaciers
    835807!------------------------
    836     if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
    837                                                             global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
    838     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)
     808    if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
     809                                                            ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
     810    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)
    839811
    840812!------------------------
    … …  
    844816! II_d.1 Update Tsurf
    845817    write(*,*) "Updating the new Tsurf"
    846     bool_sublim = .false.
    847     Tsurfavg_before_saved = tsurf_avg
    848818    do ig = 1,ngrid
    849819        do islope = 1,nslope
    850             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
     820            ! CO2 ice disappeared so we look for the closest point without CO2 ice
     821            if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then
    851822                co2ice_disappeared(ig,islope) = .true.
    852823                if (latitude_deg(ig) > 0) then
    853                     do ig_loop = ig,ngrid
     824                    outer1: do ig_loop = ig,ngrid
    854825                        do islope_loop = islope,iflat,-1
    855826                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
    856827                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
    857                                 bool_sublim = .true.
    858                                 exit
     828                                exit outer1
    859829                            endif
    860830                        enddo
    861                         if (bool_sublim) exit
    862                     enddo
     831                    enddo outer1
    863832                else
    864                     do ig_loop = ig,1,-1
     833                    outer2: do ig_loop = ig,1,-1
    865834                        do islope_loop = islope,iflat
    866835                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
    867836                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
    868                                 bool_sublim = .true.
    869                                 exit
     837                                exit outer2
    870838                            endif
    871839                        enddo
    872                         if (bool_sublim) exit
    873                     enddo
     840                    enddo outer2
    874841                endif
    875                 if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
     842                if (h2o_ice(ig,islope) > frost_albedo_threshold) then
    876843                    albedo(ig,1,islope) = albedo_h2o_frost
    877844                    albedo(ig,2,islope) = albedo_h2o_frost
    … …  
    881848                    emis(ig,islope) = emissiv
    882849                endif
    883             else if ((co2_ice(ig,islope) > 1.e-3) .and. (d_co2ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
    884                 ave = 0.
    885                 do t = 1,timelen
    886                     ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
    887                 enddo
    888                 tsurf_avg(ig,islope) = ave/timelen
     850            else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as tcond CO2
     851                call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_ini,ps_avg_global_new,tsurf_avg)
    889852            endif
    890             !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBLIMATION??? !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    891         enddo
    892     enddo
    893 
    894     do t = 1,timelen
    895         tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_avg - Tsurfavg_before_saved
    896     enddo
    897     ! for the start
    898     do ig = 1,ngrid
    899         do islope = 1,nslope
    900             tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_avg(ig,islope))
    901853        enddo
    902854    enddo
    … …  
    908860! II_d.3 Update soil temperature
    909861        write(*,*)"Updating soil temperature"
    910         allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
     862        allocate(tsoil_avg_old(ngrid,nsoilmx_PEM))
    911863        do islope = 1,nslope
     864            tsoil_avg_old = tsoil_PEM(:,:,islope)
    912865            call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
    913866            call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
    914867
    915868            do t = 1,timelen
    916                 Tsoil_locslope(:,1:nsoilmx) = tsoil_PEM(:,1:nsoilmx,islope) + tsoil_anom(:,:,islope,t)
    917                 Tsoil_locslope(:,nsoilmx + 1:) = tsoil_PEM(:,nsoilmx + 1:,islope)
    918869                do ig = 1,ngrid
    919870                    do isoil = 1,nsoilmx_PEM
    920                         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)
     871            ! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
     872            tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
     873            ! Update of watersoil density
     874                        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)
    921875                        if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
    922876                    enddo
    … …  
    925879        enddo
    926880        watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
    927 
     881        deallocate(tsoil_avg_old)
    928882        write(*,*) "Update of soil temperature done"
    929 
    930         deallocate(Tsoil_locslope)
    931883
    932884! II_d.4 Update the ice table
    … …  
    942894            do ig = 1,ngrid
    943895                do islope = 1,nslope
    944                     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)
     896                    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)
    945897                    icetable_depth(ig,islope) = ssi_depth
    946898                    ice_porefilling(ig,:,islope) = porefill
    … …  
    952904
    953905! II_d.5 Update the soil thermal properties
    954         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)
     906        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)
    955907
    956908! II_d.6 Update the mass of the regolith adsorbed
    957909        if (adsorption_pem) then
    958             call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                          &
    959                                      h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
     910            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
     911                                     tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,      &
    960912                                     h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
    961913
    … …  
    988940!    II_e Outputs
    989941!------------------------
    990     call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/global_avg_press_new/))
     942    call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
    991943    do islope = 1,nslope
    992944        write(str2(1:2),'(i2.2)') islope
    … …  
    996948        call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
    997949        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
    998         call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
     950        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
    999951        if (icetable_equilibrium) then
    1000952            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
    … …  
    1005957
    1006958        if (soil_pem) then
    1007             call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
    1008             call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
     959            call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
     960            call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
    1009961            if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
    1010962            if (adsorption_pem) then
    … …  
    1019971!    II_f Update the tendencies
    1020972!------------------------
    1021     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)
     973    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)
    1022974
    1023975!------------------------
    … …  
    1026978!------------------------
    1027979    write(*,*) "Checking the stopping criteria..."
    1028     call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
    1029     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)
     980    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
     981    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)
    1030982    i_myear_leg = i_myear_leg + dt
    1031983    i_myear = i_myear + dt
    … …  
    10611013        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
    10621014    endif
    1063 
    1064     global_avg_press_old = global_avg_press_new
    1065 
    10661015enddo ! big time iteration loop of the pem
    10671016!------------------------------ END RUN --------------------------------
    … …  
    10721021!    III_a Update surface value for the PCM start files
    10731022!------------------------
    1074 ! III_a.1 Ice update (for startfi)
     1023! III_a.1 Ice update for start file
    10751024watercap = 0.
    10761025perennial_co2ice = co2_ice
    … …  
    11001049enddo
    11011050
    1102 ! III_a.2 Tsoil update (for startfi)
     1051! III.a.2. Tsurf update for start file
     1052tsurf = tsurf_avg + tsurf_dev
     1053
     1054! III_a.3 Tsoil update for start file
    11031055if (soil_pem) then
    11041056    inertiesoil = TI_PEM(:,:nsoilmx,:)
    1105     tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen)
     1057    ! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value
     1058    do isoil = 1,nsoilmx
     1059        tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:)
     1060    enddo
     1061    tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev
    11061062#ifndef CPP_STD
    11071063    flux_geo = fluxgeo
    11081064#endif
    11091065endif
    1110 deallocate(tsoil_anom)
    1111 
    1112 ! III_a.4 Pressure (for start)
    1113 ps = ps*global_avg_press_new/global_avg_press_PCM
    1114 ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
    1115 
    1116 ! III_a.5 Tracer (for start)
    1117 allocate(zplev_new(ngrid,nlayer + 1))
    1118 
     1066
     1067! III_a.4 Pressure update for start file
     1068ps_start = ps_avg + ps_dev
     1069
     1070! III_a.5 Tracers update for start file
     1071allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
    11191072do l = 1,nlayer + 1
    1120     zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
     1073    zplev_start0(:,l) = ap(l) + bp(l)*ps_start0
     1074    zplev_new(:,l) = ap(l) + bp(l)*ps_start
    11211075enddo
    11221076
    … …  
    11251079        do l = 1,llm - 1
    11261080            do ig = 1,ngrid
    1127                 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))
     1081                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))
    11281082            enddo
    11291083            q(:,llm,nnq) = q(:,llm - 1,nnq)
    … …  
    11321086        do l = 1,llm - 1
    11331087            do ig = 1,ngrid
    1134                 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)) &
    1135                               + ((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))
     1088                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)) &
     1089                              + ((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))
    11361090            enddo
    11371091            q(:,llm,nnq) = q(:,llm - 1,nnq)
    … …  
    11401094enddo
    11411095
    1142 ! Conserving the tracers mass for PCM start files
     1096! Conserving the tracers mass for start file
    11431097do nnq = 1,nqtot
    11441098    do ig = 1,ngrid
    … …  
    11541108    enddo
    11551109enddo
     1110deallocate(zplev_start0,zplev_new)
    11561111
    11571112if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
    … …  
    11711126    do l = 1,nlayer
    11721127        do i = 1,ip1jmp1
    1173             teta(i,l)= teta(i,l)*(ps0(i)/ps(i))**kappa
     1128            teta(i,l)= teta(i,l)*(ps_start0(i)/ps_start(i))**kappa
    11741129        enddo
    11751130    enddo
    11761131    ! Correction on atmospheric pressure
    1177     call pression(ip1jmp1,ap,bp,ps,p)
     1132    call pression(ip1jmp1,ap,bp,ps_start,p)
    11781133    ! Correction on the mass of atmosphere
    11791134    call massdair(p,masse)
    11801135    call dynredem0("restart.nc",day_ini,phis)
    1181     call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps)
     1136    call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps_start)
    11821137    write(*,*) "restart.nc has been written"
    11831138#else
    1184     call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
     1139    call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
    11851140    write(*,*) "restart1D.txt has been written"
    11861141#endif
    … …  
    12311186    deallocate(new_str,new_lag1,new_lag2,current1,current2)
    12321187endif
    1233 deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
    1234 deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
    1235 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
    1236 deallocate(delta_co2_adsorbed,delta_h2o_adsorbed,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
     1188deallocate(ps_start,ps_start0,ps_timeseries,ps_avg,ps_dev)
     1189deallocate(tsurf_avg,tsurf_dev,tsurf_avg_old)
     1190deallocate(tsoil_PEM,tsoil_dev,tsoil_PEM_timeseries)
     1191deallocate(vmr_co2_PCM,vmr_co2_PEM_phys,q_co2_PEM_phys,q_h2o_PEM_phys)
     1192deallocate(watersurf_density_avg,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
     1193deallocate(delta_co2_adsorbed,delta_h2o_adsorbed,delta_h2o_icetablesublim)
    12371194deallocate(icetable_thickness_old,ice_porefilling_old,zshift_surf,zlag)
    1238 deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
     1195deallocate(is_co2ice_ini,co2ice_disappeared,is_co2ice_sublim_ini,is_h2oice_sublim_ini,stratif)
    12391196!----------------------------- END OUTPUT ------------------------------
    12401197
Note: See TracChangeset for help on using the changeset viewer.