Context Navigation

pem.F90 @ 3327

Last change on this file since 3327 was 3327, checked in by jbclement, 14 months ago

PEM:

Update of some default parameters;
Correction of a bug when the ice table was too low for the PEM subsurface discretization;
Correction of the computation for the change of sublimating ice necessary to the stopping criteria (the mistake was introduced in r3149) + simplification of the algorithm + renaming of variables with more explicit names;
Cleaning of "soil_thermalproperties_mod.F90".

JBC

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

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

Download in other formats:

Original Format