Context Navigation

pem.F90 @ 3790

Last change on this file since 3790 was 3789, checked in by jbclement, 6 weeks ago
PEM: Few corrections for the layering algorithm, in particular when a dust lag layer is created. JBC
File size: 71.4 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
[3778]	672	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
	673	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
[3770]	674	enddo
	675	enddo
[3782]	676	! We put the sublimating tendency coming from subsurface ice into the overall tendency
	677	where (zdqsdif_ssi_tot < 0.)
	678	d_h2oice = zdqsdif_ssi_tot
	679	end where
[3770]	680	endif
[3384]	681	do i = 1,ngrid
[3779]	682	totmass_atmco2_ini = totmass_atmco2_ini + cell_area(i)*ps_avg(i)/g
[3384]	683	do islope = 1,nslope
[3779]	684	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]	685	if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true.
[3498]	686	if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then
[3571]	687	is_co2ice_sublim_ini(i,islope) = .true.
	688	co2ice_sublim_surf_ini = co2ice_sublim_surf_ini + cell_area(i)*subslope_dist(i,islope)
[3384]	689	endif
[3782]	690	if (d_h2oice(i,islope) < 0.) then
	691	if (h2o_ice(i,islope) > 0.) then
	692	is_h2oice_sublim_ini(i,islope) = .true.
	693	h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope)
	694	else if (h2o_ice_depth(i,islope) > 0.) then
	695	is_h2oice_sublim_ini(i,islope) = .true.
	696	endif
[3384]	697	endif
	698	enddo
	699	enddo
[3603]	700	write(,) "Total initial surface of CO2 ice sublimating =", co2ice_sublim_surf_ini
	701	write(,) "Total initial surface of H2O ice sublimating =", h2oice_ini_surf
[3384]	702
[3779]	703	totmass_adsco2_ini = 0.
	704	totmass_adsh2o = 0.
[3028]	705	if (adsorption_pem) then
	706	do ig = 1,ngrid
[3070]	707	do islope = 1,nslope
[3028]	708	do l = 1,nsoilmx_PEM - 1
[3571]	709	if (l == 1) then
[3779]	710	totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
[3028]	711	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
[3779]	712	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
[3028]	713	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
[3264]	714	else
[3779]	715	totmass_adsco2_ini = totmass_adsco2_ini + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l-1)) &
[3264]	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) - layer_PEM(l-1)) &
[3264]	718	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	719	endif
[3028]	720	enddo
[2961]	721	enddo
[3028]	722	enddo
[3779]	723	totmass_adsco2 = totmass_adsco2_ini
	724	write(,) "Tot mass of CO2 in the regolith =", totmass_adsco2
	725	write(,) "Tot mass of H2O in the regolith =", totmass_adsh2o
[3028]	726	endif ! adsorption
[2794]	727
[2835]	728	!------------------------
[3028]	729	! I Initialization
[2835]	730	! I_i Compute orbit criterion
[3028]	731	!------------------------
[2842]	732	#ifndef CPP_STD
[3050]	733	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
[2842]	734	#else
[3050]	735	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
[2842]	736	#endif
[2794]	737
[3498]	738	n_myear_leg = Max_iter_pem
	739	if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg)
[3403]	740
[3028]	741	!-------------------------- END INITIALIZATION -------------------------
[2794]	742
[3028]	743	!-------------------------------- RUN ----------------------------------
[2794]	744	!------------------------
	745	! II Run
[3065]	746	! II_a Update pressure, ice and tracers
[2794]	747	!------------------------
[3616]	748	write(,)
[3673]	749	write(,) '******* PEM cycle *******'
[3498]	750	i_myear_leg = 0
[3149]	751	stopPEM = 0
[3319]	752	if (layering_algo) then
[3784]	753	allocate(h2o_ice_depth_old(ngrid,nslope),new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope))
[3319]	754	new_str = .true.
[3770]	755	new_lag = .true.
[3319]	756	do islope = 1,nslope
	757	do ig = 1,ngrid
