Context Navigation

pem.F90 @ 3047

Last change on this file since 3047 was 3042, checked in by jbclement, 2 years ago
Mars PEM: The date and time are redefined at the end of PEM to write the file "restarfi_evol.nc". Like so, the next GCM runs can restart at the date and time given at the beginning of the simulation since the PEM is only progressing by step of one Martian year. JBC
File size: 56.8 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
[3028]	6	! I_d READ GCM data and convert to the physical grid
	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
	10	! I_h Read the PEMstart
	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
[3028]	18	! II_e Update the tendencies
[2835]	19	! II_f Checking the stopping criterion
[2779]	20
	21	! III Output
[2835]	22	! III_a Update surface value for the PCM start files
[3028]	23	! III_b Write restart_evol.nc and restartfi_evol.nc
	24	! III_c Write restartfi_PEM.nc
[2779]	25	!------------------------
	26
	27	PROGRAM pem
	28
[3028]	29	use phyetat0_mod, only: phyetat0
	30	use phyredem, only: physdem0, physdem1
	31	use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close
	32	use turb_mod, only: q2, wstar
	33	use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, major_slope, ini_comslope_h
	34	use logic_mod, only: iflag_phys
	35	use mod_const_mpi, only: COMM_LMDZ
	36	use comconst_mod, only: rad, g, cpp, pi
	37	use infotrac
	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
	50	fluxgeo, & ! Geothermal flux for the PEM and GCM
	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
[3039]	55	use time_evol_mod, only: dt_pem, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
[3028]	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, &
[3031]	59	ini_ice_table_porefilling, computeice_table_equilibrium,compute_massh2o_exchange_ssi
[3028]	60	use soil_thermalproperties_mod, only: update_soil_thermalproperties
[3039]	61	use time_phylmdz_mod, only: daysec, dtphys, day_end
[2985]	62
[2842]	63	#ifndef CPP_STD
[3028]	64	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil
	65	use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
	66	albedodat, zmea, zstd, zsig, zgam, zthe, &
	67	hmons, summit, base,albedo_h2o_frost, &
[3032]	68	frost_albedo_threshold, emissiv, watercaptag, perenial_co2ice, &
	69	emisice, albedice
[3028]	70	use dimradmars_mod, only: totcloudfrac, albedo
	71	use dust_param_mod, only: tauscaling
	72	use tracer_mod, only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
	73	use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
	74	use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit
	75	use comcstfi_h, only: r, mugaz
[3032]	76	use paleoclimate_mod, only: albedo_perenialco2
[2842]	77	#else
[3028]	78	use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names
	79	use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, &
	80	PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, &
	81	ALBEDO_CO2_ICE_SPECTV, phys_state_var_init
[3039]	82	use aerosol_mod, only: iniaerosol
	83	use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit
	84	use comcstfi_mod, only: r, mugaz
[2842]	85	#endif
[2985]	86
[3028]	87	#ifndef CPP_1D
[3039]	88	use iniphysiq_mod, only: iniphysiq
	89	use control_mod, only: iphysiq, day_step, nsplit_phys
[3019]	90	#else
[3039]	91	use time_phylmdz_mod, only: iphysiq, day_step
[3019]	92	#endif
[3039]	93
[2980]	94	#ifdef CPP_1D
[3028]	95	use regular_lonlat_mod, only: init_regular_lonlat
	96	use physics_distribution_mod, only: init_physics_distribution
	97	use mod_grid_phy_lmdz, only: regular_lonlat
	98	use init_phys_1d_mod, only: init_phys_1d
[2980]	99	#endif
[2835]	100
[3028]	101	IMPLICIT NONE
[2980]	102
[3028]	103	include "dimensions.h"
	104	include "paramet.h"
	105	include "comgeom.h"
	106	include "iniprint.h"
[3039]	107	include "callkeys.h"
[2779]	108
[3028]	109	integer ngridmx
	110	parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
[2794]	111
[3028]	112	! Same variable names as in the GCM
	113	integer :: ngrid ! Number of physical grid points
	114	integer :: nlayer ! Number of vertical layer
	115	integer :: nq ! Number of tracer
	116	integer :: day_ini ! First day of the simulation
	117	real :: pday ! Physical day
	118	real :: time_phys ! Same as GCM
	119	real :: ptimestep ! Same as GCM
	120	real :: ztime_fin ! Same as GCM
[2794]	121
[3028]	122	! Variables to read start.nc
	123	character(len = *), parameter :: FILE_NAME_start = "start_evol.nc" ! Name of the file used for initialsing the PEM
[2779]	124
[3028]	125	! Dynamic variables
	126	real :: vcov(ip1jm,llm), ucov(ip1jmp1,llm) ! vents covariants
	127	real :: teta(ip1jmp1,llm) ! temperature potentielle
	128	real, dimension(:,:,:), allocatable :: q ! champs advectes
	129	real :: ps(ip1jmp1) ! pression au sol
	130	real, dimension(:), allocatable :: ps_start_GCM !(ngrid) pression au sol
	131	real, dimension(:,:), allocatable :: ps_timeseries !(ngrid x timelen) ! pression au sol instantannées
	132	real :: masse(ip1jmp1,llm) ! masse d'air
	133	real :: phis(ip1jmp1) ! geopotentiel au sol
	134	real :: time_0
