Context Navigation

pem.F90 @ 3116

Last change on this file since 3116 was 3114, checked in by jbclement, 21 months ago
PEM: Following r3113, addition of 'nqsoil' and 'qsoil' in the arguments of the subroutines 'phyetat0' and 'physdem1' to be able to compile. JBC
File size: 57.4 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
[3114]	75	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, volcapa, inertiesoil, nqsoil, qsoil
[3028]	76	use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, &
	77	albedodat, zmea, zstd, zsig, zgam, zthe, &
	78	hmons, summit, base,albedo_h2o_frost, &
[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
[3114]	369	call phyetat0(FILE_NAME,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, &
	370	tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, &
[3028]	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))
[3114]	385	allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope))
[3028]	386	allocate(emis_read_generic(ngrid))
	387	allocate(tsurf(ngrid,1))
	388	allocate(qsurf(ngrid,nq,1))
	389	allocate(tsoil(ngrid,nsoilmx,1))
	390	allocate(emis(ngrid,1))
	391	allocate(watercap(ngrid,1))
	392	allocate(watercaptag(ngrid))
	393	allocate(albedo_read_generic(ngrid,2))
	394	allocate(albedo(ngrid,2,1))
	395	allocate(inertiesoil(ngrid,nsoilmx,1))
[3114]	396	call phyetat0(.true.,ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, &
	397	tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, &
	398	qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, &
	399	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
[3065]	400	call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
[2842]	401
[3028]	402	nslope = 1
	403	call ini_comslope_h(ngrid,1)
[2842]	404
[3028]	405	qsurf(:,:,1) = qsurf_read_generic(:,:)
	406	tsurf(:,1) = tsurf_read_generic(:)
	407	tsoil(:,:,1) = tsoil_read_generic(:,:)
	408	emis(:,1) = emis_read_generic(:)
	409	watercap(:,1) = 0.
	410	watercaptag(:) = .false.
	411	albedo(:,1,1) = albedo_read_generic(:,1)
	412	albedo(:,2,1) = albedo_read_generic(:,2)
	413	inertiesoil(:,:,1) = inertiedat(:,:)
[2842]	414
[3028]	415	if (nslope == 1) then
	416	def_slope(1) = 0
	417	def_slope(2) = 0
	418	def_slope_mean = 0
	419	subslope_dist(:,1) = 1.
	420	endif
[2842]	421
[3070]	422	! Remove unphysical values of surface tracer
	423	qsurf(:,:,1) = qsurf_read_generic(:,:)
	424	where (qsurf < 0.) qsurf = 0.
[2842]	425	#endif
	426
[3028]	427	do nnq = 1,nqtot ! Why not using ini_tracer ?
	428	if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq
	429	if (noms(nnq) == "h2o_vap") then
	430	igcm_h2o_vap = nnq
	431	mmol(igcm_h2o_vap)=18.
	432	endif
	433	if (noms(nnq) == "co2") igcm_co2 = nnq
[3065]	434	enddo
[3039]	435	r = 8.314511*1000./mugaz
[3028]	436
[2835]	437	!------------------------
[3028]	438	! I Initialization
[2835]	439	! I_c Subslope parametrisation
	440	!------------------------
[3028]	441	! Define some slope statistics
	442	iflat = 1
	443	do islope = 2,nslope
	444	if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope
	445	enddo
[2794]	446
[3028]	447	write(,) 'Flat slope for islope = ',iflat
	448	write(,) 'corresponding criterium = ',def_slope_mean(iflat)
[2794]	449
[3028]	450	allocate(flag_co2flow(ngrid,nslope))
	451	allocate(flag_co2flow_mesh(ngrid))
	452	allocate(flag_h2oflow(ngrid,nslope))
	453	allocate(flag_h2oflow_mesh(ngrid))
[2835]	454
[3065]	455	flag_co2flow(:,:) = 0
[3028]	456	flag_co2flow_mesh(:) = 0
[3065]	457	flag_h2oflow(:,:) = 0
[3028]	458	flag_h2oflow_mesh(:) = 0
[2835]	459
[2794]	460	!------------------------
[3028]	461	! I Initialization
[3096]	462	! I_d READ PCM data and convert to the physical grid
[3028]	463	!------------------------
[3096]	464	! 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
	465	call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
[2794]	466
[3028]	467	allocate(tsoil_ave(ngrid,nsoilmx,nslope))
	468	allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
	469	allocate(vmr_co2_gcm(ngrid,timelen))
	470	allocate(ps_timeseries(ngrid,timelen))
	471	allocate(min_co2_ice_1(ngrid,nslope))
	472	allocate(min_h2o_ice_1(ngrid,nslope))
	473	allocate(min_co2_ice_2(ngrid,nslope))
	474	allocate(min_h2o_ice_2(ngrid,nslope))
	475	allocate(tsurf_ave_yr1(ngrid,nslope))
	476	allocate(tsurf_ave(ngrid,nslope))
	477	allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
	478	allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
	479	allocate(q_co2_PEM_phys(ngrid,timelen))
	480	allocate(q_h2o_PEM_phys(ngrid,timelen))
	481	allocate(co2_ice_GCM(ngrid,nslope,timelen))
	482	allocate(watersurf_density_ave(ngrid,nslope))
	483	allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
	484	allocate(Tsurfave_before_saved(ngrid,nslope))
	485	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
	486	allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
	487	allocate(delta_co2_adsorbded(ngrid))
[3031]	488	allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
	489	allocate(delta_h2o_icetablesublim(ngrid))
[3028]	490	allocate(delta_h2o_adsorbded(ngrid))
	491	allocate(vmr_co2_pem_phys(ngrid,timelen))
[2794]	492
[3028]	493	write(,) "Downloading data Y1..."
[3096]	494	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]	495	tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
	496	co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
	497	write(,) "Downloading data Y1 done"
