Context Navigation

pem.F90 @ 3563

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

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

Download in other formats:

Original Format