source: trunk/LMDZ.COMMON/libf/evolution/pem.F90 @ 3486

Last change on this file since 3486 was 3486, checked in by evos, 4 weeks ago

we added the option to use NS dynamical subsurface ice in the model to more realisticly calculate the amount of ice in the subsurface and therfore the subsurface thermal inertia
File size: 61.7 KB
RevLine 
[2779]1!------------------------
[3028]2! I   Initialization
[3161]3!    I_a Read the "run.def"
[3317]4!    I_b Read the "start.nc" and "startfi.nc"
[2835]5!    I_c Subslope parametrisation
[3161]6!    I_d Read the PCM data and convert them to the physical grid
[3028]7!    I_e Initialization of the PEM variable and soil
[3143]8!    I_f Compute tendencies
[3384]9!    I_g Save the initial situation
10!    I_h Read the "startpem.nc"
[2835]11!    I_i Compute orbit criterion
[2779]12
13! II  Run
[3028]14!    II_a Update pressure, ice and tracers
[3149]15!    II_b Evolution of ice
16!    II_c Flow of glaciers
[2835]17!    II_d Update surface and soil temperatures
[3088]18!    II_e Outputs
19!    II_f Update the tendencies
20!    II_g Checking the stopping criterion
[2779]21
22! III Output
[2835]23!    III_a Update surface value for the PCM start files
[3317]24!    III_b Write the "restart.nc" and "restartfi.nc"
[3161]25!    III_c Write the "restartpem.nc"
[2779]26!------------------------
27
28PROGRAM pem
29
[3161]30use phyetat0_mod,               only: phyetat0
31use phyredem,                   only: physdem0, physdem1
32use netcdf,                     only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
33use turb_mod,                   only: q2, wstar
[3206]34use comslope_mod,               only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h
[3161]35use logic_mod,                  only: iflag_phys
36use mod_const_mpi,              only: COMM_LMDZ
[3028]37use infotrac
[3161]38use geometry_mod,               only: latitude_deg
39use conf_pem_mod,               only: conf_pem
40use pemredem,                   only: pemdem0, pemdem1
41use glaciers_mod,               only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, &
42                                      metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice
43use stopping_crit_mod,          only: stopping_crit_h2o_ice, stopping_crit_co2
44use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
45use evol_ice_mod,               only: evol_co2_ice, evol_h2o_ice
46use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
[3330]47                                      TI_PEM,               & ! soil thermal inertia
48                                      tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
49                                      fluxgeo                 ! Geothermal flux for the PEM and PCM
[3161]50use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
51                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
52                                      co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
53use time_evol_mod,              only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
54use orbit_param_criterion_mod,  only: orbit_param_criterion
55use recomp_orb_param_mod,       only: recomp_orb_param
56use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
[3486]57                                      ini_ice_table_porefilling, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium,compute_massh2o_exchange_ssi
[3161]58use soil_thermalproperties_mod, only: update_soil_thermalproperties
[3206]59use time_phylmdz_mod,           only: daysec, dtphys
[3161]60use abort_pem_mod,              only: abort_pem
61use soil_settings_PEM_mod,      only: soil_settings_PEM
62use compute_tend_mod,           only: compute_tend
63use info_PEM_mod,               only: info_PEM
64use interpol_TI_PEM2PCM_mod,    only: interpol_TI_PEM2PCM
65use nb_time_step_PCM_mod,       only: nb_time_step_PCM
66use pemetat0_mod,               only: pemetat0
67use read_data_PCM_mod,          only: read_data_PCM
68use recomp_tend_co2_slope_mod,  only: recomp_tend_co2_slope
[3178]69use compute_soiltemp_mod,       only: compute_tsoil_pem
[3181]70use writediagpem_mod,           only: writediagpem, writediagsoilpem
[3207]71use co2condens_mod,             only: CO2cond_ps
[3319]72use layering_mod,               only: tend_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo
[2985]73
[2842]74#ifndef CPP_STD
[3206]75    use comsoil_h,          only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil
[3028]76    use surfdat_h,          only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
77                                  albedodat, zmea, zstd, zsig, zgam, zthe,     &
[3206]78                                  albedo_h2o_frost,frost_albedo_threshold,     &
79                                  emissiv, watercaptag, perennial_co2ice, emisice, albedice
[3028]80    use dimradmars_mod,     only: totcloudfrac, albedo
81    use dust_param_mod,     only: tauscaling
[3149]82    use tracer_mod,         only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
[3028]83    use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
[3096]84    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
[3442]85    use comcstfi_h,         only: pi, rad, g, mugaz, r
[3143]86    use surfini_mod,        only: surfini
[3442]87    use comconst_mod,       only: kappa, cpp
[2842]88#else
[3028]89    use tracer_h,           only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
90    use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT,   &
91                                  PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
92                                  ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
[3039]93    use aerosol_mod,        only: iniaerosol
94    use planete_mod,        only: apoastr, periastr, year_day, peri_day, obliquit
[3442]95    use comcstfi_mod,       only: pi, rad, g, mugaz, r
[2842]96#endif
[2985]97
[3028]98#ifndef CPP_1D
[3076]99    use iniphysiq_mod,            only: iniphysiq
100    use control_mod,              only: iphysiq, day_step, nsplit_phys
[3019]101#else
[3386]102    use time_phylmdz_mod,         only: iphysiq, steps_per_sol
[3028]103    use regular_lonlat_mod,       only: init_regular_lonlat
104    use physics_distribution_mod, only: init_physics_distribution
105    use mod_grid_phy_lmdz,        only: regular_lonlat
[3065]106    use init_testphys1d_mod,      only: init_testphys1d
107    use comvert_mod,              only: ap, bp
[3076]108    use writerestart1D_mod,       only: writerestart1D
[2980]109#endif
[2835]110
[3076]111implicit none
[2980]112
[3028]113include "dimensions.h"
114include "paramet.h"
115include "comgeom.h"
116include "iniprint.h"
[3039]117include "callkeys.h"
[2779]118
[3028]119integer ngridmx
120parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
[2794]121
[3096]122! Same variable names as in the PCM
[3065]123integer, parameter :: nlayer = llm ! Number of vertical layer
124integer            :: ngrid        ! Number of physical grid points
125integer            :: nq           ! Number of tracer
126integer            :: day_ini      ! First day of the simulation
127real               :: pday         ! Physical day
[3149]128real               :: time_phys    ! Same as in PCM
129real               :: ptimestep    ! Same as in PCM
[3297]130real               :: ztime_fin     ! Same as in PCM
[2794]131
[3028]132! Variables to read start.nc
[3317]133character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM
[2779]134
[3028]135! Dynamic variables
[3065]136real, dimension(ip1jm,llm)          :: vcov          ! vents covariants
137real, dimension(ip1jmp1,llm)        :: ucov          ! vents covariants
138real, dimension(ip1jmp1,llm)        :: teta          ! temperature potentielle
139real, dimension(:,:,:), allocatable :: q             ! champs advectes
[3143]140real, dimension(ip1jmp1)            :: ps            ! pression au sol
[3442]141real, dimension(ip1jmp1)            :: ps0           ! pression au sol initiale
[3149]142real, dimension(:),     allocatable :: ps_start_PCM  ! (ngrid) surface pressure
143real, dimension(:,:),   allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure
[3065]144real, dimension(ip1jmp1,llm)        :: masse         ! masse d'air
145real, dimension(ip1jmp1)            :: phis          ! geopotentiel au sol
[3028]146real                                :: time_0
[2779]147
[3028]148! Variables to read starfi.nc
[3317]149character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM
[3143]150character(2)            :: str2
[3206]151integer                 :: ncid, status                           ! Variable for handling opening of files
152integer                 :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files
153integer                 :: apvarid, bpvarid                       ! Variable ID for Netcdf files
[2794]154
[3028]155! Variables to read starfi.nc and write restartfi.nc
[3143]156real, dimension(:), allocatable :: longitude     ! Longitude read in startfi_name and written in restartfi
157real, dimension(:), allocatable :: latitude      ! Latitude read in startfi_name and written in restartfi
158real, dimension(:), allocatable :: cell_area     ! Cell_area read in startfi_name and written in restartfi
[3028]159real                            :: Total_surface ! Total surface of the planet
[2897]160
[3028]161! Variables for h2o_ice evolution
[3327]162real, dimension(:,:),    allocatable  :: h2o_ice              ! h2o ice in the PEM
163real, dimension(:,:,:),  allocatable  :: min_h2o_ice          ! Minima of h2o ice at each point for the PCM years [kg/m^2]
164real                                  :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
165logical, dimension(:,:), allocatable  :: ini_h2oice_sublim    ! Logical array to know if h2o ice is sublimating
166real                                  :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa]
167real                                  :: global_avg_press_old ! constant: Global average pressure of initial/previous time step
168real                                  :: global_avg_press_new ! constant: Global average pressure of current time step
169real,   dimension(:,:),   allocatable :: zplev_new            ! Physical x Atmospheric field: mass of the atmospheric  layers in the pem at current time step [kg/m^2]
170real,   dimension(:,:),   allocatable :: zplev_PCM            ! same but retrieved from the PCM [kg/m^2]
171real,   dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
172real,   dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
173integer                               :: 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
174real, save                            :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
175real,   dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
176integer                               :: timelen              ! # time samples
177real                                  :: ave                  ! intermediate varibale to compute average
178real,   dimension(:,:),   allocatable :: p                    ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
179real                                  :: extra_mass           ! Intermediate variables Extra mass of a tracer if it is greater than 1
[2779]180
[3130]181! Variables for co2_ice evolution
[3330]182real,    dimension(:,:), allocatable :: co2_ice           ! co2 ice in the PEM
183logical, dimension(:,:), allocatable :: is_co2ice_ini     ! Was there co2 ice initially in the PEM?
[3327]184real,  dimension(:,:,:), allocatable :: min_co2_ice       ! Minimum of co2 ice at each point for the first year [kg/m^2]
185real                                 :: co2ice_ini_surf   ! Initial surface of sublimating co2 ice
186logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating
[3330]187real,    dimension(:,:), allocatable :: vmr_co2_PCM       ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
188real,    dimension(:,:), allocatable :: vmr_co2_PEM_phys  ! Physics x Times co2 volume mixing ratio used in the PEM
189real,    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]
[3130]190
[3297]191! Variables for stratification (layering) evolution
192type(layering), dimension(:,:), allocatable :: stratif                     ! Layering (linked list of "stratum") for each grid point and slope
193type(ptrarray), dimension(:,:), allocatable :: current1, current2          ! Current active stratum in the layering
194logical,        dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm
195
[3028]196! Variables for slopes
[3149]197real, dimension(:,:),   allocatable :: tend_co2_ice       ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
[3339]198real, dimension(:,:),   allocatable :: tend_co2_ice_ini   ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
199real, dimension(:,:),   allocatable :: tend_h2o_ice       ! physical point x slope field: Tendency of evolution of perennial h2o ice
[3149]200real, dimension(:,:),   allocatable :: flag_co2flow       ! (ngrid,nslope): Flag where there is a CO2 glacier flow
201real, dimension(:),     allocatable :: flag_co2flow_mesh  ! (ngrid)       : Flag where there is a CO2 glacier flow
202real, dimension(:,:),   allocatable :: flag_h2oflow       ! (ngrid,nslope): Flag where there is a H2O glacier flow
203real, dimension(:),     allocatable :: flag_h2oflow_mesh  ! (ngrid)       : Flag where there is a H2O glacier flow
[2779]204
[3028]205! Variables for surface and soil
[3367]206real, dimension(:,:),     allocatable :: tsurf_avg                          ! Physic x SLOPE field: Averaged Surface Temperature [K]
207real, dimension(:,:,:),   allocatable :: tsoil_avg                          ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
[3339]208real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries               ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
[3426]209real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K]
210real, dimension(:,:,:,:), allocatable :: tsoil_anom                         ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K]
[3149]211real, dimension(:,:),     allocatable :: tsurf_avg_yr1                      ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
[3339]212real, dimension(:,:),     allocatable :: Tsoil_locslope                     ! Physic x Soil: Intermediate when computing Tsoil [K]
[3149]213real, dimension(:),       allocatable :: Tsurf_locslope                     ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
[3065]214real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries       ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
[3367]215real, dimension(:,:),     allocatable :: watersurf_density_avg              ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
[3065]216real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries   ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
[3367]217real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg          ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
[3149]218real, dimension(:,:),     allocatable :: Tsurfavg_before_saved              ! Surface temperature saved from previous time step [K]
[3065]219real, dimension(:),       allocatable :: delta_co2_adsorbded                ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
220real, dimension(:),       allocatable :: delta_h2o_adsorbded                ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
221real                                  :: totmassco2_adsorbded               ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
222real                                  :: totmassh2o_adsorbded               ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
223logical                               :: bool_sublim                        ! logical to check if there is sublimation or not
[3330]224logical, dimension(:,:),  allocatable :: co2ice_disappeared                 ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
[3065]225real, dimension(:,:),     allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
[3181]226real, dimension(:),       allocatable :: delta_h2o_icetablesublim           ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
[3065]227
[3486]228! Dynamic ice table
229real, dimension(:,:),     allocatable :: ss_ice_pf                          ! the amout of porefilling in each layer in each grid [m^3/m^3]
230real, dimension(:,:),     allocatable :: ss_ice_pf_new                      ! the amout of porefilling in each layer in each grid after the ss module[m^3/m^3]
231real, dimension(:,:),     allocatable :: porefillingice_depth_new           ! new pure ice table depth
[3028]232! Some variables for the PEM run
[3363]233real, parameter :: year_step = 1   ! Timestep for the pem
234integer         :: year_iter       ! Number of iteration
235integer         :: year_iter_max   ! Maximum number of iterations before stopping
236integer         :: i_myear         ! Global number of Martian years of the chained simulations
237integer         :: n_myear         ! Maximum number of Martian years of the chained simulations
238real            :: timestep        ! Timestep [s]
239character(20)   :: job_id          ! Job id provided as argument passed on the command line when the program was invoked
[3394]240logical         :: timewall        ! Flag to use the time limit stopping criterion in case of a PEM job
[3363]241integer(kind=8) :: cr              ! Number of clock ticks per second (count rate)
242integer(kind=8) :: c1, c2          ! Counts of processor clock
243character(100)  :: chtimelimit     ! Time limit for the PEM job outputted by the SLURM command
244real            :: timelimit       ! Time limit for the PEM job in seconds
[3365]245real, parameter :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default
[3363]246integer         :: cstat, days, hours, minutes, seconds
247character(1)    :: sep
[2779]248
[2842]249#ifdef CPP_STD
[3065]250    real                                :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
251    real                                :: albedo_h2o_frost              ! albedo of h2o frost
[3143]252    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
253    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
254    real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
255    real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimension when reading form generic
256    real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimension when reading form generic
[3065]257    real, dimension(:,:),   allocatable :: tsurf                         ! Subslope variable, only needed in the GENERIC case
258    real, dimension(:,:,:), allocatable :: qsurf                         ! Subslope variable, only needed in the GENERIC case
259    real, dimension(:,:,:), allocatable :: tsoil                         ! Subslope variable, only needed in the GENERIC case
260    real, dimension(:,:),   allocatable :: emis                          ! Subslope variable, only needed in the GENERIC case
261    real, dimension(:,:),   allocatable :: watercap                      ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
[3068]262    logical, dimension(:),  allocatable :: watercaptag                   ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
[3065]263    real, dimension(:,:,:), allocatable :: albedo                        ! Subslope variable, only needed in the GENERIC case
264    real, dimension(:,:,:), allocatable :: inertiesoil                   ! Subslope variable, only needed in the GENERIC case
[2842]265#endif
266
[2980]267#ifdef CPP_1D
[3143]268    integer            :: nsplit_phys
269    integer, parameter :: jjm_value = jjm - 1
[3386]270    integer            :: day_step
[3065]271
272    ! Dummy variables to use the subroutine 'init_testphys1d'
[3129]273    logical                             :: therestart1D, therestartfi
[3068]274    integer                             :: ndt, day0
275    real                                :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau
276    real, dimension(:),     allocatable :: zqsat
277    real, dimension(:,:,:), allocatable :: dq, dqdyn
278    real, dimension(nlayer)             :: play, w
279    real, dimension(nlayer + 1)         :: plev
[2980]280#else
[3143]281    integer, parameter              :: jjm_value = jjm
282    real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
283    real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
[2980]284#endif
285
[3028]286! Loop variables
[3206]287integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
[2779]288
[3363]289! Elapsed time with system clock
290call system_clock(count_rate = cr)
291call system_clock(c1)
[3394]292timewall = .true.
[3363]293timelimit = 86400 ! 86400 seconds = 24 h by default
294if (command_argument_count() > 0) then
295    ! Read the job id passed as argument to the program
296    call get_command_argument(1,job_id)
297    ! Execute the system command
298    call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
[3403]299    if (cstat /= 0) then
300        call execute_command_line('qstat -f '//trim(job_id)//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt', cmdstat = cstat)
301        if (cstat > 0) then
302            error stop 'pem: command execution failed!'
303        else if (cstat < 0) then
304            error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!'
305        endif
[3363]306    endif
307    ! Read the output
308    open(1,file = 'tmp_cmdout.txt',status = 'old')
309    read(1,'(a)') chtimelimit
310    close(1)
311    chtimelimit = trim(chtimelimit)
312    call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat)
313    if (cstat > 0) then
314        error stop 'pem: command execution failed!'
315    else if (cstat < 0) then
316        error stop 'pem: command execution not supported!'
317    endif
[3387]318    if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS"
[3388]319        read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds
[3363]320        timelimit = days*86400 + hours*3600 + minutes*60 + seconds
[3387]321    else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS"
[3388]322        read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds
[3363]323        timelimit = hours*3600 + minutes*60 + seconds
[3387]324    else ! 'chtimelimit' format is "MM:SS"
325        read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds
326        timelimit = minutes*60 + seconds
[3363]327    endif
[3394]328else
329    timewall = .false.
[3363]330endif
331
[3028]332! Parallel variables
[2842]333#ifndef CPP_STD
[3028]334    is_sequential = .true.
335    is_parallel = .false.
336    is_mpi_root = .true.
337    is_omp_root = .true.
338    is_master = .true.
[2842]339#endif
[2779]340
[3065]341! Some constants
[3028]342day_ini = 0    ! test
343time_phys = 0. ! test
344ngrid = ngridmx
345A = (1/m_co2 - 1/m_noco2)
346B = 1/m_noco2
347year_day = 669
348daysec = 88775.
[3426]349timestep = year_day*daysec*year_step
[2794]350
[3028]351!----------------------------- INITIALIZATION --------------------------
[2779]352!------------------------
[3028]353! I   Initialization
[3161]354!    I_a Read the "run.def"
[2779]355!------------------------
[2980]356#ifndef CPP_1D
[3383]357    dtphys = daysec/48. ! Dummy value (overwritten in phyetat0)
[3028]358    call conf_gcm(99,.true.)
359    call infotrac_init
360    nq = nqtot
361    allocate(q(ip1jmp1,llm,nqtot))
[3065]362    allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
[2980]363#else
[3068]364    allocate(q(1,llm,nqtot))
[3065]365    allocate(longitude(1),latitude(1),cell_area(1))
[3129]366
[3143]367    therestart1D = .false. ! Default value
[3363]368    inquire(file = 'start1D.txt',exist = therestart1D)
[3129]369    if (.not. therestart1D) then
[3363]370        write(*,*) 'There is no "start1D.txt" file!'
[3129]371        error stop 'Initialization cannot be done for the 1D PEM.'
372    endif
[3143]373    therestartfi = .false. ! Default value
[3317]374    inquire(file = 'startfi.nc',exist = therestartfi)
[3129]375    if (.not. therestartfi) then
[3317]376        write(*,*) 'There is no "startfi.nc" file!'
[3129]377        error stop 'Initialization cannot be done for the 1D PEM.'
378    endif
379
[3363]380    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,        &
381                         time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
[3207]382                         play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
[3065]383    ps(2) = ps(1)
[3028]384    nsplit_phys = 1
[3399]385    day_step = steps_per_sol
[2980]386#endif
[2779]387
[3039]388call conf_pem(i_myear,n_myear)
[2779]389
[2835]390!------------------------
[3028]391! I   Initialization
[3317]392!    I_b Read of the "start.nc" and starfi_evol.nc
[3028]393!------------------------
[3317]394! I_b.1 Read "start.nc"
[3149]395allocate(ps_start_PCM(ngrid))
[2980]396#ifndef CPP_1D
[3143]397    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
[2779]398
[3149]399    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
[2897]400
[3028]401    call iniconst !new
402    call inigeom
[2980]403
[3028]404    allocate(ap(nlayer + 1))
405    allocate(bp(nlayer + 1))
[3143]406    status = nf90_open(start_name,NF90_NOWRITE,ncid)
[3028]407    status = nf90_inq_varid(ncid,"ap",apvarid)
408    status = nf90_get_var(ncid,apvarid,ap)
409    status = nf90_inq_varid(ncid,"bp",bpvarid)
410    status = nf90_get_var(ncid,bpvarid,bp)
411    status = nf90_close(ncid)
[2779]412
[3319]413    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)
[2980]414#else
[3149]415    ps_start_PCM(1) = ps(1)
[2980]416#endif
417
[3442]418! Save initial surface pressure
419ps0 = ps
420
[3096]421! In the PCM, these values are given to the physic by the dynamic.
[3317]422! Here we simply read them in the "startfi.nc" file
[3143]423status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
[2963]424
[3028]425status = nf90_inq_varid(ncid,"longitude",lonvarid)
426status = nf90_get_var(ncid,lonvarid,longitude)
[2963]427
[3028]428status = nf90_inq_varid(ncid,"latitude",latvarid)
429status = nf90_get_var(ncid,latvarid,latitude)
[2963]430
[3028]431status = nf90_inq_varid(ncid,"area",areavarid)
432status = nf90_get_var(ncid,areavarid,cell_area)
[2963]433
[3028]434status = nf90_inq_varid(ncid,"soildepth",sdvarid)
435status = nf90_get_var(ncid,sdvarid,mlayer)
[2963]436
[3028]437status = nf90_close(ncid)
[2963]438
[3317]439! I_b.2 Read the "startfi.nc"
[2779]440! First we read the initial state (starfi.nc)
[2842]441#ifndef CPP_STD
[3143]442    call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
[3149]443                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,          &
[3130]444                  watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
[2779]445
[3070]446    ! Remove unphysical values of surface tracer
447    where (qsurf < 0.) qsurf = 0.
[2885]448
[3143]449    call surfini(ngrid,nslope,qsurf)
[2842]450#else
[3028]451    call phys_state_var_init(nq)
452    if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
453    call initracer(ngrid,nq)
454    call iniaerosol()
455    allocate(tsurf_read_generic(ngrid))
456    allocate(qsurf_read_generic(ngrid,nq))
457    allocate(tsoil_read_generic(ngrid,nsoilmx))
[3114]458    allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
[3028]459    allocate(emis_read_generic(ngrid))
460    allocate(tsurf(ngrid,1))
461    allocate(qsurf(ngrid,nq,1))
462    allocate(tsoil(ngrid,nsoilmx,1))
463    allocate(emis(ngrid,1))
464    allocate(watercap(ngrid,1))
465    allocate(watercaptag(ngrid))
466    allocate(albedo_read_generic(ngrid,2))
467    allocate(albedo(ngrid,2,1))
468    allocate(inertiesoil(ngrid,nsoilmx,1))
[3143]469    call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
[3149]470                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,               &
471                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,        &
[3114]472                  rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
[3065]473    call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
[2842]474
[3028]475    nslope = 1
476    call ini_comslope_h(ngrid,1)
[2842]477
[3149]478    qsurf(:,:,1) = qsurf_read_generic
479    tsurf(:,1) = tsurf_read_generic
480    tsoil(:,:,1) = tsoil_read_generic
481    emis(:,1) = emis_read_generic
[3028]482    watercap(:,1) = 0.
483    watercaptag(:) = .false.
484    albedo(:,1,1) = albedo_read_generic(:,1)
485    albedo(:,2,1) = albedo_read_generic(:,2)
[3149]486    inertiesoil(:,:,1) = inertiedat
[2842]487
[3028]488    if (nslope == 1) then
489        def_slope(1) = 0
490        def_slope(2) = 0
491        def_slope_mean = 0
492        subslope_dist(:,1) = 1.
493    endif
[2842]494
[3070]495    ! Remove unphysical values of surface tracer
[3149]496    qsurf(:,:,1) = qsurf_read_generic
[3070]497    where (qsurf < 0.) qsurf = 0.
[2842]498#endif
499
[3028]500do nnq = 1,nqtot  ! Why not using ini_tracer ?
501    if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
502    if (noms(nnq) == "h2o_vap") then
503        igcm_h2o_vap = nnq
[3143]504        mmol(igcm_h2o_vap) = 18.
[3028]505    endif
506    if (noms(nnq) == "co2") igcm_co2 = nnq
[3065]507enddo
[3039]508r = 8.314511*1000./mugaz
[3028]509
[2835]510!------------------------
[3028]511! I   Initialization
[2835]512!    I_c Subslope parametrisation
513!------------------------
[3028]514! Define some slope statistics
515iflat = 1
516do islope = 2,nslope
517    if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
518enddo
[2794]519
[3028]520write(*,*) 'Flat slope for islope = ',iflat
521write(*,*) 'corresponding criterium = ',def_slope_mean(iflat)
[2794]522
[3028]523allocate(flag_co2flow(ngrid,nslope))
524allocate(flag_co2flow_mesh(ngrid))
525allocate(flag_h2oflow(ngrid,nslope))
526allocate(flag_h2oflow_mesh(ngrid))
[2835]527
[3149]528flag_co2flow = 0
529flag_co2flow_mesh = 0
530flag_h2oflow = 0
531flag_h2oflow_mesh = 0
[2835]532
[2794]533!------------------------
[3028]534! I   Initialization
[3161]535!    I_d Read the PCM data and convert them to the physical grid
[3028]536!------------------------
[3096]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
538call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
[2794]539
[3367]540allocate(tsoil_avg(ngrid,nsoilmx,nslope))
541allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope))
[3149]542allocate(vmr_co2_PCM(ngrid,timelen))
[3028]543allocate(ps_timeseries(ngrid,timelen))
[3149]544allocate(min_co2_ice(ngrid,nslope,2))
545allocate(min_h2o_ice(ngrid,nslope,2))
546allocate(tsurf_avg_yr1(ngrid,nslope))
[3367]547allocate(tsurf_avg(ngrid,nslope))
[3149]548allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
[3426]549allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen))
[3028]550allocate(q_co2_PEM_phys(ngrid,timelen))
551allocate(q_h2o_PEM_phys(ngrid,timelen))
[3367]552allocate(watersurf_density_avg(ngrid,nslope))
[3028]553allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
[3149]554allocate(Tsurfavg_before_saved(ngrid,nslope))
[3028]555allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
556allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
557allocate(delta_co2_adsorbded(ngrid))
[3330]558allocate(co2ice_disappeared(ngrid,nslope))
[3031]559allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
560allocate(delta_h2o_icetablesublim(ngrid))
[3028]561allocate(delta_h2o_adsorbded(ngrid))
[3149]562allocate(vmr_co2_PEM_phys(ngrid,timelen))
[2794]563
[3028]564write(*,*) "Downloading data Y1..."
[3199]565call 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), &
[3426]566                   tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys,                     &
[3367]567                   watersurf_density_avg,watersoil_density_timeseries)
[3199]568write(*,*) "Downloading data Y1 done!"
[2985]569
[3096]570! 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
[3199]571write(*,*) "Downloading data Y2..."
[3149]572call 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), &
[3426]573                   tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys,                         &
[3367]574                   watersurf_density_avg,watersoil_density_timeseries)
[3199]575write(*,*) "Downloading data Y2 done!"
[2794]576
[2835]577!------------------------
[3028]578! I   Initialization
579!    I_e Initialization of the PEM variables and soil
[2835]580!------------------------
[3028]581call end_comsoil_h_PEM
582call ini_comsoil_h_PEM(ngrid,nslope)
583call end_adsorption_h_PEM
584call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
585call end_ice_table_porefilling
586call ini_ice_table_porefilling(ngrid,nslope)
[2794]587
[3028]588if (soil_pem) then
[3426]589    do t = 1,timelen
590        tsoil_anom(:,:,:,t) = tsoil_anom(:,:,:,t) - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - <Tsoil> to recompute properly tsoil in the restart
591    enddo
[3028]592    call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
[3367]593    tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg
[3426]594    tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_anom
[3149]595    watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries
[3070]596    do l = nsoilmx + 1,nsoilmx_PEM
[3367]597        tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:)
[3070]598        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
[3028]599    enddo
[3367]600    watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
[3028]601endif !soil_pem
[3426]602deallocate(tsoil_avg)
[2794]603
[2779]604!------------------------
[3028]605! I   Initialization
[3143]606!    I_f Compute tendencies
[3028]607!------------------------
[3149]608allocate(tend_co2_ice(ngrid,nslope),tend_h2o_ice(ngrid,nslope))
609allocate(tend_co2_ice_ini(ngrid,nslope))
[2779]610
[3028]611! Compute the tendencies of the evolution of ice over the years
[3149]612call compute_tend(ngrid,nslope,min_co2_ice,tend_co2_ice)
613call compute_tend(ngrid,nslope,min_h2o_ice,tend_h2o_ice)
614tend_co2_ice_ini = tend_co2_ice
[3365]615deallocate(min_co2_ice,min_h2o_ice)
[2895]616
[2835]617!------------------------
[3028]618! I   Initialization
[3384]619!    I_g Save the initial situation
[3028]620!------------------------
[3384]621allocate(zplev_PCM(ngrid,nlayer + 1))
622Total_surface = 0.
623do ig = 1,ngrid
624    Total_surface = Total_surface + cell_area(ig)
625    zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig)
626enddo
627global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface
628global_avg_press_PCM = global_avg_press_old
629global_avg_press_new = global_avg_press_old
630write(*,*) "Total surface of the planet =", Total_surface
631write(*,*) "Initial global average pressure =", global_avg_press_PCM
632
633!------------------------
634! I   Initialization
635!    I_h Read the "startpem.nc"
636!------------------------
[3149]637allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope))
638
[3297]639allocate(stratif(ngrid,nslope))
[3319]640if (layering_algo) then
641    do islope = 1,nslope
642        do i = 1,ngrid
643            call ini_layering(stratif(i,islope))
644        enddo
[3297]645    enddo
[3319]646endif
[3297]647
[3088]648call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
[3367]649              porefillingice_thickness,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,          &
[3149]650              tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM,              &
[3367]651              watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,delta_co2_adsorbded,                &
[3297]652              h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
[2779]653
[3384]654! We save the places where h2o ice is sublimating
655! We compute the surface of h2o ice sublimating
656allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope),is_co2ice_ini(ngrid,nslope))
657co2ice_ini_surf = 0.
658h2oice_ini_surf = 0.
659ini_co2ice_sublim = .false.
660ini_h2oice_sublim = .false.
661is_co2ice_ini = .false.
662co2ice_disappeared = .false.
663do i = 1,ngrid
664    do islope = 1,nslope
665        if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
666        if (tend_co2_ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
667            ini_co2ice_sublim(i,islope) = .true.
668            co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
669        endif
670        if (tend_h2o_ice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
671            ini_h2oice_sublim(i,islope) = .true.
672            h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
673        endif
674    enddo
675enddo
676write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf
677write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf
678
[3149]679delta_h2o_icetablesublim = 0.
[3130]680
[3028]681if (adsorption_pem) then
682    totmassco2_adsorbded = 0.
683    totmassh2o_adsorbded = 0.
684    do ig = 1,ngrid
[3070]685        do islope = 1,nslope
[3028]686            do l = 1,nsoilmx_PEM - 1
[3264]687                if (l==1) then
688                   totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
[3028]689                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
[3264]690                   totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
[3028]691                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
[3264]692                else
693                   totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
694                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
695                   totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
696                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
697                endif
[3028]698            enddo
[2961]699        enddo
[3028]700    enddo
[3143]701    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
702    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
[3028]703endif ! adsorption
[3149]704deallocate(tsurf_avg_yr1)
[2794]705
[2835]706!------------------------
[3028]707! I   Initialization
[2835]708!    I_i Compute orbit criterion
[3028]709!------------------------
[2842]710#ifndef CPP_STD
[3050]711    call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
[2842]712#else
[3050]713    call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
[2842]714#endif
[2794]715
[3403]716year_iter_max = Max_iter_pem
717if (evol_orbit_pem) call orbit_param_criterion(i_myear,year_iter_max)
718
[3028]719!-------------------------- END INITIALIZATION -------------------------
[2794]720
[3028]721!-------------------------------- RUN ----------------------------------
[2794]722!------------------------
723! II  Run
[3065]724!    II_a Update pressure, ice and tracers
[2794]725!------------------------
[3028]726year_iter = 0
[3149]727stopPEM = 0
[3319]728if (layering_algo) then
729    allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope))
730    new_str = .true.
731    new_lag1 = .true.
732    new_lag2 = .true.
733    do islope = 1,nslope
734        do ig = 1,ngrid
735            current1(ig,islope)%p => stratif(ig,islope)%top
736            current2(ig,islope)%p => stratif(ig,islope)%top
737        enddo
[3297]738    enddo
[3319]739endif
[2794]740
[3388]741do while (year_iter < year_iter_max .and. i_myear < n_myear)
[2835]742! II.a.1. Compute updated global pressure
[3028]743    write(*,*) "Recomputing the new pressure..."
744    do i = 1,ngrid
745        do islope = 1,nslope
[3188]746            global_avg_press_new = global_avg_press_new - CO2cond_ps*g*cell_area(i)*tend_co2_ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
[3028]747        enddo
748    enddo
[3065]749
[3028]750    if (adsorption_pem) then
751        do i = 1,ngrid
[3149]752            global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
[3050]753        enddo
[3028]754    endif
[3149]755    write(*,*) 'Global average pressure old time step',global_avg_press_old
756    write(*,*) 'Global average pressure new time step',global_avg_press_new
[2835]757
[3161]758! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption)
[3070]759    allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
[3028]760    write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
761    do l = 1,nlayer + 1
762        do ig = 1,ngrid
763            zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
764        enddo
765    enddo
[2779]766
[2835]767! II.a.3. Surface pressure timeseries
[3028]768    write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..."
769    do ig = 1,ngrid
[3149]770        ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old
[3028]771    enddo
[2779]772
[2835]773! II.a.4. New pressure levels timeseries
[3149]774    allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen))
[3028]775    write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
776    do l = 1,nlayer + 1
777        do ig = 1,ngrid
778            zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
779        enddo
780    enddo
[2779]781
[2835]782! II.a.5. Tracers timeseries
[3028]783    write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..."
[2794]784
[3028]785    l = 1
786    do ig = 1,ngrid
787        do t = 1,timelen
788            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))/ &
789                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
[3143]790            if (q_h2o_PEM_phys(ig,t) < 0) then
791                q_h2o_PEM_phys(ig,t) = 1.e-30
792            else if (q_h2o_PEM_phys(ig,t) > 1) then
793                q_h2o_PEM_phys(ig,t) = 1.
794            endif
[3028]795        enddo
796    enddo
[2794]797
[3028]798    do ig = 1,ngrid
[3065]799        do t = 1,timelen
[3028]800            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))/ &
[3122]801                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))                       &
802                                + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))                       &
803                                -  (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/                     &
[3028]804                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
805            if (q_co2_PEM_phys(ig,t) < 0) then
806                q_co2_PEM_phys(ig,t) = 1.e-30
[3143]807            else if (q_co2_PEM_phys(ig,t) > 1) then
[3028]808                q_co2_PEM_phys(ig,t) = 1.
809            endif
810            mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
[3149]811            vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
[3028]812        enddo
813    enddo
[2794]814
[3028]815    deallocate(zplev_new_timeseries,zplev_old_timeseries)
816
817!------------------------
[2835]818! II  Run
[3149]819!    II_b Evolution of ice
[3028]820!------------------------
[3149]821    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,tend_h2o_ice,stopPEM)
822    call evol_co2_ice(ngrid,nslope,co2_ice,tend_co2_ice)
[3319]823    if (layering_algo) then
824        do islope = 1,nslope
825            do ig = 1,ngrid
826                call make_layering(stratif(ig,islope),tend_co2_ice(ig,islope),tend_h2o_ice(ig,islope),tend_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
827            enddo
[3297]828        enddo
[3319]829    endif
[2794]830
831!------------------------
832! II  Run
[3149]833!    II_c Flow of glaciers
[2794]834!------------------------
[3181]835    if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, &
[3149]836                                            global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh)
[3367]837    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)
[3065]838
[2794]839!------------------------
840! II  Run
[2835]841!    II_d Update surface and soil temperatures
[2794]842!------------------------
[2835]843! II_d.1 Update Tsurf
[3028]844    write(*,*) "Updating the new Tsurf"
845    bool_sublim = .false.
[3367]846    Tsurfavg_before_saved = tsurf_avg
[3028]847    do ig = 1,ngrid
848        do islope = 1,nslope
[3330]849            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
[3331]850                co2ice_disappeared(ig,islope) = .true.
[3028]851                if (latitude_deg(ig) > 0) then
852                    do ig_loop = ig,ngrid
853                        do islope_loop = islope,iflat,-1
[3330]854                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
[3367]855                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
[3028]856                                bool_sublim = .true.
857                                exit
858                            endif
859                        enddo
860                        if (bool_sublim) exit
861                    enddo
862                else
863                    do ig_loop = ig,1,-1
864                        do islope_loop = islope,iflat
[3330]865                            if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
[3367]866                                tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
[3028]867                                bool_sublim = .true.
868                                exit
869                            endif
870                        enddo
871                        if (bool_sublim) exit
872                    enddo
[2835]873                endif
[3149]874                if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then
[3028]875                    albedo(ig,1,islope) = albedo_h2o_frost
876                    albedo(ig,2,islope) = albedo_h2o_frost
877                else
878                    albedo(ig,1,islope) = albedodat(ig)
[3065]879                    albedo(ig,2,islope) = albedodat(ig)
[3028]880                    emis(ig,islope) = emissiv
881                endif
[3149]882            else if ((co2_ice(ig,islope) > 1.e-3) .and. (tend_co2_ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2
[3028]883                ave = 0.
884                do t = 1,timelen
[3367]885                    ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.))
[2794]886                enddo
[3367]887                tsurf_avg(ig,islope) = ave/timelen
[2835]888            endif
889        enddo
[3028]890    enddo
[2794]891
[3028]892    do t = 1,timelen
[3367]893        tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_avg - Tsurfavg_before_saved
[3028]894    enddo
895    ! for the start
896    do ig = 1,ngrid
[2835]897        do islope = 1,nslope
[3367]898            tsurf(ig,islope) =  tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_avg(ig,islope))
[2794]899        enddo
[3028]900    enddo
[2794]901
[3028]902    if (soil_pem) then
[2794]903
[2835]904! II_d.2 Update soil temperature
[3426]905        write(*,*)"Updating soil temperature"
[3028]906        allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
[3426]907        do islope = 1,nslope
908            call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
909            call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
[2794]910
[3028]911            do t = 1,timelen
[3426]912                Tsoil_locslope(:,1:nsoilmx) = tsoil_PEM(:,1:nsoilmx,islope) + tsoil_anom(:,:,islope,t)
913                Tsoil_locslope(:,nsoilmx + 1:) = tsoil_PEM(:,nsoilmx + 1:,islope)
[3028]914                do ig = 1,ngrid
915                    do isoil = 1,nsoilmx_PEM
916                        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)
[3426]917                        if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
[3028]918                    enddo
919                enddo
920            enddo
921        enddo
[3367]922        watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
[2794]923
[3028]924        write(*,*) "Update of soil temperature done"
[2888]925
[3426]926        deallocate(Tsoil_locslope)
[2849]927
[2835]928! II_d.3 Update the ice table
[3170]929        if (icetable_equilibrium) then
930            write(*,*) "Compute ice table"
931            porefillingice_thickness_prev_iter = porefillingice_thickness
[3367]932            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
933            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
[3170]934        endif
[3486]935        ! II_d.3 Update the ice table
936        if (icetable_dynamic) then
937            write(*,*) "Compute ice table"
938            porefillingice_thickness_prev_iter = porefillingice_thickness
939            call dyn_ss_ice_m_wrapper(ngrid,nsoilmx,TI_PEM,ps,mmol(igcm_h2o_vap),tsoil_PEM,porefillingice_depth,ss_ice_pf,ss_ice_pf_new,porefillingice_depth_new)
940           
941            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
942        endif
943
[3122]944! II_d.4 Update the soil thermal properties
[3149]945        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
[2794]946
[3143]947! II_d.5 Update the mass of the regolith adsorbed
[3028]948        if (adsorption_pem) then
[3159]949            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tend_h2o_ice,tend_co2_ice,                   &
[3149]950                                     h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
951                                     h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
[2794]952
[3028]953            totmassco2_adsorbded = 0.
954            totmassh2o_adsorbded = 0.
955            do ig = 1,ngrid
956                do islope =1, nslope
[3264]957                    do l = 1,nsoilmx_PEM
958                        if (l==1) then
959                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
960                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
961                            totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
962                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
963                        else
964                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
965                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
966                            totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
967                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
968                        endif
[3028]969                    enddo
970                enddo
971            enddo
972            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
973            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
974        endif
975    endif !soil_pem
976
[2794]977!------------------------
978! II  Run
[3088]979!    II_e Outputs
[2794]980!------------------------
[3367]981    call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/global_avg_press_new/))
[3088]982    do islope = 1,nslope
983        write(str2(1:2),'(i2.2)') islope
[3181]984        call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
985        call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
986        call writediagpem(ngrid,'tend_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tend_h2o_ice(:,islope))
987        call writediagpem(ngrid,'tend_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tend_co2_ice(:,islope))
988        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
989        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
[3339]990        if (icetable_equilibrium) then
[3181]991            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
992            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
[3171]993        endif
[3486]994        if (icetable_dynamic) then
995            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
996            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
997        endif
998
[3339]999        if (soil_pem) then
[3171]1000            call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
1001            call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
1002            if (adsorption_pem) then
1003                call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
1004                call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
1005            endif                       
1006        endif
[3088]1007    enddo
1008
1009!------------------------
1010! II  Run
1011!    II_f Update the tendencies
1012!------------------------
[3149]1013    call recomp_tend_co2_slope(ngrid,nslope,timelen,tend_co2_ice,tend_co2_ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries, &
1014                               global_avg_press_PCM,global_avg_press_new)
[2794]1015
[2835]1016!------------------------
1017! II  Run
[3088]1018!    II_g Checking the stopping criterion
[2835]1019!------------------------
[3389]1020
1021    write(*,*) "Checking the stopping criteria..."
[3327]1022    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
1023    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)
[3028]1024    year_iter = year_iter + dt_pem
[3039]1025    i_myear = i_myear + dt_pem
[3430]1026    if (stopPEM <= 0 .and. year_iter >= year_iter_max) stopPEM = 5
1027    if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
[3389]1028    call system_clock(c2)
[3430]1029    if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
[3432]1030    if (stopPEM > 0) then
[3149]1031        select case (stopPEM)
1032            case(1)
[3159]1033                write(*,*) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
[3149]1034            case(2)
[3159]1035                write(*,*) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
[3149]1036            case(3)
1037                write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
1038            case(4)
1039                write(*,*) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
1040            case(5)
[3419]1041                write(*,*) "STOPPING because maximum number of iterations is reached (possibly due to orbital parameters):", stopPEM
[3149]1042            case(6)
1043                write(*,*) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
[3363]1044            case(7)
1045                write(*,*) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM
[3430]1046            case(8)
1047                write(*,*) "STOPPING because the layering algorithm met an hasty end:", stopPEM
[3149]1048            case default
1049                write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM
1050        end select
[2779]1051        exit
[3028]1052    else
[3143]1053        write(*,*) "The PEM can continue!"
[3039]1054        write(*,*) "Number of iterations of the PEM: year_iter =", year_iter
1055        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
[3028]1056    endif
[2779]1057
[3149]1058    global_avg_press_old = global_avg_press_new
[2779]1059
[3149]1060enddo ! big time iteration loop of the pem
[3028]1061!------------------------------ END RUN --------------------------------
[2779]1062
[3028]1063!------------------------------- OUTPUT --------------------------------
[2794]1064!------------------------
1065! III Output
[2835]1066!    III_a Update surface value for the PCM start files
[2794]1067!------------------------
[2835]1068! III_a.1 Ice update (for startfi)
[2779]1069
[3149]1070watercap = 0.
[3159]1071perennial_co2ice = co2_ice
[3028]1072do ig = 1,ngrid
[3159]1073    ! H2O ice metamorphism
[3161]1074    if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
[3308]1075        h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1076        qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
[3159]1077    endif
1078
1079    ! Is H2O ice still considered as an infinite reservoir for the PCM?
[3149]1080    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
[3159]1081        ! There is enough ice to be considered as an infinite reservoir
[3149]1082        watercaptag(ig) = .true.
1083    else
[3159]1084        ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
[3149]1085        watercaptag(ig) = .false.
1086        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1087        h2o_ice(ig,:) = 0.
[3028]1088    endif
[3159]1089
1090    ! CO2 ice metamorphism
[3161]1091    if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
[3308]1092        perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1093        qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
[3159]1094    endif
[3028]1095enddo
[2888]1096
[3446]1097! III_a.2 Tsoil update (for startfi)
[3028]1098if (soil_pem) then
[3149]1099    call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
[3426]1100    tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen)
[3172]1101#ifndef CPP_STD
[3181]1102    flux_geo = fluxgeo
[3172]1103#endif
[3028]1104endif
[3426]1105deallocate(tsoil_anom)
[2779]1106
[2835]1107! III_a.4 Pressure (for start)
[3149]1108ps = ps*global_avg_press_new/global_avg_press_PCM
1109ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
[2794]1110
[2835]1111! III_a.5 Tracer (for start)
[3028]1112allocate(zplev_new(ngrid,nlayer + 1))
[2835]1113
[3028]1114do l = 1,nlayer + 1
[3149]1115    zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
[3028]1116enddo
[2835]1117
[3028]1118do nnq = 1,nqtot
1119    if (noms(nnq) /= "co2") then
1120        do l = 1,llm - 1
1121            do ig = 1,ngrid
[3149]1122                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))
[3028]1123            enddo
1124            q(:,llm,nnq) = q(:,llm - 1,nnq)
1125        enddo
1126    else
1127        do l = 1,llm - 1
1128            do ig = 1,ngrid
[3149]1129                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)) &
1130                              + ((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))
[3028]1131            enddo
1132            q(:,llm,nnq) = q(:,llm - 1,nnq)
1133        enddo
1134    endif
1135enddo
[2835]1136
[3096]1137! Conserving the tracers mass for PCM start files
[3028]1138do nnq = 1,nqtot
1139    do ig = 1,ngrid
1140        do l = 1,llm - 1
1141            if (q(ig,l,nnq) > 1 .and. (noms(nnq) /= "dust_number") .and. (noms(nnq) /= "ccn_number") .and. (noms(nnq) /= "stormdust_number") .and. (noms(nnq) /= "topdust_number")) then
[3065]1142                extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1143                q(ig,l,nnq) = 1.
1144                q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
[3028]1145                write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
[2835]1146           endif
[3028]1147            if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1148        enddo
1149    enddo
1150enddo
[2779]1151
[3039]1152if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
[2779]1153
1154!------------------------
[3028]1155! III Output
[3317]1156!    III_b Write "restart.nc" and "restartfi.nc"
[3028]1157!------------------------
[3317]1158! III_b.1 Write "restart.nc"
[3042]1159ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
[3028]1160pday = day_ini
[3042]1161ztime_fin = time_phys
[2779]1162
[3028]1163allocate(p(ip1jmp1,nlayer + 1))
[2980]1164#ifndef CPP_1D
[3442]1165    ! Correction on teta due to surface pressure changes
1166    do l = 1,nlayer
1167        do i = 1,ip1jmp1
1168            teta(i,l)= teta(i,l)*(1 + (ps(i) - ps0(i))*year_iter/ps(i))**kappa
1169        enddo
1170    enddo
1171    ! Correction on atmospheric pressure
1172    call pression(ip1jmp1,ap,bp,ps,p)
1173    ! Correction on the mass of atmosphere
[3028]1174    call massdair(p,masse)
[3317]1175    call dynredem0("restart.nc",day_ini,phis)
1176    call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps)
1177    write(*,*) "restart.nc has been written"
[2980]1178#else
[3363]1179    call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1180    write(*,*) "restart1D.txt has been written"
[2980]1181#endif
1182
[3317]1183! III_b.2 Write the "restartfi.nc"
[2842]1184#ifndef CPP_STD
[3317]1185    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
[3327]1186                  nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
[3028]1187                  inertiedat,def_slope,subslope_dist)
[3327]1188    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,      &
[3114]1189                  ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,        &
1190                  albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
[3130]1191                  wstar,watercap,perennial_co2ice)
[2842]1192#else
[3317]1193    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
[3327]1194                  nlayer,nq,ptimestep,pday,time_phys,cell_area,    &
[3028]1195                  albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[3327]1196    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,       &
[3114]1197                  ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1198                  cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab,   &
1199                  tsea_ice,sea_ice)
[2842]1200#endif
[3317]1201write(*,*) "restartfi.nc has been written"
[2842]1202
[2794]1203!------------------------
1204! III Output
[3161]1205!    III_c Write the "restartpem.nc"
[2794]1206!------------------------
[3319]1207if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
[3206]1208call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
[3088]1209call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
[3297]1210             TI_PEM,porefillingice_depth,porefillingice_thickness,       &
1211             co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
[3088]1212write(*,*) "restartpem.nc has been written"
[2779]1213
[3149]1214call info_PEM(year_iter,stopPEM,i_myear,n_myear)
1215
[3039]1216write(*,*) "The PEM has run for", year_iter, "Martian years."
1217write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1218write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
1219write(*,*) "LL & RV & JBC: so far, so good!"
[2794]1220
[3319]1221if (layering_algo) then
1222    do islope = 1,nslope
1223        do i = 1,ngrid
1224            call del_layering(stratif(i,islope))
1225        enddo
[3297]1226    enddo
[3319]1227    deallocate(new_str,new_lag1,new_lag2,current1,current2)
1228endif
[3149]1229deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
[3367]1230deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
1231deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
[3149]1232deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
[3330]1233deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
[3028]1234!----------------------------- END OUTPUT ------------------------------
[2897]1235
[2779]1236END PROGRAM pem
Note: See TracBrowser for help on using the repository browser.