[2985]	498
[3096]	499	! 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]	500	write(,) "Downloading data Y2"
[3096]	501	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]	502	tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
	503	co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
	504	write(,) "Downloading data Y2 done"
[2794]	505
[2835]	506	!------------------------
[3028]	507	! I Initialization
	508	! I_e Initialization of the PEM variables and soil
[2835]	509	!------------------------
[3028]	510	call end_comsoil_h_PEM
	511	call ini_comsoil_h_PEM(ngrid,nslope)
	512	call end_adsorption_h_PEM
	513	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
	514	call end_ice_table_porefilling
	515	call ini_ice_table_porefilling(ngrid,nslope)
[2794]	516
[3028]	517	if (soil_pem) then
	518	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
[3070]	519	do l = 1,nsoilmx
	520	tsoil_PEM(:,l,:) = tsoil_ave(:,l,:)
	521	tsoil_phys_PEM_timeseries(:,l,:,:) = tsoil_GCM_timeseries(:,l,:,:)
	522	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,l,:,:)
[3028]	523	enddo
[3070]	524	do l = nsoilmx + 1,nsoilmx_PEM
	525	tsoil_PEM(:,l,:) = tsoil_ave(:,nsoilmx,:)
	526	watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
[3028]	527	enddo
[3070]	528	watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
[3028]	529	endif !soil_pem
	530	deallocate(tsoil_ave,tsoil_GCM_timeseries)
[2794]	531
[2779]	532	!------------------------
[3028]	533	! I Initialization
[2835]	534	! I_f Compute tendencies & Save initial situation
[3028]	535	!------------------------
	536	allocate(tendencies_co2_ice(ngrid,nslope))
	537	allocate(tendencies_co2_ice_ini(ngrid,nslope))
	538	allocate(tendencies_h2o_ice(ngrid,nslope))
[2779]	539
[3028]	540	! Compute the tendencies of the evolution of ice over the years
	541	call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,min_co2_ice_2,tendencies_co2_ice)
	542	tendencies_co2_ice_ini(:,:) = tendencies_co2_ice(:,:)
	543	call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,min_h2o_ice_2,tendencies_h2o_ice)
