Context Navigation

pem.F90 @ 3296

Last change on this file since 3296 was 3287, checked in by jbclement, 8 months ago

PEM:

Small correction to make the 3D PEM be able to compile.
Improvement of "launch_pem.sh": a file "kill_launch_pem.sh" is now automatically created which allows the user to kill the process of the launching script in case.

JBC

File size: 55.5 KB

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

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

Download in other formats:

Original Format