Context Navigation

pem.F90 @ 2909

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

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

Download in other formats:

Original Format