[3778]	758	current(ig,islope)%p => layerings_map(ig,islope)%top
[3319]	759	enddo
[3297]	760	enddo
[3319]	761	endif
[2794]	762
[3498]	763	do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear)
[2835]	764	! II.a.1. Compute updated global pressure
[3673]	765	write(,'(a,f10.2)') ' *** Iteration of the PEM leg (Martian years): ', i_myear_leg + 1
[3603]	766	write(,) "> Updating the surface pressure"
[3571]	767	ps_avg_global_old = ps_avg_global_new
[3028]	768	do i = 1,ngrid
	769	do islope = 1,nslope
[3611]	770	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]	771	enddo
	772	enddo
	773	if (adsorption_pem) then
	774	do i = 1,ngrid
[3611]	775	ps_avg_global_new = ps_avg_global_new - gcell_area(i)delta_co2_adsorbed(i)/total_surface
[3050]	776	enddo
[3028]	777	endif
[3680]	778	ps_avg = ps_avg*ps_avg_global_new/ps_avg_global_old
[3603]	779	write(,) 'Global average pressure old time step:',ps_avg_global_old
	780	write(,) 'Global average pressure new time step:',ps_avg_global_new
[2835]	781
[3571]	782	! II.a.2. Pressure timeseries (the values are deleted when unused because of big memory consumption)
[3603]	783	write(,) "> Updating the surface pressure timeseries for the new pressure"
[3571]	784	allocate(zplev_timeseries_old(ngrid,nlayer + 1,timelen))
[3028]	785	do l = 1,nlayer + 1
	786	do ig = 1,ngrid
[3571]	787	zplev_timeseries_old(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
[3028]	788	enddo
	789	enddo
[3673]	790	ps_timeseries(:,:) = ps_timeseries(:,:)*ps_avg_global_new/ps_avg_global_old
[3603]	791	write(,) "> Updating the pressure levels timeseries for the new pressure"
[3571]	792	allocate(zplev_timeseries_new(ngrid,nlayer + 1,timelen))
[3028]	793	do l = 1,nlayer + 1
	794	do ig = 1,ngrid
[3571]	795	zplev_timeseries_new(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
[3028]	796	enddo
	797	enddo
[2779]	798
[3571]	799	! II.a.3. Tracers timeseries
[3603]	800	write(,) "> Updating the tracer VMR timeseries for the new pressure"
[3591]	801	allocate(vmr_co2_PEM_phys(ngrid,timelen))
[3028]	802	l = 1
	803	do ig = 1,ngrid
	804	do t = 1,timelen
[3571]	805	! H2O
	806	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))/ &
	807	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
[3143]	808	if (q_h2o_PEM_phys(ig,t) < 0) then
	809	q_h2o_PEM_phys(ig,t) = 1.e-30
	810	else if (q_h2o_PEM_phys(ig,t) > 1) then
	811	q_h2o_PEM_phys(ig,t) = 1.
	812	endif
[3571]	813	! CO2
	814	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))/ &
	815	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
	816	+ ((zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t)) &
	817	- (zplev_timeseries_old(ig,l,t) - zplev_timeseries_old(ig,l + 1,t)))/ &
	818	(zplev_timeseries_new(ig,l,t) - zplev_timeseries_new(ig,l + 1,t))
[3028]	819	if (q_co2_PEM_phys(ig,t) < 0) then
	820	q_co2_PEM_phys(ig,t) = 1.e-30
[3143]	821	else if (q_co2_PEM_phys(ig,t) > 1) then
[3028]	822	q_co2_PEM_phys(ig,t) = 1.
	823	endif
[3571]	824	mmean = 1./(A*q_co2_PEM_phys(ig,t) + B)
[3149]	825	vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
[3028]	826	enddo
	827	enddo
[3571]	828	deallocate(zplev_timeseries_new,zplev_timeseries_old)
[2794]	829
[3028]	830	!------------------------
[2835]	831	! II Run
[3149]	832	! II_b Evolution of ice
[3028]	833	!------------------------
[3591]	834	allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
[3319]	835	if (layering_algo) then
[3784]	836	h2o_ice_depth_old = h2o_ice_depth
[3319]	837	do islope = 1,nslope
	838	do ig = 1,ngrid
[3778]	839	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]	840	co2_ice(ig,islope) = 0.
	841	h2o_ice(ig,islope) = 0.
	842	h2o_ice_depth(ig,islope) = 0.
