Context Navigation

pem.F90 @ 2863

Last change on this file since 2863 was 2863, checked in by llange, 23 months ago
PEM H2O adsorption is now computed. Total mass of H2o adsorbded is written in the restart_PEM. +Minor fixes and edit LL & RV
File size: 59.1 KB

Line
1
2	!------------------------
3
4	! I Initialisation
5	! I_a READ run.def
6	! I_b READ of start_evol.nc and starfi_evol.nc
7	! I_c Subslope parametrisation
8	! I_d READ GCM data and convert to the physical grid
9	! I_e Initialisation of the PEM variable and soil
10	! I_f Compute tendencies & Save initial situation
11	! I_g Save initial PCM situation
12	! I_h Read the PEMstart
13	! I_i Compute orbit criterion
14
15	! II Run
16	! II_a update pressure, ice and tracers
17	! II_b Evolution of the ice
18	! II_c CO2 glaciers flows
19	! II_d Update surface and soil temperatures
20	! II_e Update the tendencies
21	! II_f Checking the stopping criterion
22
23	! III Output
24	! III_a Update surface value for the PCM start files
25	! III_b Write start and starfi.nc
26	! III_c Write start_pem
27
28	!------------------------
29
30	PROGRAM pem
31
32	!module needed for INITIALISATION
33	#ifndef CPP_STD
34	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
35	use surfdat_h, only: tsurf, co2ice, emis,&
36	qsurf,watercap, ini_surfdat_h, &
37	albedodat, zmea, zstd, zsig, zgam, zthe, &
38	hmons, summit, base,albedo_h2o_frost, &
39	frost_albedo_threshold,emissiv
40	use dimradmars_mod, only: totcloudfrac, albedo
41	use dust_param_mod, only: tauscaling
42	use tracer_mod, only: noms,igcm_h2o_ice ! tracer names
43	#else
44	use comsoil_h, only: nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
45	use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, &
46	emissiv
47	use tracer_h, only: noms,igcm_h2o_ice,igcm_co2_ice ! tracer names
48	use phys_state_var_mod
49	#endif
50	use phyetat0_mod, only: phyetat0
51	use phyredem, only: physdem0, physdem1
52	use turb_mod, only: q2, wstar
53	use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
54	nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
55	nf90_inquire_dimension,nf90_close
56
57	! For phyredem :
58	USE control_mod, ONLY: iphysiq, day_step,nsplit_phys
59	USE iniphysiq_mod, ONLY: iniphysiq
60	USE logic_mod, ONLY: iflag_phys
61	#ifndef CPP_STD
62	use mod_phys_lmdz_para, only: is_parallel, is_sequential, &
63	is_mpi_root, is_omp_root, &
64	is_master
65	use planete_h, only: aphelie, periheli, year_day, peri_day, &
66	obliquit
67	#else
68	! USE comcstfi_mod, ONLY: rad,g,r,cpp,pi
69	! USE inifis_mod, ONLY: inifis
70	use planete_mod, only: apoastr, periastr, year_day, peri_day, &
71	obliquit
72	#endif
73	USE mod_const_mpi, ONLY: COMM_LMDZ
74	USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, &
75	subslope_dist,co2ice_slope, &
76	tsurf_slope,tsoil_slope,fluxgrd_slope,&
77	fluxrad_sky_slope,sky_slope,callsubslope,&
78	co2iceflow, beta_slope, capcal_slope,&
79	albedo_slope,emiss_slope,qsurf_slope,&
80	iflat,major_slope,ini_comslope_h
81	use time_phylmdz_mod, only: daysec,dtphys
82	USE comconst_mod, ONLY: rad,g,r,cpp,pi
83	USE infotrac
84	USE geometry_mod, only: latitude_deg
85	use conf_pem_mod, only: conf_pem
86	use pemredem, only: pemdem1
87	use co2glaciers_mod,only: co2glaciers_evol
88	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
89	use comsoil_h_PEM, only: soil_pem,ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
90	TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM, & ! soil thermal inertia
91	tsoil_PEM, mlayer_PEM,layer_PEM, & !number of subsurface layers, soil mid layer depths
92	fluxgeo, co2_adsorbded_phys, h2o_adsorbded_phys ! geothermal flux, mass of co2 & h2o in the regolith
93	use adsorption_mod, only : regolith_adsorption
94
95	!!! For orbit parameters
96	USE temps_mod_evol, ONLY: dt_pem, evol_orbit_pem, Max_iter_pem
97	use orbit_param_criterion_mod, only : orbit_param_criterion
98	use recomp_orb_param_mod, only: recomp_orb_param
99
100
101	IMPLICIT NONE
102
103	include "dimensions.h"
104	include "paramet.h"
105
106	INTEGER ngridmx
107	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
108
109	include "comdissnew.h"
110	include "comgeom.h"
111	include "iniprint.h"
112	! Same variable's name as in the GCM
113
114	INTEGER :: ngrid !Number of physical grid points
115	INTEGER :: nlayer !Number of vertical layer
116	INTEGER :: nq !Number of tracer
117	INTEGER :: day_ini !First day of the simulation
118	REAL :: pday !Physical day
119	REAL :: time_phys !Same as GCM
120	REAL :: ptimestep !Same as GCM
121	REAL :: ztime_fin !Same as GCM
122
123	! Variable for reading start.nc
124	character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
125	! variables dynamiques
126	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
127	REAL teta(ip1jmp1,llm) ! temperature potentielle
128	REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
129	REAL ps(ip1jmp1) ! pression au sol
130	REAL, dimension(:),allocatable :: ps_phys !(ngrid) ! pression au sol
131	REAL, dimension(:,:),allocatable :: ps_phys_timeseries !(ngrid x timelen) ! pression au sol instantannées
132	REAL, dimension(:,:),allocatable :: ps_phys_timeseries_yr1 !(ngrid x timelen) ! pression au sol instantannées for the first year of the gcm
133
134	REAL masse(ip1jmp1,llm) ! masse d'air
135	REAL phis(ip1jmp1) ! geopotentiel au sol
136	REAL time_0
137
138	! Variable for reading starfi.nc
139
140	character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
141	integer :: ncid, varid,status !Variable for handling opening of files
142	integer :: phydimid, subdimid, nlayerdimid, nqdimid !Variable ID for Netcdf files
143	integer :: lonvarid, latvarid, areavarid,sdvarid !Variable ID for Netcdf files
144	integer :: apvarid,bpvarid !Variable ID for Netcdf files
145
146	! Variable for reading starfi.nc and writting restartfi.nc
147
148	REAL, dimension(:),allocatable :: longitude !Longitude read in FILE_NAME and written in restartfi
149	REAL, dimension(:),allocatable :: latitude !Latitude read in FILE_NAME and written in restartfi
150	REAL, dimension(:),allocatable :: ap !Coefficient ap read in FILE_NAME_start and written in restart
151	REAL, dimension(:),allocatable :: bp !Coefficient bp read in FILE_NAME_start and written in restart
152	REAL, dimension(:),allocatable :: cell_area !Cell_area read in FILE_NAME and written in restartfi
153	REAL :: Total_surface !Total surface of the planet
154
155	! Variable for h2o_ice evolution
156
157	REAL , dimension(:,:), allocatable :: tendencies_h2o_ice ! LON x LAT field : Tendency of evolution of perenial ice over a year
158	REAL, dimension(:),allocatable :: tendencies_h2o_ice_phys ! physical point field : Tendency of evolution of perenial ice over a year
159
160	REAL , dimension(:,:), allocatable :: tendencies_co2_ice ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
161	REAL, dimension(:),allocatable :: tendencies_co2_ice_phys ! physical point field : Tendency of evolution of perenial co2 ice over a year
162
163	REAL :: ini_surf ! Initial surface of sublimating water ice
164	REAL :: ini_surf_h2o ! Initial surface of sublimating water ice
165	REAL, dimension(:),allocatable :: initial_h2o_ice ! physical point field : Logical array indicating sublimating point
166
167	REAL :: ini_surf_co2 ! Initial surface of sublimating co2 ice
168	REAL, dimension(:),allocatable :: initial_co2_ice ! physical point field : Logical array indicating sublimating point of co2 ice
169
170	REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1 ! LON x LAT field : minimum of water ice at each point for the first year
171	REAL , dimension(:,:), allocatable :: min_h2o_ice_s_2 ! LON x LAT field : minimum of water ice at each point for the second year
172
173	REAL , dimension(:,:), allocatable :: min_co2_ice_s_1 ! LON x LAT field : minimum of water ice at each point for the first year
174	REAL , dimension(:,:), allocatable :: min_co2_ice_s_2 ! LON x LAT field : minimum of water ice at each point for the second year
175
176	REAL :: global_ave_press_GCM ! constant: global average pressure retrieved in the GCM [Pa]
177	REAL :: global_ave_press_old ! constant: Global average pressure of initial/previous time step
178	REAL :: global_ave_press_new ! constant: Global average pressure of current time step
179
180	REAL , dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field : mass of the atmospheric layers in the pem at current time step [kg/m^2]
181	REAL , dimension(:,:), allocatable :: zplev_gcm ! same but retrieved from the gcm [kg/m^2]
182	REAL , dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
183	REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
184
185	LOGICAL :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
186	LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
187	LOGICAL :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
188	LOGICAL :: STOPPING_1_co2 ! Logical : is there still water ice to sublimate?
189
190	REAL, dimension(:,:,:),allocatable :: q_co2_GCM ! Lon x Lat x Time : mass mixing ratio of co2 in the first layer [kg/kg]
191
192	real,save :: m_co2, m_noco2, A , B, mmean ! Molar mass of co2, no co2 (Ar, ...), intermediate A, B for computations, mean molar mass of the layer [mol/kg]
193	real ,allocatable :: vmr_co2_gcm_phys(:,:) ! Physics x Times co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
194	real ,allocatable :: vmr_co2_pem_phys(:,:) ! Physics x Times co2 volume mixing ratio used in the PEM
195	real ,allocatable :: q_h2o_GCM_phys(:,:) ! Physics x Times h2o mass mixing ratio in the first layer from the GCM [kg/kg]
196	real ,allocatable :: q_co2_GCM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer from the GCM [kg/kg]
197	real ,allocatable :: q_co2_PEM_phys(:,:) ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM [kg/kg]
198	REAL, ALLOCATABLE :: ps_GCM(:,:,:) ! Lon x Lat x Times: surface pressure from the GCM [Pa]
199	REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:) ! Lon x Lat x Times: surface pressure from the 1st year of the GCM [Pa]
200	REAL, ALLOCATABLE :: vmr_co2_gcm(:,:,:) ! Lon x Lat x Times: co2 volumemixing ratio retrieve from the gcm [m^3/m^3]
201	REAL, ALLOCATABLE :: q_h2o_GCM(:,:,:) ! Lon x Lat x Times: h2o volume mixing ratio retrieved from the GCM
202	REAL ,allocatable :: q_h2o_PEM_phys(:,:) ! Physics x Times: h2o mass mixing ratio computed in the PEM
203	integer :: timelen ! # time samples
204	REAL :: ave ! intermediate varibale to compute average
205
206	REAL, ALLOCATABLE :: p(:,:) ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1)
207	REAL :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1
208	REAL :: beta_clap_co2 = 3182.48 ! clapeyron's law for CO2
209	REAL :: alpha_clap_co2 = 23.3494 ! Clapeyron's law for CO2
210
211
212	!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
213	REAL ,allocatable :: watercap_slope(:,:) ! Physics x Nslope: watercap per slope
214	REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1 ! LON x LAT field : minimum of co2 ice at each point for the first year [kg/m^2]
215	REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2 ! LON x LAT field : minimum of co2 ice at each point for the second year [kg/m^2]
216	REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1 ! LON x LAT field : minimum of water ice at each point for the first year [kg/m^2]
217	REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2 ! LON x LAT field : minimum of water ice at each point for the second year [kg/m^2]
218	REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM [kg/m^2]
219	REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope ! Physics x NSLOPE x Times field : co2 ice given by the GCM [kg/m^2]
220	REAL, dimension(:,:),allocatable :: initial_co2_ice_sublim_slope ! physical point field : Logical array indicating sublimating point of co2 ice
221	REAL, dimension(:,:),allocatable :: initial_h2o_ice_slope ! physical point field : Logical array indicating if there is water ice at initial state
222	REAL, dimension(:,:),allocatable :: initial_co2_ice_slope ! physical point field : Logical array indicating if there is co2 ice at initial state
223	REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope ! LON x LAT x nslope field : Tendency of evolution of perenial co2 ice over a year
224	REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope ! LON x LAT x slope field : Tendency of evolution of perenial water ice over a year
225	REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope ! physical point xslope field : Tendency of evolution of perenial co2 ice over a year
226	REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
227	REAL, dimension(:,:),allocatable :: tendencies_h2o_ice_phys_slope ! physical pointx slope field : Tendency of evolution of perenial co2 ice
228	REAL , dimension(:,:), allocatable :: flag_co2flow(:,:) !(ngrid,nslope) ! To flag where there is a glacier flow
229	REAL , dimension(:), allocatable :: flag_co2flow_mesh(:) !(ngrid) ! To flag where there is a glacier flow
230
231	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL
232
233	REAL, ALLOCATABLE :: tsurf_ave(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
234	REAL, ALLOCATABLE :: tsurf_ave_phys(:,:) ! Physic x LAT x SLOPE field : Averaged Surface Temperature [K]
235	REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
236	REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature during 1st year of the GCM [K]
237	REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) ! Physics x SLOPE field : Averaged Soil Temperature during 1st year [K]
238	REAL, ALLOCATABLE :: TI_GCM(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Thermal Inertia [SI]
239	REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:) ! LON X LAT x SLOPE XTULES field : Surface Temperature in timeseries [K]
240	REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) ! Physic x SLOPE XTULES field : NOn averaged Surf Temperature [K]
241
242	REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
243	REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
244	REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
245	REAL, ALLOCATABLE :: tsurf_ave_phys_yr1(:,:) ! Physic SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
246	REAL, ALLOCATABLE :: inertiesoil(:,:) !Physic x Depth Thermal inertia of the mesh for restart [SI]
247
248	REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:) ! Physic x Depth x Slope Averaged GCM Thermal Inertia per slope [SI]
249	REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Same but for the start
250
251	REAL,ALLOCATABLE :: ice_depth(:,:) ! Physic x SLope: Ice table depth [m]
252	REAL,ALLOCATABLE :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI]
253	REAL,ALLOCATABLE :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
254	REAL,ALLOCATABLE :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
255	REAL,ALLOCATABLE :: alph_locslope(:,:) ! Physic x Soil: Intermediate to compute Tsoil [1]
256	REAL,ALLOCATABLE :: beta_locslope(:,:) ! Physic x Soil : Intermediate tocompute Tsoil [K]
257	REAL,ALLOCATABLE :: watersurf_density_timeseries(:,:,:,:) ! Lon x Lat x Slope x Times: water surface density, time series [kg/m^3]
258	REAL,ALLOCATABLE :: watersoil_density_timeseries(:,:,:,:,:) ! Lon x Lat x Soil x Slope x Times water soil density, time series [kg /m^3]
259	REAL,ALLOCATABLE :: watersurf_density_phys_timeseries(:,:,:) ! Physic x Slope x Times, water surface density, time series [kg/m^3]
260	REAL,ALLOCATABLE :: watersurf_density_phys_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
261	REAL,ALLOCATABLE :: watersoil_density_phys_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
262	REAL,ALLOCATABLE :: watersoil_density_phys_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
263	REAL,ALLOCATABLE :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K]
264	REAL, ALLOCATABLE :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
265	REAL, ALLOCATABLE :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
266	REAL :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
267	REAL :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
268	REAL :: alpha_clap_h2o = -6143.7 ! coeffcient to compute psat, from Murphie et Kood 2005 [K]
269	REAL :: beta_clap_h2o = 28.9074 ! coefficient to compute psat, from Murphie et Kood 2005 [1]
270	LOGICAL :: bool_sublim ! logical to check if there is sublimation or not
271
272	!! Some parameters for the PEM run
273	REAL, PARAMETER :: year_step = 1 ! timestep for the pem
274	INTEGER :: year_iter ! number of iteration
275	INTEGER :: year_iter_max ! maximum number of iterations before stopping
276	REAL :: timestep ! timestep [s]
277	#ifdef CPP_STD
278	! INTEGER :: nsplit_phys=1
279	! LOGICAL :: iflag_phys=.true.
280	REAL :: frost_albedo_threshold=0.05 ! frost albedo threeshold to convert fresh frost to old ice
281	REAL :: albedo_h2o_frost ! albedo of h2o frost
282	REAL,ALLOCATABLE :: co2ice(:) ! Physics: co2 ice mesh averaged [kg/m^2]
283	#endif
284
285	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
286
287	! Loop variable
288	INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil
289	#ifndef CPP_STD
290	! Parallel variables
291	is_sequential=.true.
292	is_parallel=.false.
293	is_mpi_root=.true.
294	is_omp_root=.true.
295	is_master=.true.
296	#endif
297
298	day_ini=0 !test
299	time_phys=0. !test
300
301	! Some constants
302
303	ngrid=ngridmx
304	nlayer=llm
305
306	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
307	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
308	A =(1/m_co2 - 1/m_noco2)
309	B=1/m_noco2
310
311	year_day=669
312	daysec=88775.
313	dtphys=0
314	timestep=year_day*daysec/year_step
315
316	!------------------------
317
318	! I Initialisation
319	! I_a READ run.def
320
321	!------------------------
322
323	!----------------------------READ run.def ---------------------
324	CALL conf_gcm( 99, .TRUE. )
325	CALL conf_pem
326
327	call infotrac_init
328	nq=nqtot
329
330	!------------------------
331
332	! I Initialisation
333	! I_a READ run.def
334	! I_b READ of start_evol.nc and starfi_evol.nc
335
336	!------------------------
337
338	!----------------------------Initialisation : READ some constant of startfi_evol.nc ---------------------
339
340	! In the gcm, these values are given to the physic by the dynamic.
341	! Here we simply read them in the startfi_evol.nc file
342	status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
343
344	allocate(longitude(ngrid))
345	allocate(latitude(ngrid))
346	allocate(cell_area(ngrid))
347
348	status = nf90_inq_varid(ncid, "longitude", lonvarid)
349	status = nf90_get_var(ncid, lonvarid, longitude)
350
351	status = nf90_inq_varid(ncid, "latitude", latvarid)
352	status = nf90_get_var(ncid, latvarid, latitude)
353
354	status = nf90_inq_varid(ncid, "area", areavarid)
355	status = nf90_get_var(ncid, areavarid, cell_area)
356
357	call ini_comsoil_h(ngrid)
358
359	status = nf90_inq_varid(ncid, "soildepth", sdvarid)
360	status = nf90_get_var(ncid, sdvarid, mlayer)
361
362	status =nf90_close(ncid)
363
364	!----------------------------READ start.nc ---------------------
365
366	allocate(q(ip1jmp1,llm,nqtot))
367	CALL dynetat0(FILE_NAME_start,vcov,ucov, &
368	teta,q,masse,ps,phis, time_0)
369
370	CALL iniconst !new
371	CALL inigeom
372	allocate(ap(nlayer+1))
373	allocate(bp(nlayer+1))
374	status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
375	status = nf90_inq_varid(ncid, "ap", apvarid)
376	status = nf90_get_var(ncid, apvarid, ap)
377	status = nf90_inq_varid(ncid, "bp", bpvarid)
378	status = nf90_get_var(ncid, bpvarid, bp)
379	status =nf90_close(ncid)
380
381	CALL iniphysiq(iim,jjm,llm, &
382	(jjm-1)*iim+2,comm_lmdz, &
383	daysec,day_ini,dtphys/nsplit_phys, &
384	rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
385	iflag_phys)
386
387	!----------------------------READ startfi.nc ---------------------
388
389	! First we read the initial state (starfi.nc)
390
391	allocate(watercap_slope(ngrid,nslope))
392	allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
393	allocate(inertiesoil(ngrid,nsoilmx))
394
395	#ifndef CPP_STD
396	CALL phyetat0 (FILE_NAME,0,0, &
397	nsoilmx,ngrid,nlayer,nq, &
398	day_ini,time_phys, &
399	tsurf,tsoil,albedo,emis, &
400	q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar, &
401	watercap,inertiesoil,nslope,tsurf_slope, &
402	tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
403	subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM_start, &
404	qsurf_slope,watercap_slope)
405
406	if(soil_pem) then
407	deallocate(TI_GCM_start) !not used then
408	endif
409
410	! Remove unphysical values of surface tracer
411	DO i=1,ngrid
412	DO nnq=1,nqtot
413	DO islope=1,nslope
414	if(qsurf_slope(i,nnq,islope).LT.0) then
415	qsurf_slope(i,nnq,islope)=0.
416	endif
417	enddo
418	enddo
419	enddo
420
421	#else
422	call phys_state_var_init(nq)
423	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
424	call initracer(ngrid,nq)
425	call iniaerosol()
426	call phyetat0(.true., &
427	ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq, &
428	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
429	cloudfrac,totcloudfrac,hice, &
430	rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
431	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
432
433	call ini_comslope_h(ngrid,nsoilmx,nq)
434
435	allocate(co2ice(ngrid))
436	co2ice(:)=qsurf(:,igcm_co2_ice)
437	co2ice_slope(:,1)=co2ice(:)
438	tsurf_slope(:,1)=tsurf(:)
439
440	if (nslope.eq.1) then
441	def_slope(1) = 0
442	def_slope(2) = 0
443	def_slope_mean=0
444	subslope_dist(:,1) = 1.
445	endif
446
447	! Remove unphysical values of surface tracer
448	DO i=1,ngrid
449	DO nnq=1,nqtot
450	qsurf_slope(i,nnq,1)=qsurf(i,nnq)
451	if(qsurf(i,nnq).LT.0) then
452	qsurf(i,nnq)=0.
453	endif
454	enddo
455	enddo
456	#endif
457
458	DO nnq=1,nqtot
459	if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
460	ENDDO
461
462	!------------------------
463
464	! I Initialisation
465	! I_a READ run.def
466	! I_b READ of start_evol.nc and starfi_evol.nc
467	! I_c Subslope parametrisation
468
469	!------------------------
470
471	!----------------------------Subslope parametrisation definition ---------------------
472
473	! Define some slope statistics
474	iflat=1
475	DO islope=2,nslope
476	IF(abs(def_slope_mean(islope)).lt. &
477	abs(def_slope_mean(iflat))) THEN
478	iflat = islope
479	ENDIF
480	ENDDO
481
482	PRINT*,'Flat slope for islope = ',iflat
483	PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
484
485
486	allocate(flag_co2flow(ngrid,nslope))
487	allocate(flag_co2flow_mesh(ngrid))
488
489	flag_co2flow(:,:) = 0.
490	flag_co2flow_mesh(:) = 0.
491
492
493	!---------------------------- READ GCM data ---------------------
494
495	! I Initialisation
496	! I_a READ run.def
497	! I_b READ of start_evol.nc and starfi_evol.nc
498	! I_c Subslope parametrisation
499	! I_d READ GCM data and convert to the physical grid
500
501	!------------------------
502
503	! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value
504
505	call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
506
507	allocate(min_h2o_ice_s_1(iim+1,jjm+1))
508	allocate(min_co2_ice_s_1(iim+1,jjm+1))
509	allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
510	allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
511	allocate(q_co2_GCM(iim+1,jjm+1,timelen))
512	allocate(ps_GCM(iim+1,jjm+1,timelen))
513	allocate(ps_GCM_yr1(iim+1,jjm+1,timelen))
514	allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope))
515	allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope))
516	allocate(tsurf_ave(iim+1,jjm+1,nslope))
517	allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
518	allocate(tsurf_ave_phys(ngrid,nslope))
519	allocate(tsurf_ave_phys_yr1(ngrid,nslope))
520	allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
521	allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
522	allocate(co2_ice_GCM_phys_slope(ngrid,nslope,timelen))
523	allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
524	allocate(Tsurfave_before_saved(ngrid,nslope))
525	allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
526	allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
527	allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
528	allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
529	allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
530	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
531	allocate(delta_co2_adsorbded(ngrid))
532	allocate(delta_h2o_adsorbded(ngrid))
533	allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen))
534	allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
535	allocate(watersurf_density_phys_timeseries(ngrid,nslope,timelen))
536	allocate(watersurf_density_phys_ave(ngrid,nslope))
537	allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
538	allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope))
539	print *, "Downloading data Y1..."
540
541	call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm, min_h2o_ice_s_1,min_co2_ice_s_1,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,&
542	nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &
543	watersurf_density_timeseries,watersoil_density_timeseries)
544
545	! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
546
547	print *, "Downloading data Y1 done"
548
549	allocate(min_h2o_ice_s_2(iim+1,jjm+1))
550	allocate(min_co2_ice_s_2(iim+1,jjm+1))
551	allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope))
552	allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
553
554	print *, "Downloading data Y2"
555
556	call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm ,min_h2o_ice_s_2,min_co2_ice_s_2,vmr_co2_gcm,ps_GCM,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
557	nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &
558	watersurf_density_timeseries,watersoil_density_timeseries)
559
560	print *, "Downloading data Y2 done"
561
562	! The variables in the dynamic grid are transfered to the physical grid
563
564	allocate(vmr_co2_gcm_phys(ngrid,timelen))
565	allocate(vmr_co2_pem_phys(ngrid,timelen))
566	allocate(q_h2o_GCM_phys(ngrid,timelen))
567	allocate(q_h2o_PEM_phys(ngrid,timelen))
568	allocate(q_co2_GCM_phys(ngrid,timelen))
569	allocate(q_co2_PEM_phys(ngrid,timelen))
570	allocate(ps_phys(ngrid))
571	allocate(ps_phys_timeseries(ngrid,timelen))
572	allocate(ps_phys_timeseries_yr1(ngrid,timelen))
573
574	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
575	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
576	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
577	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
578	call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
579	call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
580	CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
581	CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
582
583	deallocate(vmr_co2_gcm)
584	deallocate(q_h2o_GCM)
585	deallocate(q_co2_GCM)
586	deallocate(ps_GCM)
587	deallocate(ps_GCM_yr1)
588	deallocate(tsurf_ave)
589	deallocate(tsurf_ave_yr1)
590
591	q_co2_PEM_phys(:,:)= q_co2_GCM_phys(:,:)
592	q_h2o_PEM_phys(:,:)= q_h2o_GCM_phys(:,:)
593
594	!------------------------
595
596	! I Initialisation
597	! I_a READ run.def
598	! I_b READ of start_evol.nc and starfi_evol.nc
599	! I_c Subslope parametrisation
600	! I_d READ GCM data and convert to the physical grid
601	! I_e Initialisation of the PEM variable and soil
602
603	!------------------------
604
605	!---------------------------- Initialisation of the PEM soil and values ---------------------
606
607	call end_comsoil_h_PEM
608	call ini_comsoil_h_PEM(ngrid,nslope)
609
610	allocate(ice_depth(ngrid,nslope))
611	ice_depth(:,:) = 0.
612	allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
613
614	DO islope = 1,nslope
615	if(soil_pem) then
616	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
617	endif !soil_pem
618	DO t=1,timelen
619	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
620	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
621	enddo
622	ENDDO
623
624	deallocate(co2_ice_GCM_slope)
625	deallocate(TI_GCM)
626	deallocate(tsurf_GCM_timeseries)
627
628
629	if(soil_pem) then
630	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
631	DO islope = 1,nslope
632	DO t=1,timelen
633	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersurf_density_timeseries(:,:,islope,t),watersurf_density_phys_timeseries(:,islope,t))
634	ENDDO
635	DO l=1,nsoilmx
636	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
637	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
638	DO t=1,timelen
639	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
640	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_phys_PEM_timeseries(:,l,islope,t))
641	ENDDO
642
643	ENDDO
644	DO l=nsoilmx+1,nsoilmx_PEM
645	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
646	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
647	DO t=1,timelen
648	watersoil_density_phys_PEM_timeseries(:,l,islope,t) = watersoil_density_phys_PEM_timeseries(:,nsoilmx,islope,t)
649	ENDDO
650	ENDDO
651	ENDDO
652	watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
653	watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen
654	deallocate(watersurf_density_timeseries)
655	deallocate(watersurf_density_phys_timeseries)
656	deallocate(watersoil_density_timeseries)
657	deallocate(tsoil_ave_yr1)
658	deallocate(tsoil_ave)
659	deallocate(tsoil_GCM_timeseries)
660	endif !soil_pem
661
662	!------------------------
663
664	! I Initialisation
665	! I_a READ run.def
666	! I_b READ of start_evol.nc and starfi_evol.nc
667	! I_c Subslope parametrisation
668	! I_d READ GCM data and convert to the physical grid
669	! I_e Initialisation of the PEM variable and soil
670	! I_f Compute tendencies & Save initial situation
671
672	!----- Compute tendencies from the PCM run
673
674	allocate(tendencies_h2o_ice(iim+1,jjm+1))
675	allocate(tendencies_h2o_ice_phys(ngrid))
676	allocate(tendencies_co2_ice(iim+1,jjm+1))
677	allocate(tendencies_co2_ice_phys(ngrid))
678	allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
679	allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
680	allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
681	allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
682	allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope))
683
684	! Compute the tendencies of the evolution of ice over the years
685
686	call compute_tendencies(tendencies_h2o_ice,min_h2o_ice_s_1,&
687	min_h2o_ice_s_2,iim,jjm,ngrid,tendencies_h2o_ice_phys)
688
689	call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,&
690	min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys)
691
692	call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,&
693	min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
694
695	tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
696
697	call compute_tendencies_slope(tendencies_h2o_ice_slope,min_h2o_ice_slope_1,&
698	min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
699
700	!------------------------
701
702	! I Initialisation
703	! I_a READ run.def
704	! I_b READ of start_evol.nc and starfi_evol.nc
705	! I_c Subslope parametrisation
706	! I_d READ GCM data and convert to the physical grid
707	! I_e Initialisation of the PEM variable and soil
708	! I_f Compute tendencies & Save initial situation
709	! I_g Save initial PCM situation
710
711	!---------------------------- Save initial PCM situation ---------------------
712
713	allocate(initial_h2o_ice(ngrid))
714	allocate(initial_co2_ice(ngrid))
715	allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
716	allocate(initial_co2_ice_slope(ngrid,nslope))
717	allocate(initial_h2o_ice_slope(ngrid,nslope))
718
719	! We save the places where water ice is sublimating
720	! We compute the surface of water ice sublimating
721	ini_surf=0.
722	ini_surf_co2=0.
723	ini_surf_h2o=0.
724	Total_surface=0.
725	do i=1,ngrid
726	Total_surface=Total_surface+cell_area(i)
727	if (tendencies_h2o_ice_phys(i).LT.0) then
728	initial_h2o_ice(i)=1.
729	ini_surf=ini_surf+cell_area(i)
730	else
731	initial_h2o_ice(i)=0.
732	endif
733	do islope=1,nslope
734	if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
735	initial_co2_ice_sublim_slope(i,islope)=1.
736	ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
737	else
738	initial_co2_ice_sublim_slope(i,islope)=0.
739	endif
740	if (co2ice_slope(i,islope).GT.0) then
741	initial_co2_ice_slope(i,islope)=1.
742	else
743	initial_co2_ice_slope(i,islope)=0.
744	endif
745	if (tendencies_h2o_ice_phys_slope(i,islope).LT.0) then
746	initial_h2o_ice_slope(i,islope)=1.
747	ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
748	else
749	initial_h2o_ice_slope(i,islope)=0.
750	endif
751	enddo
752	enddo
753
754	print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
755	print *, "Total initial surface of h2o ice sublimating=", ini_surf
756	print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
757	print *, "Total surface of the planet=", Total_surface
758
759	allocate(zplev_gcm(ngrid,nlayer+1))
760
761	DO l=1,nlayer+1
762	DO ig=1,ngrid
763	zplev_gcm(ig,l) = ap(l) + bp(l)*ps_phys(ig)
764	ENDDO
765	ENDDO
766
767	global_ave_press_old=0.
768	do i=1,ngrid
769	global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
770	enddo
771
772	global_ave_press_GCM=global_ave_press_old
773	global_ave_press_new=global_ave_press_old
774	print *, "Initial global average pressure=", global_ave_press_GCM
775
776	!------------------------
777
778	! I Initialisation
779	! I_a READ run.def
780	! I_b READ of start_evol.nc and starfi_evol.nc
781	! I_c Subslope parametrisation
782	! I_d READ GCM data and convert to the physical grid
783	! I_e Initialisation of the PEM variable and soil
784	! I_f Compute tendencies & Save initial situation
785	! I_g Save initial PCM situation
786	! I_h Read the PEMstart
787
788	!---------------------------- Read the PEMstart ---------------------
789
790	call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
791	tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsoil_phys_PEM_timeseries,&
792	tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
793	watersurf_density_phys_ave,watersoil_density_phys_PEM_ave, &
794	co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded)
795
796	if(soil_pem) then
797	totmassco2_adsorbded = 0.
798	totmassh2o_adsorbded = 0.
799	do ig = 1,ngrid
800	do islope =1, nslope
801	do l = 1,nsoilmx_PEM - 1
802	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
803	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
804	cell_area(ig)
805	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
806	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
807	cell_area(ig)
808	enddo
809	enddo
810	enddo
811
812	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
813	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
814	deallocate(tsoil_ave_phys_yr1)
815	endif !soil_pem
816	deallocate(tsurf_ave_phys_yr1)
817	deallocate(ps_phys_timeseries_yr1)
818
819	!------------------------
820
821	! I Initialisation
822	! I_a READ run.def
823	! I_b READ of start_evol.nc and starfi_evol.nc
824	! I_c Subslope parametrisation
825	! I_d READ GCM data and convert to the physical grid
826	! I_e Initialisation of the PEM variable and soil
827	! I_f Compute tendencies & Save initial situation
828	! I_g Save initial PCM situation
829	! I_h Read the PEMstar
830	! I_i Compute orbit criterion
831	#ifndef CPP_STD
832	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
833	#else
834	CALL iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
835	#endif
836
837	if(evol_orbit_pem) then
838	call orbit_param_criterion(year_iter_max)
839	else
840	year_iter_max=Max_iter_pem
841	endif
842
843	!--------------------------- END INITIALISATION ---------------------
844
845	!---------------------------- RUN ---------------------
846
847	!------------------------
848
849	! II Run
850	! II_a update pressure,ice and tracers
851
852	!------------------------
853	year_iter=0
854	do while (year_iter.LT.year_iter_max)
855
856	! II.a.1. Compute updated global pressure
857	print *, "Recomputing the new pressure..."
858	do i=1,ngrid
859	do islope=1,nslope
860	global_ave_press_new=global_ave_press_new-gcell_area(i)tendencies_co2_ice_phys_slope(i,islope)subslope_dist(i,islope)/cos(pidef_slope_mean(islope)/180.)/Total_surface
861	enddo
862	enddo
863	print *, 'Global average pressure old time step',global_ave_press_old
864	print *, 'Global average pressure new time step',global_ave_press_new
865
866	if(soil_pem) then
867	do i=1,ngrid
868	global_ave_press_new = global_ave_press_new -gcell_area(i)delta_co2_adsorbded(i)/Total_surface
869	enddo
870	endif
871	print *, 'Global average pressure old time step',global_ave_press_old
872	print *, 'Global average pressure new time step',global_ave_press_new
873
874	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
875	allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
876	print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
877	DO l=1,nlayer+1
878	DO ig=1,ngrid
879	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:)
880	ENDDO
881	ENDDO
882
883	! II.a.3. Surface pressure timeseries
884	print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
885	do i = 1,ngrid
886	ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
887	enddo
888
889	! II.a.4. New pressure levels timeseries
890	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
891	print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
892	do l=1,nlayer+1
893	do ig=1,ngrid
894	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:)
895	enddo
896	enddo
897
898	! II.a.5. Tracers timeseries
899	print *, "Recomputing of tracer VMR timeseries for the new pressure..."
900
901	l=1
902	DO ig=1,ngrid
903	DO t = 1, timelen
904	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))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
905	if(q_h2o_PEM_phys(ig,t).LT.0) then
906	q_h2o_PEM_phys(ig,t)=1E-30
907	endif
908	if(q_h2o_PEM_phys(ig,t).GT.1) then
909	q_h2o_PEM_phys(ig,t)=1.
910	endif
911	enddo
912	enddo
913
914	DO ig=1,ngrid
915	DO t = 1, timelen
916	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))/(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) &
917	+ ( (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) - &
918	(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t)) ) / &
919	(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
920	if (q_co2_PEM_phys(ig,t).LT.0) then
921	q_co2_PEM_phys(ig,t)=1E-30
922	elseif (q_co2_PEM_phys(ig,t).GT.1) then
923	q_co2_PEM_phys(ig,t)=1.
924	endif
925	mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
926	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
927	ENDDO
928	ENDDO
929
930	deallocate(zplev_new_timeseries)
931	deallocate(zplev_old_timeseries)
932
933	! II Run
934	! II_a update pressure, ice and tracers
935	! II_b Evolution of the ice
936
937	! II.b. Evolution of the ice
938	print *, "Evolution of h2o ice"
939	call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
940
941
942	print *, "Evolution of co2 ice"
943	call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
944
945	!------------------------
946
947	! II Run
948	! II_a update pressure, ice and tracers
949	! II_b Evolution of the ice
950	! II_c CO2 glaciers flows
951
952	!------------------------
953
954	print *, "Co2 glacier flows"
955
956
957
958	call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_phys_timeseries,&
959	global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
960
961
962
963
964
965	!------------------------
966
967	! II Run
968	! II_a update pressure, ice and tracers
969	! II_b Evolution of the ice
970	! II_c CO2 glaciers flows
971	! II_d Update surface and soil temperatures
972
973	!------------------------
974
975	! II_d.1 Update Tsurf
976
977	print *, "Updating the new Tsurf"
978	bool_sublim=0
979	Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
980	DO ig = 1,ngrid
981	DO islope = 1,nslope
982	if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice
983	if(latitude_deg(ig).gt.0) then
984	do ig_loop=ig,ngrid
985	DO islope_loop=islope,iflat,-1
986	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
987	tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
988	bool_sublim=1
989	exit
990	endif
991	enddo
992	if (bool_sublim.eq.1) then
993	exit
994	endif
995	enddo
996	else
997	do ig_loop=ig,1,-1
998	DO islope_loop=islope,iflat
999	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
1000	tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
1001	bool_sublim=1
1002	exit
1003	endif
1004	enddo
1005	if (bool_sublim.eq.1) then
1006	exit
1007	endif
1008	enddo
1009	endif
1010	initial_co2_ice_slope(ig,islope)=0
1011	if ((co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then
1012	albedo_slope(ig,1,islope) = albedo_h2o_frost
1013	albedo_slope(ig,2,islope) = albedo_h2o_frost
1014	else
1015	albedo_slope(ig,1,islope) = albedodat(ig)
1016	albedo_slope(ig,2,islope) = albedodat(ig)
1017	emiss_slope(ig,islope) = emissiv
1018	endif
1019	elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_slope(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
1020	ave=0.
1021	do t=1,timelen
1022	if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then
1023	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
1024	else
1025	ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
1026	endif
1027	enddo
1028	tsurf_ave_phys(ig,islope)=ave/timelen
1029	endif
1030	enddo
1031	enddo
1032
1033	do t = 1,timelen
1034	tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:))
1035	enddo
1036	! for the start
1037	do ig = 1,ngrid
1038	do islope = 1,nslope
1039	tsurf_slope(ig,islope) = tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
1040	enddo
1041	enddo
1042
1043
1044	if(soil_pem) then
1045
1046	! II_d.2 Update soil temperature
1047
1048	allocate(TI_locslope(ngrid,nsoilmx_PEM))
1049	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
1050	allocate(Tsurf_locslope(ngrid))
1051	allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
1052	allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
1053
1054	print *,"Updating soil temperature"
1055
1056	! Soil averaged
1057	do islope = 1,nslope
1058	TI_locslope(:,:) = TI_PEM(:,:,islope)
1059	alph_locslope(:,:) = alph_PEM(:,:,islope)
1060	beta_locslope(:,:) = beta_PEM(:,:,islope)
1061	do t = 1,timelen
1062	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
1063	Tsurf_locslope(:) = tsurf_phys_GCM_timeseries(:,islope,t)
1064
1065	call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
1066	call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
1067
1068
1069	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
1070	do ig = 1,ngrid
1071	do isoil = 1,nsoilmx_PEM
1072	watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
1073	if(isnan(Tsoil_locslope(ig,isoil))) then
1074	write(,)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil)
1075	stop
1076	endif
1077	enddo
1078	enddo
1079	enddo
1080	alph_PEM(:,:,islope) = alph_locslope(:,:)
1081	beta_PEM(:,:,islope) = beta_locslope(:,:)
1082	enddo
1083
1084
1085	tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
1086	watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
1087	print *, "Update of soil temperature done"
1088
1089	deallocate(TI_locslope)
1090	deallocate(Tsoil_locslope)
1091	deallocate(Tsurf_locslope)
1092	deallocate(alph_locslope)
1093	deallocate(beta_locslope)
1094
1095	write(,) "Compute ice table"
1096
1097	! II_d.3 Update the ice table
1098	call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
1099
1100	print *, "Update soil propreties"
1101	! II_d.4 Update the soil thermal properties
1102	call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
1103	ice_depth,TI_PEM)
1104
1105	! II_d.5 Update the mass of the regolith adsorbded
1106
1107	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, &
1108	tsoil_PEM,TI_PEM,ps_phys_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
1109	h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
1110
1111	endif !soil_pem
1112
1113	!------------------------
1114
1115	! II Run
1116	! II_a update pressure, ice and tracers
1117	! II_b Evolution of the ice
1118	! II_c CO2 glaciers flows
1119	! II_d Update surface and soil temperatures
1120	! II_e Update the tendencies
1121
1122	!------------------------
1123
1124	print *, "Adaptation of the new co2 tendencies to the current pressure"
1125	call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&
1126	global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
1127
1128	!------------------------
1129
1130	! II Run
1131	! II_a update pressure, ice and tracers
1132	! II_b Evolution of the ice
1133	! II_c CO2 glaciers flows
1134	! II_d Update surface and soil temperatures
1135	! II_e Update the tendencies
1136	! II_f Checking the stopping criterion
1137
1138	!------------------------
1139	call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
1140
1141	call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
1142
1143	year_iter=year_iter+dt_pem
1144
1145	print *, "Checking all the stopping criterion."
1146	if (STOPPING_water) then
1147	print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
1148	endif
1149	if (year_iter.ge.year_iter_max) then
1150	print *, "STOPPING because maximum number of iterations reached"
1151	endif
1152	if (STOPPING_1_water) then
1153	print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
1154	endif
1155	if (STOPPING_co2) then
1156	print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2
1157	endif
1158	if (STOPPING_1_co2) then
1159	print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
1160	endif
1161
1162
1163
1164	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2) then
1165	exit
1166	else
1167	print *, "We continue!"
1168	print *, "Number of iteration of the PEM : year_iter=", year_iter
1169	endif
1170
1171	global_ave_press_old=global_ave_press_new
1172
1173	enddo ! big time iteration loop of the pem
1174
1175
1176	!---------------------------- END RUN PEM ---------------------
1177
1178	!---------------------------- OUTPUT ---------------------
1179
1180	!------------------------
1181
1182	! III Output
1183	! III_a Update surface value for the PCM start files
1184
1185	!------------------------
1186
1187	! III_a.1 Ice update (for startfi)
1188	! Co2 ice
1189	DO ig = 1,ngrid
1190	co2ice(ig) = 0.
1191	DO islope = 1,nslope
1192	co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
1193	* subslope_dist(ig,islope) / &
1194	cos(pi*def_slope_mean(islope)/180.)
1195	ENDDO
1196	#ifdef CPP_STD
1197	qsurf(ig,igcm_co2_ice)=co2ice(ig)
1198	#endif
1199	ENDDO ! of DO ig=1,ngrid
1200	! H2o ice
1201	DO ig = 1,ngrid
1202	qsurf(ig,igcm_h2o_ice) = 0.
1203	DO islope = 1,nslope
1204	qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
1205	* subslope_dist(ig,islope) / &
1206	cos(pi*def_slope_mean(islope)/180.)
1207	ENDDO
1208	ENDDO ! of DO ig=1,ngrid
1209
1210	! III_a.2 Tsoil update (for startfi)
1211
1212	if(soil_pem) then
1213	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
1214	tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
1215	else
1216	TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:)
1217	endif !soil_pem
1218
1219
1220	#ifndef CPP_STD
1221	DO ig = 1,ngrid
1222	DO iloop = 1,nsoilmx
1223	tsoil(ig,iloop) = 0.
1224	inertiesoil(ig,iloop) = 0.
1225	DO islope = 1,nslope
1226	tsoil(ig,iloop) = tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
1227	* subslope_dist(ig,islope)
1228	inertiesoil(ig,iloop) = inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
1229	* subslope_dist(ig,islope)
1230	ENDDO
1231	ENDDO
1232	ENDDO ! of DO ig=1,ngrid
1233
1234	! III_a.3 Surface optical properties (for startfi)
1235
1236	DO ig = 1,ngrid
1237	DO l = 1,2
1238	albedo(ig,l) =0.
1239	DO islope = 1,nslope
1240	albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) &
1241	*subslope_dist(ig,islope)
1242	ENDDO
1243	ENDDO
1244	ENDDO
1245
1246	DO ig = 1,ngrid
1247	emis(ig) =0.
1248	DO islope = 1,nslope
1249	emis(ig)= emis(ig)+emiss_slope(ig,islope) &
1250	*subslope_dist(ig,islope)
1251	ENDDO
1252	ENDDO
1253
1254
1255	DO ig = 1,ngrid
1256	tsurf(ig) = 0.
1257	DO islope = 1,nslope
1258	tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)tsurf_slope(ig,islope))*4 &
1259	* subslope_dist(ig,islope)
1260	ENDDO
1261	tsurf(ig) = tsurf(ig)**(0.25)/emis(ig)
1262	ENDDO ! of DO ig=1,ngrid
1263
1264	#endif
1265
1266	! III_a.4 Pressure (for start)
1267	do i=1,ip1jmp1
1268	ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
1269	enddo
1270
1271	do i = 1,ngrid
1272	ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_GCM
1273	enddo
1274
1275	! III_a.5 Tracer (for start)
1276	allocate(zplev_new(ngrid,nlayer+1))
1277
1278	do l=1,nlayer+1
1279	do ig=1,ngrid
1280	zplev_new(ig,l) = ap(l) + bp(l)*ps_phys(ig)
1281	enddo
1282	enddo
1283
1284	DO nnq=1,nqtot
1285	if (noms(nnq).NE."co2") then
1286	DO l=1,llm-1
1287	DO ig=1,ngrid
1288	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))
1289	ENDDO
1290	q(:,llm,nnq)=q(:,llm-1,nnq)
1291	ENDDO
1292	else
1293	DO l=1,llm-1
1294	DO ig=1,ngrid
1295	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)) &
1296	+ ( (zplev_new(ig,l)-zplev_new(ig,l+1)) - &
1297	(zplev_gcm(ig,l)-zplev_gcm(ig,l+1)) ) / &
1298	(zplev_new(ig,l)-zplev_new(ig,l+1))
1299	ENDDO
1300	q(:,llm,nnq)=q(:,llm-1,nnq)
1301	ENDDO
1302	endif
1303	ENDDO
1304
1305	! Conserving the tracers's mass for GCM start files
1306	DO nnq=1,nqtot
1307	DO ig=1,ngrid
1308	DO l=1,llm-1
1309	if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") ) then
1310	extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
1311	q(ig,l,nnq)=1.
1312	q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
1313	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number"
1314	endif
1315	if(q(ig,l,nnq).LT.0) then
1316	q(ig,l,nnq)=1E-30
1317	endif
1318	ENDDO
1319	enddo
1320	enddo
1321
1322	!------------------------
1323	if(evol_orbit_pem) then
1324	call recomp_orb_param(year_iter)
1325	endif
1326
1327	! III Output
1328	! III_a Update surface value for the PCM start files
1329	! III_b Write start and starfi.nc
1330
1331	!------------------------
1332
1333	! III_b.1 WRITE restart.nc
1334
1335	ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
1336	pday=day_ini
1337	ztime_fin=0.
1338
1339	allocate(p(ip1jmp1,nlayer+1))
1340	CALL pression (ip1jmp1,ap,bp,ps,p)
1341	CALL massdair(p,masse)
1342
1343	CALL dynredem0("restart_evol.nc", day_ini, phis)
1344
1345	CALL dynredem1("restart_evol.nc", &
1346	time_0,vcov,ucov,teta,q,masse,ps)
1347	print *, "restart_evol.nc has been written"
1348
1349	! III_b.2 WRITE restartfi.nc
1350	#ifndef CPP_STD
1351	call physdem0("restartfi_evol.nc",longitude,latitude, &
1352	nsoilmx,ngrid,nlayer,nq, &
1353	ptimestep,pday,0.,cell_area, &
1354	albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
1355	hmons,summit,base,nslope,def_slope, &
1356	subslope_dist)
1357
1358	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
1359	ptimestep,ztime_fin, &
1360	tsurf,tsoil,co2ice,albedo,emis, &
1361	q2,qsurf,tauscaling,totcloudfrac,wstar, &
1362	watercap,inertiesoil,nslope,co2ice_slope, &
1363	tsurf_slope,tsoil_slope, albedo_slope, &
1364	emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
1365	#else
1366	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
1367	ptimestep,pday,time_phys,cell_area, &
1368	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1369
1370	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
1371	ptimestep,ztime_fin, &
1372	tsurf,tsoil,emis,q2,qsurf, &
1373	cloudfrac,totcloudfrac,hice, &
1374	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
1375	#endif
1376
1377	print *, "restartfi_evol.nc has been written"
1378	!------------------------
1379
1380	! III Output
1381	! III_a Update surface value for the PCM start files
1382	! III_b Write start and starfi.nc
1383	! III_c Write start_pem
1384
1385	!------------------------
1386
1387	call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
1388	tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys)
1389
1390	print *, "restartfi_PEM.nc has been written"
1391	print *, "The PEM had run for ", year_iter, " years."
1392	print *, "LL & RV : So far so good"
1393
1394	END PROGRAM pem

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

Download in other formats:

Original Format