Context Navigation

pem.F90 @ 3831

Last change on this file since 3831 was 3809, checked in by jbclement, 3 weeks ago

PEM:

Correction to detect whether the stratum below the dust lag is made of ice or not.
Correction for the initialization of surface ice and 'h2o_ice_depth' according to the layerings map.
Correction to handle flow of glaciers with the layering algorithm.
Simplification of "recomp_tend_h2o", making it more robust.
Making the removing of a stratum more robust.
Plotting the orbital variations over the simulated time in "visu_evol_layering.py".
Lowering the value of dust tendency.
Addition of "_co2" in the variables name for the CO2 flux correction.

JBC

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

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

Download in other formats:

Original Format