Context Navigation

pem.F90 @ 2950

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

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

Download in other formats:

Original Format