[2779]	135
[3028]	136	! Variables to read starfi.nc
	137	character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
	138	character*2 str2
	139	integer :: ncid, varid,status ! Variable for handling opening of files
	140	integer :: phydimid, subdimid, nlayerdimid, nqdimid ! Variable ID for Netcdf files
	141	integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files
	142	integer :: apvarid, bpvarid ! Variable ID for Netcdf files
[2794]	143
[3028]	144	! Variables to read starfi.nc and write restartfi.nc
	145	real, dimension(:), allocatable :: longitude ! Longitude read in FILE_NAME and written in restartfi
	146	real, dimension(:), allocatable :: latitude ! Latitude read in FILE_NAME and written in restartfi
	147	real, dimension(:), allocatable :: ap ! Coefficient ap read in FILE_NAME_start and written in restart
	148	real, dimension(:), allocatable :: bp ! Coefficient bp read in FILE_NAME_start and written in restart
	149	real, dimension(:), allocatable :: cell_area ! Cell_area read in FILE_NAME and written in restartfi
	150	real :: Total_surface ! Total surface of the planet
[2897]	151
[3028]	152	! Variables for h2o_ice evolution
	153	real :: ini_surf_h2o ! Initial surface of sublimating h2o ice
	154	real :: ini_surf_co2 ! Initial surface of sublimating co2 ice
	155	real :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa]
	156	real :: global_ave_press_old ! constant: Global average pressure of initial/previous time step
	157	real :: global_ave_press_new ! constant: Global average pressure of current time step
[3039]	158	real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2]
	159	real, dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2]
	160	real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
[3028]	161	real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
	162	logical :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
	163	logical :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
	164	logical :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
	165	logical :: STOPPING_pressure ! Logical : is the criterion (% of change in the surface pressure) reached?
	166	integer :: criterion_stop ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
	167	real, save :: A , B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg]
	168	real, allocatable :: vmr_co2_gcm(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
	169	real, allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM
	170	real, allocatable :: q_co2_PEM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]
	171	real, allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
	172	integer :: timelen ! # time samples
	173	real :: ave ! intermediate varibale to compute average
	174	real, allocatable :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
	175	real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1
[2779]	176
[3028]	177	! Variables for slopes
	178	real, dimension(:,:), allocatable :: min_co2_ice_1 ! ngrid field : minimum of co2 ice at each point for the first year [kg/m^2]
	179	real, dimension(:,:), allocatable :: min_co2_ice_2 ! ngrid field : minimum of co2 ice at each point for the second year [kg/m^2]
	180	real, dimension(:,:), allocatable :: min_h2o_ice_1 ! ngrid field : minimum of water ice at each point for the first year [kg/m^2]
	181	real, dimension(:,:), allocatable :: min_h2o_ice_2 ! ngrid field : minimum of water ice at each point for the second year [kg/m^2]
	182	real, dimension(:,:,:), allocatable :: co2_ice_GCM ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2]
	183	real, dimension(:,:), allocatable :: initial_co2_ice_sublim ! physical point field : Logical array indicating sublimating point of co2 ice
	184	real, dimension(:,:), allocatable :: initial_h2o_ice ! physical point field : Logical array indicating if there is water ice at initial state
	185	real, dimension(:,:), allocatable :: initial_co2_ice ! physical point field : Logical array indicating if there is co2 ice at initial state
	186	real, dimension(:,:), allocatable :: tendencies_co2_ice ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
	187	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 GCM
	188	real, dimension(:,:), allocatable :: tendencies_h2o_ice ! physical pointx slope field : Tendency of evolution of perenial h2o ice
[3039]	189	real, dimension(:,:), allocatable :: flag_co2flow(:,:) ! (ngrid,nslope): Flag where there is a CO2 glacier flow
	190	real, dimension(:), allocatable :: flag_co2flow_mesh(:) ! (ngrid) : Flag where there is a CO2 glacier flow
	191	real, dimension(:,:), allocatable :: flag_h2oflow(:,:) ! (ngrid,nslope): Flag where there is a H2O glacier flow
	192	real, dimension(:), allocatable :: flag_h2oflow_mesh(:) ! (ngrid) : Flag where there is a H2O glacier flow
[2779]	193
[3028]	194	! Variables for surface and soil
	195	real, allocatable :: tsurf_ave(:,:) ! Physic x SLOPE field : Averaged Surface Temperature [K]
	196	real, allocatable :: tsoil_ave(:,:,:) ! Physic x SOIL x SLOPE field : Averaged Soil Temperature [K]
	197	real, allocatable :: tsurf_GCM_timeseries(:,:,:) ! ngrid x SLOPE XTULES field : Surface Temperature in timeseries [K]
[3039]	198	real, allocatable :: tsoil_phys_PEM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
	199	real, allocatable :: tsoil_GCM_timeseries(:,:,:,:) ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
[3028]	200	real, allocatable :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
	201	real, allocatable :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI]
	202	real, allocatable :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
	203	real, allocatable :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
	204	real, allocatable :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
	205	real, allocatable :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
	206	real, allocatable :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
	207	real, allocatable :: watersoil_density_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
	208	real, allocatable :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K]
	209	real, allocatable :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
	210	real, allocatable :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
	211	real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
	212	real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
	213	logical :: bool_sublim ! logical to check if there is sublimation or not
[3031]	214	real,allocatable :: porefillingice_thickness_prev_iter(:,:) ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
	215	real,allocatable :: delta_h2o_icetablesublim(:) ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
