Context Navigation

pem.F90 @ 3130

Last change on this file since 3130 was 3130, checked in by jbclement, 21 months ago

PEM:
The perennial co2 ice is now taken into account with co2 frost (qsurf) to compute the tendency and to make the update + Rework of how co2 frost is converted to perennial co2 ice at the end of the PEM run + Correction of the value of 'threshold_co2_frost2perennial' to correspond to 10 m + Perennial co2 ice is now handled outside 'paleoclimate' in "phyetat0_mod.F90" of the Mars PCM + Some cleanings.

/!\ Commit for the PEM management of co2 ice before a rework of ice management in the PEM!
JBC

File size: 59.6 KB

Rev	Line
[2779]	1	!------------------------
[3028]	2	! I Initialization
	3	! I_a READ run.def
[2835]	4	! I_b READ of start_evol.nc and starfi_evol.nc
	5	! I_c Subslope parametrisation
[3096]	6	! I_d READ PCM data and convert to the physical grid
[3028]	7	! I_e Initialization of the PEM variable and soil
[2835]	8	! I_f Compute tendencies & Save initial situation
	9	! I_g Save initial PCM situation
[3088]	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
[2835]	15	! II_b Evolution of the ice
[3028]	16	! II_c CO2 & H2O glaciers flows
[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
[3028]	24	! III_b Write restart_evol.nc and restartfi_evol.nc
[3088]	25	! III_c Write restartpem.nc
[2779]	26	!------------------------
	27
	28	PROGRAM pem
	29
[3076]	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
	34	use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
	35	use logic_mod, only: iflag_phys
	36	use mod_const_mpi, only: COMM_LMDZ
[3028]	37	use infotrac
[3076]	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: co2glaciers_evol, h2oglaciers_evol, co2glaciersflow, h2oglaciersflow
	42	use criterion_pem_stop_mod, only: criterion_waterice_stop, criterion_co2_stop
	43	use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, &
[3130]	44	m_noco2, threshold_water_frost2perennial, threshold_co2_frost2perennial
[3076]	45	use evol_co2_ice_s_mod, only: evol_co2_ice_s
	46	use evol_h2o_ice_s_mod, only: evol_h2o_ice_s
	47	use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &
	48	TI_PEM, inertiedat_PEM, & ! soil thermal inertia
	49	tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
[3096]	50	fluxgeo, & ! Geothermal flux for the PEM and PCM
[3076]	51	water_reservoir ! Water ressources
	52	use adsorption_mod, only: regolith_adsorption, adsorption_pem, & ! Bool to check if adsorption, main subroutine
	53	ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
	54	co2_adsorbded_phys, h2o_adsorbded_phys ! Mass of co2 and h2O adsorbded
	55	use time_evol_mod, only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
	56	use orbit_param_criterion_mod, only: orbit_param_criterion
	57	use recomp_orb_param_mod, only: recomp_orb_param
	58	use ice_table_mod, only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
	59	ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
	60	use soil_thermalproperties_mod, only: update_soil_thermalproperties
	61	use time_phylmdz_mod, only: daysec, dtphys, day_end
	62	use abort_pem_mod, only: abort_pem
	63	use soil_settings_PEM_mod, only: soil_settings_PEM
	64	use compute_tendencies_slope_mod, only: compute_tendencies_slope
[3096]	65	use info_PEM_mod, only: info_PEM
[3076]	66	use interpolate_TIPEM_TIGCM_mod, only: interpolate_TIPEM_TIGCM
[3096]	67	use nb_time_step_PCM_mod, only: nb_time_step_PCM
[3076]	68	use pemetat0_mod, only: pemetat0
[3096]	69	use read_data_PCM_mod, only: read_data_PCM
[3076]	70	use recomp_tend_co2_slope_mod, only: recomp_tend_co2_slope
	71	use soil_pem_compute_mod, only: soil_pem_compute
[3088]	72	use writediagpem_mod, only: writediagpem
[2985]	73
[2842]	74	#ifndef CPP_STD
[3114]	75	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil, nqsoil, qsoil
[3028]	76	use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
	77	albedodat, zmea, zstd, zsig, zgam, zthe, &
	78	hmons, summit, base,albedo_h2o_frost, &
[3130]	79	frost_albedo_threshold, emissiv, watercaptag, perennial_co2ice, &
[3032]	80	emisice, albedice
[3028]	81	use dimradmars_mod, only: totcloudfrac, albedo
	82	use dust_param_mod, only: tauscaling
	83	use tracer_mod, only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
	84	use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
[3096]	85	use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
[3130]	86	use paleoclimate_mod, only: albedo_perennialco2
[3082]	87	use comcstfi_h, only: pi, rad, g, cpp, mugaz, r
[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
[3082]	95	use comcstfi_mod, only: pi, rad, g, cpp, 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
[3076]	102	use time_phylmdz_mod, only: iphysiq, day_step
[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
[3096]	128	real :: time_phys ! Same as PCM
	129	real :: ptimestep ! Same as PCM
	130	real :: ztime_fin ! Same as PCM
[2794]	131
[3028]	132	! Variables to read start.nc
	133	character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing 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
	140	real, dimension(ip1jmp1) :: ps ! pression au sol
	141	real, dimension(:), allocatable :: ps_start_GCM ! (ngrid) pression au sol
	142	real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) ! pression au sol instantannées
	143	real, dimension(ip1jmp1,llm) :: masse ! masse d'air
	144	real, dimension(ip1jmp1) :: phis ! geopotentiel au sol
[3028]	145	real :: time_0
[2779]	146
[3028]	147	! Variables to read starfi.nc
[3065]	148	character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" ! Name of the file used for initialsing the PEM
[3068]	149	character(2) :: str2
	150	integer :: ncid, varid, status ! Variable for handling opening of files
	151	integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf 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
[3065]	156	real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi
	157	real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi
	158	real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi
[3028]	159	real :: Total_surface ! Total surface of the planet
[2897]	160
[3028]	161	! Variables for h2o_ice evolution
[3065]	162	real :: ini_surf_h2o ! Initial surface of sublimating h2o ice
	163	real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa]
	164	real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step
	165	real :: global_ave_press_new ! constant: Global average pressure of current time step
	166	real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2]