[3778]	843	if (is_co2ice_str(layerings_map(ig,islope)%top)) then
	844	co2_ice(ig,islope) = layerings_map(ig,islope)%top%h_co2ice
	845	else if (is_h2oice_str(layerings_map(ig,islope)%top)) then
	846	h2o_ice(ig,islope) = layerings_map(ig,islope)%top%h_h2oice
	847	else
[3789]	848	call subsurface_ice_layering(layerings_map(ig,islope),h2o_ice_depth(ig,islope),h2o_ice(ig,islope))
[3778]	849	endif
[3319]	850	enddo
[3297]	851	enddo
[3789]	852	call print_layering(layerings_map(1,1))
[3553]	853	else
	854	zlag = 0.
[3770]	855	call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
	856	call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
[3319]	857	endif
[2794]	858
	859	!------------------------
	860	! II Run
[3149]	861	! II_c Flow of glaciers
[2794]	862	!------------------------
[3591]	863	allocate(flag_co2flow(ngrid,nslope),flag_h2oflow(ngrid,nslope))
[3571]	864	if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys, &
[3591]	865	ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow)
	866	if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow)
[3770]	867	if (layering_algo) then
	868	do islope = 1,nslope
	869	do ig = 1,ngrid
[3778]	870	layerings_map(ig,islope)%top%h_co2ice = co2_ice(ig,islope)
	871	layerings_map(ig,islope)%top%h_h2oice = h2o_ice(ig,islope)
[3770]	872	enddo
	873	enddo
	874	endif
[3065]	875
[2794]	876	!------------------------
	877	! II Run
[2835]	878	! II_d Update surface and soil temperatures
[2794]	879	!------------------------
[2835]	880	! II_d.1 Update Tsurf
[3603]	881	write(,) "> Updating surface temperature"
[3028]	882	do ig = 1,ngrid
	883	do islope = 1,nslope
[3571]	884	! CO2 ice disappeared so we look for the closest point without CO2 ice
	885	if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then
[3331]	886	co2ice_disappeared(ig,islope) = .true.
[3616]	887	if (latitude_deg(ig) > 0.) then ! North hemisphere
	888	outer1: do ig_loop = ig,ngrid ! Go towards equator
	889	do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be equator-ward)
[3330]	890	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
[3367]	891	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
[3571]	892	exit outer1
[3028]	893	endif
	894	enddo
[3616]	895	do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be pole-ward)
	896	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
	897	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
	898	exit outer1
	899	endif
	900	enddo
[3571]	901	enddo outer1
[3616]	902	else ! South hemisphere
	903	outer2: do ig_loop = ig,1,-1 ! Go towards equator
	904	do islope_loop = islope + 1,nslope ! Go over the slopes (forward numbering - should be equator-ward)
[3330]	905	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
[3367]	906	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
[3571]	907	exit outer2
[3028]	908	endif
	909	enddo
[3616]	910	do islope_loop = islope - 1,1,-1 ! Go over the slopes (backward numbering - should be pole-ward)
	911	if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then
	912	tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop)
	913	exit outer2
	914	endif
	915	enddo
[3571]	916	enddo outer2
[2835]	917	endif
[3616]	918	else if (co2_ice(ig,islope) > 1.e-10 .and. d_co2ice(ig,islope) > 1.e-10) then ! Put tsurf as Tcond CO2
[3571]	919	call computeTcondCO2(timelen,ngrid,nslope,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_ini,ps_avg_global_new,tsurf_avg)
[2835]	920	endif
	921	enddo
[3028]	922	enddo
[2794]	923
[3028]	924	if (soil_pem) then
[3553]	925	! II_d.2 Shifting soil temperature to surface
[3770]	926	call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM)
[2794]	927
[3553]	928	! II_d.3 Update soil temperature
[3603]	929	write(,)"> Updating soil temperature profile"
[3610]	930	allocate(tsoil_avg_old(ngrid,nsoilmx_PEM),tsoil_PEM_timeseries_old(ngrid,nsoilmx_PEM,nslope,timelen))
	931	tsoil_PEM_timeseries_old = tsoil_PEM_timeseries
[3426]	932	do islope = 1,nslope
[3571]	933	tsoil_avg_old = tsoil_PEM(:,:,islope)
[3426]	934	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
	935	call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope))