[2895]	544
[3028]	545	deallocate(min_co2_ice_1,min_co2_ice_2,min_h2o_ice_1,min_h2o_ice_2)
[2779]	546
[2835]	547	!------------------------
[3028]	548	! I Initialization
[2835]	549	! I_g Save initial PCM situation
[3028]	550	!------------------------
	551	allocate(initial_co2_ice_sublim(ngrid,nslope))
	552	allocate(initial_co2_ice(ngrid,nslope))
	553	allocate(initial_h2o_ice(ngrid,nslope))
[2835]	554
[2794]	555	! We save the places where water ice is sublimating
[2835]	556	! We compute the surface of water ice sublimating
[3028]	557	ini_surf_co2 = 0.
	558	ini_surf_h2o = 0.
	559	Total_surface = 0.
	560	do i = 1,ngrid
[3070]	561	Total_surface = Total_surface + cell_area(i)
[3028]	562	do islope = 1,nslope
	563	if (tendencies_co2_ice(i,islope) < 0) then
	564	initial_co2_ice_sublim(i,islope) = 1.
[3070]	565	ini_surf_co2 = ini_surf_co2 + cell_area(i)*subslope_dist(i,islope)
[3028]	566	else
	567	initial_co2_ice_sublim(i,islope) = 0.
	568	endif
	569	if (qsurf(i,igcm_co2,islope) > 0) then
	570	initial_co2_ice(i,islope) = 1.
	571	else
	572	initial_co2_ice(i,islope) = 0.
	573	endif
	574	if (tendencies_h2o_ice(i,islope) < 0) then
	575	initial_h2o_ice(i,islope) = 1.
	576	ini_surf_h2o=ini_surf_h2o + cell_area(i)*subslope_dist(i,islope)
	577	else
	578	initial_h2o_ice(i,islope) = 0.
	579	endif
[2779]	580	enddo
[3028]	581	enddo
[2779]	582
[3028]	583	write(,) "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
	584	write(,) "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
	585	write(,) "Total surface of the planet=", Total_surface
	586	allocate(zplev_gcm(ngrid,nlayer + 1))
[2779]	587
[3070]	588	do ig = 1,ngrid
	589	zplev_gcm(ig,:) = ap(:) + bp(:)*ps_start_GCM(ig)
[3028]	590	enddo
[2779]	591
[3028]	592	global_ave_press_old = 0.
	593	do i = 1,ngrid
[3065]	594	global_ave_press_old = global_ave_press_old + cell_area(i)*ps_start_GCM(i)/Total_surface
[3028]	595	enddo
[2779]	596
[3028]	597	global_ave_press_GCM = global_ave_press_old
	598	global_ave_press_new = global_ave_press_old
	599	write(,) "Initial global average pressure=", global_ave_press_GCM
[2779]	600
	601	!------------------------
[3028]	602	! I Initialization
[3088]	603	! I_h Read the startpem.nc
[3028]	604	!------------------------
[3088]	605	call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
[3028]	606	porefillingice_thickness,tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries, &
	607	tsoil_phys_PEM_timeseries,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_co2,:), &
	608	qsurf(:,igcm_h2o_ice,:),global_ave_press_GCM,watersurf_density_ave,watersoil_density_PEM_ave, &
	609	co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
[2779]	610
[3031]	611	delta_h2o_icetablesublim(:) = 0.
	612
[3028]	613	do ig = 1,ngrid
	614	do islope = 1,nslope
	615	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
	616	enddo
	617	enddo
[2779]	618
[3028]	619	if (adsorption_pem) then
	620	totmassco2_adsorbded = 0.
	621	totmassh2o_adsorbded = 0.
	622	do ig = 1,ngrid
[3070]	623	do islope = 1,nslope
[3028]	624	do l = 1,nsoilmx_PEM - 1
	625	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	626	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	627	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
	628	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
	629	enddo
[2961]	630	enddo
[3028]	631	enddo
[2885]	632
[3028]	633	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
	634	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
	635	endif ! adsorption