[3096]	167	real, dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the PCM [kg/m^2]
[3065]	168	real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
	169	real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
	170	logical :: STOPPING_water ! Logical: is the criterion (% of change in the surface of sublimating water ice) reached?
	171	logical :: STOPPING_1_water ! Logical: is there still water ice to sublimate?
	172	logical :: STOPPING_pressure ! Logical: is the criterion (% of change in the surface pressure) reached?
	173	integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
[3068]	174	real, save :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg]
[3096]	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]
[3065]	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
	182	real :: ini_surf_co2 ! Initial surface of sublimating co2 ice
	183	logical :: STOPPING_co2 ! Logical: is the criterion (% of change in the surface of sublimating co2 ice) reached?
	184	real, dimension(:,:), allocatable :: vmr_co2_gcm ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
	185	real, dimension(:,:), allocatable :: vmr_co2_pem_phys ! Physics x Times co2 volume mixing ratio used in the PEM
	186	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]
	187	real, dimension(:,:), allocatable :: co2frost_ini ! Initial amount of co2 frost (at the start of the PEM run)
	188	real, dimension(:,:), allocatable :: perennial_co2ice_ini ! Initial amoun of perennial co2 ice (at the start of the PEM run)
	189
	190
	191
[3028]	192	! Variables for slopes
[3065]	193	real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field: minimum of co2 ice at each point for the first year [kg/m^2]
	194	real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field: minimum of co2 ice at each point for the second year [kg/m^2]
	195	real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field: minimum of water ice at each point for the first year [kg/m^2]
	196	real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field: minimum of water ice at each point for the second year [kg/m^2]
[3096]	197	real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field: co2 ice given by the PCM [kg/m^2]
[3065]	198	real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
	199	real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field: Logical array indicating if there is water ice at initial state
	200	real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field: Logical array indicating if there is co2 ice at initial state
[3130]	201	real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year
	202	real, dimension(:,:), allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perennial co2 ice over a year in the PCM
	203	real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical point x slope field: Tendency of evolution of perennial h2o ice
[3065]	204	real, dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow
	205	real, dimension(:), allocatable :: flag_co2flow_mesh ! (ngrid) : Flag where there is a CO2 glacier flow
	206	real, dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow
	207	real, dimension(:), allocatable :: flag_h2oflow_mesh ! (ngrid) : Flag where there is a H2O glacier flow
[2779]	208
[3028]	209	! Variables for surface and soil
[3065]	210	real, dimension(:,:), allocatable :: tsurf_ave ! Physic x SLOPE field : Averaged Surface Temperature [K]
	211	real, dimension(:,:,:), allocatable :: tsoil_ave ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
	212	real, dimension(:,:,:), allocatable :: tsurf_GCM_timeseries ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
	213	real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
	214	real, dimension(:,:,:,:), allocatable :: tsoil_GCM_timeseries ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
[3096]	215	real, dimension(:,:), allocatable :: tsurf_ave_yr1 ! Physic x SLOPE field : Averaged Surface Temperature of first call of the PCM [K]
[3065]	216	real, dimension(:,:), allocatable :: TI_locslope ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI]
	217	real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: intermediate when computing Tsoil [K]
	218	real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
	219	real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
	220	real, dimension(:,:), allocatable :: watersurf_density_ave ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
	221	real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
	222	real, dimension(:,:,:), allocatable :: watersoil_density_PEM_ave ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
	223	real, dimension(:,:), allocatable :: Tsurfave_before_saved ! Surface temperature saved from previous time step [K]
	224	real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
	225	real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
	226	real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
	227	real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
	228	logical :: bool_sublim ! logical to check if there is sublimation or not
	229	real, dimension(:,:), allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
	230	real, dimension(:), allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
	231
[3028]	232	! Some variables for the PEM run
	233	real, parameter :: year_step = 1 ! timestep for the pem
	234	integer :: year_iter ! number of iteration
	235	integer :: year_iter_max ! maximum number of iterations before stopping
[3039]	236	integer :: i_myear ! Global number of Martian years of the chained simulations
	237	integer :: n_myear ! Maximum number of Martian years of the chained simulations
[3028]	238	real :: timestep ! timestep [s]
[3065]	239	real :: watercap_sum ! total mass of water cap [kg/m^2]
	240	real :: water_sum ! total mass of water in the mesh [kg/m^2]
[2779]	241
[2842]	242	#ifdef CPP_STD
[3065]	243	real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
	244	real :: albedo_h2o_frost ! albedo of h2o frost
	245	real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	246	real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	247	real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	248	real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	249	real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	250	real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case
	251	real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case
	252	real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case
	253	real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case
	254	real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
[3068]	255	logical, dimension(:), allocatable :: watercaptag ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
[3065]	256	real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case
	257	real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case
[2842]	258	#endif
	259
[2980]	260	#ifdef CPP_1D
[3065]	261	integer :: nsplit_phys
	262	integer, parameter :: jjm_value = jjm - 1
	263	integer :: ierr
	264
	265	! Dummy variables to use the subroutine 'init_testphys1d'
[3129]	266	logical :: therestart1D, therestartfi
[3068]	267	integer :: ndt, day0
	268	real :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau
	269	real, dimension(:), allocatable :: zqsat
	270	real, dimension(:,:,:), allocatable :: dq, dqdyn
	271	real, dimension(nlayer) :: play, w
	272	real, dimension(nlayer + 1) :: plev
[2980]	273	#else
[3065]	274	integer, parameter :: jjm_value = jjm
	275	real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart
	276	real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart
[2980]	277	#endif
	278
[3028]	279	! Loop variables
[3039]	280	integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil, icap
[2779]	281
[3028]	282	! Parallel variables
[2842]	283	#ifndef CPP_STD
[3028]	284	is_sequential = .true.
	285	is_parallel = .false.
	286	is_mpi_root = .true.
	287	is_omp_root = .true.
	288	is_master = .true.
