Context Navigation

pem.F90 @ 3599

Last change on this file since 3599 was 3599, checked in by jbclement, 34 hours ago
PEM: Few small optimizations. JBC
File size: 63.1 KB

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

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

Download in other formats:

Original Format