Context Navigation

pem.F90 @ 3026

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

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

Download in other formats:

Original Format