[2842]	289	#endif
[2779]	290
[3065]	291	! Some constants
[3028]	292	day_ini = 0 ! test
	293	time_phys = 0. ! test
	294	ngrid = ngridmx
	295	A = (1/m_co2 - 1/m_noco2)
	296	B = 1/m_noco2
	297	year_day = 669
	298	daysec = 88775.
	299	timestep = year_day*daysec/year_step
[2794]	300
[3028]	301	!----------------------------- INITIALIZATION --------------------------
[2779]	302	!------------------------
[3028]	303	! I Initialization
	304	! I_a READ run.def
[2779]	305	!------------------------
[2980]	306	#ifndef CPP_1D
[3028]	307	dtphys = 0
	308	call conf_gcm(99,.true.)
	309	call infotrac_init
	310	nq = nqtot
	311	allocate(q(ip1jmp1,llm,nqtot))
[3065]	312	allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid))
[2980]	313	#else
[3068]	314	allocate(q(1,llm,nqtot))
[3065]	315	allocate(longitude(1),latitude(1),cell_area(1))
[3129]	316
	317	therestart1D = .false.
	318	inquire(file = 'start1D_evol.txt',exist = therestart1D)
	319	if (.not. therestart1D) then
	320	write(,) 'There is no "start1D_evol.txt" file!'
	321	error stop 'Initialization cannot be done for the 1D PEM.'
	322	endif
	323	therestartfi = .false.
	324	inquire(file = 'startfi_evol.nc',exist = therestartfi)
	325	if (.not. therestartfi) then
	326	write(,) 'There is no "startfi_evol.nc" file!'
	327	error stop 'Initialization cannot be done for the 1D PEM.'
	328	endif
	329
	330	call init_testphys1d('start1D_evol.txt','startfi_evol.nc',.true.,therestart1D,therestartfi,ngrid,nlayer,610., &
	331	nq,q,time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
[3068]	332	play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
[3065]	333	ps(2) = ps(1)
[3028]	334	nsplit_phys = 1
[2980]	335	#endif
[2779]	336
[3039]	337	call conf_pem(i_myear,n_myear)
[2779]	338
[2835]	339	!------------------------
[3028]	340	! I Initialization
	341	! I_b READ of start_evol.nc and starfi_evol.nc
	342	!------------------------
	343	! I_b.1 READ start_evol.nc
	344	allocate(ps_start_GCM(ngrid))
[2980]	345	#ifndef CPP_1D
[3028]	346	call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0)
[2779]	347
[3028]	348	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
[2897]	349
[3028]	350	call iniconst !new
	351	call inigeom
[2980]	352
[3028]	353	allocate(ap(nlayer + 1))
	354	allocate(bp(nlayer + 1))
	355	status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid)
	356	status = nf90_inq_varid(ncid,"ap",apvarid)
	357	status = nf90_get_var(ncid,apvarid,ap)
	358	status = nf90_inq_varid(ncid,"bp",bpvarid)
	359	status = nf90_get_var(ncid,bpvarid,bp)
	360	status = nf90_close(ncid)
[2779]	361
[3028]	362	call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, &
	363	rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
[2980]	364	#else
[3028]	365	ps_start_GCM(1) = ps(1)
[2980]	366	#endif
	367
[3096]	368	! In the PCM, these values are given to the physic by the dynamic.
[2963]	369	! Here we simply read them in the startfi_evol.nc file
[3070]	370	status = nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
[2963]	371
[3028]	372	status = nf90_inq_varid(ncid,"longitude",lonvarid)
	373	status = nf90_get_var(ncid,lonvarid,longitude)
[2963]	374
[3028]	375	status = nf90_inq_varid(ncid,"latitude",latvarid)
	376	status = nf90_get_var(ncid,latvarid,latitude)
[2963]	377
[3028]	378	status = nf90_inq_varid(ncid,"area",areavarid)
	379	status = nf90_get_var(ncid,areavarid,cell_area)
[2963]	380
[3028]	381	status = nf90_inq_varid(ncid,"soildepth",sdvarid)
	382	status = nf90_get_var(ncid,sdvarid,mlayer)
[2963]	383
[3028]	384	status = nf90_close(ncid)
[2963]	385
[3028]	386	! I_b.2 READ startfi_evol.nc
[2779]	387	! First we read the initial state (starfi.nc)
[2842]	388	#ifndef CPP_STD
[3114]	389	call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
	390	tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, &
[3130]	391	watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist)
[2779]	392
[3070]	393	! Remove unphysical values of surface tracer
	394	where (qsurf < 0.) qsurf = 0.
[2885]	395
[3028]	396	call surfini(ngrid,qsurf)
[2842]	397	#else
[3028]	398	call phys_state_var_init(nq)
	399	if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
	400	call initracer(ngrid,nq)
	401	call iniaerosol()
	402	allocate(tsurf_read_generic(ngrid))
	403	allocate(qsurf_read_generic(ngrid,nq))
	404	allocate(tsoil_read_generic(ngrid,nsoilmx))
[3114]	405	allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
[3028]	406	allocate(emis_read_generic(ngrid))
	407	allocate(tsurf(ngrid,1))
	408	allocate(qsurf(ngrid,nq,1))
	409	allocate(tsoil(ngrid,nsoilmx,1))
	410	allocate(emis(ngrid,1))
	411	allocate(watercap(ngrid,1))
	412	allocate(watercaptag(ngrid))
	413	allocate(albedo_read_generic(ngrid,2))
	414	allocate(albedo(ngrid,2,1))
	415	allocate(inertiesoil(ngrid,nsoilmx,1))
[3114]	416	call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
	417	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
	418	qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, &
	419	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
[3065]	420	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
[2842]	421
[3028]	422	nslope = 1
	423	call ini_comslope_h(ngrid,1)
[2842]	424
[3028]	425	qsurf(:,:,1) = qsurf_read_generic(:,:)
	426	tsurf(:,1) = tsurf_read_generic(:)
	427	tsoil(:,:,1) = tsoil_read_generic(:,:)
	428	emis(:,1) = emis_read_generic(:)
	429	watercap(:,1) = 0.
	430	watercaptag(:) = .false.
	431	albedo(:,1,1) = albedo_read_generic(:,1)
	432	albedo(:,2,1) = albedo_read_generic(:,2)
	433	inertiesoil(:,:,1) = inertiedat(:,:)
