Context Navigation

pem.F90 @ 3015

Last change on this file since 3015 was 3002, checked in by llange, 2 years ago

MARS PEM

Fix a bug in conf_pem ("Flux_geo" was expected in the .def while it should be 'fluxgeo', corrected)
Fix a bug in glaciers modules (the length for 'name_ice' was too large)
Fix a bug when writing tsoil in the PEM (ngrid x nsoil_PEM x nslope was given, while it was expecting a ngrid x nsoil_GCM x nslope)

PEM runs correctly now
LL

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	call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
842
843	if(adsorption_pem) then
844	do i=1,ngrid
845	global_ave_press_new = global_ave_press_new -gcell_area(i)delta_co2_adsorbded(i)/Total_surface
846	enddo
847	print *, 'Global average pressure old time step',global_ave_press_old
848	print *, 'Global average pressure new time step',global_ave_press_new
849	endif
850
851	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
852	allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
853	print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
854	DO l=1,nlayer+1
855	DO ig=1,ngrid
856	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
857	ENDDO
858	ENDDO
859
860	! II.a.3. Surface pressure timeseries
861	print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
862	do ig = 1,ngrid
863	ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_ave_press_new/global_ave_press_old
864	enddo
865
866	! II.a.4. New pressure levels timeseries
867	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
868	print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
869	do l=1,nlayer+1
870	do ig=1,ngrid
871	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:)
872	enddo
873	enddo
874
875	! II.a.5. Tracers timeseries
876	print *, "Recomputing of tracer VMR timeseries for the new pressure..."
877
878	l=1
879	DO ig=1,ngrid
880	DO t = 1, timelen
881	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))
882	if(q_h2o_PEM_phys(ig,t).LT.0) then
883	q_h2o_PEM_phys(ig,t)=1E-30
884	endif
885	if(q_h2o_PEM_phys(ig,t).GT.1) then
886	q_h2o_PEM_phys(ig,t)=1.
887	endif
888	enddo
889	enddo
890
891	DO ig=1,ngrid
892	DO t = 1, timelen
893	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)) &
894	+ ( (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) - &
895	(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t)) ) / &
896	(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
897	if (q_co2_PEM_phys(ig,t).LT.0) then
898	q_co2_PEM_phys(ig,t)=1E-30
899	elseif (q_co2_PEM_phys(ig,t).GT.1) then
900	q_co2_PEM_phys(ig,t)=1.
901	endif
902	mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
903	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
904	ENDDO
905	ENDDO
906
907	deallocate(zplev_new_timeseries)
908	deallocate(zplev_old_timeseries)
909
910	! II Run
911	! II_a update pressure, ice and tracers
912	! II_b Evolution of the ice
913
914	! II.b. Evolution of the ice
915	print *, "Evolution of h2o ice"
916	call evol_h2o_ice_s(qsurf(:,igcm_h2o_ice,:),tendencies_h2o_ice,iim,jjm_value,ngrid,cell_area,STOPPING_1_water,nslope)
917
918	print *, "Evolution of co2 ice"
919	call evol_co2_ice_s(qsurf(:,igcm_co2,:),tendencies_co2_ice,iim,jjm_value,ngrid,cell_area,nslope)
920
921	DO islope=1, nslope
922	write(str2(1:2),'(i2.2)') islope
923	call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf(:,igcm_h2o_ice,islope))
924	call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice(:,islope))
925	call WRITEDIAGFI(ngrid,'c2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
926	call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
927	ENDDO
928
929	!------------------------
930
931	! II Run
932	! II_a update pressure, ice and tracers
933	! II_b Evolution of the ice
934	! II_c CO2 & H2O glaciers flows
935
936	!------------------------
937
938	print *, "CO2 glacier flows"
939
940	if (co2glaciersflow) then
941	call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_timeseries,&
942	global_ave_press_GCM,global_ave_press_new,qsurf(:,igcm_co2,:),flag_co2flow,flag_co2flow_mesh)
943	endif
944
945	print *, "H2O glacier flows"
946
947	if (h2oglaciersflow) then
948	call h2oglaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_ave,qsurf(:,igcm_h2o_ice,:),flag_h2oflow,flag_h2oflow_mesh)
949	endif
950	DO islope=1, nslope
951	write(str2(1:2),'(i2.2)') islope
952	call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,qsurf(:,igcm_co2,islope))
953	call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice(:,islope))
954	call WRITEDIAGFI(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope))
955	ENDDO
956
957	!------------------------
958
959	! II Run
960	! II_a update pressure, ice and tracers
961	! II_b Evolution of the ice
962	! II_c CO2 glaciers flows
963	! II_d Update surface and soil temperatures
964
965	!------------------------
966
967	! II_d.1 Update Tsurf
968
969	print *, "Updating the new Tsurf"
970	bool_sublim=.false.
971	Tsurfave_before_saved(:,:) = tsurf_ave(:,:)
972	DO ig = 1,ngrid
973	DO islope = 1,nslope
974	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
975	if(latitude_deg(ig).gt.0) then
976	do ig_loop=ig,ngrid
977	DO islope_loop=islope,iflat,-1
978	if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then
979	tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
980	bool_sublim=.true.
981	exit
982	endif
983	enddo
984	if (bool_sublim.eqv. .true.) then
985	exit
986	endif
987	enddo
988	else
989	do ig_loop=ig,1,-1
990	DO islope_loop=islope,iflat
991	if(initial_co2_ice(ig_loop,islope_loop).lt.0.5 .and. qsurf(ig_loop,igcm_co2,islope_loop).LT. 1E-10) then
992	tsurf_ave(ig,islope)=tsurf_ave(ig_loop,islope_loop)
993	bool_sublim=.true.
994	exit
995	endif
996	enddo
997	if (bool_sublim.eqv. .true.) then
998	exit
999	endif
1000	enddo
1001	endif
1002	initial_co2_ice(ig,islope)=0
1003	if ((qsurf(ig,igcm_co2,islope).lt.1e-10).and. (qsurf(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then
1004	albedo(ig,1,islope) = albedo_h2o_frost
1005	albedo(ig,2,islope) = albedo_h2o_frost
1006	else
1007	albedo(ig,1,islope) = albedodat(ig)
1008	albedo(ig,2,islope) = albedodat(ig)
1009	emis(ig,islope) = emissiv
1010	endif
1011	elseif( (qsurf(ig,igcm_co2,islope).GT. 1E-3).and.(tendencies_co2_ice(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
1012	ave=0.
1013	do t=1,timelen
1014	if(co2_ice_GCM(ig,islope,t).gt.1e-3) then
1015	ave = ave + beta_clap_co2/(alpha_clap_co2-log(vmr_co2_pem_phys(ig,t)*ps_timeseries(ig,t)/100.))
1016	else
1017	ave = ave + tsurf_GCM_timeseries(ig,islope,t)
1018	endif
1019	enddo
1020	tsurf_ave(ig,islope)=ave/timelen
1021	endif
1022	enddo
1023	enddo
1024
1025	do t = 1,timelen
1026	tsurf_GCM_timeseries(:,:,t) = tsurf_GCM_timeseries(:,:,t) +( tsurf_ave(:,:) -Tsurfave_before_saved(:,:))
1027	enddo
1028	! for the start
1029	do ig = 1,ngrid
1030	do islope = 1,nslope
1031	tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave(ig,islope))
1032	enddo
1033	enddo
1034
1035	DO islope=1, nslope
1036	write(str2(1:2),'(i2.2)') islope
1037	call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
1038	ENDDO
1039
1040	if(soil_pem) then
1041
1042	! II_d.2 Update soil temperature
1043
1044	allocate(TI_locslope(ngrid,nsoilmx_PEM))
1045	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
1046	allocate(Tsurf_locslope(ngrid))
1047	print *,"Updating soil temperature"
1048
1049	! Soil averaged
1050	do islope = 1,nslope
1051	TI_locslope(:,:) = TI_PEM(:,:,islope)
1052	do t = 1,timelen
1053	Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t)
1054	Tsurf_locslope(:) = tsurf_GCM_timeseries(:,islope,t)
1055	call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
1056	call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope)
1057	tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:)
1058	do ig = 1,ngrid
1059	do isoil = 1,nsoilmx_PEM
1060	watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)
1061	if(isnan(Tsoil_locslope(ig,isoil))) then
1062	call abort_pem("PEM - Update Tsoil","NAN detected in Tsoil ",1)
1063	endif
1064	enddo
1065	enddo
1066	enddo
1067	enddo
1068	tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
1069	watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
1070
1071	print *, "Update of soil temperature done"
1072
1073	deallocate(TI_locslope)
1074	deallocate(Tsoil_locslope)
1075	deallocate(Tsurf_locslope)
1076	write(,) "Compute ice table"
1077
1078	! II_d.3 Update the ice table
1079	call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
1080	print *, "Update soil propreties"
1081
1082	! II_d.4 Update the soil thermal properties
1083	call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, &
1084	porefillingice_depth,porefillingice_thickness,TI_PEM)
1085
1086	! II_d.5 Update the mass of the regolith adsorbded
1087	if(adsorption_pem) then
1088	call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,tendencies_h2o_ice,tendencies_co2_ice,qsurf(:,igcm_h2o_ice,:),qsurf(:,igcm_co2,:), &
1089	tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
1090	h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
1091	endif
1092	endif !soil_pem
1093
1094	!------------------------
1095
1096	! II Run
1097	! II_a update pressure, ice and tracers
1098	! II_b Evolution of the ice
1099	! II_c CO2 glaciers flows
1100	! II_d Update surface and soil temperatures
1101	! II_e Update the tendencies
1102
1103	!------------------------
1104
1105	print *, "Adaptation of the new co2 tendencies to the current pressure"
1106	call recomp_tend_co2_slope(tendencies_co2_ice,tendencies_co2_ice_ini,qsurf(:,igcm_co2,:),vmr_co2_gcm,vmr_co2_pem_phys,ps_timeseries,&
1107	global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
1108
1109	!------------------------
1110
1111	! II Run
1112	! II_a update pressure, ice and tracers
1113	! II_b Evolution of the ice
1114	! II_c CO2 glaciers flows
1115	! II_d Update surface and soil temperatures
1116	! II_e Update the tendencies
1117	! II_f Checking the stopping criterion
1118
1119	!------------------------
1120	call criterion_waterice_stop(cell_area,ini_surf_h2o,qsurf(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice)
1121
1122	call criterion_co2_stop(cell_area,ini_surf_co2,qsurf(:,igcm_co2,:),STOPPING_co2,STOPPING_pressure,ngrid, &
1123	initial_co2_ice_sublim,global_ave_press_GCM,global_ave_press_new,nslope)
1124
1125	year_iter=year_iter+dt_pem
1126
1127	print *, "Checking all the stopping criterion."
1128	if (STOPPING_water) then
1129	print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
1130	criterion_stop=1
1131	endif
1132	if (STOPPING_1_water) then
1133	print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
1134	criterion_stop=1
1135	endif
1136	if (STOPPING_co2) then
1137	print *, "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2
1138	criterion_stop=2
1139	endif
1140	if (STOPPING_pressure) then
1141	print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure
1142	criterion_stop=3
1143	endif
1144	if (year_iter.ge.year_iter_max) then
1145	print *, "STOPPING because maximum number of iterations reached"
1146	criterion_stop=4
1147	endif
1148
1149	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_pressure) then
1150	exit
1151	else
1152	print *, "We continue!"
1153	print *, "Number of iteration of the PEM : year_iter=", year_iter
1154	endif
1155
1156	global_ave_press_old=global_ave_press_new
1157
1158	enddo ! big time iteration loop of the pem
1159
1160
1161	!---------------------------- END RUN PEM ---------------------
1162
1163	!---------------------------- OUTPUT ---------------------
1164
1165	!------------------------
1166
1167	! III Output
1168	! III_a Update surface value for the PCM start files
1169
1170	!------------------------
1171
1172	! III_a.1 Ice update (for startfi)
1173
1174	! H2O ice
1175	DO ig=1,ngrid
1176	if(watercaptag(ig)) then
1177	watercap_sum=0.
1178	DO islope=1,nslope
1179	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
1180	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
1181	else
1182	! 2nd case: we have sublimate ice: then let's put qsurf = 0. and add the difference in watercap
1183	watercap(ig,islope)=watercap(ig,islope)+qsurf(ig,igcm_h2o_ice,islope)-(watercap(ig,islope)+water_reservoir(ig)cos(pidef_slope_mean(islope)/180.))
1184	qsurf(ig,igcm_h2o_ice,islope)=0.
1185	endif
1186	watercap_sum=watercap_sum+watercap(ig,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
1187	watercap(ig,islope)=0.
1188	enddo
1189	water_reservoir(ig)=water_reservoir(ig)+watercap_sum
1190	endif
1191	enddo
1192
1193	DO ig=1,ngrid
1194	water_sum = 0.
1195	DO islope=1,nslope
1196	water_sum = water_sum + qsurf(ig,igcm_h2o_ice,islope)subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.)
1197	ENDDO
1198	if(watercaptag(ig).eqv..false.) then ! let's check if we have an 'infinite' source of water that has been forming.
1199	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
1200	watercaptag(ig)=.true.
1201	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
1202	DO islope = 1,nslope
1203	qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope) - threshold_water_frost2perenial/2.cos(pidef_slope_mean(islope)/180.)
1204	ENDDO
1205	endif
1206	else ! let's check that the infinite source of water disapear
1207	if((water_sum + water_reservoir(ig)).lt.threshold_water_frost2perenial) then
1208	watercaptag(ig)=.false.
1209	DO islope = 1,nslope
1210	qsurf(ig,igcm_h2o_ice,islope)=qsurf(ig,igcm_h2o_ice,islope)+water_reservoir(ig)cos(pidef_slope_mean(islope)/180.)
1211	ENDDO
1212	water_reservoir(ig) = 0.
1213	endif
1214	endif
1215	enddo
1216
1217	! CO2 ice
1218
1219	DO ig=1,ngrid
1220	DO islope=1,nslope
1221	if(qsurf(ig,igcm_co2,islope).gt.threshold_co2_frost2perenial) then
1222	perenial_co2ice(ig,islope) = 0.5*qsurf(ig,igcm_co2,islope)
1223	qsurf(ig,igcm_co2,islope) = 0.5*qsurf(ig,igcm_co2,islope)
1224	endif
1225	ENDDO
1226	ENDDO
1227
1228
1229	! III_a.2 Tsoil update (for startfi)
1230
1231	if(soil_pem) then
1232	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,inertiesoil)
1233	tsoil(:,:,:) = tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,timelen)
1234	endif !soil_pem
1235
1236	! III_a.4 Pressure (for start)
1237	do i=1,ip1jmp1
1238	ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
1239	enddo
1240
1241	do i = 1,ngrid
1242	ps_start_GCM(i)=ps_start_GCM(i)*global_ave_press_new/global_ave_press_GCM
1243	enddo
1244
1245	! III_a.5 Tracer (for start)
1246	allocate(zplev_new(ngrid,nlayer+1))
1247
1248	do l=1,nlayer+1
1249	do ig=1,ngrid
1250	zplev_new(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
1251	enddo
1252	enddo
1253
1254	DO nnq=1,nqtot
1255	if (noms(nnq).NE."co2") then
1256	DO l=1,llm-1
1257	DO ig=1,ngrid
1258	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))
1259	ENDDO
1260	q(:,llm,nnq)=q(:,llm-1,nnq)
1261	ENDDO
1262	else
1263	DO l=1,llm-1
1264	DO ig=1,ngrid
1265	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)) &
1266	+ ( (zplev_new(ig,l)-zplev_new(ig,l+1)) - &
1267	(zplev_gcm(ig,l)-zplev_gcm(ig,l+1)) ) / &
1268	(zplev_new(ig,l)-zplev_new(ig,l+1))
1269	ENDDO
1270	q(:,llm,nnq)=q(:,llm-1,nnq)
1271	ENDDO
1272	endif
1273	ENDDO
1274
1275	! Conserving the tracers's mass for GCM start files
1276	DO nnq=1,nqtot
1277	DO ig=1,ngrid
1278	DO l=1,llm-1
1279	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
1280	extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
1281	q(ig,l,nnq)=1.
1282	q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
1283	write(,) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number"
1284	endif
1285	if(q(ig,l,nnq).LT.0) then
1286	q(ig,l,nnq)=1E-30
1287	endif
1288	ENDDO
1289	enddo
1290	enddo
1291
1292	!------------------------
1293	if(evol_orbit_pem) then
1294	call recomp_orb_param(year_iter)
1295	endif
1296
1297	! III Output
1298	! III_a Update surface value for the PCM start files
1299	! III_b Write start and starfi.nc
1300
1301	!------------------------
1302
1303	! III_b.1 WRITE restart.nc
1304
1305	ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
1306	pday=day_ini
1307	ztime_fin=0.
1308
1309	allocate(p(ip1jmp1,nlayer+1))
1310	#ifndef CPP_1D
1311	CALL pression (ip1jmp1,ap,bp,ps,p)
1312	CALL massdair(p,masse)
1313
1314	CALL dynredem0("restart_evol.nc", day_ini, phis)
1315
1316	CALL dynredem1("restart_evol.nc", &
1317	time_0,vcov,ucov,teta,q,masse,ps)
1318	print *, "restart_evol.nc has been written"
1319
1320	#else
1321	DO nnq = 1, nqtot
1322	call writeprofile(nlayer,q(1,:,nnq),noms(nnq),nnq,qsurf)
1323	ENDDO
1324	#endif
1325
1326	! III_b.2 WRITE restartfi.nc
1327	#ifndef CPP_STD
1328	call physdem0("restartfi_evol.nc",longitude,latitude, &
1329	nsoilmx,ngrid,nlayer,nq, &
1330	ptimestep,pday,0.,cell_area, &
1331	albedodat,inertiedat,def_slope, &
1332	subslope_dist)
1333
1334	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
1335	ptimestep,ztime_fin, &
1336	tsurf,tsoil,inertiesoil,albedo, &
1337	emis,q2,qsurf,tauscaling,totcloudfrac,wstar,&
1338	watercap,perenial_co2ice)
1339	#else
1340	call physdem0("restartfi_evol.nc",longitude,latitude,nsoilmx,ngrid,nlayer,nq, &
1341	ptimestep,pday,time_phys,cell_area, &
1342	albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
1343
1344	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
1345	ptimestep,ztime_fin, &
1346	tsurf,tsoil,emis,q2,qsurf, &
1347	cloudfrac,totcloudfrac,hice, &
1348	rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
1349	#endif
1350
1351	print *, "restartfi_evol.nc has been written"
1352	!------------------------
1353
1354	! III Output
1355	! III_a Update surface value for the PCM start files
1356	! III_b Write start and starfi.nc
1357	! III_c Write start_pem
1358
1359	!------------------------
1360	call pemdem0("restartfi_PEM.nc",longitude,latitude,cell_area,nsoilmx_PEM,ngrid, &
1361	float(day_ini),0.,nslope,def_slope,subslope_dist)
1362
1363
1364	call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
1365	tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, &
1366	co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
1367
1368	call info_run_PEM(year_iter, criterion_stop)
1369
1370	print *, "restartfi_PEM.nc has been written"
1371	print *, "The PEM had run for ", year_iter, " years."
1372	print *, "LL & RV : So far so good"
1373
1374	deallocate(vmr_co2_gcm)
1375	deallocate(ps_timeseries)
1376	deallocate(tsurf_GCM_timeseries)
1377	deallocate(q_co2_PEM_phys)
1378	deallocate(q_h2o_PEM_phys)
1379	deallocate(co2_ice_GCM)
1380	deallocate(watersurf_density_ave)
1381	deallocate(watersoil_density_timeseries)
1382	deallocate(Tsurfave_before_saved)
1383	deallocate(tsoil_phys_PEM_timeseries)
1384	deallocate(watersoil_density_PEM_timeseries)
1385	deallocate(watersoil_density_PEM_ave)
1386	deallocate(delta_co2_adsorbded)
1387	deallocate(delta_h2o_adsorbded)
1388	deallocate(vmr_co2_pem_phys)
1389
1390	END PROGRAM pem

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

Download in other formats:

Original Format