[3028]	216	! Some variables for the PEM run
	217	real, parameter :: year_step = 1 ! timestep for the pem
	218	integer :: year_iter ! number of iteration
	219	integer :: year_iter_max ! maximum number of iterations before stopping
[3039]	220	integer :: i_myear ! Global number of Martian years of the chained simulations
	221	integer :: n_myear ! Maximum number of Martian years of the chained simulations
[3028]	222	real :: timestep ! timestep [s]
	223	real :: watercap_sum ! total mass of water cap [kg/m^2]
	224	real :: water_sum ! total mass of water in the mesh [kg/m^2]
[2779]	225
[2842]	226	#ifdef CPP_STD
[3028]	227	real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice
	228	real :: albedo_h2o_frost ! albedo of h2o frost
	229	real, allocatable :: tsurf_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	230	real, allocatable :: qsurf_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	231	real, allocatable :: tsoil_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	232	real, allocatable :: emis_read_generic(:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	233	real, allocatable :: albedo_read_generic(:,:) ! Temporary variable to do the subslope transfert dimensiion when reading form generic
	234	real, allocatable :: tsurf(:,:) ! Subslope variable, only needed in the GENERIC case
	235	real, allocatable :: qsurf(:,:,:) ! Subslope variable, only needed in the GENERIC case
	236	real, allocatable :: tsoil(:,:,:) ! Subslope variable, only needed in the GENERIC case
	237	real, allocatable :: emis(:,:) ! Subslope variable, only needed in the GENERIC case
	238	real, allocatable :: watercap(:,:) ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model
	239	logical, allocatable :: WATERCAPTAG(:) ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model
	240	real, allocatable :: albedo(:,:,:) ! Subslope variable, only needed in the GENERIC case
	241	real, allocatable :: inertiesoil(:,:,:) ! Subslope variable, only needed in the GENERIC case
[2842]	242	#endif
	243
[2980]	244	#ifdef CPP_1D
[3028]	245	integer :: nsplit_phys
	246	integer, parameter :: jjm_value = jjm - 1
	247	integer :: ierr
[2980]	248	#else
[3028]	249	integer, parameter :: jjm_value = jjm
[2980]	250	#endif
	251
[3028]	252	! Loop variables
[3039]	253	integer :: i, j, ig0, l, ig, nnq, t, islope, ig_loop, islope_loop, iloop, isoil, icap
[2779]	254
[3028]	255	! Parallel variables
[2842]	256	#ifndef CPP_STD
[3028]	257	is_sequential = .true.
	258	is_parallel = .false.
	259	is_mpi_root = .true.
	260	is_omp_root = .true.
	261	is_master = .true.
[2842]	262	#endif
[2779]	263
[3028]	264	day_ini = 0 ! test
	265	time_phys = 0. ! test
[2779]	266
[2835]	267	! Some constants
[3028]	268	ngrid = ngridmx
	269	nlayer = llm
[2794]	270
[3028]	271	A = (1/m_co2 - 1/m_noco2)
	272	B = 1/m_noco2
[2794]	273
[3028]	274	year_day = 669
	275	daysec = 88775.
	276	timestep = year_day*daysec/year_step
[2794]	277
[3028]	278	!----------------------------- INITIALIZATION --------------------------
[2779]	279	!------------------------
[3028]	280	! I Initialization
	281	! I_a READ run.def
[2779]	282	!------------------------
[2980]	283	#ifndef CPP_1D
[3028]	284	dtphys = 0
	285	call conf_gcm(99,.true.)
	286	call infotrac_init
	287	nq = nqtot
	288	allocate(q(ip1jmp1,llm,nqtot))
[2980]	289	#else
[3028]	290	! load tracer names from file 'traceur.def'
	291	open(90,file = 'traceur.def',status = 'old',form = 'formatted',iostat = ierr)
	292	if (ierr /= 0) then
	293	write(,) 'Cannot find required file "traceur.def"'
	294	write(,) ' If you want to run with tracers, I need it'
	295	write(,) ' ... might as well stop here ...'
	296	stop
	297	else
	298	write(,) "pem1d: Reading file traceur.def"
	299	! read number of tracers:
	300	read(90,*,iostat = ierr) nq
	301	write(,) "nq",nq
	302	nqtot = nq ! set value of nqtot (in infotrac module) as nq
	303	if (ierr /= 0) then
[2980]	304	write(,) "pem1d: error reading number of tracers"
	305	write(,) " (first line of traceur.def) "
	306	stop
[3028]	307	endif
	308	if (nq < 1) then
[2980]	309	write(,) "pem1d: error number of tracers"
	310	write(,) "is nq=",nq," but must be >=1!"
	311	stop
	312	endif
[3028]	313	endif
	314	nq = nqtot
	315	allocate(q(ip1jmp1,llm,nqtot))
	316	allocate(ap(nlayer + 1))
	317	allocate(bp(nlayer + 1))
	318	call init_phys_1d(llm,nqtot,vcov,ucov,teta,q,ps,time_0,ap,bp)
	319	pi = 2.*asin(1.)
	320	g = 3.72
	321	nsplit_phys = 1
[2980]	322	#endif
[2779]	323
[3039]	324	call conf_pem(i_myear,n_myear)
[2779]	325
[2835]	326	!------------------------
[3028]	327	! I Initialization
	328	! I_b READ of start_evol.nc and starfi_evol.nc
	329	!------------------------
	330	! I_b.1 READ start_evol.nc
	331	allocate(ps_start_GCM(ngrid))
[2980]	332	#ifndef CPP_1D
[3028]	333	call dynetat0(FILE_NAME_start,vcov,ucov,teta,q,masse,ps,phis,time_0)
[2779]	334
[3028]	335	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
[2897]	336
[3028]	337	call iniconst !new
	338	call inigeom
[2980]	339
[3028]	340	allocate(ap(nlayer + 1))
	341	allocate(bp(nlayer + 1))
	342	status = nf90_open(FILE_NAME_start,NF90_NOWRITE,ncid)
	343	status = nf90_inq_varid(ncid,"ap",apvarid)
	344	status = nf90_get_var(ncid,apvarid,ap)
	345	status = nf90_inq_varid(ncid,"bp",bpvarid)
	346	status = nf90_get_var(ncid,bpvarid,bp)
	347	status = nf90_close(ncid)
[2779]	348
[3028]	349	call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys, &
	350	rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys)
