Context Navigation

pem.F90 @ 3586

Last change on this file since 3586 was 3584, checked in by jbclement, 6 months ago

PEM:

Albedo is now updated only at the end of the PEM (and not at every iteration) + Correct way to set it taking into account CO2/H2O ice and frost.
Cosmetic cleanings.

JBC

File size: 62.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format