[3065]	636	deallocate(tsurf_ave_yr1)
[2794]	637
[2835]	638	!------------------------
[3028]	639	! I Initialization
[2835]	640	! I_i Compute orbit criterion
[3028]	641	!------------------------
[2842]	642	#ifndef CPP_STD
[3050]	643	call iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
[2842]	644	#else
[3050]	645	call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
[2842]	646	#endif
[2794]	647
[3028]	648	if (evol_orbit_pem) then
[3039]	649	call orbit_param_criterion(i_myear,year_iter_max)
[3028]	650	else
	651	year_iter_max = Max_iter_pem
	652	endif
	653	!-------------------------- END INITIALIZATION -------------------------
[2794]	654
[3028]	655	!-------------------------------- RUN ----------------------------------
[2794]	656	!------------------------
	657	! II Run
[3065]	658	! II_a Update pressure, ice and tracers
[2794]	659	!------------------------
[3028]	660	year_iter = 0
[2794]	661
[3039]	662	do while (year_iter < year_iter_max .and. i_myear < n_myear)
[2835]	663	! II.a.1. Compute updated global pressure
[3028]	664	write(,) "Recomputing the new pressure..."
	665	do i = 1,ngrid
	666	do islope = 1,nslope
[3065]	667	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]	668	enddo
	669	enddo
[3065]	670
[3028]	671	if (adsorption_pem) then
	672	do i = 1,ngrid
	673	global_ave_press_new = global_ave_press_new - gcell_area(i)delta_co2_adsorbded(i)/Total_surface
[3050]	674	enddo
[3028]	675	endif
[3050]	676	write(,) 'Global average pressure old time step',global_ave_press_old
	677	write(,) 'Global average pressure new time step',global_ave_press_new
[2835]	678
	679	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
[3070]	680	allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
[3028]	681	write(,) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
	682	do l = 1,nlayer + 1
	683	do ig = 1,ngrid
	684	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
	685	enddo
	686	enddo
[2779]	687
[2835]	688	! II.a.3. Surface pressure timeseries
[3028]	689	write(,) "Recomputing the surface pressure timeserie adapted to the new pressure..."
	690	do ig = 1,ngrid
	691	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
	692	enddo
[2779]	693
[2835]	694	! II.a.4. New pressure levels timeseries
[3028]	695	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
	696	write(,) "Recomputing the new pressure levels timeserie adapted to the new pressure..."
	697	do l = 1,nlayer + 1
	698	do ig = 1,ngrid
	699	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
	700	enddo
	701	enddo
[2779]	702
[2835]	703	! II.a.5. Tracers timeseries
[3028]	704	write(,) "Recomputing of tracer VMR timeseries for the new pressure..."
[2794]	705
[3028]	706	l = 1
	707	do ig = 1,ngrid
	708	do t = 1,timelen
	709	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))/ &
	710	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
	711	if (q_h2o_PEM_phys(ig,t) < 0) q_h2o_PEM_phys(ig,t) = 1.e-30
	712	if (q_h2o_PEM_phys(ig,t) > 1) q_h2o_PEM_phys(ig,t) = 1.
	713	enddo
	714	enddo
[2794]	715
[3028]	716	do ig = 1,ngrid
[3065]	717	do t = 1,timelen
[3028]	718	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]	719	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
	720	+ ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) &
[3028]	721	- (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
	722	(zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
	723	if (q_co2_PEM_phys(ig,t) < 0) then
	724	q_co2_PEM_phys(ig,t) = 1.e-30
	725	elseif (q_co2_PEM_phys(ig,t) > 1) then
	726	q_co2_PEM_phys(ig,t) = 1.
	727	endif
	728	mmean=1/(A*q_co2_PEM_phys(ig,t) + B)
	729	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
	730	enddo
	731	enddo
[2794]	732
[3028]	733	deallocate(zplev_new_timeseries,zplev_old_timeseries)
	734
	735	!------------------------
[2835]	736	! II Run
	737	! II_b Evolution of the ice
[3028]	738	!------------------------
	739	write(,) "Evolution of h2o ice"
[3031]	740	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]	741
[3028]	742	write(,) "Evolution of co2 ice"
	743	call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
[2998]	744
[2794]	745	!------------------------
	746	! II Run
[3028]	747	! II_c CO2 & H2O glaciers flows
[2794]	748	!------------------------
[3028]	749	write(,) "CO2 glacier flows"
	750	if (co2glaciersflow) call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries, &
	751	global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
[3065]	752
[3028]	753	write(,) "H2O glacier flows"
	754	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]	755