[2794]	936
[3028]	937	do t = 1,timelen
	938	do ig = 1,ngrid
	939	do isoil = 1,nsoilmx_PEM
[3609]	940	! Update of soil temperature timeseries which is needed to compute the water soil density timeseries
	941	tsoil_PEM_timeseries(ig,isoil,islope,t) = tsoil_PEM_timeseries(ig,isoil,islope,t)*tsoil_PEM(ig,isoil,islope)/tsoil_avg_old(ig,isoil)
	942	! Update of watersoil density
[3571]	943	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]	944	if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1)
[3028]	945	enddo
	946	enddo
	947	enddo
	948	enddo
[3367]	949	watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen
[3571]	950	deallocate(tsoil_avg_old)
[2888]	951
[3553]	952	! II_d.4 Update the ice table
[3610]	953	allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope),icetable_depth_old(ngrid,nslope))
[3170]	954	if (icetable_equilibrium) then
[3603]	955	write(,) "> Updating ice table (equilibrium method)"
[3525]	956	icetable_thickness_old = icetable_thickness
[3493]	957	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
[3525]	958	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]	959	else if (icetable_dynamic) then
[3603]	960	write(,) "> Updating ice table (dynamic method)"
[3525]	961	ice_porefilling_old = ice_porefilling
[3610]	962	icetable_depth_old = icetable_depth
[3493]	963	allocate(porefill(nsoilmx_PEM))
	964	do ig = 1,ngrid
	965	do islope = 1,nslope
[3571]	966	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]	967	icetable_depth(ig,islope) = ssi_depth
	968	ice_porefilling(ig,:,islope) = porefill
	969	enddo
	970	enddo
	971	deallocate(porefill)
[3673]	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
[3486]	973	endif
[3591]	974	deallocate(icetable_thickness_old,ice_porefilling_old)
[3486]	975
[3553]	976	! II_d.5 Update the soil thermal properties
[3571]	977	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]	978
[3553]	979	! II_d.6 Update the mass of the regolith adsorbed
[3779]	980	totmass_adsco2 = 0.
	981	totmass_adsh2o = 0.
[3028]	982	if (adsorption_pem) then
[3571]	983	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,h2o_ice,co2_ice, &
	984	tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
[3553]	985	h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
[3028]	986	do ig = 1,ngrid
[3493]	987	do islope = 1,nslope
[3532]	988	do l = 1,nsoilmx_PEM
[3493]	989	if (l == 1) then
[3779]	990	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
[3264]	991	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
[3779]	992	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l)) &
[3264]	993	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	994	else
[3779]	995	totmass_adsco2 = totmass_adsco2 + co2_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
[3264]	996	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
[3779]	997	totmass_adsh2o = totmass_adsh2o + h2o_adsorbed_phys(ig,l,islope)(layer_PEM(l) - layer_PEM(l - 1)) &
[3264]	998	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	999	endif
[3028]	1000	enddo
	1001	enddo
	1002	enddo
[3779]	1003	write(,) "Total mass of CO2 in the regolith =", totmass_adsco2
	1004	write(,) "Total mass of H2O in the regolith =", totmass_adsh2o
[3028]	1005	endif
	1006	endif !soil_pem
[3666]	1007	deallocate(zshift_surf,zlag)
[3028]	1008
[2794]	1009	!------------------------
	1010	! II Run
[3088]	1011	! II_e Outputs
[2794]	1012	!------------------------
[3571]	1013	call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/ps_avg_global_new/))
[3088]	1014	do islope = 1,nslope
	1015	write(str2(1:2),'(i2.2)') islope
[3181]	1016	call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope))
	1017	call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope))
[3498]	1018	call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope))
	1019	call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope))
[3591]	1020	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,real(flag_co2flow(:,islope)))
	1021	call writediagpem(ngrid,'Flow_h2oice_slope'//str2,'H2O ice flow','Boolean',2,real(flag_h2oflow(:,islope)))
[3571]	1022	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_avg(:,islope))
[3339]	1023	if (icetable_equilibrium) then
[3493]	1024	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
[3537]	1025	call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope))
[3490]	1026	else if (icetable_dynamic) then
[3493]	1027	call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
[3486]	1028	endif
	1029