[2842]	434
[3028]	435	if (nslope == 1) then
	436	def_slope(1) = 0
	437	def_slope(2) = 0
	438	def_slope_mean = 0
	439	subslope_dist(:,1) = 1.
	440	endif
[2842]	441
[3070]	442	! Remove unphysical values of surface tracer
	443	qsurf(:,:,1) = qsurf_read_generic(:,:)
	444	where (qsurf < 0.) qsurf = 0.
[2842]	445	#endif
	446
[3028]	447	do nnq = 1,nqtot ! Why not using ini_tracer ?
	448	if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
	449	if (noms(nnq) == "h2o_vap") then
	450	igcm_h2o_vap = nnq
	451	mmol(igcm_h2o_vap)=18.
	452	endif
	453	if (noms(nnq) == "co2") igcm_co2 = nnq
[3065]	454	enddo
[3039]	455	r = 8.314511*1000./mugaz
[3028]	456
[2835]	457	!------------------------
[3028]	458	! I Initialization
[2835]	459	! I_c Subslope parametrisation
	460	!------------------------
[3028]	461	! Define some slope statistics
	462	iflat = 1
	463	do islope = 2,nslope
	464	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
	465	enddo
[2794]	466
[3028]	467	write(,) 'Flat slope for islope = ',iflat
	468	write(,) 'corresponding criterium = ',def_slope_mean(iflat)
[2794]	469
[3028]	470	allocate(flag_co2flow(ngrid,nslope))
	471	allocate(flag_co2flow_mesh(ngrid))
	472	allocate(flag_h2oflow(ngrid,nslope))
	473	allocate(flag_h2oflow_mesh(ngrid))
[2835]	474
[3065]	475	flag_co2flow(:,:) = 0
[3028]	476	flag_co2flow_mesh(:) = 0
[3065]	477	flag_h2oflow(:,:) = 0
[3028]	478	flag_h2oflow_mesh(:) = 0
[2835]	479
[2794]	480	!------------------------
[3028]	481	! I Initialization
[3096]	482	! I_d READ PCM data and convert to the physical grid
[3028]	483	!------------------------
[3096]	484	! 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
	485	call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
[2794]	486
[3028]	487	allocate(tsoil_ave(ngrid,nsoilmx,nslope))
	488	allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
	489	allocate(vmr_co2_gcm(ngrid,timelen))
	490	allocate(ps_timeseries(ngrid,timelen))
	491	allocate(min_co2_ice_1(ngrid,nslope))
	492	allocate(min_h2o_ice_1(ngrid,nslope))
	493	allocate(min_co2_ice_2(ngrid,nslope))
	494	allocate(min_h2o_ice_2(ngrid,nslope))
	495	allocate(tsurf_ave_yr1(ngrid,nslope))
	496	allocate(tsurf_ave(ngrid,nslope))
	497	allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
	498	allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
	499	allocate(q_co2_PEM_phys(ngrid,timelen))
	500	allocate(q_h2o_PEM_phys(ngrid,timelen))
	501	allocate(co2_ice_GCM(ngrid,nslope,timelen))
	502	allocate(watersurf_density_ave(ngrid,nslope))
	503	allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
	504	allocate(Tsurfave_before_saved(ngrid,nslope))
	505	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
	506	allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
	507	allocate(delta_co2_adsorbded(ngrid))
[3031]	508	allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
	509	allocate(delta_h2o_icetablesublim(ngrid))
[3028]	510	allocate(delta_h2o_adsorbded(ngrid))
	511	allocate(vmr_co2_pem_phys(ngrid,timelen))
[2794]	512
[3028]	513	write(,) "Downloading data Y1..."
[3096]	514	call read_data_PCM("data_PCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
[3028]	515	tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
	516	co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
	517	write(,) "Downloading data Y1 done"
[2985]	518
[3096]	519	! 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
[3028]	520	write(,) "Downloading data Y2"
[3096]	521	call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
[3028]	522	tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
	523	co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
	524	write(,) "Downloading data Y2 done"
[2794]	525
[2835]	526	!------------------------
[3028]	527	! I Initialization
	528	! I_e Initialization of the PEM variables and soil
[2835]	529	!------------------------
[3028]	530	call end_comsoil_h_PEM
	531	call ini_comsoil_h_PEM(ngrid,nslope)
	532	call end_adsorption_h_PEM
	533	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
	534	call end_ice_table_porefilling
	535	call ini_ice_table_porefilling(ngrid,nslope)
[2794]	536
[3028]	537	if (soil_pem) then
	538	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
[3070]	539	do l = 1,nsoilmx
	540	tsoil_PEM(:,l,:) = tsoil_ave(:,l,:)
	541	tsoil_phys_PEM_timeseries(:,l,:,:) = tsoil_GCM_timeseries(:,l,:,:)
	542	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,l,:,:)
[3028]	543	enddo
[3070]	544	do l = nsoilmx + 1,nsoilmx_PEM
	545	tsoil_PEM(:,l,:) = tsoil_ave(:,nsoilmx,:)
	546	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
[3028]	547	enddo
[3070]	548	watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
[3028]	549	endif !soil_pem
	550	deallocate(tsoil_ave,tsoil_GCM_timeseries)
[2794]	551
[2779]	552	!------------------------
[3028]	553	! I Initialization
[2835]	554	! I_f Compute tendencies & Save initial situation
[3028]	555	!------------------------
	556	allocate(tendencies_co2_ice(ngrid,nslope))
	557	allocate(tendencies_co2_ice_ini(ngrid,nslope))
	558	allocate(tendencies_h2o_ice(ngrid,nslope))
[2779]	559
[3028]	560	! Compute the tendencies of the evolution of ice over the years
	561	call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice)
	562	tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:)
	563	call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice)
[2895]	564
[3028]	565	deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2)
[2779]	566
[2835]	567	!------------------------
[3028]	568	! I Initialization
[2835]	569	! I_g Save initial PCM situation
[3028]	570	!------------------------
	571	allocate(initial_co2_ice_sublim(ngrid,nslope))
	572	allocate(initial_co2_ice(ngrid,nslope))
	573	allocate(initial_h2o_ice(ngrid,nslope))
