Context Navigation

pem.F90 @ 3000

Last change on this file since 3000 was 2998, checked in by llange, 2 years ago

MARS PEM

Fix some bug in water conservations (transfer from watercap to waterreservoir was wrong)
Add the H2O glaciers in the main program.
Add the distinction between co2 frost (qsurf) and perenial co2 ice. Will be developed in the Mars PCM in a future commit
Thickness thresolhd to pass from frost to glaciers is now not harcoded but a constant fixed in the mars_pem_constants

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

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

Download in other formats:

Original Format