[2980]	351	#else
[3028]	352	ps_start_GCM(1) = ps(1)
[2980]	353	#endif
	354
[2963]	355	! In the gcm, these values are given to the physic by the dynamic.
	356	! Here we simply read them in the startfi_evol.nc file
[3028]	357	status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
[2963]	358
[3028]	359	allocate(longitude(ngrid))
	360	allocate(latitude(ngrid))
	361	allocate(cell_area(ngrid))
[2963]	362
[3028]	363	status = nf90_inq_varid(ncid,"longitude",lonvarid)
	364	status = nf90_get_var(ncid,lonvarid,longitude)
[2963]	365
[3028]	366	status = nf90_inq_varid(ncid,"latitude",latvarid)
	367	status = nf90_get_var(ncid,latvarid,latitude)
[2963]	368
[3028]	369	status = nf90_inq_varid(ncid,"area",areavarid)
	370	status = nf90_get_var(ncid,areavarid,cell_area)
[2963]	371
[3028]	372	status = nf90_inq_varid(ncid,"soildepth",sdvarid)
	373	status = nf90_get_var(ncid,sdvarid,mlayer)
[2963]	374
[3028]	375	status = nf90_close(ncid)
[2963]	376
[3028]	377	! I_b.2 READ startfi_evol.nc
[2779]	378	! First we read the initial state (starfi.nc)
[2842]	379	#ifndef CPP_STD
[3028]	380	call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,day_ini,time_phys,tsurf, &
	381	tsoil,albedo,emis,q2,qsurf,tauscaling,totcloudfrac,wstar, &
	382	watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
[2779]	383
[2835]	384	! Remove unphysical values of surface tracer
[3028]	385	do i = 1,ngrid
	386	do nnq = 1,nqtot
	387	do islope = 1,nslope
	388	if (qsurf(i,nnq,islope) < 0) qsurf(i,nnq,islope) = 0.
	389	enddo
	390	enddo
	391	enddo
[2885]	392
[3028]	393	call surfini(ngrid,qsurf)
[2842]	394	#else
[3028]	395	call phys_state_var_init(nq)
	396	if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
	397	call initracer(ngrid,nq)
	398	call iniaerosol()
	399	allocate(tsurf_read_generic(ngrid))
	400	allocate(qsurf_read_generic(ngrid,nq))
	401	allocate(tsoil_read_generic(ngrid,nsoilmx))
	402	allocate(emis_read_generic(ngrid))
	403	allocate(tsurf(ngrid,1))
	404	allocate(qsurf(ngrid,nq,1))
	405	allocate(tsoil(ngrid,nsoilmx,1))
	406	allocate(emis(ngrid,1))
	407	allocate(watercap(ngrid,1))
	408	allocate(watercaptag(ngrid))
	409	allocate(albedo_read_generic(ngrid,2))
	410	allocate(albedo(ngrid,2,1))
	411	allocate(inertiesoil(ngrid,nsoilmx,1))
	412	call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,day_ini,time_phys, &
	413	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
	414	qsurf_read_generic,cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, &
	415	tslab, tsea_ice,sea_ice)
	416	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
[2842]	417
[3028]	418	nslope = 1
	419	call ini_comslope_h(ngrid,1)
[2842]	420
[3028]	421	qsurf(:,:,1) = qsurf_read_generic(:,:)
	422	tsurf(:,1) = tsurf_read_generic(:)
	423	tsoil(:,:,1) = tsoil_read_generic(:,:)
	424	emis(:,1) = emis_read_generic(:)
	425	watercap(:,1) = 0.
	426	watercaptag(:) = .false.
	427	albedo(:,1,1) = albedo_read_generic(:,1)
	428	albedo(:,2,1) = albedo_read_generic(:,2)
	429	inertiesoil(:,:,1) = inertiedat(:,:)
[2842]	430
[3028]	431	if (nslope == 1) then
	432	def_slope(1) = 0
	433	def_slope(2) = 0
	434	def_slope_mean = 0
	435	subslope_dist(:,1) = 1.
	436	endif
[2842]	437
	438	! Remove unphysical values of surface tracer
