Context Navigation

pem.F90 @ 2937

Last change on this file since 2937 was 2937, checked in by llange, 23 months ago
PEM Thermal Inertia is now only read in the start and not read in the xios file Fixing few bug (problem of allocation and index for update_soil and some variables in the pem) LL
File size: 55.0 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,co2glaciersflow
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 :: 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(:,:,:) ! Same but for the start
226
227	REAL,ALLOCATABLE :: ice_depth(:,:) ! Physic x SLope: Ice table depth [m]
228	REAL,ALLOCATABLE :: TI_locslope(:,:) ! Physic x Soil: Intermediate thermal inertia to compute Tsoil [SI]
229	REAL,ALLOCATABLE :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
230	REAL,ALLOCATABLE :: Tsurf_locslope(:) ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
231	REAL,ALLOCATABLE :: watersoil_density_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
232	REAL,ALLOCATABLE :: watersurf_density_ave(:,:) ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
233	REAL,ALLOCATABLE :: watersoil_density_PEM_timeseries(:,:,:,:) ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
234	REAL,ALLOCATABLE :: watersoil_density_PEM_ave(:,:,:) ! Physic x Soil x SLopes, water soil density, yearly averaged [kg/m^3]
235	REAL,ALLOCATABLE :: Tsurfave_before_saved(:,:) ! Surface temperature saved from previous time step [K]
236	REAL, ALLOCATABLE :: delta_co2_adsorbded(:) ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
237	REAL, ALLOCATABLE :: delta_h2o_adsorbded(:) ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
238	REAL :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
239	REAL :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
240	REAL :: alpha_clap_h2o = -6143.7 ! coeffcient to compute psat, from Murphie et Kood 2005 [K]
241	REAL :: beta_clap_h2o = 28.9074 ! coefficient to compute psat, from Murphie et Kood 2005 [1]
242	LOGICAL :: bool_sublim ! logical to check if there is sublimation or not
243
244
245	!! Some parameters for the PEM run
246	REAL, PARAMETER :: year_step = 1 ! timestep for the pem
247	INTEGER :: year_iter ! number of iteration
248	INTEGER :: year_iter_max ! maximum number of iterations before stopping
249	REAL :: timestep ! timestep [s]
250	REAL :: watercap_sum ! total mass of water cap [kg/m^2]
251
252	#ifdef CPP_STD
253	! INTEGER :: nsplit_phys=1
254	! LOGICAL :: iflag_phys=.true.
255	REAL :: frost_albedo_threshold=0.05 ! frost albedo threeshold to convert fresh frost to old ice
256	REAL :: albedo_h2o_frost ! albedo of h2o frost
257	REAL,ALLOCATABLE :: co2ice(:) ! Physics: co2 ice mesh averaged [kg/m^2]
258	#endif
259
260	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
261
262	! Loop variable
263	INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil
264	#ifndef CPP_STD
265	! Parallel variables
266	is_sequential=.true.
267	is_parallel=.false.
268	is_mpi_root=.true.
269	is_omp_root=.true.
270	is_master=.true.
271	#endif
272
273	day_ini=0 !test
274	time_phys=0. !test
275
276	! Some constants
277
278	ngrid=ngridmx
279	nlayer=llm
280
281	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
282	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
283	A =(1/m_co2 - 1/m_noco2)
284	B=1/m_noco2
285
286	year_day=669
287	daysec=88775.
288	dtphys=0
289	timestep=year_day*daysec/year_step
290
291	!------------------------
292
293	! I Initialisation
294	! I_a READ run.def
295
296	!------------------------
297
298	!----------------------------READ run.def ---------------------
299	CALL conf_gcm( 99, .TRUE. )
300	CALL conf_pem
301
302	call infotrac_init
303	nq=nqtot
304
305	!------------------------
306
307	! I Initialisation
308	! I_a READ run.def
309	! I_b READ of start_evol.nc and starfi_evol.nc
310
311	!------------------------
312
313	!----------------------------Initialisation : READ some constant of startfi_evol.nc ---------------------
314
315	! In the gcm, these values are given to the physic by the dynamic.
316	! Here we simply read them in the startfi_evol.nc file
317	status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
318
319	allocate(longitude(ngrid))
320	allocate(latitude(ngrid))
321	allocate(cell_area(ngrid))
322
323	status = nf90_inq_varid(ncid, "longitude", lonvarid)
324	status = nf90_get_var(ncid, lonvarid, longitude)
325
326	status = nf90_inq_varid(ncid, "latitude", latvarid)
327	status = nf90_get_var(ncid, latvarid, latitude)
328
329	status = nf90_inq_varid(ncid, "area", areavarid)
330	status = nf90_get_var(ncid, areavarid, cell_area)
331
332	call ini_comsoil_h(ngrid)
333
334	status = nf90_inq_varid(ncid, "soildepth", sdvarid)
335	status = nf90_get_var(ncid, sdvarid, mlayer)
336
337	status =nf90_close(ncid)
338
339	!----------------------------READ start.nc ---------------------
340
341	allocate(q(ip1jmp1,llm,nqtot))
342	CALL dynetat0(FILE_NAME_start,vcov,ucov, &
343	teta,q,masse,ps,phis, time_0)
344
345	allocate(ps_start_GCM(ngrid))
346	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_GCM)
347
348	CALL iniconst !new
349	CALL inigeom
350	allocate(ap(nlayer+1))
351	allocate(bp(nlayer+1))
352	status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
353	status = nf90_inq_varid(ncid, "ap", apvarid)
354	status = nf90_get_var(ncid, apvarid, ap)
355	status = nf90_inq_varid(ncid, "bp", bpvarid)
356	status = nf90_get_var(ncid, bpvarid, bp)
357	status =nf90_close(ncid)
358
359	CALL iniphysiq(iim,jjm,llm, &
360	(jjm-1)*iim+2,comm_lmdz, &
361	daysec,day_ini,dtphys/nsplit_phys, &
362	rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
363	iflag_phys)
364
365	!----------------------------READ startfi.nc ---------------------
366
367	! First we read the initial state (starfi.nc)
368
369	allocate(watercap_slope(ngrid,nslope))
370	allocate(TI_GCM(ngrid,nsoilmx,nslope))
371	allocate(inertiesoil(ngrid,nsoilmx))
372
373	#ifndef CPP_STD
374	CALL phyetat0 (FILE_NAME,0,0, &
375	nsoilmx,ngrid,nlayer,nq, &
376	day_ini,time_phys, &
377	tsurf,tsoil,albedo,emis, &
378	q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar, &
379	watercap,inertiesoil,nslope,tsurf_slope, &
380	tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
381	subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM, &
382	qsurf_slope,watercap_slope)
383
384
385
386	! Remove unphysical values of surface tracer
387	DO i=1,ngrid
388	DO nnq=1,nqtot
389	DO islope=1,nslope
390	if(qsurf_slope(i,nnq,islope).LT.0) then
391	qsurf_slope(i,nnq,islope)=0.
392	endif
393	enddo
394	enddo
395	enddo
396
397	call surfini(ngrid,qsurf)
398
399	#else
400	call phys_state_var_init(nq)
401	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
402	call initracer(ngrid,nq)
403	call iniaerosol()
404	call phyetat0(.true., &
405	ngrid,nlayer,FILE_NAME,0,0,nsoilmx,nq, &
406	day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf, &
407	cloudfrac,totcloudfrac,hice, &
408	rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
409	call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
410
411	call ini_comslope_h(ngrid,nsoilmx,nq)
412
413	allocate(co2ice(ngrid))
414	co2ice(:)=qsurf(:,igcm_co2_ice)
415	co2ice_slope(:,1)=co2ice(:)
416	tsurf_slope(:,1)=tsurf(:)
417
418	if (nslope.eq.1) then
419	def_slope(1) = 0
420	def_slope(2) = 0
421	def_slope_mean=0
422	subslope_dist(:,1) = 1.
423	endif
424
425	! Remove unphysical values of surface tracer
426	DO i=1,ngrid
427	DO nnq=1,nqtot
428	qsurf_slope(i,nnq,1)=qsurf(i,nnq)
429	if(qsurf(i,nnq).LT.0) then
430	qsurf(i,nnq)=0.
431	endif
432	enddo
433	enddo
434	#endif
435
436	DO nnq=1,nqtot
437	if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
438	ENDDO
439
440	!------------------------
441
442	! I Initialisation
443	! I_a READ run.def
444	! I_b READ of start_evol.nc and starfi_evol.nc
445	! I_c Subslope parametrisation
446
447	!------------------------
448
449	!----------------------------Subslope parametrisation definition ---------------------
450
451	! Define some slope statistics
452	iflat=1
453	DO islope=2,nslope
454	IF(abs(def_slope_mean(islope)).lt. &
455	abs(def_slope_mean(iflat))) THEN
456	iflat = islope
457	ENDIF
458	ENDDO
459
460	PRINT*,'Flat slope for islope = ',iflat
461	PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
462
463	allocate(flag_co2flow(ngrid,nslope))
464	allocate(flag_co2flow_mesh(ngrid))
465
466	flag_co2flow(:,:) = 0.
467	flag_co2flow_mesh(:) = 0.
468
469	!---------------------------- READ GCM data ---------------------
470
471	! I Initialisation
472	! I_a READ run.def
473	! I_b READ of start_evol.nc and starfi_evol.nc
474	! I_c Subslope parametrisation
475	! I_d READ GCM data
476
477	!------------------------
478
479	! 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
480
481	call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
482
483	allocate(vmr_co2_gcm(ngrid,timelen))
484	allocate(ps_timeseries(ngrid,timelen))
485	allocate(min_co2_ice_1(ngrid,nslope))
486	allocate(min_h2o_ice_1(ngrid,nslope))
487	allocate(min_co2_ice_2(ngrid,nslope))
488	allocate(min_h2o_ice_2(ngrid,nslope))
489	allocate(tsurf_ave_yr1(ngrid,nslope))
490	allocate(tsurf_ave(ngrid,nslope))
491	allocate(tsoil_ave_yr1(ngrid,nsoilmx,nslope))
492	allocate(tsoil_ave(ngrid,nsoilmx,nslope))
493	allocate(tsurf_GCM_timeseries(ngrid,nslope,timelen))
494	allocate(tsoil_GCM_timeseries(ngrid,nsoilmx,nslope,timelen))
495	allocate(q_co2_PEM_phys(ngrid,timelen))
496	allocate(q_h2o_PEM_phys(ngrid,timelen))
497	allocate(co2_ice_GCM_slope(ngrid,nslope,timelen))
498	allocate(watersurf_density_ave(ngrid,nslope))
499	allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
500
501	allocate(tsoil_ave_PEM_yr1(ngrid,nsoilmx_PEM,nslope))
502	allocate(Tsurfave_before_saved(ngrid,nslope))
503	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
504	allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
505	allocate(watersoil_density_PEM_ave(ngrid,nsoilmx_PEM,nslope))
506	allocate(delta_co2_adsorbded(ngrid))
507	allocate(delta_h2o_adsorbded(ngrid))
508	allocate(vmr_co2_pem_phys(ngrid,timelen))
509
510	print *, "Downloading data Y1..."
511
512	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,&
513	tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM_slope, &
514	watersurf_density_ave,watersoil_density_timeseries)
515
516	! 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
517
518	print *, "Downloading data Y1 done"
519	print *, "Downloading data Y2"
520
521	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, &
522	tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,co2_ice_GCM_slope, &
523	watersurf_density_ave,watersoil_density_timeseries)
524
525	print *, "Downloading data Y2 done"
526
527	!------------------------
528
529	! I Initialisation
530	! I_a READ run.def
531	! I_b READ of start_evol.nc and starfi_evol.nc
532	! I_c Subslope parametrisation
533	! I_d READ GCM data and convert to the physical grid
534	! I_e Initialisation of the PEM variable and soil
535
536	!------------------------
537
538	!---------------------------- Initialisation of the PEM soil and values ---------------------
539
540	call end_comsoil_h_PEM
541	call ini_comsoil_h_PEM(ngrid,nslope)
542	call end_adsorption_h_PEM
543	call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
544
545	allocate(ice_depth(ngrid,nslope))
546	ice_depth(:,:) = 0.
547
548	if(soil_pem) then
549	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM,TI_PEM)
550	DO l=1,nsoilmx
551	tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,l,:)
552	tsoil_PEM(:,l,:)=tsoil_ave(:,l,:)
553	tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:)
554	watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:)
555	ENDDO
556	DO l=nsoilmx+1,nsoilmx_PEM
557	tsoil_ave_PEM_yr1(:,l,:)=tsoil_ave_yr1(:,nsoilmx,:)
558	tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:)
559	watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:)
560	ENDDO
561	watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
562	deallocate(tsoil_ave_yr1)
563	deallocate(tsoil_ave)
564	deallocate(tsoil_GCM_timeseries)
565	endif !soil_pem
566
567	!------------------------
568
569	! I Initialisation
570	! I_a READ run.def
571	! I_b READ of start_evol.nc and starfi_evol.nc
572	! I_c Subslope parametrisation
573	! I_d READ GCM data and convert to the physical grid
574	! I_e Initialisation of the PEM variable and soil
575	! I_f Compute tendencies & Save initial situation
576
577	!----- Compute tendencies from the PCM run
578
579	allocate(tendencies_co2_ice(ngrid,nslope))
580	allocate(tendencies_co2_ice_ini(ngrid,nslope))
581	allocate(tendencies_h2o_ice(ngrid,nslope))
582
583	! Compute the tendencies of the evolution of ice over the years
584
585	call compute_tendencies_slope(ngrid,nslope,min_co2_ice_1,&
586	min_co2_ice_2,tendencies_co2_ice)
587
588	tendencies_co2_ice_ini(:,:)=tendencies_co2_ice(:,:)
589
590	call compute_tendencies_slope(ngrid,nslope,min_h2o_ice_1,&
591	min_h2o_ice_2,tendencies_h2o_ice)
592
593	deallocate(min_co2_ice_1)
594	deallocate(min_co2_ice_2)
595	deallocate(min_h2o_ice_1)
596	deallocate(min_h2o_ice_2)
597
598	!------------------------
599
600	! I Initialisation
601	! I_a READ run.def
602	! I_b READ of start_evol.nc and starfi_evol.nc
603	! I_c Subslope parametrisation
604	! I_d READ GCM data and convert to the physical grid
605	! I_e Initialisation of the PEM variable and soil
606	! I_f Compute tendencies & Save initial situation
607	! I_g Save initial PCM situation
608
609	!---------------------------- Save initial PCM situation ---------------------
610
611	allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
612	allocate(initial_co2_ice_slope(ngrid,nslope))
613	allocate(initial_h2o_ice_slope(ngrid,nslope))
614
615	! We save the places where water ice is sublimating
616	! We compute the surface of water ice sublimating
617	ini_surf_co2=0.
618	ini_surf_h2o=0.
619	Total_surface=0.
620	do i=1,ngrid
621	Total_surface=Total_surface+cell_area(i)
622	do islope=1,nslope
623	if (tendencies_co2_ice(i,islope).LT.0) then
624	initial_co2_ice_sublim_slope(i,islope)=1.
625	ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
626	else
627	initial_co2_ice_sublim_slope(i,islope)=0.
628	endif
629	if (co2ice_slope(i,islope).GT.0) then
630	initial_co2_ice_slope(i,islope)=1.
631	else
632	initial_co2_ice_slope(i,islope)=0.
633	endif
634	if (tendencies_h2o_ice(i,islope).LT.0) then
635	initial_h2o_ice_slope(i,islope)=1.
636	ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
637	else
638	initial_h2o_ice_slope(i,islope)=0.
639	endif
640	enddo
641	enddo
642
643	print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
644	print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
645	print *, "Total surface of the planet=", Total_surface
646
647	allocate(zplev_gcm(ngrid,nlayer+1))
648
649	DO l=1,nlayer+1
650	DO ig=1,ngrid
651	zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
652	ENDDO
653	ENDDO
654
655	global_ave_press_old=0.
656	do i=1,ngrid
657	global_ave_press_old=global_ave_press_old+cell_area(i)*ps_start_GCM(i)/Total_surface
658	enddo
659
660	global_ave_press_GCM=global_ave_press_old
661	global_ave_press_new=global_ave_press_old
662	print *, "Initial global average pressure=", global_ave_press_GCM
663
664	!------------------------
665
666	! I Initialisation
667	! I_a READ run.def
668	! I_b READ of start_evol.nc and starfi_evol.nc
669	! I_c Subslope parametrisation
670	! I_d READ GCM data and convert to the physical grid
671	! I_e Initialisation of the PEM variable and soil
672	! I_f Compute tendencies & Save initial situation
673	! I_g Save initial PCM situation
674	! I_h Read the PEMstart
675
676	!---------------------------- Read the PEMstart ---------------------
677
678	call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_PEM_yr1,tsoil_PEM,ice_depth, &
679	tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,&
680	tendencies_h2o_ice,tendencies_co2_ice,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
681	watersurf_density_ave,watersoil_density_PEM_ave, &
682	co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
683
684	do ig = 1,ngrid
685	do islope = 1,nslope
686	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.)
687	enddo
688	enddo
689
690	if(adsorption_pem) then
691	totmassco2_adsorbded = 0.
692	totmassh2o_adsorbded = 0.
693	do ig = 1,ngrid
694	do islope =1, nslope
695	do l = 1,nsoilmx_PEM - 1
696	totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
697	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
698	cell_area(ig)
699	totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
700	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
701	cell_area(ig)
702	enddo
703	enddo
704	enddo
705
706	write(,) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
707	write(,) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
708	endif ! adsorption
709	deallocate(tsoil_ave_PEM_yr1)
710	deallocate(tsurf_ave_yr1)
711
712	!------------------------
713
714	! I Initialisation
715	! I_a READ run.def
716	! I_b READ of start_evol.nc and starfi_evol.nc
717	! I_c Subslope parametrisation
718	! I_d READ GCM data and convert to the physical grid
719	! I_e Initialisation of the PEM variable and soil
720	! I_f Compute tendencies & Save initial situation
721	! I_g Save initial PCM situation
722	! I_h Read the PEMstar
723	! I_i Compute orbit criterion
724	#ifndef CPP_STD
725	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
726	#else
727	CALL iniorbit(apoastr, periastr, year_day, peri_day,obliquit)
728	#endif
729
730	if(evol_orbit_pem) then
731	call orbit_param_criterion(year_iter_max)
732	else
733	year_iter_max=Max_iter_pem
734	endif
735
736	!--------------------------- END INITIALISATION ---------------------
737
738	!---------------------------- RUN ---------------------
739
740	!------------------------
741
742	! II Run
743	! II_a update pressure,ice and tracers
744
745	!------------------------
746	year_iter=0
747	do while (year_iter.LT.year_iter_max)
748
749	! II.a.1. Compute updated global pressure
750	print *, "Recomputing the new pressure..."
751	do i=1,ngrid
752	do islope=1,nslope
753	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
754	enddo
755	enddo
756	print *, 'Global average pressure old time step',global_ave_press_old
757
758	call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
759
760	if(adsorption_pem) then
761	do i=1,ngrid
762	global_ave_press_new = global_ave_press_new -gcell_area(i)delta_co2_adsorbded(i)/Total_surface
763	enddo
764	print *, 'Global average pressure old time step',global_ave_press_old
765	print *, 'Global average pressure new time step',global_ave_press_new
766	endif
767
768	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
769	allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
770	print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
771	DO l=1,nlayer+1
772	DO ig=1,ngrid
773	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
774	ENDDO
775	ENDDO
776
777	! II.a.3. Surface pressure timeseries
778	print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
779	do ig = 1,ngrid
780	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
781	enddo
782
783	! II.a.4. New pressure levels timeseries
784	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
785	print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
786	do l=1,nlayer+1
787	do ig=1,ngrid
788	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
789	enddo
790	enddo
791
792	! II.a.5. Tracers timeseries
793	print *, "Recomputing of tracer VMR timeseries for the new pressure..."
794
795	l=1
796	DO ig=1,ngrid
797	DO t = 1, timelen
798	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))
799	if(q_h2o_PEM_phys(ig,t).LT.0) then
800	q_h2o_PEM_phys(ig,t)=1E-30
801	endif
802	if(q_h2o_PEM_phys(ig,t).GT.1) then
803	q_h2o_PEM_phys(ig,t)=1.
804	endif
805	enddo
806	enddo
807
808	DO ig=1,ngrid
809	DO t = 1, timelen
810	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)) &
811	+ ( (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) - &
812	(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t)) ) / &
813	(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
814	if (q_co2_PEM_phys(ig,t).LT.0) then
815	q_co2_PEM_phys(ig,t)=1E-30
816	elseif (q_co2_PEM_phys(ig,t).GT.1) then
817	q_co2_PEM_phys(ig,t)=1.
818	endif
819	mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
820	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
821	ENDDO
822	ENDDO
823
824	deallocate(zplev_new_timeseries)
825	deallocate(zplev_old_timeseries)
826
827	! II Run
828	! II_a update pressure, ice and tracers
829	! II_b Evolution of the ice
830
831	! II.b. Evolution of the ice
832	print *, "Evolution of h2o ice"
833
834	call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
835
836	DO islope=1, nslope
837	write(str2(1:2),'(i2.2)') islope
838	call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf_slope(:,igcm_h2o_ice,islope))
839	call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
840	ENDDO
841
842	print *, "Evolution of co2 ice"
843	call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
844
845	!------------------------
846
847	! II Run
848	! II_a update pressure, ice and tracers
849	! II_b Evolution of the ice
850	! II_c CO2 glaciers flows
851
852	!------------------------
853
854	print *, "CO2 glacier flows"
855
856	if (co2glaciersflow) then
857	call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,&
858	global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
859	endif
860
861
862
863	DO islope=1, nslope
864	write(str2(1:2),'(i2.2)') islope
865	call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,co2ice_slope(:,islope))
866	call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
867	ENDDO
868
869	!------------------------
870
871	! II Run
872	! II_a update pressure, ice and tracers
873	! II_b Evolution of the ice
874	! II_c CO2 glaciers flows
875	! II_d Update surface and soil temperatures
876
877	!------------------------
878
879	! II_d.1 Update Tsurf
880
881	print *, "Updating the new Tsurf"
882	bool_sublim=0
883	Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
884	DO ig = 1,ngrid
885	DO islope = 1,nslope
886	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
887	if(latitude_deg(ig).gt.0) then
888	do ig_loop=ig,ngrid
889	DO islope_loop=islope,iflat,-1
890	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
891	tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
892	bool_sublim=1
893	exit
894	endif
895	enddo
896	if (bool_sublim.eq.1) then
897	exit
898	endif
899	enddo
900	else
901	do ig_loop=ig,1,-1
902	DO islope_loop=islope,iflat
903	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
904	tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
905	bool_sublim=1
906	exit
907	endif
908	enddo
909	if (bool_sublim.eq.1) then
910	exit
911	endif
912	enddo
913	endif
914	initial_co2_ice_slope(ig,islope)=0
915	if ((co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then
916	albedo_slope(ig,1,islope) = albedo_h2o_frost
917	albedo_slope(ig,2,islope) = albedo_h2o_frost
918	else
919	albedo_slope(ig,1,islope) = albedodat(ig)
920	albedo_slope(ig,2,islope) = albedodat(ig)
921	emiss_slope(ig,islope) = emissiv
922	endif
923	elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
924	ave=0.
925	do t=1,timelen
926	if(co2_ice_GCM_slope(ig,islope,t).gt.1e-3) then
927	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
928	else
929	ave = ave + tsurf_GCM_timeseries(ig,islope,t)
930	endif
931	enddo
932	tsurf_ave(ig,islope)=ave/timelen
933	endif
934	enddo
935	enddo
936
937	do t = 1,timelen
938	tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) +( tsurf_ave(:,:) -Tsurfave_before_saved(:,:))
939	enddo
940	! for the start
941	do ig = 1,ngrid
942	do islope = 1,nslope
943	tsurf_slope(ig,islope) = tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope))
944	enddo
945	enddo
946
947	DO islope=1, nslope
948	write(str2(1:2),'(i2.2)') islope
949	call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_slope(:,islope))
950	ENDDO
951
952	if(soil_pem) then
953
954	! II_d.2 Update soil temperature
955
956	allocate(TI_locslope(ngrid,nsoilmx_PEM))
957	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
958	allocate(Tsurf_locslope(ngrid))
959	print *,"Updating soil temperature"
960
961	! Soil averaged
962	do islope = 1,nslope
963	TI_locslope(:,:) = TI_PEM(:,:,islope)
964	do t = 1,timelen
965	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
966	Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
967	call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
968	call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
969	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
970	do ig = 1,ngrid
971	do isoil = 1,nsoilmx_PEM
972	watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil)
973	if(isnan(Tsoil_locslope(ig,isoil))) then
974	call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1)
975	endif
976	enddo
977	enddo
978	enddo
979	enddo
980	tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
981	watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
982
983	print *, "Update of soil temperature done"
984
985	deallocate(TI_locslope)
986	deallocate(Tsoil_locslope)
987	deallocate(Tsurf_locslope)
988	write(,) "Compute ice table"
989
990	! II_d.3 Update the ice table
991	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,ice_depth)
992
993	print *, "Update soil propreties"
994
995	! II_d.4 Update the soil thermal properties
996	call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
997	ice_depth,TI_PEM)
998
999	! II_d.5 Update the mass of the regolith adsorbded
1000	if(adsorption_pem) then
1001	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf_slope(:,igcm_h2o_ice,:),co2ice_slope, &
1002	tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
1003	h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
1004	endif
1005	endif !soil_pem
1006
1007	!------------------------
1008
1009	! II Run
1010	! II_a update pressure, ice and tracers
1011	! II_b Evolution of the ice
1012	! II_c CO2 glaciers flows
1013	! II_d Update surface and soil temperatures
1014	! II_e Update the tendencies
1015
1016	!------------------------
1017
1018	print *, "Adaptation of the new co2 tendencies to the current pressure"
1019	call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,co2ice_slope,vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,&
1020	global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
1021
1022	!------------------------
1023
1024	! II Run
1025	! II_a update pressure, ice and tracers
1026	! II_b Evolution of the ice
1027	! II_c CO2 glaciers flows
1028	! II_d Update surface and soil temperatures
1029	! II_e Update the tendencies
1030	! II_f Checking the stopping criterion
1031
1032	!------------------------
1033	call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
1034
1035	call criterion_co2_stop(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid, &
1036	initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
1037
1038	year_iter=year_iter+dt_pem
1039
1040	print *, "Checking all the stopping criterion."
1041	if (STOPPING_water) then
1042	print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
1043	criterion_stop=1
1044	endif
1045	if (STOPPING_1_water) then
1046	print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
1047	criterion_stop=1
1048	endif
1049	if (STOPPING_co2) then
1050	print *, "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
1051	criterion_stop=2
1052	endif
1053	if (STOPPING_1_co2) then
1054	print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
1055	criterion_stop=2
1056	endif
1057	if (STOPPING_pressure) then
1058	print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
1059	criterion_stop=3
1060	endif
1061	if (year_iter.ge.year_iter_max) then
1062	print *, "STOPPING because maximum number of iterations reached"
1063	criterion_stop=4
1064	endif
1065
1066	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2 .or. STOPPING_pressure) then
1067	exit
1068	else
1069	print *, "We continue!"
1070	print *, "Number of iteration of the PEM : year_iter=", year_iter
1071	endif
1072
1073	global_ave_press_old=global_ave_press_new
1074
1075	enddo ! big time iteration loop of the pem
1076
1077
1078	!---------------------------- END RUN PEM ---------------------
1079
1080	!---------------------------- OUTPUT ---------------------
1081
1082	!------------------------
1083
1084	! III Output
1085	! III_a Update surface value for the PCM start files
1086
1087	!------------------------
1088
1089	! III_a.1 Ice update (for startfi)
1090	! Co2 ice
1091	DO ig = 1,ngrid
1092	co2ice(ig) = 0.
1093	DO islope = 1,nslope
1094	co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
1095	* subslope_dist(ig,islope) / &
1096	cos(pi*def_slope_mean(islope)/180.)
1097	ENDDO
1098	#ifdef CPP_STD
1099	qsurf(ig,igcm_co2_ice)=co2ice(ig)
1100	#endif
1101	ENDDO ! of DO ig=1,ngrid
1102
1103	! H2O ice
1104	DO ig=1,ngrid
1105	if(watercaptag(ig)) then
1106	watercap_sum=0.
1107	DO islope=1,nslope
1108	watercap_slope_saved = watercap_slope(ig,islope)
1109	if(qsurf_slope(ig,igcm_h2o_ice,islope).GT. (watercap_slope(ig,islope)+water_reservoir(ig)cos(pidef_slope_mean(islope)/180.))) then
1110	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.))
1111	else
1112	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.))
1113	qsurf_slope(ig,igcm_h2o_ice,islope)=0.
1114	endif
1115	watercap_sum=watercap_sum+(watercap_slope(ig,islope)-watercap_slope_saved)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
1116	watercap_slope(ig,islope)=0.
1117	enddo
1118	water_reservoir(ig)=water_reservoir(ig)+watercap_sum
1119	endif
1120	enddo
1121
1122	DO ig = 1,ngrid
1123	qsurf(ig,igcm_h2o_ice) = 0.
1124	DO islope = 1,nslope
1125	qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
1126	* subslope_dist(ig,islope) / &
1127	cos(pi*def_slope_mean(islope)/180.)
1128	ENDDO
1129	ENDDO ! of DO ig=1,ngrid
1130
1131	DO ig=1,ngrid
1132	DO islope=1,nslope
1133	if(qsurf_slope(ig,igcm_h2o_ice,islope).GT.500) then
1134	watercaptag(ig)=.true.
1135	qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250
1136	water_reservoir(ig)=water_reservoir(ig)+250
1137	endif
1138	enddo
1139	enddo
1140
1141	DO ig=1,ngrid
1142	if(water_reservoir(ig).LT. 10) then
1143	watercaptag(ig)=.false.
1144	qsurf(ig,igcm_h2o_ice)=qsurf(ig,igcm_h2o_ice)+water_reservoir(ig)
1145	DO islope=1,nslope
1146	qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+water_reservoir(ig)
1147	ENDDO
1148	endif
1149	enddo
1150
1151	DO ig = 1,ngrid
1152	watercap(ig) = 0.
1153	DO islope = 1,nslope
1154	watercap(ig) = watercap(ig) + watercap_slope(ig,islope) &
1155	* subslope_dist(ig,islope) / &
1156	cos(pi*def_slope_mean(islope)/180.)
1157	ENDDO
1158	ENDDO ! of DO ig=1,ngrid
1159
1160	! III_a.2 Tsoil update (for startfi)
1161
1162	if(soil_pem) then
1163	fluxgeo_slope(:,:) = fluxgeo
1164	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM)
1165	tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
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