Context Navigation

pem.F90 @ 2857

Last change on this file since 2857 was 2857, checked in by llange, 2 years ago
PEM Fixing bug when updating tracer mass mixing ratio (ccn_number and dust_number were previously set to 1) Tiny edits in evol_h2o_ice and criterion_ice_stop_mod_water LL
File size: 58.3 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
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 ! geothermal flux, mass of co2 in the regolith
93	use adsorption_mod, only : regolith_co2adsorption
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 :: totmass_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
266	REAL :: alpha_clap_h2o = -6143.7 ! coeffcient to compute psat, from Murphie et Kood 2005 [K]
267	REAL :: beta_clap_h2o = 28.9074 ! coefficient to compute psat, from Murphie et Kood 2005 [1]
268	LOGICAL :: bool_sublim ! logical to check if there is sublimation or not
269
270	!! Some parameters for the PEM run
271	REAL, PARAMETER :: year_step = 1 ! timestep for the pem
272	INTEGER :: year_iter ! number of iteration
273	INTEGER :: year_iter_max ! maximum number of iterations before stopping
274	REAL :: timestep ! timestep [s]
275	#ifdef CPP_STD
276	! INTEGER :: nsplit_phys=1
277	! LOGICAL :: iflag_phys=.true.
278	REAL :: frost_albedo_threshold=0.05 ! frost albedo threeshold to convert fresh frost to old ice
279	REAL :: albedo_h2o_frost ! albedo of h2o frost
280	REAL,ALLOCATABLE :: co2ice(:) ! Physics: co2 ice mesh averaged [kg/m^2]
281	#endif
282
283	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
284
285	! Loop variable
286	INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil
287	#ifndef CPP_STD
288	! Parallel variables
289	is_sequential=.true.
290	is_parallel=.false.
291	is_mpi_root=.true.
292	is_omp_root=.true.
293	is_master=.true.
294	#endif
295
296	day_ini=0 !test
297	time_phys=0. !test
298
299	! Some constants
300
301	ngrid=ngridmx
302	nlayer=llm
303
304	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
305	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
306	A =(1/m_co2 - 1/m_noco2)
307	B=1/m_noco2
308
309	year_day=669
310	daysec=88775.
311	dtphys=0
312	timestep=year_day*daysec/year_step
313
314	!------------------------
315
316	! I Initialisation
317	! I_a READ run.def
318
319	!------------------------
320
321	!----------------------------READ run.def ---------------------
322	CALL conf_gcm( 99, .TRUE. )
323
324	call infotrac_init
325	nq=nqtot
326
327	!------------------------
328
329	! I Initialisation
330	! I_a READ run.def
331	! I_b READ of start_evol.nc and starfi_evol.nc
332
333	!------------------------
334
335	!----------------------------Initialisation : READ some constant of startfi_evol.nc ---------------------
336
337	! In the gcm, these values are given to the physic by the dynamic.
338	! Here we simply read them in the startfi_evol.nc file
339	status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
340
341	allocate(longitude(ngrid))
342	allocate(latitude(ngrid))
343	allocate(cell_area(ngrid))
344
345	status = nf90_inq_varid(ncid, "longitude", lonvarid)
346	status = nf90_get_var(ncid, lonvarid, longitude)
347
348	status = nf90_inq_varid(ncid, "latitude", latvarid)
349	status = nf90_get_var(ncid, latvarid, latitude)
350
351	status = nf90_inq_varid(ncid, "area", areavarid)
352	status = nf90_get_var(ncid, areavarid, cell_area)
353
354	call ini_comsoil_h(ngrid)
355
356	status = nf90_inq_varid(ncid, "soildepth", sdvarid)
357	status = nf90_get_var(ncid, sdvarid, mlayer)
358
359	status =nf90_close(ncid)
360
361	!----------------------------READ start.nc ---------------------
362
363	allocate(q(ip1jmp1,llm,nqtot))
364	CALL dynetat0(FILE_NAME_start,vcov,ucov, &
365	teta,q,masse,ps,phis, time_0)
366
367	CALL iniconst !new
368	CALL inigeom
369	allocate(ap(nlayer+1))
370	allocate(bp(nlayer+1))
371	status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
372	status = nf90_inq_varid(ncid, "ap", apvarid)
373	status = nf90_get_var(ncid, apvarid, ap)
374	status = nf90_inq_varid(ncid, "bp", bpvarid)
375	status = nf90_get_var(ncid, bpvarid, bp)
376	status =nf90_close(ncid)
377
378	CALL iniphysiq(iim,jjm,llm, &
379	(jjm-1)*iim+2,comm_lmdz, &
380	daysec,day_ini,dtphys/nsplit_phys, &
381	rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
382	iflag_phys)
383
384	!----------------------------READ startfi.nc ---------------------
385
386	! First we read the initial state (starfi.nc)
387
388	allocate(watercap_slope(ngrid,nslope))
389	allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
390	allocate(inertiesoil(ngrid,nsoilmx))
391
392	#ifndef CPP_STD
393	CALL phyetat0 (FILE_NAME,0,0, &
394	nsoilmx,ngrid,nlayer,nq, &
395	day_ini,time_phys, &
396	tsurf,tsoil,albedo,emis, &
397	q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar, &
398	watercap,inertiesoil,nslope,tsurf_slope, &
399	tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
400	subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM_start, &
401	qsurf_slope,watercap_slope)
402
403	if(soil_pem) then
404	deallocate(TI_GCM_start) !not used then
405	endif
406
407	! Remove unphysical values of surface tracer
408	DO i=1,ngrid
409	DO nnq=1,nqtot
410	DO islope=1,nslope
411	if(qsurf_slope(i,nnq,islope).LT.0) then
412	qsurf_slope(i,nnq,islope)=0.
413	endif
414	enddo
415	enddo
416	enddo
417
418	#else
419	call phys_state_var_init(nq)
420	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
421	call initracer(ngrid,nq)
422	call iniaerosol()
423	call phyetat0(.true., &
424	ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq, &
425	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
426	cloudfrac,totcloudfrac,hice, &
427	rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
428	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
429
430	call ini_comslope_h(ngrid,nsoilmx,nq)
431
432	allocate(co2ice(ngrid))
433	co2ice(:)=qsurf(:,igcm_co2_ice)
434	co2ice_slope(:,1)=co2ice(:)
435	tsurf_slope(:,1)=tsurf(:)
436
437	if (nslope.eq.1) then
438	def_slope(1) = 0
439	def_slope(2) = 0
440	def_slope_mean=0
441	subslope_dist(:,1) = 1.
442	endif
443
444	! Remove unphysical values of surface tracer
445	DO i=1,ngrid
446	DO nnq=1,nqtot
447	qsurf_slope(i,nnq,1)=qsurf(i,nnq)
448	if(qsurf(i,nnq).LT.0) then
449	qsurf(i,nnq)=0.
450	endif
451	enddo
452	enddo
453	#endif
454
455	DO nnq=1,nqtot
456	if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
457	ENDDO
458
459	!------------------------
460
461	! I Initialisation
462	! I_a READ run.def
463	! I_b READ of start_evol.nc and starfi_evol.nc
464	! I_c Subslope parametrisation
465
466	!------------------------
467
468	!----------------------------Subslope parametrisation definition ---------------------
469
470	! Define some slope statistics
471	iflat=1
472	DO islope=2,nslope
473	IF(abs(def_slope_mean(islope)).lt. &
474	abs(def_slope_mean(iflat))) THEN
475	iflat = islope
476	ENDIF
477	ENDDO
478
479	PRINT*,'Flat slope for islope = ',iflat
480	PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
481
482
483	allocate(flag_co2flow(ngrid,nslope))
484	allocate(flag_co2flow_mesh(ngrid))
485
486	flag_co2flow(:,:) = 0.
487	flag_co2flow_mesh(:) = 0.
488
489
490	!---------------------------- READ GCM data ---------------------
491
492	! I Initialisation
493	! I_a READ run.def
494	! I_b READ of start_evol.nc and starfi_evol.nc
495	! I_c Subslope parametrisation
496	! I_d READ GCM data and convert to the physical grid
497
498	!------------------------
499
500	! 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
501
502	call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
503
504	allocate(min_h2o_ice_s_1(iim+1,jjm+1))
505	allocate(min_co2_ice_s_1(iim+1,jjm+1))
506	allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
507	allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
508	allocate(q_co2_GCM(iim+1,jjm+1,timelen))
509	allocate(ps_GCM(iim+1,jjm+1,timelen))
510	allocate(ps_GCM_yr1(iim+1,jjm+1,timelen))
511	allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope))
512	allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope))
513	allocate(tsurf_ave(iim+1,jjm+1,nslope))
514	allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
515	allocate(tsurf_ave_phys(ngrid,nslope))
516	allocate(tsurf_ave_phys_yr1(ngrid,nslope))
517	allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
518	allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
519	allocate(co2_ice_GCM_phys_slope(ngrid,nslope,timelen))
520	allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
521	allocate(Tsurfave_before_saved(ngrid,nslope))
522	allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
523	allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
524	allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
525	allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
526	allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
527	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
528	allocate(delta_co2_adsorbded(ngrid))
529	allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen))
530	allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
531	allocate(watersurf_density_phys_timeseries(ngrid,nslope,timelen))
532	allocate(watersurf_density_phys_ave(ngrid,nslope))
533	allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
534	allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope))
535	print *, "Downloading data Y1..."
536
537	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,&
538	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, &
539	watersurf_density_timeseries,watersoil_density_timeseries)
540
541	! 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
542
543	print *, "Downloading data Y1 done"
544
545	allocate(min_h2o_ice_s_2(iim+1,jjm+1))
546	allocate(min_co2_ice_s_2(iim+1,jjm+1))
547	allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope))
548	allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
549
550	print *, "Downloading data Y2"
551
552	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, &
553	nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &
554	watersurf_density_timeseries,watersoil_density_timeseries)
555
556	print *, "Downloading data Y2 done"
557
558	! The variables in the dynamic grid are transfered to the physical grid
559
560	allocate(vmr_co2_gcm_phys(ngrid,timelen))
561	allocate(vmr_co2_pem_phys(ngrid,timelen))
562	allocate(q_h2o_GCM_phys(ngrid,timelen))
563	allocate(q_h2o_PEM_phys(ngrid,timelen))
564	allocate(q_co2_GCM_phys(ngrid,timelen))
565	allocate(q_co2_PEM_phys(ngrid,timelen))
566	allocate(ps_phys(ngrid))
567	allocate(ps_phys_timeseries(ngrid,timelen))
568	allocate(ps_phys_timeseries_yr1(ngrid,timelen))
569
570	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
571	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
572	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
573	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
574	call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
575	call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
576	CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
577	CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
578
579	deallocate(vmr_co2_gcm)
580	deallocate(q_h2o_GCM)
581	deallocate(q_co2_GCM)
582	deallocate(ps_GCM)
583	deallocate(ps_GCM_yr1)
584	deallocate(tsurf_ave)
585	deallocate(tsurf_ave_yr1)
586
587	q_co2_PEM_phys(:,:)= q_co2_GCM_phys(:,:)
588	q_h2o_PEM_phys(:,:)= q_h2o_GCM_phys(:,:)
589
590	!------------------------
591
592	! I Initialisation
593	! I_a READ run.def
594	! I_b READ of start_evol.nc and starfi_evol.nc
595	! I_c Subslope parametrisation
596	! I_d READ GCM data and convert to the physical grid
597	! I_e Initialisation of the PEM variable and soil
598
599	!------------------------
600
601	!---------------------------- Initialisation of the PEM soil and values ---------------------
602
603	call end_comsoil_h_PEM
604	call ini_comsoil_h_PEM(ngrid,nslope)
605
606	allocate(ice_depth(ngrid,nslope))
607	ice_depth(:,:) = 0.
608	allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
609
610	DO islope = 1,nslope
611	if(soil_pem) then
612	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
613	endif !soil_pem
614	DO t=1,timelen
615	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
616	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
617	enddo
618	ENDDO
619
620	deallocate(co2_ice_GCM_slope)
621	deallocate(TI_GCM)
622	deallocate(tsurf_GCM_timeseries)
623
624
625	if(soil_pem) then
626	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
627	DO islope = 1,nslope
628	DO t=1,timelen
629	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersurf_density_timeseries(:,:,islope,t),watersurf_density_phys_timeseries(:,islope,t))
630	ENDDO
631	DO l=1,nsoilmx
632	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
633	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
634	DO t=1,timelen
635	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
636	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_phys_PEM_timeseries(:,l,islope,t))
637	ENDDO
638
639	ENDDO
640	DO l=nsoilmx+1,nsoilmx_PEM
641	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
642	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
643	DO t=1,timelen
644	watersoil_density_phys_PEM_timeseries(:,l,islope,t) = watersoil_density_phys_PEM_timeseries(:,nsoilmx,islope,t)
645	ENDDO
646	ENDDO
647	ENDDO
648	watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
649	watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen
650	deallocate(watersurf_density_timeseries)
651	deallocate(watersurf_density_phys_timeseries)
652	deallocate(watersoil_density_timeseries)
653	deallocate(tsoil_ave_yr1)
654	deallocate(tsoil_ave)
655	deallocate(tsoil_GCM_timeseries)
656	endif !soil_pem
657
658	!------------------------
659
660	! I Initialisation
661	! I_a READ run.def
662	! I_b READ of start_evol.nc and starfi_evol.nc
663	! I_c Subslope parametrisation
664	! I_d READ GCM data and convert to the physical grid
665	! I_e Initialisation of the PEM variable and soil
666	! I_f Compute tendencies & Save initial situation
667
668	!----- Compute tendencies from the PCM run
669
670	allocate(tendencies_h2o_ice(iim+1,jjm+1))
671	allocate(tendencies_h2o_ice_phys(ngrid))
672	allocate(tendencies_co2_ice(iim+1,jjm+1))
673	allocate(tendencies_co2_ice_phys(ngrid))
674	allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
675	allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
676	allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
677	allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
678	allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope))
679
680	! Compute the tendencies of the evolution of ice over the years
681
682	call compute_tendencies(tendencies_h2o_ice,min_h2o_ice_s_1,&
683	min_h2o_ice_s_2,iim,jjm,ngrid,tendencies_h2o_ice_phys)
684
685	call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,&
686	min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys)
687
688	call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,&
689	min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
690
691	tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
692
693	call compute_tendencies_slope(tendencies_h2o_ice_slope,min_h2o_ice_slope_1,&
694	min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
695
696	!------------------------
697
698	! I Initialisation
699	! I_a READ run.def
700	! I_b READ of start_evol.nc and starfi_evol.nc
701	! I_c Subslope parametrisation
702	! I_d READ GCM data and convert to the physical grid
703	! I_e Initialisation of the PEM variable and soil
704	! I_f Compute tendencies & Save initial situation
705	! I_g Save initial PCM situation
706
707	!---------------------------- Save initial PCM situation ---------------------
708
709	allocate(initial_h2o_ice(ngrid))
710	allocate(initial_co2_ice(ngrid))
711	allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
712	allocate(initial_co2_ice_slope(ngrid,nslope))
713	allocate(initial_h2o_ice_slope(ngrid,nslope))
714
715	! We save the places where water ice is sublimating
716	! We compute the surface of water ice sublimating
717	ini_surf=0.
718	ini_surf_co2=0.
719	ini_surf_h2o=0.
720	Total_surface=0.
721	do i=1,ngrid
722	Total_surface=Total_surface+cell_area(i)
723	if (tendencies_h2o_ice_phys(i).LT.0) then
724	initial_h2o_ice(i)=1.
725	ini_surf=ini_surf+cell_area(i)
726	else
727	initial_h2o_ice(i)=0.
728	endif
729	do islope=1,nslope
730	if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
731	initial_co2_ice_sublim_slope(i,islope)=1.
732	ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
733	else
734	initial_co2_ice_sublim_slope(i,islope)=0.
735	endif
736	if (co2ice_slope(i,islope).GT.0) then
737	initial_co2_ice_slope(i,islope)=1.
738	else
739	initial_co2_ice_slope(i,islope)=0.
740	endif
741	if (tendencies_h2o_ice_phys_slope(i,islope).LT.0) then
742	initial_h2o_ice_slope(i,islope)=1.
743	ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
744	else
745	initial_h2o_ice_slope(i,islope)=0.
746	endif
747	enddo
748	enddo
749
750	print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
751	print *, "Total initial surface of h2o ice sublimating=", ini_surf
752	print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
753	print *, "Total surface of the planet=", Total_surface
754
755	allocate(zplev_gcm(ngrid,nlayer+1))
756
757	DO l=1,nlayer+1
758	DO ig=1,ngrid
759	zplev_gcm(ig,l) = ap(l) + bp(l)*ps_phys(ig)
760	ENDDO
761	ENDDO
762
763	global_ave_press_old=0.
764	do i=1,ngrid
765	global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
766	enddo
767
768	global_ave_press_GCM=global_ave_press_old
769	global_ave_press_new=global_ave_press_old
770	print *, "Initial global average pressure=", global_ave_press_GCM
771
772	!------------------------
773
774	! I Initialisation
775	! I_a READ run.def
776	! I_b READ of start_evol.nc and starfi_evol.nc
777	! I_c Subslope parametrisation
778	! I_d READ GCM data and convert to the physical grid
779	! I_e Initialisation of the PEM variable and soil
780	! I_f Compute tendencies & Save initial situation
781	! I_g Save initial PCM situation
782	! I_h Read the PEMstart
783
784	!---------------------------- Read the PEMstart ---------------------
785
786	call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
787	tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsoil_phys_PEM_timeseries,&
788	tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM, co2_adsorbded_phys,delta_co2_adsorbded,&
789	watersurf_density_phys_ave,watersoil_density_phys_PEM_ave)
790
791	if(soil_pem) then
792	totmass_adsorbded = 0.
793
794	do ig = 1,ngrid
795	do islope =1, nslope
796	do l = 1,nsoilmx_PEM - 1
797	totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
798	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
799	cell_area(ig)
800	enddo
801	enddo
802	enddo
803
804	write(,) "Tot mass in the regolith=", totmass_adsorbded
805	deallocate(tsoil_ave_phys_yr1)
806	endif !soil_pem
807	deallocate(tsurf_ave_phys_yr1)
808	deallocate(ps_phys_timeseries_yr1)
809
810	!------------------------
811
812	! I Initialisation
813	! I_a READ run.def
814	! I_b READ of start_evol.nc and starfi_evol.nc
815	! I_c Subslope parametrisation
816	! I_d READ GCM data and convert to the physical grid
817	! I_e Initialisation of the PEM variable and soil
818	! I_f Compute tendencies & Save initial situation
819	! I_g Save initial PCM situation
820	! I_h Read the PEMstar
821	! I_i Compute orbit criterion
822	#ifndef CPP_STD
823	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
824	#else
825	CALL iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
826	#endif
827
828	if(evol_orbit_pem) then
829	call orbit_param_criterion(year_iter_max)
830	else
831	year_iter_max=Max_iter_pem
832	endif
833
834	!--------------------------- END INITIALISATION ---------------------
835
836	!---------------------------- RUN ---------------------
837
838	!------------------------
839
840	! II Run
841	! II_a update pressure,ice and tracers
842
843	!------------------------
844	year_iter=0
845	do while (year_iter.LT.year_iter_max)
846
847	! II.a.1. Compute updated global pressure
848	print *, "Recomputing the new pressure..."
849	do i=1,ngrid
850	do islope=1,nslope
851	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
852	enddo
853	enddo
854	print *, 'Global average pressure old time step',global_ave_press_old
855	print *, 'Global average pressure new time step',global_ave_press_new
856
857	do i=1,ngrid
858	if(soil_pem) then
859	global_ave_press_new = global_ave_press_new -gcell_area(i)delta_co2_adsorbded(i)/Total_surface
860	endif
861	enddo
862	print *, 'Global average pressure old time step',global_ave_press_old
863	print *, 'Global average pressure new time step',global_ave_press_new
864
865	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
866	allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
867	print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
868	DO l=1,nlayer+1
869	DO ig=1,ngrid
870	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:)
871	ENDDO
872	ENDDO
873
874	! II.a.3. Surface pressure timeseries
875	print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
876	do i = 1,ngrid
877	ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
878	enddo
879
880	! II.a.4. New pressure levels timeseries
881	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
882	print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
883	do l=1,nlayer+1
884	do ig=1,ngrid
885	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:)
886	enddo
887	enddo
888
889	! II.a.5. Tracers timeseries
890	print *, "Recomputing of tracer VMR timeseries for the new pressure..."
891
892	l=1
893	DO ig=1,ngrid
894	DO t = 1, timelen
895	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))
896	if(q_h2o_PEM_phys(ig,t).LT.0) then
897	q_h2o_PEM_phys(ig,t)=1E-30
898	endif
899	if(q_h2o_PEM_phys(ig,t).GT.1) then
900	q_h2o_PEM_phys(ig,t)=1.
901	endif
902	enddo
903	enddo
904
905	DO ig=1,ngrid
906	DO t = 1, timelen
907	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)) &
908	+ ( (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) - &
909	(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t)) ) / &
910	(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
911	if (q_co2_PEM_phys(ig,t).LT.0) then
912	q_co2_PEM_phys(ig,t)=1E-30
913	elseif (q_co2_PEM_phys(ig,t).GT.1) then
914	q_co2_PEM_phys(ig,t)=1.
915	endif
916	mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
917	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
918	ENDDO
919	ENDDO
920
921	deallocate(zplev_new_timeseries)
922	deallocate(zplev_old_timeseries)
923
924	! II Run
925	! II_a update pressure, ice and tracers
926	! II_b Evolution of the ice
927
928	! II.b. Evolution of the ice
929	print *, "Evolution of h2o ice"
930	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)
931
932
933	print *, "Evolution of co2 ice"
934	call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
935
936	!------------------------
937
938	! II Run
939	! II_a update pressure, ice and tracers
940	! II_b Evolution of the ice
941	! II_c CO2 glaciers flows
942
943	!------------------------
944
945	print *, "Co2 glacier flows"
946
947
948
949	call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_phys_timeseries,&
950	global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
951
952
953
954
955
956	!------------------------
957
958	! II Run
959	! II_a update pressure, ice and tracers
960	! II_b Evolution of the ice
961	! II_c CO2 glaciers flows
962	! II_d Update surface and soil temperatures
963
964	!------------------------
965
966	! II_d.1 Update Tsurf
967
968	print *, "Updating the new Tsurf"
969	bool_sublim=0
970	Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
971	DO ig = 1,ngrid
972	DO islope = 1,nslope
973	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
974	if(latitude_deg(ig).gt.0) then
975	do ig_loop=ig,ngrid
976	DO islope_loop=islope,iflat,-1
977	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
978	tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
979	bool_sublim=1
980	exit
981	endif
982	enddo
983	if (bool_sublim.eq.1) then
984	exit
985	endif
986	enddo
987	else
988	do ig_loop=ig,1,-1
989	DO islope_loop=islope,iflat
990	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
991	tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
992	bool_sublim=1
993	exit
994	endif
995	enddo
996	if (bool_sublim.eq.1) then
997	exit
998	endif
999	enddo
1000	endif
1001	initial_co2_ice_slope(ig,islope)=0
1002	if ((co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then
1003	albedo_slope(ig,1,islope) = albedo_h2o_frost
1004	albedo_slope(ig,2,islope) = albedo_h2o_frost
1005	else
1006	albedo_slope(ig,1,islope) = albedodat(ig)
1007	albedo_slope(ig,2,islope) = albedodat(ig)
1008	emiss_slope(ig,islope) = emissiv
1009	endif
1010	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
1011	ave=0.
1012	do t=1,timelen
1013	if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then
1014	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
1015	else
1016	ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
1017	endif
1018	enddo
1019	tsurf_ave_phys(ig,islope)=ave/timelen
1020	endif
1021	enddo
1022	enddo
1023
1024	do t = 1,timelen
1025	tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:))
1026	enddo
1027	! for the start
1028	do ig = 1,ngrid
1029	do islope = 1,nslope
1030	tsurf_slope(ig,islope) = tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
1031	enddo
1032	enddo
1033
1034
1035	if(soil_pem) then
1036
1037	! II_d.2 Update soil temperature
1038
1039	allocate(TI_locslope(ngrid,nsoilmx_PEM))
1040	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
1041	allocate(Tsurf_locslope(ngrid))
1042	allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
1043	allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
1044
1045	print *,"Updating soil temperature"
1046
1047	! Soil averaged
1048	do islope = 1,nslope
1049	TI_locslope(:,:) = TI_PEM(:,:,islope)
1050	alph_locslope(:,:) = alph_PEM(:,:,islope)
1051	beta_locslope(:,:) = beta_PEM(:,:,islope)
1052	do t = 1,timelen
1053	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
1054	Tsurf_locslope(:) = tsurf_phys_GCM_timeseries(:,islope,t)
1055
1056	call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
1057	call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
1058
1059
1060	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
1061	do ig = 1,ngrid
1062	do isoil = 1,nsoilmx_PEM
1063	watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
1064	if(isnan(Tsoil_locslope(ig,isoil))) then
1065	write(,)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil)
1066	stop
1067	endif
1068	enddo
1069	enddo
1070	enddo
1071	alph_PEM(:,:,islope) = alph_locslope(:,:)
1072	beta_PEM(:,:,islope) = beta_locslope(:,:)
1073	enddo
1074
1075
1076	tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
1077	watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
1078	print *, "Update of soil temperature done"
1079
1080	deallocate(TI_locslope)
1081	deallocate(Tsoil_locslope)
1082	deallocate(Tsurf_locslope)
1083	deallocate(alph_locslope)
1084	deallocate(beta_locslope)
1085
1086	write(,) "Compute ice table"
1087
1088	! II_d.3 Update the ice table
1089	call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth)
1090
1091	print *, "Update soil propreties"
1092	! II_d.4 Update the soil thermal properties
1093	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, &
1094	ice_depth,TI_PEM)
1095
1096	! II_d.5 Update the mass of the regolith adsorbded
1097	call regolith_co2adsorption(ngrid,nslope,nsoilmx_PEM,timelen,ps_phys_timeseries,tsoil_PEM,TI_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope, &
1098	co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded)
1099
1100	endif !soil_pem
1101
1102	!------------------------
1103
1104	! II Run
1105	! II_a update pressure, ice and tracers
1106	! II_b Evolution of the ice
1107	! II_c CO2 glaciers flows
1108	! II_d Update surface and soil temperatures
1109	! II_e Update the tendencies
1110
1111	!------------------------
1112
1113	print *, "Adaptation of the new co2 tendencies to the current pressure"
1114	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,&
1115	global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
1116
1117	!------------------------
1118
1119	! II Run
1120	! II_a update pressure, ice and tracers
1121	! II_b Evolution of the ice
1122	! II_c CO2 glaciers flows
1123	! II_d Update surface and soil temperatures
1124	! II_e Update the tendencies
1125	! II_f Checking the stopping criterion
1126
1127	!------------------------
1128	call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
1129
1130	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)
1131
1132	year_iter=year_iter+dt_pem
1133
1134	print *, "Checking all the stopping criterion."
1135	if (STOPPING_water) then
1136	print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
1137	endif
1138	if (year_iter.ge.year_iter_max) then
1139	print *, "STOPPING because maximum number of iterations reached"
1140	endif
1141	if (STOPPING_1_water) then
1142	print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
1143	endif
1144	if (STOPPING_co2) then
1145	print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2
1146	endif
1147	if (STOPPING_1_co2) then
1148	print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
1149	endif
1150
1151
1152
1153	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2) then
1154	exit
1155	else
1156	print *, "We continue!"
1157	print *, "Number of iteration of the PEM : year_iter=", year_iter
1158	endif
1159
1160	global_ave_press_old=global_ave_press_new
1161
1162	enddo ! big time iteration loop of the pem
1163
1164
1165	!---------------------------- END RUN PEM ---------------------
1166
1167	!---------------------------- OUTPUT ---------------------
1168
1169	!------------------------
1170
1171	! III Output
1172	! III_a Update surface value for the PCM start files
1173
1174	!------------------------
1175
1176	! III_a.1 Ice update (for startfi)
1177	! Co2 ice
1178	DO ig = 1,ngrid
1179	co2ice(ig) = 0.
1180	DO islope = 1,nslope
1181	co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
1182	* subslope_dist(ig,islope) / &
1183	cos(pi*def_slope_mean(islope)/180.)
1184	ENDDO
1185	#ifdef CPP_STD
1186	qsurf(ig,igcm_co2_ice)=co2ice(ig)
1187	#endif
1188	ENDDO ! of DO ig=1,ngrid
1189	! H2o ice
1190	DO ig = 1,ngrid
1191	qsurf(ig,igcm_h2o_ice) = 0.
1192	DO islope = 1,nslope
1193	qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
1194	* subslope_dist(ig,islope) / &
1195	cos(pi*def_slope_mean(islope)/180.)
1196	ENDDO
1197	ENDDO ! of DO ig=1,ngrid
1198
1199	! III_a.2 Tsoil update (for startfi)
1200
1201	if(soil_pem) then
1202	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
1203	tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
1204	else
1205	TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:)
1206	endif !soil_pem
1207
1208
1209	#ifndef CPP_STD
1210	DO ig = 1,ngrid
1211	DO iloop = 1,nsoilmx
1212	tsoil(ig,iloop) = 0.
1213	inertiesoil(ig,iloop) = 0.
1214	DO islope = 1,nslope
1215	tsoil(ig,iloop) = tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
1216	* subslope_dist(ig,islope)
1217	inertiesoil(ig,iloop) = inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
1218	* subslope_dist(ig,islope)
1219	ENDDO
1220	ENDDO
1221	ENDDO ! of DO ig=1,ngrid
1222
1223	! III_a.3 Surface optical properties (for startfi)
1224
1225	DO ig = 1,ngrid
1226	DO l = 1,2
1227	albedo(ig,l) =0.
1228	DO islope = 1,nslope
1229	albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) &
1230	*subslope_dist(ig,islope)
1231	ENDDO
1232	ENDDO
1233	ENDDO
1234
1235	DO ig = 1,ngrid
1236	emis(ig) =0.
1237	DO islope = 1,nslope
1238	emis(ig)= emis(ig)+emiss_slope(ig,islope) &
1239	*subslope_dist(ig,islope)
1240	ENDDO
1241	ENDDO
1242
1243
1244	DO ig = 1,ngrid
1245	tsurf(ig) = 0.
1246	DO islope = 1,nslope
1247	tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)tsurf_slope(ig,islope))*4 &
1248	* subslope_dist(ig,islope)
1249	ENDDO
1250	tsurf(ig) = tsurf(ig)**(0.25)/emis(ig)
1251	ENDDO ! of DO ig=1,ngrid
1252
1253	#endif
1254
1255	! III_a.4 Pressure (for start)
1256	do i=1,ip1jmp1
1257	ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
1258	enddo
1259
1260	do i = 1,ngrid
1261	ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_GCM
1262	enddo
1263
1264	write(,) 'rapport ps',global_ave_press_new/global_ave_press_GCM
1265
1266	! III_a.5 Tracer (for start)
1267	allocate(zplev_new(ngrid,nlayer+1))
1268
1269	do l=1,nlayer+1
1270	do ig=1,ngrid
1271	zplev_new(ig,l) = ap(l) + bp(l)*ps_phys(ig)
1272	enddo
1273	enddo
1274
1275	DO nnq=1,nqtot
1276	if (noms(nnq).NE."co2") then
1277	DO l=1,llm-1
1278	DO ig=1,ngrid
1279	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))
1280	ENDDO
1281	q(:,llm,nnq)=q(:,llm-1,nnq)
1282	ENDDO
1283	else
1284	DO l=1,llm-1
1285	DO ig=1,ngrid
1286	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)) &
1287	+ ( (zplev_new(ig,l)-zplev_new(ig,l+1)) - &
1288	(zplev_gcm(ig,l)-zplev_gcm(ig,l+1)) ) / &
1289	(zplev_new(ig,l)-zplev_new(ig,l+1))
1290	ENDDO
1291	q(:,llm,nnq)=q(:,llm-1,nnq)
1292	ENDDO
1293	endif
1294	ENDDO
1295
1296	! Conserving the tracers's mass for GCM start files
1297	DO nnq=1,nqtot
1298	DO ig=1,ngrid
1299	DO l=1,llm-1
1300	if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") ) then
1301	extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
1302	q(ig,l,nnq)=1.
1303	q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
1304	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number"
1305	endif
1306	if(q(ig,l,nnq).LT.0) then
1307	q(ig,l,nnq)=1E-30
1308	endif
1309	ENDDO
1310	enddo
1311	enddo
1312
1313	!------------------------
1314	if(evol_orbit_pem) then
1315	call recomp_orb_param(year_iter)
1316	endif
1317
1318	! III Output
1319	! III_a Update surface value for the PCM start files
1320	! III_b Write start and starfi.nc
1321
1322	!------------------------
1323
1324	! III_b.1 WRITE restart.nc
1325
1326	ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
1327	pday=day_ini
1328	ztime_fin=0.
1329
1330	allocate(p(ip1jmp1,nlayer+1))
1331	CALL pression (ip1jmp1,ap,bp,ps,p)
1332	CALL massdair(p,masse)
1333
1334	CALL dynredem0("restart_evol.nc", day_ini, phis)
1335
1336	CALL dynredem1("restart_evol.nc", &
1337	time_0,vcov,ucov,teta,q,masse,ps)
1338	print *, "restart_evol.nc has been written"
1339
1340	! III_b.2 WRITE restartfi.nc
1341	#ifndef CPP_STD
1342	call physdem0("restartfi_evol.nc",longitude,latitude, &
1343	nsoilmx,ngrid,nlayer,nq, &
1344	ptimestep,pday,0.,cell_area, &
1345	albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
1346	hmons,summit,base,nslope,def_slope, &
1347	subslope_dist)
1348
1349	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
1350	ptimestep,ztime_fin, &
1351	tsurf,tsoil,co2ice,albedo,emis, &
1352	q2,qsurf,tauscaling,totcloudfrac,wstar, &
1353	watercap,inertiesoil,nslope,co2ice_slope, &
1354	tsurf_slope,tsoil_slope, albedo_slope, &
1355	emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
1356	#else
1357	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
1358	ptimestep,pday,time_phys,cell_area, &
1359	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1360
1361	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
1362	ptimestep,ztime_fin, &
1363	tsurf,tsoil,emis,q2,qsurf, &
1364	cloudfrac,totcloudfrac,hice, &
1365	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
1366	#endif
1367
1368	print *, "restartfi_evol.nc has been written"
1369	!------------------------
1370
1371	! III Output
1372	! III_a Update surface value for the PCM start files
1373	! III_b Write start and starfi.nc
1374	! III_c Write start_pem
1375
1376	!------------------------
1377
1378	call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
1379	tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys)
1380
1381	print *, "restartfi_PEM.nc has been written"
1382	print *, "The PEM had run for ", year_iter, " years."
1383	print *, "LL & RV : So far so good"
1384
1385	END PROGRAM pem

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

Download in other formats:

Original Format