[3339]	1030	if (soil_pem) then
[3571]	1031	call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil','K',3,tsoil_PEM(:,:,islope))
	1032	call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI','K',3,TI_PEM(:,:,islope))
[3543]	1033	if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
[3171]	1034	if (adsorption_pem) then
[3553]	1035	call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope))
	1036	call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope))
[3532]	1037	endif
[3171]	1038	endif
[3088]	1039	enddo
[3591]	1040	deallocate(flag_co2flow,flag_h2oflow)
[3088]	1041
[3779]	1042	! Checking mass balance for CO2
	1043	totmass_co2ice = 0.
	1044	totmass_atmco2 = 0.
	1045	do ig = 1,ngrid
	1046	totmass_atmco2 = totmass_atmco2 + cell_area(ig)*ps_avg(ig)/g
	1047	do islope = 1,nslope
	1048	totmass_co2ice = totmass_co2ice + co2_ice(ig,islope)cell_area(ig)subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
	1049	enddo
	1050	enddo
	1051	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), ' %'
	1052	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
	1053	write(,) ' /!\ Warning: mass balance is not conseved!'
	1054	write(,'(a,f8.3,a)') ' Atmospheric CO2 mass balance = ', 100.(totmass_atmco2 - totmass_atmco2_ini)/totmass_atmco2_ini, ' %'
	1055	write(,'(a,f8.3,a)') ' CO2 ice mass balance = ', 100.(totmass_co2ice - totmass_co2ice_ini)/totmass_co2ice_ini, ' %'
	1056	write(,'(a,f8.3,a)') ' Adsorbed CO2 mass balance = ', 100.(totmass_adsco2 - totmass_adsco2_ini)/totmass_adsco2_ini, ' %'
	1057	endif
	1058
[3088]	1059	!------------------------
	1060	! II Run
	1061	! II_f Update the tendencies
	1062	!------------------------
[3571]	1063	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]	1064	write(,) "> Updating the H2O sub-surface ice tendency due to lag layer"
[3770]	1065	if (layering_algo) then
	1066	do ig = 1,ngrid
	1067	do islope = 1,nslope
[3789]	1068	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]	1069	enddo
[3707]	1070	enddo
[3770]	1071	!~ else
	1072	!~ do ig = 1,ngrid
	1073	!~ do islope = 1,nslope
	1074	!~ 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))
	1075	!~ enddo
	1076	!~ enddo
	1077	endif
[3666]	1078	if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old)
	1079	deallocate(vmr_co2_PEM_phys)
[3711]	1080	write(,) "> Updating the H2O sub-surface ice depth"
[2794]	1081
[2835]	1082	!------------------------
	1083	! II Run
[3088]	1084	! II_g Checking the stopping criterion
[2835]	1085	!------------------------
[3603]	1086	write(,) "> Checking the stopping criteria"
[3571]	1087	call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
	1088	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]	1089	i_myear_leg = i_myear_leg + dt
	1090	i_myear = i_myear + dt
	1091	if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5
[3430]	1092	if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6
[3389]	1093	call system_clock(c2)
[3430]	1094	if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7
[3432]	1095	if (stopPEM > 0) then
[3149]	1096	select case (stopPEM)
	1097	case(1)
[3673]	1098	write(,'(a,i0,a)') " *** STOPPING because surface of h2o ice sublimating is too low: ", stopPEM, ". See message above."
[3149]	1099	case(2)
[3673]	1100	write(,'(a,i0,a)') " *** STOPPING because tendencies on h2o ice = 0: ", stopPEM, ". See message above."
[3149]	1101	case(3)
[3673]	1102	write(,'(a,i0,a)') " *** STOPPING because surface of co2 ice sublimating is too low: ", stopPEM, ". See message above."
[3149]	1103	case(4)
[3673]	1104	write(,'(a,i0,a)') " *** STOPPING because surface global pressure changed too much: ", stopPEM, ". See message above."
[3149]	1105	case(5)
[3673]	1106	write(,'(a,i0)') " *** STOPPING because maximum number of iterations is reached (possibly due to orbital parameters): ", stopPEM
[3149]	1107	case(6)
[3673]	1108	write(,'(a,i0)') " *** STOPPING because maximum number of Martian years to be simulated is reached: ", stopPEM
[3363]	1109	case(7)
[3673]	1110	write(,'(a,i0)') " *** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM
[3149]	1111	case default
[3673]	1112	write(,'(a,i0)') " *** STOPPING with unknown stopping criterion code: ", stopPEM
[3149]	1113	end select
[2779]	1114	exit
[3028]	1115	else
[3673]	1116	write(,'(a,f10.2,a)') ' *** The chained simulation has run for ',i_myear,' Martian years.'
	1117	write(,) '**** The PEM can continue!'
	1118	write(,) '****'