[2835]	574
[2794]	575	! We save the places where water ice is sublimating
[2835]	576	! We compute the surface of water ice sublimating
[3028]	577	ini_surf_co2 = 0.
	578	ini_surf_h2o = 0.
	579	Total_surface = 0.
	580	do i = 1,ngrid
[3070]	581	Total_surface = Total_surface + cell_area(i)
[3028]	582	do islope = 1,nslope
	583	if (tendencies_co2_ice(i,islope) < 0) then
	584	initial_co2_ice_sublim(i,islope) = 1.
[3070]	585	ini_surf_co2 = ini_surf_co2 + cell_area(i)*subslope_dist(i,islope)
[3028]	586	else
	587	initial_co2_ice_sublim(i,islope) = 0.
	588	endif
	589	if (qsurf(i,igcm_co2,islope) > 0) then
	590	initial_co2_ice(i,islope) = 1.
	591	else
	592	initial_co2_ice(i,islope) = 0.
	593	endif
	594	if (tendencies_h2o_ice(i,islope) < 0) then
	595	initial_h2o_ice(i,islope) = 1.
[3130]	596	ini_surf_h2o = ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
[3028]	597	else
	598	initial_h2o_ice(i,islope) = 0.
	599	endif
[2779]	600	enddo
[3028]	601	enddo
[2779]	602
[3130]	603	write(,) "Total initial surface of co2 ice sublimating (slope) =", ini_surf_co2
	604	write(,) "Total initial surface of h2o ice sublimating (slope) =", ini_surf_h2o
[3028]	605	write(,) "Total surface of the planet=", Total_surface
	606	allocate(zplev_gcm(ngrid,nlayer + 1))
[2779]	607
[3070]	608	do ig = 1,ngrid
	609	zplev_gcm(ig,:) = ap(:) + bp(:)*ps_start_GCM(ig)
[3028]	610	enddo
[2779]	611
[3028]	612	global_ave_press_old = 0.
	613	do i = 1,ngrid
[3065]	614	global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
[3028]	615	enddo
[2779]	616
[3028]	617	global_ave_press_GCM = global_ave_press_old
	618	global_ave_press_new = global_ave_press_old
	619	write(,) "Initial global average pressure=", global_ave_press_GCM
[2779]	620
	621	!------------------------
[3028]	622	! I Initialization
[3088]	623	! I_h Read the startpem.nc
[3028]	624	!------------------------
[3088]	625	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
[3122]	626	porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries, &
	627	tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:), &
	628	qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave, &
[3028]	629	co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
[2779]	630
[3031]	631	delta_h2o_icetablesublim(:) = 0.
	632
[3130]	633	! We save the initial values for the co2 frost and perennial ice
	634	allocate(co2frost_ini(ngrid,nslope),perennial_co2ice_ini(ngrid,nslope))
	635	co2frost_ini = qsurf(:,igcm_co2,:)
	636	perennial_co2ice_ini = perennial_co2ice
	637
[3028]	638	do ig = 1,ngrid
	639	do islope = 1,nslope
	640	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
[3130]	641	qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) + perennial_co2ice(ig,islope)
[3028]	642	enddo
	643	enddo
[2779]	644
[3028]	645	if (adsorption_pem) then
	646	totmassco2_adsorbded = 0.
	647	totmassh2o_adsorbded = 0.
	648	do ig = 1,ngrid
[3070]	649	do islope = 1,nslope
[3028]	650	do l = 1,nsoilmx_PEM - 1
	651	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	652	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	653	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	654	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	655	enddo
[2961]	656	enddo
[3028]	657	enddo
[2885]	658
[3028]	659	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
	660	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
	661	endif ! adsorption
[3065]	662	deallocate(tsurf_ave_yr1)
[2794]	663
[2835]	664	!------------------------
[3028]	665	! I Initialization
[2835]	666	! I_i Compute orbit criterion
[3028]	667	!------------------------
[2842]	668	#ifndef CPP_STD
[3050]	669	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
[2842]	670	#else
[3050]	671	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
[2842]	672	#endif
[2794]	673
[3028]	674	if (evol_orbit_pem) then
[3039]	675	call orbit_param_criterion(i_myear,year_iter_max)
[3028]	676	else
	677	year_iter_max = Max_iter_pem
	678	endif
	679	!-------------------------- END INITIALIZATION -------------------------
[2794]	680
[3028]	681	!-------------------------------- RUN ----------------------------------
[2794]	682	!------------------------
	683	! II Run
[3065]	684	! II_a Update pressure, ice and tracers
[2794]	685	!------------------------
[3028]	686	year_iter = 0
[2794]	687
[3039]	688	do while (year_iter < year_iter_max .and. i_myear < n_myear)
[2835]	689	! II.a.1. Compute updated global pressure
[3028]	690	write(,) "Recomputing the new pressure..."
	691	do i = 1,ngrid
	692	do islope = 1,nslope
[3065]	693	global_ave_press_new = global_ave_press_new - gcell_area(i)tendencies_co2_ice(i,islope)subslope_dist(i,islope)/cos(pidef_slope_mean(islope)/180.)/Total_surface
[3028]	694	enddo
	695	enddo
[3065]	696
[3028]	697	if (adsorption_pem) then
	698	do i = 1,ngrid
	699	global_ave_press_new = global_ave_press_new - gcell_area(i)delta_co2_adsorbded(i)/Total_surface
[3050]	700	enddo
[3028]	701	endif
[3050]	702	write(,) 'Global average pressure old time step',global_ave_press_old
	703	write(,) 'Global average pressure new time step',global_ave_press_new
[2835]	704
	705	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
[3070]	706	allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
[3028]	707	write(,) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
	708	do l = 1,nlayer + 1
	709	do ig = 1,ngrid
	710	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
	711	enddo
	712	enddo
[2779]	713
[2835]	714	! II.a.3. Surface pressure timeseries
[3028]	715	write(,) "Recomputing the surface pressure timeserie adapted to the new pressure..."
	716	do ig = 1,ngrid
	717	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
	718	enddo