[2794]	756	!------------------------
	757	! II Run
[2835]	758	! II_d Update surface and soil temperatures
[2794]	759	!------------------------
[2835]	760	! II_d.1 Update Tsurf
[3028]	761	write(,) "Updating the new Tsurf"
	762	bool_sublim = .false.
	763	Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
	764	do ig = 1,ngrid
	765	do islope = 1,nslope
	766	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
	767	if (latitude_deg(ig) > 0) then
	768	do ig_loop = ig,ngrid
	769	do islope_loop = islope,iflat,-1
	770	if (initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
	771	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
	772	bool_sublim = .true.
	773	exit
	774	endif
	775	enddo
	776	if (bool_sublim) exit
	777	enddo
	778	else
	779	do ig_loop = ig,1,-1
	780	do islope_loop = islope,iflat
	781	if(initial_co2_ice(ig_loop,islope_loop) < 0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop) < 1.e-10) then
	782	tsurf_ave(ig,islope) = tsurf_ave(ig_loop,islope_loop)
	783	bool_sublim = .true.
	784	exit
	785	endif
	786	enddo
	787	if (bool_sublim) exit
	788	enddo
[2835]	789	endif
[3028]	790	initial_co2_ice(ig,islope) = 0
	791	if ((qsurf(ig,igcm_co2,islope) < 1.e-10) .and. (qsurf(ig,igcm_h2o_ice,islope) > frost_albedo_threshold)) then
	792	albedo(ig,1,islope) = albedo_h2o_frost
	793	albedo(ig,2,islope) = albedo_h2o_frost
	794	else
	795	albedo(ig,1,islope) = albedodat(ig)
[3065]	796	albedo(ig,2,islope) = albedodat(ig)
[3028]	797	emis(ig,islope) = emissiv
	798	endif
	799	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
	800	ave = 0.
	801	do t = 1,timelen
	802	if (co2_ice_GCM(ig,islope,t) > 1.e-3) then
	803	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
	804	else
	805	ave = ave + tsurf_GCM_timeseries(ig,islope,t)
	806	endif
[2794]	807	enddo
[3028]	808	tsurf_ave(ig,islope) = ave/timelen
[3032]	809	! set the surface albedo to be the ice albedo
	810	if (latitude_deg(ig) > 0) then
	811	icap = 1
	812	else
	813	icap = 2
	814	endif
	815	albedo(ig,1,islope) = albedice(icap)
	816	albedo(ig,2,islope) = albedice(icap)
	817	emis(ig,islope) = emisice(icap)
[2835]	818	endif
	819	enddo
[3028]	820	enddo
[2794]	821
[3028]	822	do t = 1,timelen
[3065]	823	tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) + (tsurf_ave(:,:) - Tsurfave_before_saved(:,:))
[3028]	824	enddo
	825	! for the start
	826	do ig = 1,ngrid
[2835]	827	do islope = 1,nslope
[3028]	828	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope) - tsurf_ave(ig,islope))
[2794]	829	enddo
[3028]	830	enddo
[2794]	831
[3028]	832	if (soil_pem) then
[2794]	833
[2835]	834	! II_d.2 Update soil temperature
[3028]	835	allocate(TI_locslope(ngrid,nsoilmx_PEM))
	836	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
	837	allocate(Tsurf_locslope(ngrid))
	838	write(,)"Updating soil temperature"
[2794]	839
[3028]	840	! Soil averaged
	841	do islope = 1,nslope
	842	TI_locslope(:,:) = TI_PEM(:,:,islope)
	843	do t = 1,timelen
	844	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
	845	Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