[3028]	1119	endif
[3149]	1120	enddo ! big time iteration loop of the pem
[3591]	1121	deallocate(vmr_co2_PCM,q_co2_PEM_phys,q_h2o_PEM_phys,delta_co2_adsorbed)
[3598]	1122	deallocate(watersoil_density_PEM_avg,watersurf_density_avg)
[3591]	1123	deallocate(ps_timeseries,tsoil_PEM_timeseries,watersoil_density_PEM_timeseries)
	1124	deallocate(co2ice_disappeared,delta_h2o_adsorbed,delta_h2o_icetablesublim)
	1125	deallocate(d_co2ice,d_co2ice_ini,d_h2oice)
	1126	deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini)
[3784]	1127	if (layering_algo) deallocate(h2o_ice_depth_old,new_str,new_lag,current)
[3028]	1128	!------------------------------ END RUN --------------------------------
[2779]	1129
[3028]	1130	!------------------------------- OUTPUT --------------------------------
[2794]	1131	!------------------------
	1132	! III Output
[3770]	1133	! III_a Update surface values for the PCM start files
[2794]	1134	!------------------------
[3616]	1135	write(,)
[3673]	1136	write(,) '******* PEM finalization *******'
[3778]	1137	! III_a.1 Ice update for start file
[3603]	1138	write(,) '> Reconstructing perennial ice and frost for the PCM'
[3149]	1139	watercap = 0.
[3770]	1140	perennial_co2ice = co2_ice
[3028]	1141	do ig = 1,ngrid
[3159]	1142	! H2O ice metamorphism
[3705]	1143	!if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_h2oice_threshold) then
	1144	! h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold
	1145	! qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold
	1146	!endif
[3159]	1147
	1148	! Is H2O ice still considered as an infinite reservoir for the PCM?
[3770]	1149	if (sum(h2o_ice(ig,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > inf_h2oice_threshold) then
	1150	! There is enough ice to be considered as an infinite reservoir
[3149]	1151	watercaptag(ig) = .true.
	1152	else
[3770]	1153	! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost
[3149]	1154	watercaptag(ig) = .false.
[3770]	1155	qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:)
[3149]	1156	h2o_ice(ig,:) = 0.
[3028]	1157	endif
[3159]	1158
	1159	! CO2 ice metamorphism
[3705]	1160	!if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)subslope_dist(ig,:)/cos(pidef_slope_mean(:)/180.)) > metam_co2ice_threshold) then
	1161	! perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold
	1162	! qsurf(ig,igcm_co2,:) = metam_co2ice_threshold
	1163	!endif
[3028]	1164	enddo
[2888]	1165
[3770]	1166	! III.a.3. Tsurf update for start file
[3603]	1167	write(,) '> Reconstructing the surface temperature for the PCM'
[3571]	1168	tsurf = tsurf_avg + tsurf_dev
[3591]	1169	deallocate(tsurf_dev)
[3571]	1170
[3770]	1171	! III_a.4 Tsoil update for start file
[3028]	1172	if (soil_pem) then
[3603]	1173	write(,) '> Reconstructing the soil temperature profile for the PCM'
[3532]	1174	inertiesoil = TI_PEM(:,:nsoilmx,:)
[3571]	1175	! Tsurf has evolved and so the soil temperature profile needs to be adapted to match this new value
	1176	do isoil = 1,nsoilmx
	1177	tsoil_dev(:,isoil,:) = tsoil_dev(:,isoil,:)*(tsurf_avg(:,:) - tsoil_PEM(:,1,:))/tsoil_dev(:,1,:)
	1178	enddo
	1179	tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_dev