[2779]	719
[2835]	720	! II.a.4. New pressure levels timeseries
[3028]	721	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
	722	write(,) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
	723	do l = 1,nlayer + 1
	724	do ig = 1,ngrid
	725	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
	726	enddo
	727	enddo
[2779]	728
[2835]	729	! II.a.5. Tracers timeseries
[3028]	730	write(,) "Recomputing of tracer VMR timeseries for the new pressure..."
[2794]	731
[3028]	732	l = 1
	733	do ig = 1,ngrid
	734	do t = 1,timelen
	735	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))/ &
	736	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
	737	if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30
	738	if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1.
	739	enddo
	740	enddo
[2794]	741
[3028]	742	do ig = 1,ngrid
[3065]	743	do t = 1,timelen
[3028]	744	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]	745	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
	746	+ ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
	747	- (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
[3028]	748	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
	749	if (q_co2_PEM_phys(ig,t) < 0) then
	750	q_co2_PEM_phys(ig,t) = 1.e-30
	751	elseif (q_co2_PEM_phys(ig,t) > 1) then
	752	q_co2_PEM_phys(ig,t) = 1.
	753	endif
	754	mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
	755	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
	756	enddo
	757	enddo
[2794]	758
[3028]	759	deallocate(zplev_new_timeseries,zplev_old_timeseries)
	760
	761	!------------------------
[2835]	762	! II Run
	763	! II_b Evolution of the ice
[3028]	764	!------------------------
	765	write(,) "Evolution of h2o ice"
[3031]	766	call evol_h2o_ice_s(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,STOPPING_1_water)
[2794]	767
[3028]	768	write(,) "Evolution of co2 ice"
	769	call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
[2998]	770
[2794]	771	!------------------------
	772	! II Run
[3028]	773	! II_c CO2 & H2O glaciers flows
[2794]	774	!------------------------
[3028]	775	write(,) "CO2 glacier flows"
	776	if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
	777	global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
[3065]	778
[3028]	779	write(,) "H2O glacier flows"
	780	if (h2oglaciersflow) call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
[2856]	781
[2794]	782	!------------------------
	783	! II Run
[2835]	784	! II_d Update surface and soil temperatures
[2794]	785	!------------------------
[2835]	786	! II_d.1 Update Tsurf
[3028]	787	write(,) "Updating the new Tsurf"
	788	bool_sublim = .false.
	789	Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
	790	do ig = 1,ngrid
	791	do islope = 1,nslope
	792	if (initial_co2_ice(ig,islope) > 0.5 .and. qsurf(ig,igcm_co2,islope) < 1.e-10) then !co2ice disappeared, look for closest point without co2ice
	793	if (latitude_deg(ig) > 0) then
	794	do ig_loop = ig,ngrid
	795	do islope_loop = islope,iflat,-1
	796	if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
	797	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
	798	bool_sublim = .true.
	799	exit
	800	endif
	801	enddo
	802	if (bool_sublim) exit
	803	enddo
	804	else
	805	do ig_loop = ig,1,-1
	806	do islope_loop = islope,iflat
	807	if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
	808	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
	809	bool_sublim = .true.
	810	exit
	811	endif
	812	enddo
	813	if (bool_sublim) exit
	814	enddo
[2835]	815	endif
[3028]	816	initial_co2_ice(ig,islope) = 0
	817	if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
	818	albedo(ig,1,islope) = albedo_h2o_frost
	819	albedo(ig,2,islope) = albedo_h2o_frost
	820	else
	821	albedo(ig,1,islope) = albedodat(ig)
[3065]	822	albedo(ig,2,islope) = albedodat(ig)
[3028]	823	emis(ig,islope) = emissiv
	824	endif
	825	else if ((qsurf(ig,igcm_co2,islope) > 1.e-3) .and. (tendencies_co2_ice(ig,islope) > 1.e-10)) then !Put tsurf as tcond co2
	826	ave = 0.
	827	do t = 1,timelen
	828	if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
	829	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
	830	else
	831	ave = ave + tsurf_GCM_timeseries(ig,islope,t)
	832	endif
[2794]	833	enddo
[3028]	834	tsurf_ave(ig,islope) = ave/timelen
[3032]	835	! set the surface albedo to be the ice albedo
	836	if (latitude_deg(ig) > 0) then
	837	icap = 1
	838	else
	839	icap = 2
	840	endif
	841	albedo(ig,1,islope) = albedice(icap)
	842	albedo(ig,2,islope) = albedice(icap)
	843	emis(ig,islope) = emisice(icap)
[2835]	844	endif
	845	enddo
[3028]	846	enddo
[2794]	847
[3028]	848	do t = 1,timelen
[3065]	849	tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:))
[3028]	850	enddo
	851	! for the start
	852	do ig = 1,ngrid
[2835]	853	do islope = 1,nslope
[3028]	854	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
[2794]	855	enddo
[3028]	856	enddo
[2794]	857
[3028]	858	if (soil_pem) then
[2794]	859
[2835]	860	! II_d.2 Update soil temperature
[3028]	861	allocate(TI_locslope(ngrid,nsoilmx_PEM))
	862	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
	863	allocate(Tsurf_locslope(ngrid))
	864	write(,)"Updating soil temperature"
[2794]	865
[3028]	866	! Soil averaged
	867	do islope = 1,nslope
	868	TI_locslope(:,:) = TI_PEM(:,:,islope)
	869	do t = 1,timelen
	870	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
	871	Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
[3076]	872	call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
	873	call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
[3028]	874	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
	875	do ig = 1,ngrid
	876	do isoil = 1,nsoilmx_PEM
	877	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]	878	if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1)
[3028]	879	enddo
	880	enddo
	881	enddo
	882	enddo
[3070]	883	tsoil_PEM(:,:,:) = sum(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
	884	watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
[2794]	885
[3028]	886	write(,) "Update of soil temperature done"
[2888]	887
[3028]	888	deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
[2849]	889
[2835]	890	! II_d.3 Update the ice table
[3122]	891	write(,) "Compute ice table"
[3031]	892	porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:)
[3028]	893	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
[3070]	894	call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
[3031]	895
[3122]	896	! II_d.4 Update the soil thermal properties
[3028]	897	write(,) "Update soil propreties"
[3070]	898	call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
[2794]	899
[2835]	900	! II_d.5 Update the mass of the regolith adsorbded
[3028]	901	if (adsorption_pem) then
	902	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
	903	qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, &
	904	q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