[3076]	846	call soil_pem_compute(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
	847	call soil_pem_compute(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
[3028]	848	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
	849	do ig = 1,ngrid
	850	do isoil = 1,nsoilmx_PEM
	851	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]	852	if (isnan(Tsoil_locslope(ig,isoil))) call abort_pem("PEM - Update Tsoil","NaN detected in Tsoil ",1)
[3028]	853	enddo
	854	enddo
	855	enddo
	856	enddo
[3070]	857	tsoil_PEM(:,:,:) = sum(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
	858	watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
[2794]	859
[3028]	860	write(,) "Update of soil temperature done"
[2888]	861
[3028]	862	deallocate(TI_locslope,Tsoil_locslope,Tsurf_locslope)
	863	write(,) "Compute ice table"
[2849]	864
[2835]	865	! II_d.3 Update the ice table
[3031]	866	porefillingice_thickness_prev_iter(:,:) = porefillingice_thickness(:,:)
[3028]	867	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
[3031]	868
[3070]	869	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]	870
[3028]	871	write(,) "Update soil propreties"
[2937]	872
[2835]	873	! II_d.4 Update the soil thermal properties
[3070]	874	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]	875
[2835]	876	! II_d.5 Update the mass of the regolith adsorbded
[3028]	877	if (adsorption_pem) then
	878	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice, &
	879	qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:),tsoil_PEM,TI_PEM,ps_timeseries, &
	880	q_co2_PEM_phys,q_h2o_PEM_phys,h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
[2794]	881
[3028]	882	totmassco2_adsorbded = 0.
	883	totmassh2o_adsorbded = 0.
	884	do ig = 1,ngrid
	885	do islope =1, nslope
	886	do l = 1,nsoilmx_PEM - 1
	887	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
[3070]	888	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
[3028]	889	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
[3070]	890	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)cell_area(ig)
[3028]	891	enddo
	892	enddo
	893	enddo
	894	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
	895	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
	896	endif
	897	endif !soil_pem
	898
[2794]	899	!------------------------
	900	! II Run