[3172]	1180	#ifndef CPP_STD
[3181]	1181	flux_geo = fluxgeo
[3172]	1182	#endif
[3028]	1183	endif
[3591]	1184	deallocate(tsurf_avg,tsoil_dev)
[2779]	1185
[3770]	1186	! III_a.5 Pressure update for start file
[3603]	1187	write(,) '> Reconstructing the pressure for the PCM'
[3591]	1188	allocate(ps_start(ngrid))
[3680]	1189	! The pressure deviation is rescaled as well to avoid disproportionate oscillations in case of huge average pressure drop
	1190	ps_start = ps_avg + ps_dev*ps_avg_global_new/ps_avg_global_ini
[3591]	1191	deallocate(ps_avg,ps_dev)
[2794]	1192
[3770]	1193	! III_a.6 Tracers update for start file
[3603]	1194	write(,) '> Reconstructing the tracer VMR for the PCM'
[3571]	1195	allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1))
[3028]	1196	do l = 1,nlayer + 1
[3571]	1197	zplev_start0(:,l) = ap(l) + bp(l)*ps_start0
	1198	zplev_new(:,l) = ap(l) + bp(l)*ps_start
[3028]	1199	enddo
[2835]	1200
[3028]	1201	do nnq = 1,nqtot
	1202	if (noms(nnq) /= "co2") then
	1203	do l = 1,llm - 1
	1204	do ig = 1,ngrid
[3571]	1205	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]	1206	enddo
	1207	q(:,llm,nnq) = q(:,llm - 1,nnq)
	1208	enddo
	1209	else
	1210	do l = 1,llm - 1
	1211	do ig = 1,ngrid
[3571]	1212	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)) &
	1213	+ ((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]	1214	enddo
	1215	q(:,llm,nnq) = q(:,llm - 1,nnq)
	1216	enddo
	1217	endif
	1218	enddo
[3591]	1219	deallocate(zplev_start0)
[2835]	1220
[3571]	1221	! Conserving the tracers mass for start file
[3028]	1222	do nnq = 1,nqtot
	1223	do ig = 1,ngrid
	1224	do l = 1,llm - 1
	1225	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]	1226	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
	1227	q(ig,l,nnq) = 1.
	1228	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
[3028]	1229	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
[3770]	1230	endif
[3028]	1231	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
	1232	enddo
	1233	enddo
	1234	enddo
[3591]	1235	deallocate(zplev_new)
[2779]	1236
[3770]	1237	! III_a.7 Albedo update for start file
[3603]	1238	write(,) '> Reconstructing the albedo for the PCM'
[3584]	1239	do ig = 1,ngrid
	1240	if (latitude(ig) < 0.) then
	1241	icap = 2 ! Southern hemisphere
	1242	else
	1243	icap = 1 ! Northern hemisphere
	1244	endif
[3603]	1245	do islope = 1,nslope
[3584]	1246	! Bare ground
	1247	albedo(ig,:,islope) = albedodat(ig)
	1248	emis(ig,islope) = emissiv
	1249
	1250	! CO2 ice/frost is treated after H20 ice/frost because it is considered dominant
	1251	! H2O ice
	1252	if (h2o_ice(ig,islope) > 0.) then
	1253	albedo(ig,:,islope) = albedo_h2o_cap
	1254	emis(ig,islope) = 1.
	1255	endif
	1256	! CO2 ice
	1257	if (co2_ice(ig,islope) > 0.) then
	1258	albedo(ig,:,islope) = albedo_perennialco2(icap)
	1259	emis(ig,islope) = emisice(icap)
	1260	endif
	1261	! H2O frost
	1262	if (qsurf(ig,igcm_h2o_ice,islope) > 0.) then
	1263	albedo(ig,:,islope) = albedo_h2o_frost
	1264	emis(ig,islope) = 1.
	1265	endif
	1266	! CO2 frost
	1267	if (qsurf(ig,igcm_co2,islope) > 0.) then
	1268	albedo(ig,:,islope) = albedice(icap)
	1269	emis(ig,islope) = emisice(icap)
	1270	endif
	1271	enddo
	1272	enddo
	1273