[2794]	905
[3028]	906	totmassco2_adsorbded = 0.
	907	totmassh2o_adsorbded = 0.
	908	do ig = 1,ngrid
	909	do islope =1, nslope
	910	do l = 1,nsoilmx_PEM - 1
	911	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
[3070]	912	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
[3028]	913	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
[3070]	914	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
[3028]	915	enddo
	916	enddo
	917	enddo
	918	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
	919	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
	920	endif
	921	endif !soil_pem
	922
[2794]	923	!------------------------
	924	! II Run
[3088]	925	! II_e Outputs
[2794]	926	!------------------------
[3088]	927	call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_ave_press_new/))
	928	do islope = 1,nslope
	929	write(str2(1:2),'(i2.2)') islope
	930	call writediagpem(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
	931	call writediagpem(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
	932	call writediagpem(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
	933	call writediagpem(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
	934	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
	935	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
	936	enddo
	937
	938	!------------------------
	939	! II Run
	940	! II_f Update the tendencies
	941	!------------------------
[3028]	942	write(,) "Adaptation of the new co2 tendencies to the current pressure"
[3076]	943	call recomp_tend_co2_slope(ngrid,nslope,timelen,tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),emis,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
	944	global_ave_press_GCM,global_ave_press_new)
[2794]	945
[2835]	946	!------------------------
	947	! II Run
[3088]	948	! II_g Checking the stopping criterion
[2835]	949	!------------------------
[3028]	950	call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
[2779]	951
[3028]	952	call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
	953	initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
[2794]	954
[3028]	955	year_iter = year_iter + dt_pem
[3039]	956	i_myear = i_myear + dt_pem
[2794]	957
[3122]	958	write(,) "Checking all the stopping criteria..."
[3028]	959	if (STOPPING_water) then
	960	write(,) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
	961	criterion_stop = 1
	962	endif
	963	if (STOPPING_1_water) then
	964	write(,) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
	965	criterion_stop = 1
	966	endif
	967	if (STOPPING_co2) then
	968	write(,) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
	969	criterion_stop = 2
	970	endif
	971	if (STOPPING_pressure) then
	972	write(,) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
	973	criterion_stop = 3
	974	endif
	975	if (year_iter >= year_iter_max) then
	976	write(,) "STOPPING because maximum number of iterations reached"
	977	criterion_stop = 4
	978	endif
[3039]	979	if (i_myear >= n_myear) then
	980	write(,) "STOPPING because maximum number of Martian years to be simulated reached"
	981	criterion_stop = 5
	982	endif
[2794]	983
[3028]	984	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then
[2779]	985	exit
[3028]	986	else
	987	write(,) "We continue!"
[3039]	988	write(,) "Number of iterations of the PEM: year_iter =", year_iter
	989	write(,) "Number of simulated Martian years: i_myear =", i_myear
[3028]	990	endif
[2779]	991
[3065]	992	global_ave_press_old = global_ave_press_new
[2779]	993
[3028]	994	enddo ! big time iteration loop of the pem
	995	!------------------------------ END RUN --------------------------------
[2779]	996
[3028]	997	!------------------------------- OUTPUT --------------------------------
[2794]	998	!------------------------
	999	! III Output
[2835]	1000	! III_a Update surface value for the PCM start files
[2794]	1001	!------------------------
[2835]	1002	! III_a.1 Ice update (for startfi)
[2779]	1003
[2888]	1004	! H2O ice
[3028]	1005	do ig = 1,ngrid
	1006	if (watercaptag(ig)) then
	1007	watercap_sum = 0.
	1008	do islope = 1,nslope
	1009	if (qsurf(ig,igcm_h2o_ice,islope) > (watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.))) then ! water_reservoir and water cap have not changed since PCM call: here we check if we have accumulate frost or not. 1st case we have more ice than initialy
	1010	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)) ! put back ancien frost
	1011	else
	1012	! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
	1013	watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.))
[3130]	1014	qsurf(ig,igcm_h2o_ice,islope) = 0.
[3028]	1015	endif
[3130]	1016	watercap_sum = watercap_sum + watercap(ig,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
[3028]	1017	watercap(ig,islope) = 0.
	1018	enddo
	1019	water_reservoir(ig) = water_reservoir(ig) + watercap_sum
	1020	endif
	1021	enddo
[2888]	1022
[3028]	1023	do ig = 1,ngrid
	1024	water_sum = 0.
	1025	do islope = 1,nslope
	1026	water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
	1027	enddo
[3065]	1028	if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
[3130]	1029	if (water_sum > threshold_water_frost2perennial) then ! the overall mesh can be considered as an infite source of water. No need to change the albedo: done in II.d.1
[3028]	1030	watercaptag(ig) = .true.
[3130]	1031	water_reservoir(ig) = water_reservoir(ig) + threshold_water_frost2perennial/2. ! half of the excess ices goes to the reservoir, we let the rest to be frost
[3028]	1032	do islope = 1,nslope
[3130]	1033	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perennial/2.cos(pidef_slope_mean(islope)/180.)
[3028]	1034	enddo
	1035	endif
	1036	else ! let's check that the infinite source of water disapear
[3130]	1037	if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perennial) then
[3028]	1038	watercaptag(ig) = .false.
	1039	do islope = 1,nslope
	1040	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
	1041	enddo
	1042	water_reservoir(ig) = 0.
	1043	endif
	1044	endif
	1045	enddo
[2888]	1046
[2998]	1047	! CO2 ice
[3028]	1048	do ig = 1,ngrid
	1049	do islope = 1,nslope