[3088]	901	! II_e Outputs
[2794]	902	!------------------------
[3088]	903	call writediagpem(ngrid,'ps_ave','Global average pressure','Pa',0,(/global_ave_press_new/))
	904	do islope = 1,nslope
	905	write(str2(1:2),'(i2.2)') islope
	906	call writediagpem(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
	907	call writediagpem(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
	908	call writediagpem(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
	909	call writediagpem(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
	910	call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
	911	call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
	912	enddo
	913
	914	!------------------------
	915	! II Run
	916	! II_f Update the tendencies
	917	!------------------------
[3028]	918	write(,) "Adaptation of the new co2 tendencies to the current pressure"
[3076]	919	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, &
	920	global_ave_press_GCM,global_ave_press_new)
[2794]	921
[2835]	922	!------------------------
	923	! II Run
[3088]	924	! II_g Checking the stopping criterion
[2835]	925	!------------------------
[3028]	926	call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
[2779]	927
[3028]	928	call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
	929	initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
[2794]	930
[3028]	931	year_iter = year_iter + dt_pem
[3039]	932	i_myear = i_myear + dt_pem
[2794]	933
[3028]	934	write(,) "Checking all the stopping criterion."
	935	if (STOPPING_water) then
	936	write(,) "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
	937	criterion_stop = 1
	938	endif
	939	if (STOPPING_1_water) then
	940	write(,) "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
	941	criterion_stop = 1
	942	endif
	943	if (STOPPING_co2) then
	944	write(,) "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
	945	criterion_stop = 2
	946	endif
	947	if (STOPPING_pressure) then
	948	write(,) "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
	949	criterion_stop = 3
	950	endif
	951	if (year_iter >= year_iter_max) then
	952	write(,) "STOPPING because maximum number of iterations reached"
	953	criterion_stop = 4
	954	endif
[3039]	955	if (i_myear >= n_myear) then
	956	write(,) "STOPPING because maximum number of Martian years to be simulated reached"
	957	criterion_stop = 5
	958	endif
[2794]	959
[3028]	960	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then
[2779]	961	exit
[3028]	962	else
	963	write(,) "We continue!"
[3039]	964	write(,) "Number of iterations of the PEM: year_iter =", year_iter
	965	write(,) "Number of simulated Martian years: i_myear =", i_myear
[3028]	966	endif
[2779]	967
[3065]	968	global_ave_press_old = global_ave_press_new
[2779]	969
[3028]	970	enddo ! big time iteration loop of the pem
	971	!------------------------------ END RUN --------------------------------
[2779]	972
[3028]	973	!------------------------------- OUTPUT --------------------------------
[2794]	974	!------------------------
	975	! III Output
[2835]	976	! III_a Update surface value for the PCM start files
[2794]	977	!------------------------
[2835]	978	! III_a.1 Ice update (for startfi)
[2779]	979
[2888]	980	! H2O ice
[3028]	981	do ig = 1,ngrid
	982	if (watercaptag(ig)) then
	983	watercap_sum = 0.
	984	do islope = 1,nslope
	985	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
	986	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
	987	else
	988	! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
	989	watercap(ig,islope) = watercap(ig,islope) + qsurf(ig,igcm_h2o_ice,islope) - (watercap(ig,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.))
	990	qsurf(ig,igcm_h2o_ice,islope)=0.
	991	endif
	992	watercap_sum = watercap_sum+watercap(ig,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
	993	watercap(ig,islope) = 0.
	994	enddo
	995	water_reservoir(ig) = water_reservoir(ig) + watercap_sum
	996	endif
	997	enddo
[2888]	998
[3028]	999	do ig = 1,ngrid
	1000	water_sum = 0.
	1001	do islope = 1,nslope
	1002	water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
	1003	enddo
[3065]	1004	if (.not. watercaptag(ig)) then ! let's check if we have an 'infinite' source of water that has been forming.
[3028]	1005	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
	1006	watercaptag(ig) = .true.
	1007	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
	1008	do islope = 1,nslope
	1009	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.cos(pidef_slope_mean(islope)/180.)
	1010	enddo
	1011	endif
	1012	else ! let's check that the infinite source of water disapear
	1013	if ((water_sum + water_reservoir(ig)) < threshold_water_frost2perenial) then
	1014	watercaptag(ig) = .false.
	1015	do islope = 1,nslope
	1016	qsurf(ig,igcm_h2o_ice,islope) = qsurf(ig,igcm_h2o_ice,islope) + water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
	1017	enddo
	1018	water_reservoir(ig) = 0.
	1019	endif
	1020	endif
	1021	enddo
[2888]	1022
[2998]	1023	! CO2 ice
[3028]	1024	do ig = 1,ngrid
	1025	do islope = 1,nslope
	1026	if (qsurf(ig,igcm_co2,islope) > threshold_co2_frost2perenial) then
	1027	perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
	1028	qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope)
[3032]	1029	albedo(ig,1,islope) = albedo_perenialco2
	1030	albedo(ig,2,islope) = albedo_perenialco2
[3028]	1031	endif
	1032	enddo
	1033	enddo
[2998]	1034
[2849]	1035	! III_a.2 Tsoil update (for startfi)
[3028]	1036	if (soil_pem) then
	1037	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
[3065]	1038	tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
[3028]	1039	endif
[2779]	1040
[2835]	1041	! III_a.4 Pressure (for start)
[3065]	1042	ps(:) = ps(:)*global_ave_press_new/global_ave_press_GCM
	1043	ps_start_GCM(:) = ps_start_GCM(:)*global_ave_press_new/global_ave_press_GCM
[2794]	1044
[2835]	1045	! III_a.5 Tracer (for start)
[3028]	1046	allocate(zplev_new(ngrid,nlayer + 1))
[2835]	1047
[3028]	1048	do l = 1,nlayer + 1
[3065]	1049	zplev_new(:,l) = ap(l) + bp(l)*ps_start_GCM(:)
[3028]	1050	enddo
[2835]	1051
[3028]	1052	do nnq = 1,nqtot
	1053	if (noms(nnq) /= "co2") then
	1054	do l = 1,llm - 1
	1055	do ig = 1,ngrid
	1056	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))
	1057	enddo
	1058	q(:,llm,nnq) = q(:,llm - 1,nnq)
	1059	enddo
	1060	else
	1061	do l = 1,llm - 1
	1062	do ig = 1,ngrid
	1063	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]	1064	+ ((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]	1065	enddo
	1066	q(:,llm,nnq) = q(:,llm - 1,nnq)
	1067	enddo
	1068	endif
	1069	enddo
[2835]	1070
[3096]	1071	! Conserving the tracers mass for PCM start files
[3028]	1072	do nnq = 1,nqtot
	1073	do ig = 1,ngrid
	1074	do l = 1,llm - 1
	1075	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]	1076	extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1))
	1077	q(ig,l,nnq) = 1.
	1078	q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2))
