Context Navigation

pem.F90 @ 3097

Last change on this file since 3097 was 3096, checked in by jbclement, 2 years ago

PEM:
The management of files during the chained simulation of PCM/PEM runs has been simplified:

"tmp_PEMyears.txt" and "info_run_PEM.txt" have been merged into one file called "info_PEM.txt";
"reshape_XIOS_output.F90" now creates directly the "data_PCM_Y*.nc" files needed by the PEM;
where it is relevant, 'GCM' has been replaced by 'PCM' in the files naming;
the files in deftank have been updated consequently.

Following r3095, 'iniorbit' is now a subroutine of "planete_h.F90".
JBC

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

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

Download in other formats:

Original Format