[3130]	1050	if (qsurf(ig,igcm_co2,islope) == threshold_co2_frost2perennial) then
	1051	! If co2 ice is equal to the threshold, then everything is transformed into perennial ice
	1052	perennial_co2ice(ig,islope) = threshold_co2_frost2perennial
	1053	qsurf(ig,igcm_co2,islope) = 0.
	1054	else if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perennial) then
	1055	! If co2 ice is superior to the threshold, then co2 frost is equal to the threshold (max possible value)
	1056	! and the leftover is transformed into perennial ice
	1057	perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope) - threshold_co2_frost2perennial
	1058	qsurf(ig,igcm_co2,islope) = threshold_co2_frost2perennial
	1059	else
	1060	! If co2 ice is inferior to the threshold, then we compare with the initial state
	1061	if (qsurf(ig,igcm_co2,islope) > perennial_co2ice_ini(ig,islope)) then
	1062	! If co2 ice higher than the initial perennial ice, then the change is affected only to the frost
	1063	perennial_co2ice(ig,islope) = perennial_co2ice_ini(ig,islope)
	1064	qsurf(ig,igcm_co2,islope) = qsurf(ig,igcm_co2,islope) - perennial_co2ice_ini(ig,islope)
	1065	else
	1066	! If co2 ice is lower than the initial perennial ice, then there is no frost anymore
	1067	perennial_co2ice(ig,islope) = qsurf(ig,igcm_co2,islope)
	1068	qsurf(ig,igcm_co2,islope) = 0.
	1069	endif
[3028]	1070	endif
	1071	enddo
	1072	enddo
[2998]	1073
[2849]	1074	! III_a.2 Tsoil update (for startfi)
[3028]	1075	if (soil_pem) then
	1076	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
[3065]	1077	tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
[3028]	1078	endif
[2779]	1079
[2835]	1080	! III_a.4 Pressure (for start)
[3065]	1081	ps(:) = ps(:)*global_ave_press_new/global_ave_press_GCM
	1082	ps_start_GCM(:) = ps_start_GCM(:)*global_ave_press_new/global_ave_press_GCM
[2794]	1083
[2835]	1084	! III_a.5 Tracer (for start)
[3028]	1085	allocate(zplev_new(ngrid,nlayer + 1))
[2835]	1086
[3028]	1087	do l = 1,nlayer + 1
[3065]	1088	zplev_new(:,l) = ap(l) + bp(l)*ps_start_GCM(:)
[3028]	1089	enddo
[2835]	1090
[3028]	1091	do nnq = 1,nqtot
	1092	if (noms(nnq) /= "co2") then
	1093	do l = 1,llm - 1
	1094	do ig = 1,ngrid
	1095	q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
	1096	enddo
	1097	q(:,llm,nnq) = q(:,llm - 1,nnq)
	1098	enddo
	1099	else
	1100	do l = 1,llm - 1
	1101	do ig = 1,ngrid
	1102	q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
[3070]	1103	+ ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
[3028]	1104	enddo
	1105	q(:,llm,nnq) = q(:,llm - 1,nnq)
	1106	enddo
	1107	endif
	1108	enddo
[2835]	1109
[3096]	1110	! Conserving the tracers mass for PCM start files
[3028]	1111	do nnq = 1,nqtot
	1112	do ig = 1,ngrid
	1113	do l = 1,llm - 1
	1114	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]	1115	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
	1116	q(ig,l,nnq) = 1.
	1117	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
[3028]	1118	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
[2835]	1119	endif
[3028]	1120	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
	1121	enddo
	1122	enddo
	1123	enddo
[2779]	1124
[3039]	1125	if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
[2779]	1126
	1127	!------------------------
[3028]	1128	! III Output
	1129	! III_b Write restart_evol.nc and restartfi_evol.nc
	1130	!------------------------
	1131	! III_b.1 Write restart_evol.nc
[3042]	1132	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
[3028]	1133	pday = day_ini
[3042]	1134	ztime_fin = time_phys
[2779]	1135
[3028]	1136	allocate(p(ip1jmp1,nlayer + 1))
[2980]	1137	#ifndef CPP_1D
[3028]	1138	call pression (ip1jmp1,ap,bp,ps,p)
	1139	call massdair(p,masse)
[3039]	1140	call dynredem0("restart_evol.nc",day_ini,phis)
	1141	call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
[3028]	1142	write(,) "restart_evol.nc has been written"
[2980]	1143	#else
[3069]	1144	call writerestart1D('restart1D_evol.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
[3065]	1145	write(,) "restart1D_evol.txt has been written"
[2980]	1146	#endif
	1147
[3028]	1148	! III_b.2 Write restartfi_evol.nc
[2842]	1149	#ifndef CPP_STD
[3028]	1150	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
	1151	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
	1152	inertiedat,def_slope,subslope_dist)
[3114]	1153	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
	1154	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
	1155	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
[3130]	1156	wstar,watercap,perennial_co2ice)
[2842]	1157	#else
[3028]	1158	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
	1159	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
	1160	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[3114]	1161	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
	1162	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
	1163	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
	1164	tsea_ice,sea_ice)
[2842]	1165	#endif
[3028]	1166	write(,) "restartfi_evol.nc has been written"
[2842]	1167
[2794]	1168	!------------------------
	1169	! III Output
[3088]	1170	! III_c Write restartpem.nc
[2794]	1171	!------------------------
[3088]	1172	call pemdem0("restartpem.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
[3028]	1173	float(day_ini),0.,nslope,def_slope,subslope_dist)
[3088]	1174	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
[3122]	1175	TI_PEM, porefillingice_depth,porefillingice_thickness, &
[3028]	1176	co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
[3088]	1177	write(,) "restartpem.nc has been written"
[3096]	1178	call info_PEM(year_iter,criterion_stop,i_myear,n_myear)
[2779]	1179
[3039]	1180	write(,) "The PEM has run for", year_iter, "Martian years."
	1181	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
	1182	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
	1183	write(,) "LL & RV & JBC: so far, so good!"
[2794]	1184
[3028]	1185	deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
	1186	deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
	1187	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
[3031]	1188	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
[3130]	1189	deallocate(co2frost_ini,perennial_co2ice_ini)
[3028]	1190	!----------------------------- END OUTPUT ------------------------------
[2897]	1191
[2779]	1192	END PROGRAM pem

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

Download in other formats:

Original Format