Context Navigation

pem.F90 @ 2889

Last change on this file since 2889 was 2888, checked in by llange, 3 years ago

PEM

Following r-2886-5, debugging some issues ( conservation of H2O, water_reservoir not initialized, info_pem which was not working)
Cleaning of some redundant routines (e.g., stopping criterion of the PEM) and variables renamed (e.g., alpha_criterion now transofrmed in ice_criterion)
Ice table subroutine is now in a module and recalled "compute_ice_table_equilibrium" (v.s. dynamic ice table, efforts ongoing)
Abort_pem introduced to avoid just stopping the pem with raw "stop" in the codes
conf_pem has been improved so that values are not hard coded (e.g., geothermal flux, mass of water_reservoir)
Regolith thermal properties updated in a cleaner way

(Not atomic commit, sorry)
LL & RV

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

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

Download in other formats:

Original Format