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

Last change on this file since 3525 was 3525, checked in by jbclement, 5 days ago

PEM:
Computation of <soil thermal intertia> and <H2O mass subsurface/surface exchange> according to the presence of subsurface ice provided by the (Norbert's) dynamic method + few cleanings.
JBC
File size: 61.8 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
[3498]53use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
[3161]54use orbit_param_criterion_mod,  only: orbit_param_criterion
55use recomp_orb_param_mod,       only: recomp_orb_param
[3493]56use ice_table_mod,              only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, &
57                                      ini_ice_table, 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
[3498]72use layering_mod,               only: d_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
[3498]197real, dimension(:,:),   allocatable :: d_co2ice          ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
198real, dimension(:,:),   allocatable :: d_co2ice_ini      ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM
199real, dimension(:,:),   allocatable :: d_h2oice          ! physical point x slope field: Tendency of evolution of perennial h2o ice
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
[3493]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]
208real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries             ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
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]
211real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
212real, dimension(:,:),     allocatable :: Tsoil_locslope                   ! Physic x Soil: Intermediate when computing Tsoil [K]
213real, dimension(:),       allocatable :: Tsurf_locslope                   ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
214real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
215real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
216real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
217real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
218real, dimension(:,:),     allocatable :: Tsurfavg_before_saved            ! Surface temperature saved from previous time step [K]
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
224logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
[3525]225real, dimension(:,:),     allocatable :: icetable_thickness_old           ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
226real, dimension(:,:,:),   allocatable :: ice_porefilling_old              ! ngrid x nslope: Ice pore filling at the previous iteration [m]
[3493]227real, dimension(:),       allocatable :: delta_h2o_icetablesublim         ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
228real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
229real                                  :: ssi_depth                        ! Ice table depth (output) to compute the dynamic ice table
[3065]230
[3028]231! Some variables for the PEM run
[3363]232real, parameter :: year_step = 1   ! Timestep for the pem
[3498]233real            :: i_myear_leg       ! Number of iteration
234real            :: n_myear_leg     ! Maximum number of iterations before stopping
235real            :: i_myear         ! Global number of Martian years of the chained simulations
236real            :: n_myear         ! Maximum number of Martian years of the chained simulations
[3363]237real            :: timestep        ! Timestep [s]
238character(20)   :: job_id          ! Job id provided as argument passed on the command line when the program was invoked
[3394]239logical         :: timewall        ! Flag to use the time limit stopping criterion in case of a PEM job
[3363]240integer(kind=8) :: cr              ! Number of clock ticks per second (count rate)
241integer(kind=8) :: c1, c2          ! Counts of processor clock
242character(100)  :: chtimelimit     ! Time limit for the PEM job outputted by the SLURM command
243real            :: timelimit       ! Time limit for the PEM job in seconds
[3365]244real, parameter :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default
[3363]245integer         :: cstat, days, hours, minutes, seconds
246character(1)    :: sep
[2779]247
[2842]248#ifdef CPP_STD
[3065]249    real                                :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
250    real                                :: albedo_h2o_frost              ! albedo of h2o frost
[3143]251    real, dimension(:),     allocatable :: tsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
252    real, dimension(:,:),   allocatable :: qsurf_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
253    real, dimension(:,:),   allocatable :: tsoil_read_generic            ! Temporary variable to do the subslope transfert dimension when reading form generic
254    real, dimension(:),     allocatable :: emis_read_generic             ! Temporary variable to do the subslope transfert dimension when reading form generic
255    real, dimension(:,:),   allocatable :: albedo_read_generic           ! Temporary variable to do the subslope transfert dimension when reading form generic
[3065]256    real, dimension(:,:),   allocatable :: tsurf                         ! Subslope variable, only needed in the GENERIC case
257    real, dimension(:,:,:), allocatable :: qsurf                         ! Subslope variable, only needed in the GENERIC case
258    real, dimension(:,:,:), allocatable :: tsoil                         ! Subslope variable, only needed in the GENERIC case
259    real, dimension(:,:),   allocatable :: emis                          ! Subslope variable, only needed in the GENERIC case
260    real, dimension(:,:),   allocatable :: watercap                      ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
[3068]261    logical, dimension(:),  allocatable :: watercaptag                   ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
[3065]262    real, dimension(:,:,:), allocatable :: albedo                        ! Subslope variable, only needed in the GENERIC case
263    real, dimension(:,:,:), allocatable :: inertiesoil                   ! Subslope variable, only needed in the GENERIC case
[2842]264#endif
265
[2980]266#ifdef CPP_1D
[3143]267    integer            :: nsplit_phys
268    integer, parameter :: jjm_value = jjm - 1
[3386]269    integer            :: day_step
[3065]270
271    ! Dummy variables to use the subroutine 'init_testphys1d'
[3129]272    logical                             :: therestart1D, therestartfi
[3068]273    integer                             :: ndt, day0
274    real                                :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau
275    real, dimension(:),     allocatable :: zqsat
276    real, dimension(:,:,:), allocatable :: dq, dqdyn
277    real, dimension(nlayer)             :: play, w
278    real, dimension(nlayer + 1)         :: plev
[2980]279#else
[3143]280    integer, parameter              :: jjm_value = jjm
281    real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart
282    real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart
[2980]283#endif
284
[3028]285! Loop variables
[3206]286integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap
[2779]287
[3363]288! Elapsed time with system clock
289call system_clock(count_rate = cr)
290call system_clock(c1)
[3394]291timewall = .true.
[3363]292timelimit = 86400 ! 86400 seconds = 24 h by default
293if (command_argument_count() > 0) then
294    ! Read the job id passed as argument to the program
295    call get_command_argument(1,job_id)
296    ! Execute the system command
297    call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat)
[3403]298    if (cstat /= 0) then
299        call execute_command_line('qstat -f '//trim(job_id)//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt', cmdstat = cstat)
300        if (cstat > 0) then
301            error stop 'pem: command execution failed!'
302        else if (cstat < 0) then
303            error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!'
304        endif
[3363]305    endif
306    ! Read the output
307    open(1,file = 'tmp_cmdout.txt',status = 'old')
308    read(1,'(a)') chtimelimit
309    close(1)
310    chtimelimit = trim(chtimelimit)
311    call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat)
312    if (cstat > 0) then
313        error stop 'pem: command execution failed!'
314    else if (cstat < 0) then
315        error stop 'pem: command execution not supported!'
316    endif
[3387]317    if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS"
[3388]318        read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds
[3363]319        timelimit = days*86400 + hours*3600 + minutes*60 + seconds
[3387]320    else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS"
[3388]321        read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds
[3363]322        timelimit = hours*3600 + minutes*60 + seconds
[3387]323    else ! 'chtimelimit' format is "MM:SS"
324        read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds
325        timelimit = minutes*60 + seconds
[3363]326    endif
[3394]327else
328    timewall = .false.
[3363]329endif
330
[3028]331! Parallel variables
[2842]332#ifndef CPP_STD
[3028]333    is_sequential = .true.
334    is_parallel = .false.
335    is_mpi_root = .true.
336    is_omp_root = .true.
337    is_master = .true.
[2842]338#endif
[2779]339
[3065]340! Some constants
[3028]341day_ini = 0    ! test
342time_phys = 0. ! test
343ngrid = ngridmx
344A = (1/m_co2 - 1/m_noco2)
345B = 1/m_noco2
346year_day = 669
347daysec = 88775.
[3426]348timestep = year_day*daysec*year_step
[2794]349
[3028]350!----------------------------- INITIALIZATION --------------------------
[2779]351!------------------------
[3028]352! I   Initialization
[3161]353!    I_a Read the "run.def"
[2779]354!------------------------
[2980]355#ifndef CPP_1D
[3383]356    dtphys = daysec/48. ! Dummy value (overwritten in phyetat0)
[3028]357    call conf_gcm(99,.true.)
358    call infotrac_init
359    nq = nqtot
360    allocate(q(ip1jmp1,llm,nqtot))
[3065]361    allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
[2980]362#else
[3068]363    allocate(q(1,llm,nqtot))
[3065]364    allocate(longitude(1),latitude(1),cell_area(1))
[3129]365
[3143]366    therestart1D = .false. ! Default value
[3363]367    inquire(file = 'start1D.txt',exist = therestart1D)
[3129]368    if (.not. therestart1D) then
[3363]369        write(*,*) 'There is no "start1D.txt" file!'
[3129]370        error stop 'Initialization cannot be done for the 1D PEM.'
371    endif
[3143]372    therestartfi = .false. ! Default value
[3317]373    inquire(file = 'startfi.nc',exist = therestartfi)
[3129]374    if (.not. therestartfi) then
[3317]375        write(*,*) 'There is no "startfi.nc" file!'
[3129]376        error stop 'Initialization cannot be done for the 1D PEM.'
377    endif
378
[3363]379    call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q,        &
380                         time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
[3207]381                         play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
[3065]382    ps(2) = ps(1)
[3028]383    nsplit_phys = 1
[3399]384    day_step = steps_per_sol
[2980]385#endif
[2779]386
[3039]387call conf_pem(i_myear,n_myear)
[2779]388
[2835]389!------------------------
[3028]390! I   Initialization
[3317]391!    I_b Read of the "start.nc" and starfi_evol.nc
[3028]392!------------------------
[3317]393! I_b.1 Read "start.nc"
[3149]394allocate(ps_start_PCM(ngrid))
[2980]395#ifndef CPP_1D
[3143]396    call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0)
[2779]397
[3149]398    call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM)
[2897]399
[3028]400    call iniconst !new
401    call inigeom
[2980]402
[3028]403    allocate(ap(nlayer + 1))
404    allocate(bp(nlayer + 1))
[3143]405    status = nf90_open(start_name,NF90_NOWRITE,ncid)
[3028]406    status = nf90_inq_varid(ncid,"ap",apvarid)
407    status = nf90_get_var(ncid,apvarid,ap)
408    status = nf90_inq_varid(ncid,"bp",bpvarid)
409    status = nf90_get_var(ncid,bpvarid,bp)
410    status = nf90_close(ncid)
[2779]411
[3319]412    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]413#else
[3149]414    ps_start_PCM(1) = ps(1)
[2980]415#endif
416
[3442]417! Save initial surface pressure
418ps0 = ps
419
[3096]420! In the PCM, these values are given to the physic by the dynamic.
[3317]421! Here we simply read them in the "startfi.nc" file
[3143]422status = nf90_open(startfi_name, NF90_NOWRITE, ncid)
[2963]423
[3028]424status = nf90_inq_varid(ncid,"longitude",lonvarid)
425status = nf90_get_var(ncid,lonvarid,longitude)
[2963]426
[3028]427status = nf90_inq_varid(ncid,"latitude",latvarid)
428status = nf90_get_var(ncid,latvarid,latitude)
[2963]429
[3028]430status = nf90_inq_varid(ncid,"area",areavarid)
431status = nf90_get_var(ncid,areavarid,cell_area)
[2963]432
[3028]433status = nf90_inq_varid(ncid,"soildepth",sdvarid)
434status = nf90_get_var(ncid,sdvarid,mlayer)
[2963]435
[3028]436status = nf90_close(ncid)
[2963]437
[3317]438! I_b.2 Read the "startfi.nc"
[2779]439! First we read the initial state (starfi.nc)
[2842]440#ifndef CPP_STD
[3143]441    call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
[3149]442                  tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar,          &
[3130]443                  watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
[2779]444
[3070]445    ! Remove unphysical values of surface tracer
446    where (qsurf < 0.) qsurf = 0.
[2885]447
[3143]448    call surfini(ngrid,nslope,qsurf)
[2842]449#else
[3028]450    call phys_state_var_init(nq)
451    if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
452    call initracer(ngrid,nq)
453    call iniaerosol()
454    allocate(tsurf_read_generic(ngrid))
455    allocate(qsurf_read_generic(ngrid,nq))
456    allocate(tsoil_read_generic(ngrid,nsoilmx))
[3114]457    allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
[3028]458    allocate(emis_read_generic(ngrid))
459    allocate(tsurf(ngrid,1))
460    allocate(qsurf(ngrid,nq,1))
461    allocate(tsoil(ngrid,nsoilmx,1))
462    allocate(emis(ngrid,1))
463    allocate(watercap(ngrid,1))
464    allocate(watercaptag(ngrid))
465    allocate(albedo_read_generic(ngrid,2))
466    allocate(albedo(ngrid,2,1))
467    allocate(inertiesoil(ngrid,nsoilmx,1))
[3143]468    call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
[3149]469                  tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2,               &
470                  qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice,        &
[3114]471                  rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
[3065]472    call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
[2842]473
[3028]474    nslope = 1
475    call ini_comslope_h(ngrid,1)
[2842]476
[3149]477    qsurf(:,:,1) = qsurf_read_generic
478    tsurf(:,1) = tsurf_read_generic
479    tsoil(:,:,1) = tsoil_read_generic
480    emis(:,1) = emis_read_generic
[3028]481    watercap(:,1) = 0.
482    watercaptag(:) = .false.
483    albedo(:,1,1) = albedo_read_generic(:,1)
484    albedo(:,2,1) = albedo_read_generic(:,2)
[3149]485    inertiesoil(:,:,1) = inertiedat
[2842]486
[3028]487    if (nslope == 1) then
488        def_slope(1) = 0
489        def_slope(2) = 0
490        def_slope_mean = 0
491        subslope_dist(:,1) = 1.
492    endif
[2842]493
[3070]494    ! Remove unphysical values of surface tracer
[3149]495    qsurf(:,:,1) = qsurf_read_generic
[3070]496    where (qsurf < 0.) qsurf = 0.
[2842]497#endif
498
[3028]499do nnq = 1,nqtot  ! Why not using ini_tracer ?
500    if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
501    if (noms(nnq) == "h2o_vap") then
502        igcm_h2o_vap = nnq
[3143]503        mmol(igcm_h2o_vap) = 18.
[3028]504    endif
505    if (noms(nnq) == "co2") igcm_co2 = nnq
[3065]506enddo
[3039]507r = 8.314511*1000./mugaz
[3028]508
[2835]509!------------------------
[3028]510! I   Initialization
[2835]511!    I_c Subslope parametrisation
512!------------------------
[3028]513! Define some slope statistics
514iflat = 1
515do islope = 2,nslope
516    if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
517enddo
[2794]518
[3028]519write(*,*) 'Flat slope for islope = ',iflat
520write(*,*) 'corresponding criterium = ',def_slope_mean(iflat)
[2794]521
[3028]522allocate(flag_co2flow(ngrid,nslope))
523allocate(flag_co2flow_mesh(ngrid))
524allocate(flag_h2oflow(ngrid,nslope))
525allocate(flag_h2oflow_mesh(ngrid))
[2835]526
[3149]527flag_co2flow = 0
528flag_co2flow_mesh = 0
529flag_h2oflow = 0
530flag_h2oflow_mesh = 0
[2835]531
[2794]532!------------------------
[3028]533! I   Initialization
[3161]534!    I_d Read the PCM data and convert them to the physical grid
[3028]535!------------------------
[3096]536! 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
537call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
[2794]538
[3367]539allocate(tsoil_avg(ngrid,nsoilmx,nslope))
540allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope))
[3149]541allocate(vmr_co2_PCM(ngrid,timelen))
[3028]542allocate(ps_timeseries(ngrid,timelen))
[3149]543allocate(min_co2_ice(ngrid,nslope,2))
544allocate(min_h2o_ice(ngrid,nslope,2))
545allocate(tsurf_avg_yr1(ngrid,nslope))
[3367]546allocate(tsurf_avg(ngrid,nslope))
[3149]547allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
[3426]548allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen))
[3028]549allocate(q_co2_PEM_phys(ngrid,timelen))
550allocate(q_h2o_PEM_phys(ngrid,timelen))
[3367]551allocate(watersurf_density_avg(ngrid,nslope))
[3028]552allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
[3149]553allocate(Tsurfavg_before_saved(ngrid,nslope))
[3028]554allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
555allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
556allocate(delta_co2_adsorbded(ngrid))
[3330]557allocate(co2ice_disappeared(ngrid,nslope))
[3525]558allocate(icetable_thickness_old(ngrid,nslope))
559allocate(ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
[3031]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)
[3493]585call end_ice_table
586call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
[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!------------------------
[3498]608allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope))
609allocate(d_co2ice_ini(ngrid,nslope))
[2779]610
[3028]611! Compute the tendencies of the evolution of ice over the years
[3498]612call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice)
613call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice)
614d_co2ice_ini = d_co2ice
[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
[3493]648call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
649              icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
[3498]650              ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice,                       &
651              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,         &
[3493]652              delta_co2_adsorbded,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.
[3498]666        if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
[3384]667            ini_co2ice_sublim(i,islope) = .true.
668            co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope)
669        endif
[3498]670        if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then
[3384]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
[3498]716n_myear_leg = Max_iter_pem
717if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
[3403]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!------------------------
[3498]726i_myear_leg = 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
[3498]741do while (i_myear_leg < n_myear_leg .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
[3498]746            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
[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!------------------------
[3498]821    call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
822    call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
[3319]823    if (layering_algo) then
824        do islope = 1,nslope
825            do ig = 1,ngrid
[3498]826                call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)
[3319]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, &
[3498]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
[3498]882            else if ((co2_ice(ig,islope) > 1.e-3) .and. (d_co2ice(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
[3493]930            write(*,*) "Compute ice table (equilibrium method)"
[3525]931            icetable_thickness_old = icetable_thickness
[3493]932            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
[3525]933            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
[3490]934        else if (icetable_dynamic) then
[3493]935            write(*,*) "Compute ice table (dynamic method)"
[3525]936            ice_porefilling_old = ice_porefilling
[3493]937            allocate(porefill(nsoilmx_PEM))
938            do ig = 1,ngrid
939                do islope = 1,nslope
940                    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)
941                    icetable_depth(ig,islope) = ssi_depth
942                    ice_porefilling(ig,:,islope) = porefill
943                enddo
944            enddo
945            deallocate(porefill)
[3525]946            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
[3486]947        endif
948
[3122]949! II_d.4 Update the soil thermal properties
[3525]950        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)
[2794]951
[3143]952! II_d.5 Update the mass of the regolith adsorbed
[3028]953        if (adsorption_pem) then
[3512]954            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                           &
[3149]955                                     h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
956                                     h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
[2794]957
[3028]958            totmassco2_adsorbded = 0.
959            totmassh2o_adsorbded = 0.
960            do ig = 1,ngrid
[3493]961                do islope = 1,nslope
[3264]962                    do l = 1,nsoilmx_PEM
[3493]963                        if (l == 1) then
[3264]964                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
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))* &
967                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
968                        else
[3493]969                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
[3264]970                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
[3493]971                            totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
[3264]972                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
973                        endif
[3028]974                    enddo
975                enddo
976            enddo
977            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
978            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
979        endif
980    endif !soil_pem
981
[2794]982!------------------------
983! II  Run
[3088]984!    II_e Outputs
[2794]985!------------------------
[3367]986    call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/global_avg_press_new/))
[3088]987    do islope = 1,nslope
988        write(str2(1:2),'(i2.2)') islope
[3181]989        call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
990        call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
[3498]991        call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
992        call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
[3181]993        call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
994        call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
[3339]995        if (icetable_equilibrium) then
[3493]996            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
997            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
[3490]998        else if (icetable_dynamic) then
[3493]999            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
1000            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
[3486]1001        endif
1002
[3339]1003        if (soil_pem) then
[3171]1004            call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope))
1005            call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope))
1006            if (adsorption_pem) then
1007                call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
1008                call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope))
1009            endif                       
1010        endif
[3088]1011    enddo
1012
1013!------------------------
1014! II  Run
1015!    II_f Update the tendencies
1016!------------------------
[3516]1017    call recomp_tend_co2_slope(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)
[2794]1018
[2835]1019!------------------------
1020! II  Run
[3088]1021!    II_g Checking the stopping criterion
[2835]1022!------------------------
[3389]1023
1024    write(*,*) "Checking the stopping criteria..."
[3327]1025    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
1026    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)
[3498]1027    i_myear_leg = i_myear_leg + dt
1028    i_myear = i_myear + dt
1029    if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
[3430]1030    if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
[3389]1031    call system_clock(c2)
[3430]1032    if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
[3432]1033    if (stopPEM > 0) then
[3149]1034        select case (stopPEM)
1035            case(1)
[3159]1036                write(*,*) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above."
[3149]1037            case(2)
[3159]1038                write(*,*) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above."
[3149]1039            case(3)
1040                write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above."
1041            case(4)
1042                write(*,*) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above."
1043            case(5)
[3419]1044                write(*,*) "STOPPING because maximum number of iterations is reached (possibly due to orbital parameters):", stopPEM
[3149]1045            case(6)
1046                write(*,*) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM
[3363]1047            case(7)
1048                write(*,*) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM
[3430]1049            case(8)
1050                write(*,*) "STOPPING because the layering algorithm met an hasty end:", stopPEM
[3149]1051            case default
1052                write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM
1053        end select
[2779]1054        exit
[3028]1055    else
[3143]1056        write(*,*) "The PEM can continue!"
[3498]1057        write(*,*) "Number of iterations of the PEM: i_myear_leg =", i_myear_leg
[3039]1058        write(*,*) "Number of simulated Martian years: i_myear =", i_myear
[3028]1059    endif
[2779]1060
[3149]1061    global_avg_press_old = global_avg_press_new
[2779]1062
[3149]1063enddo ! big time iteration loop of the pem
[3028]1064!------------------------------ END RUN --------------------------------
[2779]1065
[3028]1066!------------------------------- OUTPUT --------------------------------
[2794]1067!------------------------
1068! III Output
[2835]1069!    III_a Update surface value for the PCM start files
[2794]1070!------------------------
[2835]1071! III_a.1 Ice update (for startfi)
[2779]1072
[3149]1073watercap = 0.
[3159]1074perennial_co2ice = co2_ice
[3028]1075do ig = 1,ngrid
[3159]1076    ! H2O ice metamorphism
[3161]1077    if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then
[3308]1078        h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
1079        qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
[3159]1080    endif
1081
1082    ! Is H2O ice still considered as an infinite reservoir for the PCM?
[3149]1083    if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then
[3159]1084        ! There is enough ice to be considered as an infinite reservoir
[3149]1085        watercaptag(ig) = .true.
1086    else
[3159]1087        ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost
[3149]1088        watercaptag(ig) = .false.
1089        qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
1090        h2o_ice(ig,:) = 0.
[3028]1091    endif
[3159]1092
1093    ! CO2 ice metamorphism
[3161]1094    if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then
[3308]1095        perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
1096        qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
[3159]1097    endif
[3028]1098enddo
[2888]1099
[3446]1100! III_a.2 Tsoil update (for startfi)
[3028]1101if (soil_pem) then
[3149]1102    call interpol_TI_PEM2PCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
[3426]1103    tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen)
[3172]1104#ifndef CPP_STD
[3181]1105    flux_geo = fluxgeo
[3172]1106#endif
[3028]1107endif
[3426]1108deallocate(tsoil_anom)
[2779]1109
[2835]1110! III_a.4 Pressure (for start)
[3149]1111ps = ps*global_avg_press_new/global_avg_press_PCM
1112ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM
[2794]1113
[2835]1114! III_a.5 Tracer (for start)
[3028]1115allocate(zplev_new(ngrid,nlayer + 1))
[2835]1116
[3028]1117do l = 1,nlayer + 1
[3149]1118    zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM
[3028]1119enddo
[2835]1120
[3028]1121do nnq = 1,nqtot
1122    if (noms(nnq) /= "co2") then
1123        do l = 1,llm - 1
1124            do ig = 1,ngrid
[3149]1125                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]1126            enddo
1127            q(:,llm,nnq) = q(:,llm - 1,nnq)
1128        enddo
1129    else
1130        do l = 1,llm - 1
1131            do ig = 1,ngrid
[3149]1132                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)) &
1133                              + ((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]1134            enddo
1135            q(:,llm,nnq) = q(:,llm - 1,nnq)
1136        enddo
1137    endif
1138enddo
[2835]1139
[3096]1140! Conserving the tracers mass for PCM start files
[3028]1141do nnq = 1,nqtot
1142    do ig = 1,ngrid
1143        do l = 1,llm - 1
1144            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]1145                extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
1146                q(ig,l,nnq) = 1.
1147                q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
[3028]1148                write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
[2835]1149           endif
[3028]1150            if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
1151        enddo
1152    enddo
1153enddo
[2779]1154
[3498]1155if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
[2779]1156
1157!------------------------
[3028]1158! III Output
[3317]1159!    III_b Write "restart.nc" and "restartfi.nc"
[3028]1160!------------------------
[3317]1161! III_b.1 Write "restart.nc"
[3042]1162ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
[3028]1163pday = day_ini
[3042]1164ztime_fin = time_phys
[2779]1165
[3028]1166allocate(p(ip1jmp1,nlayer + 1))
[2980]1167#ifndef CPP_1D
[3442]1168    ! Correction on teta due to surface pressure changes
1169    do l = 1,nlayer
1170        do i = 1,ip1jmp1
[3492]1171            teta(i,l)= teta(i,l)*(ps0(i)/ps(i))**kappa
[3442]1172        enddo
1173    enddo
1174    ! Correction on atmospheric pressure
1175    call pression(ip1jmp1,ap,bp,ps,p)
1176    ! Correction on the mass of atmosphere
[3028]1177    call massdair(p,masse)
[3317]1178    call dynredem0("restart.nc",day_ini,phis)
1179    call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps)
1180    write(*,*) "restart.nc has been written"
[2980]1181#else
[3363]1182    call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
1183    write(*,*) "restart1D.txt has been written"
[2980]1184#endif
1185
[3317]1186! III_b.2 Write the "restartfi.nc"
[2842]1187#ifndef CPP_STD
[3317]1188    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
[3327]1189                  nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
[3028]1190                  inertiedat,def_slope,subslope_dist)
[3327]1191    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,      &
[3114]1192                  ptimestep,ztime_fin,tsurf,tsoil,inertiesoil,        &
1193                  albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
[3130]1194                  wstar,watercap,perennial_co2ice)
[2842]1195#else
[3317]1196    call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
[3327]1197                  nlayer,nq,ptimestep,pday,time_phys,cell_area,    &
[3028]1198                  albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[3327]1199    call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil,       &
[3114]1200                  ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
1201                  cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab,   &
1202                  tsea_ice,sea_ice)
[2842]1203#endif
[3317]1204write(*,*) "restartfi.nc has been written"
[2842]1205
[2794]1206!------------------------
1207! III Output
[3161]1208!    III_c Write the "restartpem.nc"
[2794]1209!------------------------
[3319]1210if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings)
[3206]1211call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
[3088]1212call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
[3493]1213             TI_PEM,icetable_depth,icetable_thickness,ice_porefilling,   &
[3297]1214             co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
[3088]1215write(*,*) "restartpem.nc has been written"
[2779]1216
[3498]1217call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
[3149]1218
[3498]1219write(*,*) "The PEM has run for", i_myear_leg, "Martian years."
[3039]1220write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
1221write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
[3498]1222write(*,*) "PEM: so far, so good!"
[2794]1223
[3319]1224if (layering_algo) then
1225    do islope = 1,nslope
1226        do i = 1,ngrid
1227            call del_layering(stratif(i,islope))
1228        enddo
[3297]1229    enddo
[3319]1230    deallocate(new_str,new_lag1,new_lag2,current1,current2)
1231endif
[3149]1232deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
[3367]1233deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
1234deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
[3525]1235deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
1236deallocate(icetable_thickness_old,ice_porefilling_old)
[3330]1237deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
[3028]1238!----------------------------- END OUTPUT ------------------------------
[2897]1239
[2779]1240END PROGRAM pem
Note: See TracBrowser for help on using the repository browser.