[3028]	439	do i = 1,ngrid
	440	do nnq = 1,nqtot
	441	qsurf(i,nnq,1) = qsurf_read_generic(i,nnq)
	442	if (qsurf(i,nnq,1) < 0) qsurf(i,nnq,1) = 0.
	443	enddo
	444	enddo
[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
	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
[3028]	475	flag_co2flow(:,:) = 0
	476	flag_co2flow_mesh(:) = 0
	477	flag_h2oflow(:,:) = 0
	478	flag_h2oflow_mesh(:) = 0
[2835]	479
[2794]	480	!------------------------
[3028]	481	! I Initialization
	482	! I_d READ GCM data and convert to the physical grid
	483	!------------------------
[2794]	484	! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value
[3028]	485	call nb_time_step_GCM("data_GCM_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..."
	514	call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
	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
[2794]	519	! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
[3028]	520	write(,) "Downloading data Y2"
	521	call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
	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)
	539	do l=1,nsoilmx
[2897]	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
	544	do l=nsoilmx+1,nsoilmx_PEM
[2897]	545	tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:)
	546	watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:)
[3028]	547	enddo
	548	watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
	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
	581	Total_surface = Total_surface+cell_area(i)
	582	do islope = 1,nslope
	583	if (tendencies_co2_ice(i,islope) < 0) then
	584	initial_co2_ice_sublim(i,islope) = 1.
	585	ini_surf_co2 = ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
	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.
	596	ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
	597	else
	598	initial_h2o_ice(i,islope) = 0.
	599	endif
[2779]	600	enddo
[3028]	601	enddo
[2779]	602
[3028]	603	write(,) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
	604	write(,) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
	605	write(,) "Total surface of the planet=", Total_surface
	606	allocate(zplev_gcm(ngrid,nlayer + 1))
[2779]	607
[3028]	608	do l = 1,nlayer + 1
	609	do ig = 1,ngrid
	610	zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
	611	enddo
	612	enddo
[2779]	613
[3028]	614	global_ave_press_old = 0.
	615	do i = 1,ngrid
	616	global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
	617	enddo
[2779]	618
[3028]	619	global_ave_press_GCM = global_ave_press_old
	620	global_ave_press_new = global_ave_press_old
	621	write(,) "Initial global average pressure=", global_ave_press_GCM
[2779]	622
	623	!------------------------
[3028]	624	! I Initialization
[2835]	625	! I_h Read the PEMstart
[3028]	626	!------------------------
	627	call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
	628	porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries, &
	629	tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:), &
	630	qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave, &
	631	co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
[2779]	632
[3031]	633	delta_h2o_icetablesublim(:) = 0.
	634
[3028]	635	do ig = 1,ngrid
	636	do islope = 1,nslope
	637	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
	638	enddo
	639	enddo
[2779]	640
[3028]	641	if (adsorption_pem) then
	642	totmassco2_adsorbded = 0.
	643	totmassh2o_adsorbded = 0.
	644	do ig = 1,ngrid
	645	do islope = 1, nslope
	646	do l = 1,nsoilmx_PEM - 1
	647	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	648	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	649	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	650	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	651	enddo
[2961]	652	enddo
[3028]	653	enddo
[2885]	654
[3028]	655	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
	656	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
	657	endif ! adsorption
	658	deallocate(tsurf_ave_yr1)
[2794]	659
[2835]	660	!------------------------
[3028]	661	! I Initialization
[2835]	662	! I_i Compute orbit criterion
[3028]	663	!------------------------
[2842]	664	#ifndef CPP_STD
[3028]	665	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
[2842]	666	#else
[3039]	667	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
[2842]	668	#endif
[2794]	669
[3028]	670	if (evol_orbit_pem) then
[3039]	671	call orbit_param_criterion(i_myear,year_iter_max)
[3028]	672	else
	673	year_iter_max = Max_iter_pem
	674	endif
	675	!-------------------------- END INITIALIZATION -------------------------
[2794]	676
[3028]	677	!-------------------------------- RUN ----------------------------------
[2794]	678	!------------------------
	679	! II Run
[3028]	680	! II_a Update pressure, ice and tracers
[2794]	681	!------------------------
[3028]	682	year_iter = 0
[2794]	683
[3039]	684	do while (year_iter < year_iter_max .and. i_myear < n_myear)
[2835]	685	! II.a.1. Compute updated global pressure
[3028]	686	write(,) "Recomputing the new pressure..."
[2939]	687
[3028]	688	do i = 1,ngrid
	689	do islope = 1,nslope
	690	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
	691	enddo
	692	enddo
	693	write(,) 'Global average pressure old time step',global_ave_press_old
	694	call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
[2863]	695
[3028]	696	if (adsorption_pem) then
	697	do i = 1,ngrid
	698	global_ave_press_new = global_ave_press_new - gcell_area(i)delta_co2_adsorbded(i)/Total_surface
	699	enddo
	700	write(,) 'Global average pressure old time step',global_ave_press_old
	701	write(,) 'Global average pressure new time step',global_ave_press_new
	702	endif
[2835]	703
	704	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
[3028]	705	allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
	706	write(,) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
	707	do l = 1,nlayer + 1
	708	do ig = 1,ngrid
	709	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
	710	enddo
	711	enddo
[2779]	712
[2835]	713	! II.a.3. Surface pressure timeseries
[3028]	714	write(,) "Recomputing the surface pressure timeserie adapted to the new pressure..."
	715	do ig = 1,ngrid
	716	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
	717	enddo