[3028]	1079	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number"
[2835]	1080	endif
[3028]	1081	if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30
	1082	enddo
	1083	enddo
	1084	enddo
[2779]	1085
[3039]	1086	if (evol_orbit_pem) call recomp_orb_param(i_myear,year_iter)
[2779]	1087
	1088	!------------------------
[3028]	1089	! III Output
	1090	! III_b Write restart_evol.nc and restartfi_evol.nc
	1091	!------------------------
	1092	! III_b.1 Write restart_evol.nc
[3042]	1093	ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys
[3028]	1094	pday = day_ini
[3042]	1095	ztime_fin = time_phys
[2779]	1096
[3028]	1097	allocate(p(ip1jmp1,nlayer + 1))
[2980]	1098	#ifndef CPP_1D
[3028]	1099	call pression (ip1jmp1,ap,bp,ps,p)
	1100	call massdair(p,masse)
[3039]	1101	call dynredem0("restart_evol.nc",day_ini,phis)
	1102	call dynredem1("restart_evol.nc",time_0,vcov,ucov,teta,q,masse,ps)
[3028]	1103	write(,) "restart_evol.nc has been written"
[2980]	1104	#else
[3069]	1105	call writerestart1D('restart1D_evol.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q)
[3065]	1106	write(,) "restart1D_evol.txt has been written"
[2980]	1107	#endif
	1108
[3028]	1109	! III_b.2 Write restartfi_evol.nc
[2842]	1110	#ifndef CPP_STD
[3028]	1111	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
	1112	nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, &
	1113	inertiedat,def_slope,subslope_dist)
[2779]	1114
[3114]	1115	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
	1116	ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, &
	1117	albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, &
[3028]	1118	wstar,watercap,perenial_co2ice)
[2842]	1119	#else
[3028]	1120	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid, &
	1121	nlayer,nq,ptimestep,pday,time_phys,cell_area, &
	1122	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
[2779]	1123
[3114]	1124	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, &
	1125	ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, &
	1126	cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, &
	1127	tsea_ice,sea_ice)
[2842]	1128	#endif
[3028]	1129	write(,) "restartfi_evol.nc has been written"
[2842]	1130
[2794]	1131	!------------------------
	1132	! III Output
[3088]	1133	! III_c Write restartpem.nc
[2794]	1134	!------------------------
[3088]	1135	call pemdem0("restartpem.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
[3028]	1136	float(day_ini),0.,nslope,def_slope,subslope_dist)
[3088]	1137	call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
[3039]	1138	TI_PEM, porefillingice_depth,porefillingice_thickness, &
[3028]	1139	co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
[3088]	1140	write(,) "restartpem.nc has been written"
[2897]	1141
[3096]	1142	call info_PEM(year_iter,criterion_stop,i_myear,n_myear)
[2779]	1143
[3039]	1144	write(,) "The PEM has run for", year_iter, "Martian years."
	1145	write(,) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years."
	1146	write(,) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years."
	1147	write(,) "LL & RV & JBC: so far, so good!"
[2794]	1148
[3028]	1149	deallocate(vmr_co2_gcm,ps_timeseries,tsurf_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys)
	1150	deallocate(co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries,Tsurfave_before_saved)
	1151	deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_ave)
[3031]	1152	deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_pem_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
[3028]	1153	!----------------------------- END OUTPUT ------------------------------
[2897]	1154
[2779]	1155	END PROGRAM pem

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

Download in other formats:

Original Format