[3770]	1274	! III_a.8 Orbital parameters update for start file
[3603]	1275	write(,) '> Setting the new orbital parameters'
[3498]	1276	if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg)
[2779]	1277
	1278	!------------------------
[3028]	1279	! III Output
[3317]	1280	! III_b Write "restart.nc" and "restartfi.nc"
[3028]	1281	!------------------------
[3317]	1282	! III_b.1 Write "restart.nc"
[3042]	1283	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
[3028]	1284	pday = day_ini
[3042]	1285	ztime_fin = time_phys
[2980]	1286	#ifndef CPP_1D
[3603]	1287	write(,) '> Writing "restart.nc"'
[3442]	1288	! Correction on teta due to surface pressure changes
[3603]	1289	allocate(pdyn(ip1jmp1))
	1290	call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start0/ps_start,pdyn)
	1291	do i = 1,ip1jmp1
[3779]	1292	teta(i,:) = teta(i,:)pdyn(i)*rcp
[3442]	1293	enddo
	1294	! Correction on atmospheric pressure
[3603]	1295	allocate(p(ip1jmp1,nlayer + 1))
	1296	call gr_fi_dyn(1,ngrid,iip1,jjp1,ps_start,pdyn)
	1297	call pression(ip1jmp1,ap,bp,pdyn,p)
[3442]	1298	! Correction on the mass of atmosphere
[3028]	1299	call massdair(p,masse)
[3317]	1300	call dynredem0("restart.nc",day_ini,phis)
[3628]	1301	call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,pdyn)
[3603]	1302	deallocate(ap,bp,p,pdyn)
[2980]	1303	#else
[3603]	1304	write(,) '> Writing "restart1D.txt"'
[3571]	1305	call writerestart1D('restart1D.txt',ps_start(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
[2980]	1306	#endif
[3591]	1307	deallocate(ps_start0,ps_start)
[2980]	1308
[3317]	1309	! III_b.2 Write the "restartfi.nc"
[3603]	1310	write(,) '> Writing "restartfi.nc"'
[2842]	1311	#ifndef CPP_STD
[3317]	1312	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
[3327]	1313	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
[3028]	1314	inertiedat,def_slope,subslope_dist)
[3327]	1315	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
[3114]	1316	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
	1317	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
[3130]	1318	wstar,watercap,perennial_co2ice)
[2842]	1319	#else
[3591]	1320	if (allocated(noms)) deallocate(noms)
	1321	deallocate(qsurf,tsurf,tsoil,emis,watercap,watercaptag,albedo,inertiesoil)
[3317]	1322	call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, &
[3327]	1323	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
[3028]	1324	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[3327]	1325	call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
[3114]	1326	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
	1327	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
	1328	tsea_ice,sea_ice)
[2842]	1329	#endif
	1330
[2794]	1331	!------------------------
	1332	! III Output
[3161]	1333	! III_c Write the "restartpem.nc"
[2794]	1334	!------------------------
[3603]	1335	write(,) '> Writing "restartpem.nc"'
[3778]	1336	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]	1337	call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
[3537]	1338	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
[3778]	1339	co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,layerings_map)
[2779]	1340
[3498]	1341	call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear)
[3149]	1342
[3616]	1343	write(,)
[3738]	1344	write(,) '**** PEM final information *****'
[3770]	1345	write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years."
	1346	write(,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myearconvert_years, " Earth years."
	1347	write(,'(a,f16.4,a)') " + The reached date is now ", (year_bp_ini + i_myear)convert_years, " Earth years."
[3616]	1348	write(,) "+ PEM: so far, so good!"
[3673]	1349	write(,) '************************************'
[2794]	1350
[3319]	1351	if (layering_algo) then
	1352	do islope = 1,nslope
	1353	do i = 1,ngrid
[3778]	1354	call del_layering(layerings_map(i,islope))
[3319]	1355	enddo
[3297]	1356	enddo
[3319]	1357	endif
[3591]	1358	deallocate(q,longitude,latitude,cell_area,tsoil_PEM)
[3778]	1359	deallocate(co2_ice,h2o_ice,layerings_map)
[3028]	1360	!----------------------------- END OUTPUT ------------------------------
[2897]	1361
[2779]	1362	END PROGRAM pem

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

Download in other formats:

Original Format