[2779]	718
[2835]	719	! II.a.4. New pressure levels timeseries
[3028]	720	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
	721	write(,) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
	722	do l = 1,nlayer + 1
	723	do ig = 1,ngrid
	724	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
	725	enddo
	726	enddo
[2779]	727
[2835]	728	! II.a.5. Tracers timeseries
[3028]	729	write(,) "Recomputing of tracer VMR timeseries for the new pressure..."
[2794]	730
[3028]	731	l = 1
	732	do ig = 1,ngrid
	733	do t = 1,timelen
	734	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))/ &
	735	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
	736	if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30
	737	if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1.
	738	enddo
	739	enddo
[2794]	740
[3028]	741	do ig = 1,ngrid
	742	do t = 1, timelen
	743	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))/ &
	744	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
	745	+ ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
	746	- (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
	747	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
	748	if (q_co2_PEM_phys(ig,t) < 0) then
	749	q_co2_PEM_phys(ig,t) = 1.e-30
	750	elseif (q_co2_PEM_phys(ig,t) > 1) then
	751	q_co2_PEM_phys(ig,t) = 1.
	752	endif
	753	mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
	754	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
	755	enddo
	756	enddo
[2794]	757
[3028]	758	deallocate(zplev_new_timeseries,zplev_old_timeseries)
	759
	760	!------------------------
[2835]	761	! II Run
	762	! II_b Evolution of the ice
[3028]	763	!------------------------
	764	write(,) "Evolution of h2o ice"
[3031]	765	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]	766
[3028]	767	write(,) "Evolution of co2 ice"
	768	call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
[2998]	769
[3028]	770	do islope=1, nslope
	771	write(str2(1:2),'(i2.2)') islope
	772	call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
	773	call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
	774	call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
	775	call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
	776	enddo
[2857]	777
[2794]	778	!------------------------
	779	! II Run
[3028]	780	! II_c CO2 & H2O glaciers flows
[2794]	781	!------------------------
[3028]	782	write(,) "CO2 glacier flows"
[2794]	783
[3028]	784	if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
	785	global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
	786
	787	write(,) "H2O glacier flows"
[2779]	788
[3028]	789	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]	790
[3028]	791	do islope = 1,nslope
	792	write(str2(1:2),'(i2.2)') islope
	793	call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
	794	call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
	795	call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
	796	enddo
[2998]	797
[2794]	798	!------------------------
	799	! II Run
[2835]	800	! II_d Update surface and soil temperatures
[2794]	801	!------------------------
[2835]	802	! II_d.1 Update Tsurf
[3028]	803	write(,) "Updating the new Tsurf"
	804	bool_sublim = .false.
	805	Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
	806	do ig = 1,ngrid
	807	do islope = 1,nslope
	808	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
	809	if (latitude_deg(ig) > 0) then
	810	do ig_loop = ig,ngrid
	811	do islope_loop = islope,iflat,-1
	812	if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
	813	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
	814	bool_sublim = .true.
	815	exit
	816	endif
	817	enddo
	818	if (bool_sublim) exit
	819	enddo
	820	else
	821	do ig_loop = ig,1,-1
	822	do islope_loop = islope,iflat
	823	if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
	824	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
	825	bool_sublim = .true.
	826	exit
	827	endif
	828	enddo
	829	if (bool_sublim) exit
	830	enddo
[2835]	831	endif
[3028]	832	initial_co2_ice(ig,islope) = 0
	833	if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
	834	albedo(ig,1,islope) = albedo_h2o_frost
	835	albedo(ig,2,islope) = albedo_h2o_frost
	836	else
	837	albedo(ig,1,islope) = albedodat(ig)
	838	albedo(ig,2,islope) = albedodat(ig)
	839	emis(ig,islope) = emissiv
	840	endif
	841	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
	842	ave = 0.
	843	do t = 1,timelen
	844	if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
	845	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
	846	else
	847	ave = ave + tsurf_GCM_timeseries(ig,islope,t)
	848	endif
[2794]	849	enddo
[3028]	850	tsurf_ave(ig,islope) = ave/timelen
[3032]	851	! set the surface albedo to be the ice albedo
	852	if (latitude_deg(ig) > 0) then
	853	icap = 1
	854	else
	855	icap = 2
	856	endif
	857	albedo(ig,1,islope) = albedice(icap)
	858	albedo(ig,2,islope) = albedice(icap)
	859	emis(ig,islope) = emisice(icap)
[2835]	860	endif
	861	enddo
[3028]	862	enddo
[2794]	863
[3028]	864	do t = 1,timelen
	865	tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:))
	866	enddo
	867	! for the start
	868	do ig = 1,ngrid
[2835]	869	do islope = 1,nslope
[3028]	870	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
[2794]	871	enddo
[3028]	872	enddo
[2794]	873
[3028]	874	do islope = 1,nslope
	875	write(str2(1:2),'(i2.2)') islope
	876	call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
	877	enddo
[2849]	878
[3028]	879	if (soil_pem) then
[2794]	880
[2835]	881	! II_d.2 Update soil temperature
[3028]	882	allocate(TI_locslope(ngrid,nsoilmx_PEM))
	883	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
	884	allocate(Tsurf_locslope(ngrid))
	885	write(,)"Updating soil temperature"
