Context Navigation

pem.F90 @ 2840

Last change on this file since 2840 was 2835, checked in by romain.vande, 3 years ago

Mars PEM:
Introduction of the possibility to follow an orbital forcing.
Introduction of new control parameters.
Cleaning of the PEM (removing unused files, add comments and new files)

A file named run_PEM.def can be added to the run.def. It contains the following variables:

_ evol_orbit_pem: Boolean. Do you want to follow an orbital forcing predefined (read in ob_ex_lsp.asc for example)? (default=false)
_ year_bp_ini: Integer. Number of year before present to start the pem run if evol_orbit_pem=.true. , default=0
_ Max_iter_pem: Integer. Maximal number of iteration if none of the stopping criterion is reached and if evol_orbit_pem=.false., default=99999999
_ dt_pem: Integer. Time step of the PEM in year, default=1
_ alpha_criterion: Real. Acceptance rate of sublimating ice surface change, default=0.2
_ soil_pem: Boolean. Do you want to run with subsurface physical processes in the PEM? default=.true.

File size: 51.6 KB

Line
1
2	!------------------------
3
4	! I Initialisation
5	! I_a READ run.def
6	! I_b READ of start_evol.nc and starfi_evol.nc
7	! I_c Subslope parametrisation
8	! I_d READ GCM data and convert to the physical grid
9	! I_e Initialisation of the PEM variable and soil
10	! I_f Compute tendencies & Save initial situation
11	! I_g Save initial PCM situation
12	! I_h Read the PEMstart
13	! I_i Compute orbit criterion
14
15	! II Run
16	! II_a update pressure, ice and tracers
17	! II_b Evolution of the ice
18	! II_c CO2 glaciers flows
19	! II_d Update surface and soil temperatures
20	! II_e Update the tendencies
21	! II_f Checking the stopping criterion
22
23	! III Output
24	! III_a Update surface value for the PCM start files
25	! III_b Write start and starfi.nc
26	! III_c Write start_pem
27
28	!------------------------
29
30	PROGRAM pem
31
32	!module needed for INITIALISATION
33	use phyetat0_mod, only: phyetat0
34	use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h,inertiedat, mlayer,volcapa
35	use surfdat_h, only: tsurf, co2ice, emis,&
36	& qsurf,watercap, ini_surfdat_h, &
37	albedodat, zmea, zstd, zsig, zgam, zthe, &
38	hmons, summit, base,albedo_h2o_frost, &
39	frost_albedo_threshold,emissiv
40	use dimradmars_mod, only: totcloudfrac, albedo
41	use turb_mod, only: q2, wstar
42	use dust_param_mod, only: tauscaling
43	use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, &
44	nf90_get_var, nf90_inq_varid, nf90_inq_dimid, &
45	nf90_inquire_dimension,nf90_close
46	use phyredem, only: physdem0, physdem1
47	use tracer_mod, only: noms,igcm_h2o_ice ! tracer names
48
49	! For phyredem :
50	USE control_mod, ONLY: iphysiq, day_step,nsplit_phys
51	use time_phylmdz_mod, only: daysec
52	use mod_phys_lmdz_para, only: is_parallel, is_sequential, &
53	is_mpi_root, is_omp_root, &
54	is_master
55	use time_phylmdz_mod, only: dtphys
56	USE mod_const_mpi, ONLY: COMM_LMDZ
57	USE comconst_mod, ONLY: rad,g,r,cpp
58	USE logic_mod, ONLY: iflag_phys
59	USE iniphysiq_mod, ONLY: iniphysiq
60	USE infotrac
61	USE comcstfi_h, only: pi
62
63	USE comslope_mod, ONLY: nslope,def_slope,def_slope_mean, &
64	subslope_dist,co2ice_slope, &
65	tsurf_slope,tsoil_slope,fluxgrd_slope,&
66	fluxrad_sky_slope,sky_slope,callsubslope,&
67	co2iceflow, beta_slope, capcal_slope,&
68	albedo_slope,emiss_slope,qsurf_slope,&
69	iflat
70
71	USE geometry_mod, only: latitude_deg
72
73	use pemredem, only: pemdem1
74
75	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SOIL
76	USE soil_evolution_mod, ONLY: soil_pem, soil_pem_CN
77	use comsoil_h_PEM, only: ini_comsoil_h_PEM,end_comsoil_h_PEM,nsoilmx_PEM, &
78	TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM, & ! soil thermal inertia
79	tsoil_PEM, mlayer_PEM,layer_PEM, & !number of subsurface layers, soil mid layer depths
80	fluxgeo, co2_adsorbded_phys ! geothermal flux, mass of co2 in the regolith
81	use adsorption_mod, only : regolith_co2adsorption
82
83	!!! For orbit parameters
84	USE temps_mod_evol, ONLY: dt_pem, evol_orbit_pem, Max_iter_pem
85	use planete_h, only: aphelie, periheli, year_day, peri_day, &
86	obliquit
87	use orbit_param_criterion_mod, only : orbit_param_criterion
88	use recomp_orb_param_mod, only: recomp_orb_param
89
90
91	IMPLICIT NONE
92
93	include "dimensions.h"
94	include "paramet.h"
95
96	INTEGER ngridmx
97	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
98
99	include "comdissnew.h"
100	include "comgeom.h"
101	include "iniprint.h"
102	! Same variable's name as in the GCM
103
104	INTEGER :: ngrid !Number of physical grid points
105	INTEGER :: nlayer !Number of vertical layer
106	INTEGER :: nq !Number of tracer
107	INTEGER :: day_ini !First day of the simulation
108	REAL :: pday !Physical day
109	REAL :: time_phys !Same as GCM
110	REAL :: ptimestep !Same as GCM
111	REAL :: ztime_fin !Same as GCM
112
113	! Variable for reading start.nc
114	character (len = *), parameter :: FILE_NAME_start = "start_evol.nc" !Name of the file used for initialsing the PEM
115	! variables dynamiques
116	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
117	REAL teta(ip1jmp1,llm) ! temperature potentielle
118	REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
119	REAL ps(ip1jmp1) ! pression au sol
120	REAL, dimension(:),allocatable :: ps_phys !(ngrid) ! pression au sol
121	REAL, dimension(:,:),allocatable :: ps_phys_timeseries !(ngrid x timelen) ! pression au sol instantannées
122	REAL, dimension(:,:),allocatable :: ps_phys_timeseries_yr1 !(ngrid x timelen) ! pression au sol instantannées for the first year of the gcm
123
124	REAL masse(ip1jmp1,llm) ! masse d'air
125	REAL phis(ip1jmp1) ! geopotentiel au sol
126	REAL time_0
127
128	! Variable for reading starfi.nc
129
130	character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
131	integer :: ncid, varid,status !Variable for handling opening of files
132	integer :: phydimid, subdimid, nlayerdimid, nqdimid !Variable ID for Netcdf files
133	integer :: lonvarid, latvarid, areavarid,sdvarid !Variable ID for Netcdf files
134	integer :: apvarid,bpvarid !Variable ID for Netcdf files
135
136	! Variable for reading starfi.nc and writting restartfi.nc
137
138	REAL, dimension(:),allocatable :: longitude !Longitude read in FILE_NAME and written in restartfi
139	REAL, dimension(:),allocatable :: latitude !Latitude read in FILE_NAME and written in restartfi
140	REAL, dimension(:),allocatable :: ap !Coefficient ap read in FILE_NAME_start and written in restart
141	REAL, dimension(:),allocatable :: bp !Coefficient bp read in FILE_NAME_start and written in restart
142	REAL, dimension(:),allocatable :: cell_area !Cell_area read in FILE_NAME and written in restartfi
143	REAL :: Total_surface !Total surface of the planet
144
145	! Variable for h2o_ice evolution
146
147	REAL , dimension(:,:), allocatable :: tendencies_h2o_ice ! LON x LAT field : Tendency of evolution of perenial ice over a year
148	REAL, dimension(:),allocatable :: tendencies_h2o_ice_phys ! physical point field : Tendency of evolution of perenial ice over a year
149
150	REAL , dimension(:,:), allocatable :: tendencies_co2_ice ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
151	REAL, dimension(:),allocatable :: tendencies_co2_ice_phys ! physical point field : Tendency of evolution of perenial co2 ice over a year
152
153	REAL :: ini_surf ! Initial surface of sublimating water ice
154	REAL :: ini_surf_h2o ! Initial surface of sublimating water ice
155	REAL, dimension(:),allocatable :: initial_h2o_ice ! physical point field : Logical array indicating sublimating point
156
157	REAL :: ini_surf_co2 ! Initial surface of sublimating co2 ice
158	REAL, dimension(:),allocatable :: initial_co2_ice ! physical point field : Logical array indicating sublimating point of co2 ice
159
160	REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1 ! LON x LAT field : minimum of water ice at each point for the first year
161	REAL , dimension(:,:), allocatable :: min_h2o_ice_s_2 ! LON x LAT field : minimum of water ice at each point for the second year
162
163	REAL , dimension(:,:), allocatable :: min_co2_ice_s_1 ! LON x LAT field : minimum of water ice at each point for the first year
164	REAL , dimension(:,:), allocatable :: min_co2_ice_s_2 ! LON x LAT field : minimum of water ice at each point for the second year
165
166	REAL :: global_ave_press_GCM
167	REAL :: global_ave_press_old ! physical point field : Global average pressure of initial/previous time step
168	REAL :: global_ave_press_new ! physical point field : Global average pressure of current time step
169
170	REAL , dimension(:,:), allocatable :: zplev_new
171	REAL , dimension(:,:), allocatable :: zplev_gcm
172	REAL , dimension(:,:,:), allocatable :: zplev_new_timeseries ! same but with the time series
173	REAL , dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series
174
175	LOGICAL :: STOPPING_water ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
176	LOGICAL :: STOPPING_1_water ! Logical : is there still water ice to sublimate?
177	LOGICAL :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached?
178	LOGICAL :: STOPPING_1_co2 ! Logical : is there still water ice to sublimate?
179
180	REAL, dimension(:,:,:),allocatable :: q_co2_GCM ! Initial amount of co2 in the first layer
181
182	real,save :: m_co2, m_noco2, A , B, mmean
183	real ,allocatable :: vmr_co2_gcm_phys(:,:) !(ngrid) ! co2 volume mixing ratio
184	real ,allocatable :: vmr_co2_pem_phys(:,:) !(ngrid) ! co2 volume mixing ratio
185	real ,allocatable :: q_h2o_GCM_phys(:,:)
186	real ,allocatable :: q_co2_GCM_phys(:,:)
187	real ,allocatable :: q_co2_PEM_phys(:,:)
188	REAL, ALLOCATABLE :: ps_GCM(:,:,:)
189	REAL, ALLOCATABLE :: ps_GCM_yr1(:,:,:)
190	REAL, ALLOCATABLE :: vmr_co2_gcm(:,:,:)
191	REAL, ALLOCATABLE :: q_h2o_GCM(:,:,:)
192	REAL ,allocatable :: q_h2o_PEM_phys(:,:)
193	integer :: timelen
194	REAL :: ave
195
196	REAL, ALLOCATABLE :: p(:,:) !(ngrid,llmp1)
197	REAL :: extra_mass ! Extra mass of a tracer if it is greater than 1
198
199	!!!!!!!!!!!!!!!!!!!!!!!! SLOPE
200	character*2 :: str2
201	REAL ,allocatable :: watercap_slope(:,:) !(ngrid,nslope)
202	REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_1 ! LON x LAT field : minimum of water ice at each point for the first year
203	REAL , dimension(:,:,:), allocatable :: min_co2_ice_slope_2 ! LON x LAT field : minimum of water ice at each point for the second year
204	REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_1 ! LON x LAT field : minimum of water ice at each point for the first year
205	REAL , dimension(:,:,:), allocatable :: min_h2o_ice_slope_2 ! LON x LAT field : minimum of water ice at each point for the second year
206
207
208	REAL , dimension(:,:,:,:), allocatable :: co2_ice_GCM_slope ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
209	REAL , dimension(:,:,:), allocatable :: co2_ice_GCM_phys_slope ! LON x LATX NSLOPE x Times field : co2 ice given by the GCM
210
211
212	REAL, dimension(:,:),allocatable :: initial_co2_ice_sublim_slope ! physical point field : Logical array indicating sublimating point of co2 ice
213	REAL, dimension(:,:),allocatable :: initial_h2o_ice_slope ! physical point field : Logical array indicating sublimating point of h2o ice
214	REAL, dimension(:,:),allocatable :: initial_co2_ice_slope ! physical point field : Logical array indicating sublimating point of co2 ice
215
216	REAL , dimension(:,:,:), allocatable :: tendencies_co2_ice_slope ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
217	REAL , dimension(:,:,:), allocatable :: tendencies_h2o_ice_slope ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
218	REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope ! physical point field : Tendency of evolution of perenial co2 ice over a year
219	REAL, dimension(:,:),allocatable :: tendencies_co2_ice_phys_slope_ini ! physical point field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
220	REAL, dimension(:,:),allocatable :: tendencies_h2o_ice_phys_slope ! physical point field : Tendency of evolution of perenial co2 ice over a year
221
222	REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: co2_hmax ! Maximum height for CO2 deposit on slopes (m)
223
224	REAL, PARAMETER :: rho_co2 = 1600 ! CO2 ice density (kg/m^3)
225	INTEGER :: iaval ! Index of the neighboord slope ()
226	REAL , dimension(:,:), allocatable :: flag_co2flow(:,:) !(ngrid,nslope) ! To flag where there is a glacier flow
227	REAL , dimension(:), allocatable :: flag_co2flow_mesh(:) !(ngrid) ! To flag where there is a glacier flow
228
229	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SURFACE/SOIL
230
231	REAL, ALLOCATABLE :: tsurf_ave(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature [K]
232	REAL, ALLOCATABLE :: tsurf_ave_phys(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature [K]
233	REAL, ALLOCATABLE :: tsoil_ave(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
234	REAL, ALLOCATABLE :: tsoil_ave_yr1(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Soil Temperature [K]
235	REAL, ALLOCATABLE :: tsoil_ave_phys_yr1(:,:,:) !IG x SLOPE field : Averaged Soil Temperature [K]
236
237	REAL, ALLOCATABLE :: TI_GCM(:,:,:,:) ! LON x LAT x SLOPE field : Averaged Thermal Inertia [SI]
238	REAL, ALLOCATABLE :: tsurf_GCM_timeseries(:,:,:,:) !LONX LAT x SLOPE XTULES field : NOn averaged Surf Temperature [K]
239	REAL, ALLOCATABLE :: tsurf_phys_GCM_timeseries(:,:,:) !IG x SLOPE XTULES field : NOn averaged Surf Temperature [K]
240
241	REAL, ALLOCATABLE :: tsoil_phys_PEM_timeseries(:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
242	REAL, ALLOCATABLE :: tsoil_GCM_timeseries(:,:,:,:,:) !IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
243
244	REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:,:) ! LON x LAT x SLOPE field : Averaged Surface Temperature of the first year of the gcm [K]
245	REAL, ALLOCATABLE :: tsurf_ave_phys_yr1(:,:) ! IG x LAT x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
246
247	REAL, PARAMETER :: year_step = 1
248	INTEGER :: year_iter !Counter for the number of PEM iteration
249	INTEGER :: year_iter_max
250	REAL :: timestep
251
252	REAL, ALLOCATABLE :: inertiesoil(:,:) ! Thermal inertia of the mesh for restart
253
254	REAL, ALLOCATABLE :: TI_GCM_phys(:,:,:) ! Averaged GCM Thermal Inertia [SI]
255	REAL, ALLOCATABLE :: TI_GCM_start(:,:,:) ! Averaged GCM Thermal Inertia [SI]
256
257	REAL,ALLOCATABLE :: ice_depth(:,:)
258	REAL,ALLOCATABLE :: TI_locslope(:,:)
259	REAL,ALLOCATABLE :: Tsoil_locslope(:,:)
260	REAL,ALLOCATABLE :: Tsurf_locslope(:)
261	REAL,ALLOCATABLE :: alph_locslope(:,:)
262	REAL,ALLOCATABLE :: beta_locslope(:,:)
263
264	REAL,ALLOCATABLE :: Tsurfave_before_saved(:,:)
265	REAL, ALLOCATABLE :: delta_co2_adsorbded(:)
266	REAL :: totmass_adsorbded
267
268	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
269
270	! Loop variable
271	LOGICAL :: bool_sublim
272	INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop
273	REAL :: beta = 3182.48
274	REAL :: alpha = 23.3494
275
276	! Parallel variables
277	is_sequential=.true.
278	is_parallel=.false.
279	is_mpi_root=.true.
280	is_omp_root=.true.
281	is_master=.true.
282
283	day_ini=0 !test
284	time_phys=0. !test
285
286	! Some constants
287
288	ngrid=ngridmx
289	nlayer=llm
290
291	m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
292	m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol)
293	A =(1/m_co2 - 1/m_noco2)
294	B=1/m_noco2
295
296	year_day=669
297	daysec=88775.
298	dtphys=0
299	timestep=year_day*daysec/year_step
300
301	!------------------------
302
303	! I Initialisation
304	! I_a READ run.def
305
306	!------------------------
307
308	!----------------------------READ run.def ---------------------
309	CALL conf_gcm( 99, .TRUE. )
310
311	call infotrac_init
312	nq=nqtot
313
314	!------------------------
315
316	! I Initialisation
317	! I_a READ run.def
318	! I_b READ of start_evol.nc and starfi_evol.nc
319
320	!------------------------
321
322	!----------------------------Initialisation : READ some constant of startfi_evol.nc ---------------------
323
324	! In the gcm, these values are given to the physic by the dynamic.
325	! Here we simply read them in the startfi_evol.nc file
326	status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
327
328	allocate(longitude(ngrid))
329	allocate(latitude(ngrid))
330	allocate(cell_area(ngrid))
331
332	status = nf90_inq_varid(ncid, "longitude", lonvarid)
333	status = nf90_get_var(ncid, lonvarid, longitude)
334
335	status = nf90_inq_varid(ncid, "latitude", latvarid)
336	status = nf90_get_var(ncid, latvarid, latitude)
337
338	status = nf90_inq_varid(ncid, "area", areavarid)
339	status = nf90_get_var(ncid, areavarid, cell_area)
340
341	call ini_comsoil_h(ngrid)
342
343	status = nf90_inq_varid(ncid, "soildepth", sdvarid)
344	status = nf90_get_var(ncid, sdvarid, mlayer)
345
346	status =nf90_close(ncid)
347
348	!----------------------------READ start.nc ---------------------
349
350	allocate(q(ip1jmp1,llm,nqtot))
351	CALL dynetat0(FILE_NAME_start,vcov,ucov, &
352	teta,q,masse,ps,phis, time_0)
353
354	CALL iniconst !new
355	CALL inigeom
356	allocate(ap(nlayer+1))
357	allocate(bp(nlayer+1))
358	status =nf90_open(FILE_NAME_start, NF90_NOWRITE, ncid)
359	status = nf90_inq_varid(ncid, "ap", apvarid)
360	status = nf90_get_var(ncid, apvarid, ap)
361	status = nf90_inq_varid(ncid, "bp", bpvarid)
362	status = nf90_get_var(ncid, bpvarid, bp)
363	status =nf90_close(ncid)
364
365	CALL iniphysiq(iim,jjm,llm, &
366	(jjm-1)*iim+2,comm_lmdz, &
367	daysec,day_ini,dtphys/nsplit_phys, &
368	rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
369	iflag_phys)
370
371	DO nnq=1,nqtot
372	if(noms(nnq).eq."h2o_ice") igcm_h2o_ice = nnq
373	ENDDO
374
375	!----------------------------READ startfi.nc ---------------------
376
377	! First we read the initial state (starfi.nc)
378
379	allocate(watercap_slope(ngrid,nslope))
380	allocate(TI_GCM_start(ngrid,nsoilmx,nslope))
381	allocate(inertiesoil(ngrid,nsoilmx))
382
383	CALL phyetat0 (FILE_NAME,0,0, &
384	nsoilmx,ngrid,nlayer,nq, &
385	day_ini,time_phys, &
386	tsurf,tsoil,albedo,emis, &
387	q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar, &
388	watercap,inertiesoil,nslope,tsurf_slope, &
389	tsoil_slope,co2ice_slope,def_slope,def_slope_mean, &
390	subslope_dist,albedo_slope,emiss_slope, TI_GCM_start, &
391	qsurf_slope,watercap_slope)
392
393	if(soil_pem) then
394	deallocate(TI_GCM_start) !not used then
395	endif
396
397	! Remove unphysical values of surface tracer
398	DO i=1,ngrid
399	DO nnq=1,nqtot
400	DO islope=1,nslope
401	if(qsurf_slope(i,nnq,islope).LT.0) then
402	qsurf_slope(i,nnq,islope)=0.
403	endif
404	enddo
405	enddo
406	enddo
407
408	!------------------------
409
410	! I Initialisation
411	! I_a READ run.def
412	! I_b READ of start_evol.nc and starfi_evol.nc
413	! I_c Subslope parametrisation
414
415	!------------------------
416
417	!----------------------------Subslope parametrisation definition ---------------------
418
419	! Define some slope statistics
420	iflat=1
421	DO islope=2,nslope
422	IF(abs(def_slope_mean(islope)).lt. &
423	abs(def_slope_mean(iflat))) THEN
424	iflat = islope
425	ENDIF
426	ENDDO
427
428	PRINT*,'Flat slope for islope = ',iflat
429	PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
430
431	! CO2 max thickness (for glaciers flows)
432	allocate(co2_hmax(nslope))
433	if(nslope.eq.7) then ! ugly way to implement that ...
434	! CF documentation that explain how the values are computed
435	co2_hmax(1) = 1.5
436	co2_hmax(7) = co2_hmax(1)
437	co2_hmax(2) = 2.4
438	co2_hmax(6) = co2_hmax(2)
439	co2_hmax(3) = 5.6
440	co2_hmax(5) = co2_hmax(3)
441	co2_hmax(4) = 1000000.
442	endif
443
444	allocate(flag_co2flow(ngrid,nslope))
445	allocate(flag_co2flow_mesh(ngrid))
446
447	flag_co2flow(:,:) = 0.
448	flag_co2flow_mesh(:) = 0.
449
450	!---------------------------- READ GCM data ---------------------
451
452	! I Initialisation
453	! I_a READ run.def
454	! I_b READ of start_evol.nc and starfi_evol.nc
455	! I_c Subslope parametrisation
456	! I_d READ GCM data and convert to the physical grid
457
458	!------------------------
459
460	! 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
461
462	call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
463
464	allocate(min_h2o_ice_s_1(iim+1,jjm+1))
465	allocate(min_co2_ice_s_1(iim+1,jjm+1))
466	allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
467	allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
468	allocate(q_co2_GCM(iim+1,jjm+1,timelen))
469	allocate(ps_GCM(iim+1,jjm+1,timelen))
470	allocate(ps_GCM_yr1(iim+1,jjm+1,timelen))
471	allocate(min_co2_ice_slope_1(iim+1,jjm+1,nslope))
472	allocate(min_h2o_ice_slope_1(iim+1,jjm+1,nslope))
473	allocate(tsurf_ave(iim+1,jjm+1,nslope))
474	allocate(tsurf_ave_yr1(iim+1,jjm+1,nslope))
475	allocate(tsurf_ave_phys(ngrid,nslope))
476	allocate(tsurf_ave_phys_yr1(ngrid,nslope))
477	allocate(tsurf_GCM_timeseries(iim+1,jjm+1,nslope,timelen))
478	allocate(tsurf_phys_GCM_timeseries(ngrid,nslope,timelen))
479	allocate(co2_ice_GCM_phys_slope(ngrid,nslope,timelen))
480	allocate(co2_ice_GCM_slope(iim+1,jjm+1,nslope,timelen))
481	allocate(Tsurfave_before_saved(ngrid,nslope))
482	allocate(tsoil_ave(iim+1,jjm+1,nsoilmx,nslope))
483	allocate(tsoil_ave_yr1(iim+1,jjm+1,nsoilmx,nslope))
484	allocate(tsoil_ave_phys_yr1(ngrid,nsoilmx_PEM,nslope))
485	allocate(TI_GCM(iim+1,jjm+1,nsoilmx,nslope))
486	allocate(tsoil_GCM_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen))
487	allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
488	allocate(delta_co2_adsorbded(ngrid))
489
490	print *, "Downloading data Y1..."
491
492	call read_data_GCM("data_GCM_Y1.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_yr1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,&
493	nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
494
495	! 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
496
497	print *, "Downloading data Y1 done"
498
499	allocate(min_h2o_ice_s_2(iim+1,jjm+1))
500	allocate(min_co2_ice_s_2(iim+1,jjm+1))
501	allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope))
502	allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
503
504	print *, "Downloading data Y2"
505
506	call read_data_GCM("data_GCM_Y2.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
507	nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope)
508
509	print *, "Downloading data Y2 done"
510
511	! The variables in the dynamic grid are transfered to the physical grid
512
513	allocate(vmr_co2_gcm_phys(ngrid,timelen))
514	allocate(vmr_co2_pem_phys(ngrid,timelen))
515	allocate(q_h2o_GCM_phys(ngrid,timelen))
516	allocate(q_h2o_PEM_phys(ngrid,timelen))
517	allocate(q_co2_GCM_phys(ngrid,timelen))
518	allocate(q_co2_PEM_phys(ngrid,timelen))
519	allocate(ps_phys(ngrid))
520	allocate(ps_phys_timeseries(ngrid,timelen))
521	allocate(ps_phys_timeseries_yr1(ngrid,timelen))
522
523	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,vmr_co2_gcm,vmr_co2_gcm_phys)
524	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_h2o_GCM,q_h2o_GCM_phys)
525	CALL gr_dyn_fi(timelen,iip1,jjp1,ngridmx,q_co2_GCM,q_co2_GCM_phys)
526	call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_phys)
527	call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM,ps_phys_timeseries)
528	call gr_dyn_fi(timelen,iip1,jjp1,ngridmx,ps_GCM_yr1,ps_phys_timeseries_yr1)
529	CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave,tsurf_ave_phys)
530	CALL gr_dyn_fi(nslope,iip1,jjp1,ngridmx,tsurf_ave_yr1,tsurf_ave_phys_yr1)
531
532	deallocate(vmr_co2_gcm)
533	deallocate(q_h2o_GCM)
534	deallocate(q_co2_GCM)
535	deallocate(ps_GCM)
536	deallocate(ps_GCM_yr1)
537	deallocate(tsurf_ave)
538	deallocate(tsurf_ave_yr1)
539
540	q_co2_PEM_phys(:,:)= q_co2_GCM_phys(:,:)
541	q_h2o_PEM_phys(:,:)= q_h2o_GCM_phys(:,:)
542
543	!------------------------
544
545	! I Initialisation
546	! I_a READ run.def
547	! I_b READ of start_evol.nc and starfi_evol.nc
548	! I_c Subslope parametrisation
549	! I_d READ GCM data and convert to the physical grid
550	! I_e Initialisation of the PEM variable and soil
551
552	!------------------------
553
554	!---------------------------- Initialisation of the PEM soil and values ---------------------
555
556	call end_comsoil_h_PEM
557	call ini_comsoil_h_PEM(ngrid,nslope)
558
559	allocate(ice_depth(ngrid,nslope))
560	allocate(TI_GCM_phys(ngrid,nsoilmx,nslope))
561
562	DO islope = 1,nslope
563	if(soil_pem) then
564	CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,TI_GCM(:,:,:,islope),TI_GCM_phys(:,:,islope))
565	endif !soil_pem
566	DO t=1,timelen
567	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsurf_GCM_timeseries(:,:,islope,t),tsurf_phys_GCM_timeseries(:,islope,t))
568	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,co2_ice_GCM_slope(:,:,islope,t),co2_ice_GCM_phys_slope(:,islope,t))
569	enddo
570	ENDDO
571
572	deallocate(co2_ice_GCM_slope)
573	deallocate(TI_GCM)
574
575	if(soil_pem) then
576
577	deallocate(tsurf_GCM_timeseries)
578	call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM)
579
580	DO islope = 1,nslope
581	DO l=1,nsoilmx
582	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope))
583	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,l,islope),tsoil_PEM(:,l,islope))
584	DO t=1,timelen
585	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t))
586	ENDDO
587	ENDDO
588	DO l=nsoilmx+1,nsoilmx_PEM
589	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope))
590	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope))
591	ENDDO
592	ENDDO
593
594	deallocate(tsoil_ave_yr1)
595	deallocate(tsoil_ave)
596	deallocate(tsoil_GCM_timeseries)
597	endif !soil_pem
598
599	!------------------------
600
601	! I Initialisation
602	! I_a READ run.def
603	! I_b READ of start_evol.nc and starfi_evol.nc
604	! I_c Subslope parametrisation
605	! I_d READ GCM data and convert to the physical grid
606	! I_e Initialisation of the PEM variable and soil
607	! I_f Compute tendencies & Save initial situation
608
609	!----- Compute tendencies from the PCM run
610
611	allocate(tendencies_h2o_ice(iim+1,jjm+1))
612	allocate(tendencies_h2o_ice_phys(ngrid))
613	allocate(tendencies_co2_ice(iim+1,jjm+1))
614	allocate(tendencies_co2_ice_phys(ngrid))
615	allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
616	allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
617	allocate(tendencies_co2_ice_phys_slope_ini(ngrid,nslope))
618	allocate(tendencies_h2o_ice_slope(iim+1,jjm+1,nslope))
619	allocate(tendencies_h2o_ice_phys_slope(ngrid,nslope))
620
621	! Compute the tendencies of the evolution of ice over the years
622
623	call compute_tendencies(tendencies_h2o_ice,min_h2o_ice_s_1,&
624	min_h2o_ice_s_2,iim,jjm,ngrid,tendencies_h2o_ice_phys)
625
626	call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,&
627	min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys)
628
629	call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,&
630	min_co2_ice_slope_2,iim,jjm,ngrid,tendencies_co2_ice_phys_slope,nslope)
631
632	tendencies_co2_ice_phys_slope_ini(:,:)=tendencies_co2_ice_phys_slope(:,:)
633
634	call compute_tendencies_slope(tendencies_h2o_ice_slope,min_h2o_ice_slope_1,&
635	min_h2o_ice_slope_2,iim,jjm,ngrid,tendencies_h2o_ice_phys_slope,nslope)
636
637	!------------------------
638
639	! I Initialisation
640	! I_a READ run.def
641	! I_b READ of start_evol.nc and starfi_evol.nc
642	! I_c Subslope parametrisation
643	! I_d READ GCM data and convert to the physical grid
644	! I_e Initialisation of the PEM variable and soil
645	! I_f Compute tendencies & Save initial situation
646	! I_g Save initial PCM situation
647
648	!---------------------------- Save initial PCM situation ---------------------
649
650	allocate(initial_h2o_ice(ngrid))
651	allocate(initial_co2_ice(ngrid))
652	allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
653	allocate(initial_co2_ice_slope(ngrid,nslope))
654	allocate(initial_h2o_ice_slope(ngrid,nslope))
655
656	! We save the places where water ice is sublimating
657	! We compute the surface of water ice sublimating
658	ini_surf=0.
659	ini_surf_co2=0.
660	ini_surf_h2o=0.
661	Total_surface=0.
662	do i=1,ngrid
663	Total_surface=Total_surface+cell_area(i)
664	if (tendencies_h2o_ice_phys(i).LT.0) then
665	initial_h2o_ice(i)=1.
666	ini_surf=ini_surf+cell_area(i)
667	else
668	initial_h2o_ice(i)=0.
669	endif
670	do islope=1,nslope
671	if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
672	initial_co2_ice_sublim_slope(i,islope)=1.
673	ini_surf_co2=ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
674	else
675	initial_co2_ice_sublim_slope(i,islope)=0.
676	endif
677	if (co2ice_slope(i,islope).GT.0) then
678	initial_co2_ice_slope(i,islope)=1.
679	else
680	initial_co2_ice_slope(i,islope)=0.
681	endif
682	if (tendencies_h2o_ice_phys_slope(i,islope).LT.0) then
683	initial_h2o_ice_slope(i,islope)=1.
684	ini_surf_h2o=ini_surf_h2o+cell_area(i)*subslope_dist(i,islope)
685	else
686	initial_h2o_ice_slope(i,islope)=0.
687	endif
688	enddo
689	enddo
690
691	print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
692	print *, "Total initial surface of h2o ice sublimating=", ini_surf
693	print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
694	print *, "Total surface of the planet=", Total_surface
695
696	allocate(zplev_gcm(ngrid,nlayer+1))
697
698	DO l=1,nlayer+1
699	DO ig=1,ngrid
700	zplev_gcm(ig,l) = ap(l) + bp(l)*ps_phys(ig)
701	ENDDO
702	ENDDO
703
704	global_ave_press_old=0.
705	do i=1,ngrid
706	global_ave_press_old=global_ave_press_old+cell_area(i)*ps_phys(i)/Total_surface
707	enddo
708
709	global_ave_press_GCM=global_ave_press_old
710	global_ave_press_new=global_ave_press_old
711	print *, "Initial global average pressure=", global_ave_press_GCM
712
713	!------------------------
714
715	! I Initialisation
716	! I_a READ run.def
717	! I_b READ of start_evol.nc and starfi_evol.nc
718	! I_c Subslope parametrisation
719	! I_d READ GCM data and convert to the physical grid
720	! I_e Initialisation of the PEM variable and soil
721	! I_f Compute tendencies & Save initial situation
722	! I_g Save initial PCM situation
723	! I_h Read the PEMstart
724
725	!---------------------------- Read the PEMstart ---------------------
726
727	call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, &
728	tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,&
729	tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded)
730
731
732	if(soil_pem) then
733	totmass_adsorbded = 0.
734
735	do ig = 1,ngrid
736	do islope =1, nslope
737	do l = 1,nsoilmx_PEM - 1
738	totmass_adsorbded = totmass_adsorbded + co2_adsorbded_phys(ig,l,islope)(layer_PEM(l+1) - layer_PEM(l)) &
739	subslope_dist(ig,islope)/cos(pidef_slope_mean(islope)/180.) &
740	cell_area(ig)
741	enddo
742	enddo
743	enddo
744
745	write(,) "Tot mass in the regolith=", totmass_adsorbded
746	deallocate(tsoil_ave_phys_yr1)
747	endif !soil_pem
748	deallocate(tsurf_ave_phys_yr1)
749	deallocate(ps_phys_timeseries_yr1)
750
751	!------------------------
752
753	! I Initialisation
754	! I_a READ run.def
755	! I_b READ of start_evol.nc and starfi_evol.nc
756	! I_c Subslope parametrisation
757	! I_d READ GCM data and convert to the physical grid
758	! I_e Initialisation of the PEM variable and soil
759	! I_f Compute tendencies & Save initial situation
760	! I_g Save initial PCM situation
761	! I_h Read the PEMstar
762	! I_i Compute orbit criterion
763
764	CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
765
766	if(evol_orbit_pem) then
767	call orbit_param_criterion(year_iter_max)
768	else
769	year_iter_max=Max_iter_pem
770	endif
771
772	!--------------------------- END INITIALISATION ---------------------
773
774	!---------------------------- RUN ---------------------
775
776	!------------------------
777
778	! II Run
779	! II_a update pressure,ice and tracers
780
781	!------------------------
782	year_iter=0
783	print *, "Max_iter_pem", Max_iter_pem
784	do while (year_iter.LT.year_iter_max)
785
786	! II.a.1. Compute updated global pressure
787	print *, "Recomputing the new pressure..."
788	do i=1,ngrid
789	do islope=1,nslope
790	global_ave_press_new=global_ave_press_new-gcell_area(i)tendencies_co2_ice_phys_slope(i,islope)subslope_dist(i,islope)/cos(pidef_slope_mean(islope)/180.)/Total_surface
791	enddo
792	if(soil_pem) then
793	global_ave_press_new = global_ave_press_new -gcell_area(i)delta_co2_adsorbded(i)/Total_surface
794	endif
795	enddo
796
797	! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
798	allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
799	print *, "Recomputing the old pressure levels timeserie adapted to the old pressure..."
800	DO l=1,nlayer+1
801	DO ig=1,ngrid
802	zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:)
803	ENDDO
804	ENDDO
805
806	! II.a.3. Surface pressure timeseries
807	print *, "Recomputing the surface pressure timeserie adapted to the new pressure..."
808	do i = 1,ngrid
809	ps_phys_timeseries(i,:) = ps_phys_timeseries(i,:)*global_ave_press_new/global_ave_press_old
810	enddo
811
812	! II.a.4. New pressure levels timeseries
813	allocate(zplev_new_timeseries(ngrid,nlayer+1,timelen))
814	print *, "Recomputing the new pressure levels timeserie adapted to the new pressure..."
815	do l=1,nlayer+1
816	do ig=1,ngrid
817	zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_phys_timeseries(ig,:)
818	enddo
819	enddo
820
821	! II.a.5. Tracers timeseries
822	print *, "Recomputing of tracer VMR timeseries for the new pressure..."
823
824	l=1
825	DO ig=1,ngrid
826	DO t = 1, timelen
827	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))
828	if(q_h2o_PEM_phys(ig,t).LT.0) then
829	q_h2o_PEM_phys(ig,t)=1E-30
830	endif
831	if(q_h2o_PEM_phys(ig,t).GT.1) then
832	q_h2o_PEM_phys(ig,t)=1.
833	endif
834	enddo
835	enddo
836
837	DO ig=1,ngrid
838	DO t = 1, timelen
839	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)) &
840	+ ( (zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t)) - &
841	(zplev_old_timeseries(ig,l,t)-zplev_old_timeseries(ig,l+1,t)) ) / &
842	(zplev_new_timeseries(ig,l,t)-zplev_new_timeseries(ig,l+1,t))
843	if (q_co2_PEM_phys(ig,t).LT.0) then
844	q_co2_PEM_phys(ig,t)=1E-30
845	elseif (q_co2_PEM_phys(ig,t).GT.1) then
846	q_co2_PEM_phys(ig,t)=1.
847	endif
848	mmean=1/(A*q_co2_PEM_phys(ig,t) +B)
849	vmr_co2_pem_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2
850	ENDDO
851	ENDDO
852
853	deallocate(zplev_new_timeseries)
854	deallocate(zplev_old_timeseries)
855
856	! II Run
857	! II_a update pressure, ice and tracers
858	! II_b Evolution of the ice
859
860	! II.b. Evolution of the ice
861	print *, "Evolution of h2o ice"
862	call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
863
864	print *, "Evolution of co2 ice"
865	call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
866
867	!------------------------
868
869	! II Run
870	! II_a update pressure, ice and tracers
871	! II_b Evolution of the ice
872	! II_c CO2 glaciers flows
873
874	!------------------------
875
876	print *, "Co2 glacier flow"
877	DO ig = 1,ngrid
878	DO islope = 1,nslope
879	IF(islope.ne.iflat) THEN ! ice can be infinite on flat ground
880	! First: check that CO2 ice must flow (excess of ice on the slope), ice can accumulate on flat ground
881	IF(co2ice_slope(ig,islope).ge.rho_co2co2_hmax(islope) &
882	cos(pi*def_slope_mean(islope)/180.)) THEN
883	! Second: determine the flatest slopes possible:
884	IF(islope.gt.iflat) THEN
885	iaval=islope-1
886	ELSE
887	iaval=islope+1
888	ENDIF
889	do while ((iaval.ne.iflat).and. &
890	(subslope_dist(ig,iaval).eq.0))
891	IF(iaval.gt.iflat) THEN
892	iaval=iaval-1
893	ELSE
894	iaval=iaval+1
895	ENDIF
896	enddo
897	co2ice_slope(ig,iaval) = co2ice_slope(ig,iaval) + &
898	(co2ice_slope(ig,islope) - rho_co2* co2_hmax(islope) * &
899	cos(pidef_slope_mean(islope)/180.)) &
900	subslope_dist(ig,islope)/subslope_dist(ig,iaval) * &
901	cos(pi*def_slope_mean(iaval)/180.) / &
902	cos(pi*def_slope_mean(islope)/180.)
903
904	co2ice_slope(ig,islope)=rho_co2co2_hmax(islope) &
905	cos(pi*def_slope_mean(islope)/180.)
906
907	flag_co2flow(ig,islope) = 1.
908	flag_co2flow_mesh(ig) = 1.
909	ENDIF ! co2ice > hmax
910	ENDIF ! iflattsoil_lope
911	ENDDO !islope
912	ENDDO !ig
913
914	!------------------------
915
916	! II Run
917	! II_a update pressure, ice and tracers
918	! II_b Evolution of the ice
919	! II_c CO2 glaciers flows
920	! II_d Update surface and soil temperatures
921
922	!------------------------
923
924	! II_d.1 Update Tsurf
925
926	print *, "Updating the new Tsurf"
927	bool_sublim=0
928	Tsurfave_before_saved(:,:) = tsurf_ave_phys(:,:)
929	DO ig = 1,ngrid
930	DO islope = 1,nslope
931	if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-5) THEN !co2ice disappeared, look for closest point without co2ice
932	if(latitude_deg(ig).gt.0) then
933	do ig_loop=ig,ngrid
934	DO islope_loop=islope,iflat,-1
935	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
936	tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
937	bool_sublim=1
938	exit
939	endif
940	enddo
941	if (bool_sublim.eq.1) then
942	exit
943	endif
944	enddo
945	else
946	do ig_loop=ig,1,-1
947	DO islope_loop=islope,iflat
948	if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
949	tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
950	bool_sublim=1
951	exit
952	endif
953	enddo
954	if (bool_sublim.eq.1) then
955	exit
956	endif
957	enddo
958	endif
959	initial_co2_ice_slope(ig,islope)=0
960	if ((co2ice_slope(ig,islope).lt.1e-5).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold)) then
961	albedo_slope(ig,1,islope) = albedo_h2o_frost
962	albedo_slope(ig,2,islope) = albedo_h2o_frost
963	else
964	albedo_slope(ig,1,islope) = albedodat(ig)
965	albedo_slope(ig,2,islope) = albedodat(ig)
966	emiss_slope(ig,islope) = emissiv
967	endif
968	elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2
969	ave=0.
970	do t=1,timelen
971	if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then
972	ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.))
973	else
974	ave = ave + tsurf_phys_GCM_timeseries(ig,islope,t)
975	endif
976	enddo
977	tsurf_ave_phys(ig,islope)=ave/timelen
978	endif
979	enddo
980	enddo
981
982	do t = 1,timelen
983	tsurf_phys_GCM_timeseries(:,:,t) = tsurf_phys_GCM_timeseries(:,:,t) +( tsurf_ave_phys(:,:) -Tsurfave_before_saved(:,:))
984	enddo
985	! for the start
986	do ig = 1,ngrid
987	do islope = 1,nslope
988	tsurf_slope(ig,islope) = tsurf_slope(ig,islope) - (Tsurfave_before_saved(ig,islope)-tsurf_ave_phys(ig,islope))
989	enddo
990	enddo
991
992	if(soil_pem) then
993
994	! II_d.2 Update soil temperature
995
996	allocate(TI_locslope(ngrid,nsoilmx_PEM))
997	allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
998	allocate(Tsurf_locslope(ngrid))
999	allocate(alph_locslope(ngrid,nsoilmx_PEM-1))
1000	allocate(beta_locslope(ngrid,nsoilmx_PEM-1))
1001
1002	print *,"Updating soil temperature"
1003
1004	! Soil averaged
1005	do islope = 1,nslope
1006	TI_locslope(:,:) = TI_PEM(:,:,islope)
1007	Tsoil_locslope(:,:) = tsoil_PEM(:,:,islope)
1008	Tsurf_locslope(:) = tsurf_ave_phys(:,islope)
1009	alph_locslope(:,:) = alph_PEM(:,:,islope)
1010	beta_locslope(:,:) = beta_PEM(:,:,islope)
1011	call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
1012	call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope)
1013	do t = 1,timelen
1014	tsoil_phys_PEM_timeseries(:,:,islope,t) = tsoil_phys_PEM_timeseries(:,:,islope,t) + (Tsoil_locslope(:,:)- tsoil_PEM(:,:,islope))
1015	enddo
1016	TI_PEM(:,:,islope) = TI_locslope(:,:)
1017	tsoil_PEM(:,:,islope) = Tsoil_locslope(:,:)
1018	tsurf_ave_phys(:,islope) = Tsurf_locslope(:)
1019	alph_PEM(:,:,islope) = alph_locslope(:,:)
1020	beta_PEM(:,:,islope) = beta_locslope(:,:)
1021	enddo
1022
1023	print *, "Update of soil temperature done"
1024
1025	! Check Nan in the time series
1026	do ig = 1,ngrid
1027	do islope = 1,nslope
1028	do iloop = 1,nsoilmx_PEM
1029	if(isnan(tsoil_PEM(ig,iloop,islope))) then
1030	write(,) "Problem : There is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope)
1031	stop
1032	endif
1033	enddo
1034	enddo
1035	enddo
1036
1037	write(,) "Compute ice table"
1038
1039	! II_d.3 Update the ice table
1040	call computeice_table(timelen,ngrid,nslope,nsoilmx,nsoilmx_PEM, tsoil_phys_PEM_timeseries,tsurf_phys_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
1041	ps_phys_timeseries, ice_depth)
1042
1043	print *, "Update soil propreties"
1044	! II_d.4 Update the soil thermal properties
1045	call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
1046	ice_depth,TI_PEM)
1047
1048	! II_d.5 Update the mass of the regolith adsorbded
1049	call regolith_co2adsorption(ngrid,nslope,nsoilmx_PEM,timelen,ps_phys_timeseries,tsoil_PEM,TI_PEM,tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope, &
1050	co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), q_co2_PEM_phys,q_h2o_PEM_phys,co2_adsorbded_phys,delta_co2_adsorbded)
1051
1052	endif !soil_pem
1053
1054	!------------------------
1055
1056	! II Run
1057	! II_a update pressure, ice and tracers
1058	! II_b Evolution of the ice
1059	! II_c CO2 glaciers flows
1060	! II_d Update surface and soil temperatures
1061	! II_e Update the tendencies
1062
1063	!------------------------
1064
1065	print *, "Adaptation of the new co2 tendencies to the current pressure"
1066	call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&
1067	global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
1068
1069	!------------------------
1070
1071	! II Run
1072	! II_a update pressure, ice and tracers
1073	! II_b Evolution of the ice
1074	! II_c CO2 glaciers flows
1075	! II_d Update surface and soil temperatures
1076	! II_e Update the tendencies
1077	! II_f Checking the stopping criterion
1078
1079	!------------------------
1080	call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope)
1081
1082	call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)
1083
1084	year_iter=year_iter+dt_pem
1085
1086	print *, "Checking all the stopping criterion."
1087	if (STOPPING_water) then
1088	print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water
1089	endif
1090	if (STOPPING_1_water) then
1091	print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water
1092	endif
1093	if (STOPPING_co2) then
1094	print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2
1095	endif
1096	if (STOPPING_1_co2) then
1097	print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
1098	endif
1099
1100	if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2) then
1101	exit
1102	else
1103	print *, "We continue!"
1104	print *, "Number of iteration of the PEM : year_iter=", year_iter
1105	endif
1106
1107	global_ave_press_old=global_ave_press_new
1108
1109	enddo ! big time iteration loop of the pem
1110
1111
1112	!---------------------------- END RUN PEM ---------------------
1113
1114	!---------------------------- OUTPUT ---------------------
1115
1116	!------------------------
1117
1118	! III Output
1119	! III_a Update surface value for the PCM start files
1120
1121	!------------------------
1122
1123	! III_a.1 Ice update (for startfi)
1124	! Co2 ice
1125	DO ig = 1,ngrid
1126	co2ice(ig) = 0.
1127	DO islope = 1,nslope
1128	co2ice(ig) = co2ice(ig) + co2ice_slope(ig,islope) &
1129	* subslope_dist(ig,islope) / &
1130	cos(pi*def_slope_mean(islope)/180.)
1131	ENDDO
1132	ENDDO ! of DO ig=1,ngrid
1133	! H2o ice
1134	DO ig = 1,ngrid
1135	qsurf(ig,igcm_h2o_ice) = 0.
1136	DO islope = 1,nslope
1137	qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &
1138	* subslope_dist(ig,islope) / &
1139	cos(pi*def_slope_mean(islope)/180.)
1140	ENDDO
1141	ENDDO ! of DO ig=1,ngrid
1142
1143	! III_a.2 Tsurf and Tsoil update (for startfi)
1144
1145	if(soil_pem) then
1146	call interpolate_TIPEM_TIGCM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_PEM,TI_GCM_phys)
1147	tsoil_slope(:,:,:) = tsoil_phys_PEM_timeseries(:,:,:,timelen)
1148	else
1149	TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:)
1150	endif !soil_pem
1151
1152	DO ig = 1,ngrid
1153	tsurf(ig) = 0.
1154	DO islope = 1,nslope
1155	tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) &
1156	* subslope_dist(ig,islope)
1157	ENDDO
1158	ENDDO ! of DO ig=1,ngrid
1159
1160	DO ig = 1,ngrid
1161	DO iloop = 1,nsoilmx
1162	tsoil(ig,iloop) = 0.
1163	inertiesoil(ig,iloop) = 0.
1164	DO islope = 1,nslope
1165	tsoil(ig,iloop) = tsoil(ig,iloop) + tsoil_slope(ig,iloop,islope) &
1166	* subslope_dist(ig,islope)
1167	inertiesoil(ig,iloop) = inertiesoil(ig,iloop) + TI_GCM_phys(ig,iloop,islope) &
1168	* subslope_dist(ig,islope)
1169	ENDDO
1170	ENDDO
1171	ENDDO ! of DO ig=1,ngrid
1172
1173	! III_a.3 Surface optical properties (for startfi)
1174
1175	DO ig = 1,ngrid
1176	DO l = 1,2
1177	albedo(ig,l) =0.
1178	DO islope = 1,nslope
1179	albedo(ig,l)= albedo(ig,l)+albedo_slope(ig,l,islope) &
1180	*subslope_dist(ig,islope)
1181	ENDDO
1182	ENDDO
1183	ENDDO
1184
1185	DO ig = 1,ngrid
1186	emis(ig) =0.
1187	DO islope = 1,nslope
1188	emis(ig)= emis(ig)+emiss_slope(ig,islope) &
1189	*subslope_dist(ig,islope)
1190	ENDDO
1191	ENDDO
1192
1193	! III_a.4 Pressure (for start)
1194	do i=1,ip1jmp1
1195	ps(i)=ps(i)*global_ave_press_new/global_ave_press_GCM
1196	enddo
1197
1198	do i = 1,ngrid
1199	ps_phys(i)=ps_phys(i)*global_ave_press_new/global_ave_press_GCM
1200	enddo
1201
1202	! III_a.5 Tracer (for start)
1203	allocate(zplev_new(ngrid,nlayer+1))
1204
1205	do l=1,nlayer+1
1206	do ig=1,ngrid
1207	zplev_new(ig,l) = ap(l) + bp(l)*ps_phys(ig)
1208	enddo
1209	enddo
1210
1211	DO nnq=1,nqtot
1212	if (noms(nnq).NE."co2") then
1213	DO l=1,llm-1
1214	DO ig=1,ngrid
1215	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))
1216	ENDDO
1217	q(:,llm,nnq)=q(:,llm-1,nnq)
1218	ENDDO
1219	else
1220	DO l=1,llm-1
1221	DO ig=1,ngrid
1222	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)) &
1223	+ ( (zplev_new(ig,l)-zplev_new(ig,l+1)) - &
1224	(zplev_gcm(ig,l)-zplev_gcm(ig,l+1)) ) / &
1225	(zplev_new(ig,l)-zplev_new(ig,l+1))
1226	ENDDO
1227	q(:,llm,nnq)=q(:,llm-1,nnq)
1228	ENDDO
1229	endif
1230	ENDDO
1231
1232	! Conserving the tracers's mass for GCM start files
1233	DO nnq=1,nqtot
1234	DO ig=1,ngrid
1235	DO l=1,llm-1
1236	if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number" .or. noms(nnq).NE."ccn_number") ) then
1237	extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
1238	q(ig,l,nnq)=1.
1239	q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
1240	endif
1241	if(q(ig,l,nnq).LT.0) then
1242	q(ig,l,nnq)=1E-30
1243	endif
1244	ENDDO
1245	enddo
1246	enddo
1247
1248	!------------------------
1249	if(evol_orbit_pem) then
1250	call recomp_orb_param(year_iter)
1251	endif
1252
1253	! III Output
1254	! III_a Update surface value for the PCM start files
1255	! III_b Write start and starfi.nc
1256
1257	!------------------------
1258
1259	! III_b.1 WRITE restart.nc
1260
1261	ptimestep=iphysiq*daysec/REAL(day_step)/nsplit_phys
1262	pday=day_ini
1263	ztime_fin=0.
1264
1265	allocate(p(ip1jmp1,nlayer+1))
1266	CALL pression (ip1jmp1,ap,bp,ps,p)
1267	CALL massdair(p,masse)
1268
1269	CALL dynredem0("restart_evol.nc", day_ini, phis)
1270
1271	CALL dynredem1("restart_evol.nc", &
1272	time_0,vcov,ucov,teta,q,masse,ps)
1273	print *, "restart_evol.nc has been written"
1274
1275	! III_b.2 WRITE restartfi.nc
1276
1277	call physdem0("restartfi_evol.nc",longitude,latitude, &
1278	nsoilmx,ngrid,nlayer,nq, &
1279	ptimestep,pday,0.,cell_area, &
1280	albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe, &
1281	hmons,summit,base,nslope,def_slope, &
1282	subslope_dist)
1283
1284	call physdem1("restartfi_evol.nc",nsoilmx,ngrid,nlayer,nq, &
1285	ptimestep,ztime_fin, &
1286	tsurf,tsoil,co2ice,albedo,emis, &
1287	q2,qsurf,tauscaling,totcloudfrac,wstar, &
1288	watercap,inertiesoil,nslope,co2ice_slope, &
1289	tsurf_slope,tsoil_slope, albedo_slope, &
1290	emiss_slope,qsurf_slope,watercap_slope, TI_GCM_phys)
1291
1292	print *, "restartfi_evol.nc has been written"
1293	!------------------------
1294
1295	! III Output
1296	! III_a Update surface value for the PCM start files
1297	! III_b Write start and starfi.nc
1298	! III_c Write start_pem
1299
1300	!------------------------
1301
1302	call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
1303	tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys)
1304
1305	print *, "restartfi_PEM.nc has been written"
1306	print *, "The PEM had run for ", year_iter, " years."
1307	print *, "LL & RV : So far so good"
1308
1309	END PROGRAM pem

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

Download in other formats:

Original Format