[2794]	886
[3028]	887	! Soil averaged
	888	do islope = 1,nslope
	889	TI_locslope(:,:) = TI_PEM(:,:,islope)
	890	do t = 1,timelen
	891	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
	892	Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
	893	call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
	894	call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
	895	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
	896	do ig = 1,ngrid
	897	do isoil = 1,nsoilmx_PEM
	898	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)
	899	if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1)
	900	enddo
	901	enddo
	902	enddo
	903	enddo
	904	tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
	905	watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
[2794]	906
[3028]	907	write(,) "Update of soil temperature done"
[2888]	908
[3028]	909	deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
	910	write(,) "Compute ice table"
[2849]	911
[2835]	912	! II_d.3 Update the ice table
[3031]	913	porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:)
[3028]	914	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
[3031]	915
	916	call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM, &
[3039]	917	delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
[3031]	918
[3028]	919	write(,) "Update soil propreties"
[2937]	920
[2835]	921	! II_d.4 Update the soil thermal properties
[3028]	922	call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, &
	923	porefillingice_depth,porefillingice_thickness,TI_PEM)
[2794]	924
[2835]	925	! II_d.5 Update the mass of the regolith adsorbded
[3028]	926	if (adsorption_pem) then
	927	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
	928	qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, &
	929	q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
[2794]	930
[3028]	931	totmassco2_adsorbded = 0.
	932	totmassh2o_adsorbded = 0.
	933	do ig = 1,ngrid
	934	do islope =1, nslope
	935	do l = 1,nsoilmx_PEM - 1
	936	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	937	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
	938	cell_area(ig)
	939	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	940	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
	941	cell_area(ig)
	942	enddo
	943	enddo
	944	enddo
	945	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
	946	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
	947	endif
	948	endif !soil_pem
	949
[2794]	950	!------------------------
	951	! II Run
[3028]	952	! II_e Update the tendencies
[2794]	953	!------------------------
[3028]	954	write(,) "Adaptation of the new co2 tendencies to the current pressure"
	955	call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries, &
[2794]	956	global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
	957
[2835]	958	!------------------------
	959	! II Run
	960	! II_f Checking the stopping criterion
	961	!------------------------
[3028]	962	call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
[2779]	963
[3028]	964	call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
	965	initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
[2794]	966
[3028]	967	year_iter = year_iter + dt_pem
[3039]	968	i_myear = i_myear + dt_pem
[2794]	969
[3028]	970	write(,) "Checking all the stopping criterion."
	971	if (STOPPING_water) then
	972	write(,) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
	973	criterion_stop = 1
	974	endif
	975	if (STOPPING_1_water) then
	976	write(,) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
	977	criterion_stop = 1
	978	endif
	979	if (STOPPING_co2) then
	980	write(,) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
	981	criterion_stop = 2
	982	endif
	983	if (STOPPING_pressure) then
	984	write(,) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
	985	criterion_stop = 3
	986	endif
	987	if (year_iter >= year_iter_max) then
	988	write(,) "STOPPING because maximum number of iterations reached"
	989	criterion_stop = 4
	990	endif
[3039]	991	if (i_myear >= n_myear) then
	992	write(,) "STOPPING because maximum number of Martian years to be simulated reached"
	993	criterion_stop = 5
	994	endif
[2794]	995
[3028]	996	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then
[2779]	997	exit
[3028]	998	else
	999	write(,) "We continue!"
[3039]	1000	write(,) "Number of iterations of the PEM: year_iter =", year_iter
	1001	write(,) "Number of simulated Martian years: i_myear =", i_myear
[3028]	1002	endif
[2779]	1003
[3028]	1004	global_ave_press_old=global_ave_press_new
[2779]	1005
[3028]	1006	enddo ! big time iteration loop of the pem
	1007	!------------------------------ END RUN --------------------------------
[2779]	1008
[3028]	1009	!------------------------------- OUTPUT --------------------------------
[2794]	1010	!------------------------
	1011	! III Output
[2835]	1012	! III_a Update surface value for the PCM start files
[2794]	1013	!------------------------
[2835]	1014	! III_a.1 Ice update (for startfi)
[2779]	1015
[2888]	1016	! H2O ice
[3028]	1017	do ig = 1,ngrid
	1018	if (watercaptag(ig)) then
	1019	watercap_sum = 0.
	1020	do islope = 1,nslope
	1021	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
	1022	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
	1023	else
	1024	! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
	1025	watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.))
	1026	qsurf(ig,igcm_h2o_ice,islope)=0.
	1027	endif
	1028	watercap_sum = watercap_sum+watercap(ig,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
	1029	watercap(ig,islope) = 0.
	1030	enddo
	1031	water_reservoir(ig) = water_reservoir(ig) + watercap_sum
	1032	endif
	1033	enddo
[2888]	1034
[3028]	1035	do ig = 1,ngrid
	1036	water_sum = 0.
	1037	do islope = 1,nslope
	1038	water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
	1039	enddo
	1040	if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
	1041	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
	1042	watercaptag(ig) = .true.
	1043	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
	1044	do islope = 1,nslope
	1045	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.cos(pidef_slope_mean(islope)/180.)
	1046	enddo
	1047	endif
	1048	else ! let's check that the infinite source of water disapear
	1049	if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then
	1050	watercaptag(ig) = .false.
	1051	do islope = 1,nslope
	1052	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
	1053	enddo
	1054	water_reservoir(ig) = 0.
	1055	endif
	1056	endif
	1057	enddo
[2888]	1058
[2998]	1059	! CO2 ice
[3028]	1060	do ig = 1,ngrid
	1061	do islope = 1,nslope
	1062	if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then
	1063	perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
	1064	qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope)
[3032]	1065	albedo(ig,1,islope) = albedo_perenialco2
	1066	albedo(ig,2,islope) = albedo_perenialco2
[3028]	1067	endif
	1068	enddo
	1069	enddo
[2998]	1070
[2849]	1071	! III_a.2 Tsoil update (for startfi)
[3028]	1072	if (soil_pem) then
	1073	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
	1074	tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
	1075	endif
[2779]	1076
[2835]	1077	! III_a.4 Pressure (for start)
[3028]	1078	do i = 1,ip1jmp1
	1079	ps(i) = ps(i)*global_ave_press_new/global_ave_press_GCM
	1080	enddo
[2794]	1081
[3028]	1082	do i = 1,ngrid
	1083	ps_start_GCM(i) = ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM
	1084	enddo
[2835]	1085
	1086	! III_a.5 Tracer (for start)
[3028]	1087	allocate(zplev_new(ngrid,nlayer + 1))
[2835]	1088
[3028]	1089	do l = 1,nlayer + 1
	1090	do ig = 1,ngrid
	1091	zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
	1092	enddo
	1093	enddo
[2835]	1094
[3028]	1095	do nnq = 1,nqtot
	1096	if (noms(nnq) /= "co2") then
	1097	do l = 1,llm - 1
	1098	do ig = 1,ngrid
	1099	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))
	1100	enddo
	1101	q(:,llm,nnq) = q(:,llm - 1,nnq)
	1102	enddo
	1103	else
	1104	do l = 1,llm - 1
	1105	do ig = 1,ngrid
	1106	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)) &
	1107	+ ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/ &
	1108	(zplev_new(ig,l) - zplev_new(ig,l + 1))
	1109	enddo
	1110	q(:,llm,nnq) = q(:,llm - 1,nnq)
	1111	enddo
	1112	endif
	1113	enddo
[2835]	1114
[3028]	1115	! Conserving the tracers mass for GCM start files
	1116	do nnq = 1,nqtot
	1117	do ig = 1,ngrid
	1118	do l = 1,llm - 1
	1119	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
	1120	extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
	1121	q(ig,l,nnq)=1.
	1122	q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
	1123	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
[2835]	1124	endif
[3028]	1125	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
	1126	enddo
	1127	enddo
	1128	enddo
[2779]	1129
[3039]	1130	if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
[2779]	1131
	1132	!------------------------
[3028]	1133	! III Output
	1134	! III_b Write restart_evol.nc and restartfi_evol.nc
	1135	!------------------------
	1136	! III_b.1 Write restart_evol.nc
[3042]	1137	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
[3028]	1138	pday = day_ini
[3042]	1139	ztime_fin = time_phys
[2779]	1140
[3028]	1141	allocate(p(ip1jmp1,nlayer + 1))
[2980]	1142	#ifndef CPP_1D
[3028]	1143	call pression (ip1jmp1,ap,bp,ps,p)
	1144	call massdair(p,masse)
[3039]	1145	call dynredem0("restart_evol.nc",day_ini,phis)
	1146	call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
[3028]	1147	write(,) "restart_evol.nc has been written"
[2980]	1148	#else
[3028]	1149	do nnq = 1, nqtot
	1150	call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf)
	1151	enddo
[2980]	1152	#endif
	1153
[3028]	1154	! III_b.2 Write restartfi_evol.nc
[2842]	1155	#ifndef CPP_STD
[3028]	1156	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
	1157	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
	1158	inertiedat,def_slope,subslope_dist)
[2779]	1159
[3028]	1160	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
	1161	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
	1162	albedo,emis,q2,qsurf,tauscaling,totcloudfrac, &
	1163	wstar,watercap,perenial_co2ice)
[2842]	1164	#else
[3028]	1165	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
	1166	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
	1167	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[2779]	1168
[3028]	1169	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
	1170	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf, &
	1171	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic, &
	1172	tslab,tsea_ice,sea_ice)
[2842]	1173	#endif
[3028]	1174	write(,) "restartfi_evol.nc has been written"
[2842]	1175
[2794]	1176	!------------------------
	1177	! III Output
[3028]	1178	! III_c Write restartfi_PEM.nc
[2794]	1179	!------------------------
[3028]	1180	call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
	1181	float(day_ini),0.,nslope,def_slope,subslope_dist)
[2779]	1182
[2888]	1183
[3039]	1184	call pemdem1("restartfi_PEM.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
	1185	TI_PEM, porefillingice_depth,porefillingice_thickness, &
[3028]	1186	co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
[3039]	1187	write(,) "restartfi_PEM.nc has been written"
[2897]	1188
[3039]	1189	call info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
[2779]	1190
[3039]	1191	write(,) "The PEM has run for", year_iter, "Martian years."
	1192	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
	1193	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
	1194	write(,) "LL & RV & JBC: so far, so good!"
[2794]	1195
[3028]	1196	deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
	1197	deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
	1198	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
[3031]	1199	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
[3028]	1200	!----------------------------- END OUTPUT ------------------------------
[2897]	1201
[2779]	1202	END PROGRAM pem
[3028]	1203

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

Download in other formats:

Original Format