1 | !------------------------ |
---|
2 | ! I Initialization |
---|
3 | ! I_a Read the "run.def" |
---|
4 | ! I_b Read the "start.nc" and "startfi.nc" |
---|
5 | ! I_c Subslope parametrisation |
---|
6 | ! I_d Read the PCM data and convert them to the physical grid |
---|
7 | ! I_e Initialization of the PEM variable and soil |
---|
8 | ! I_f Compute tendencies |
---|
9 | ! I_g Save the initial situation |
---|
10 | ! I_h Read the "startpem.nc" |
---|
11 | ! I_i Compute orbit criterion |
---|
12 | |
---|
13 | ! II Run |
---|
14 | ! II_a Update pressure, ice and tracers |
---|
15 | ! II_b Evolution of ice |
---|
16 | ! II_c Flow of glaciers |
---|
17 | ! II_d Update surface and soil temperatures |
---|
18 | ! II_e Outputs |
---|
19 | ! II_f Update the tendencies |
---|
20 | ! II_g Checking the stopping criterion |
---|
21 | |
---|
22 | ! III Output |
---|
23 | ! III_a Update surface value for the PCM start files |
---|
24 | ! III_b Write the "restart.nc" and "restartfi.nc" |
---|
25 | ! III_c Write the "restartpem.nc" |
---|
26 | !------------------------ |
---|
27 | |
---|
28 | PROGRAM pem |
---|
29 | |
---|
30 | use phyetat0_mod, only: phyetat0 |
---|
31 | use phyredem, only: physdem0, physdem1 |
---|
32 | use netcdf, only: nf90_open, NF90_NOWRITE, nf90_get_var, nf90_inq_varid, nf90_close |
---|
33 | use turb_mod, only: q2, wstar |
---|
34 | use comslope_mod, only: nslope, def_slope, def_slope_mean, subslope_dist, iflat, ini_comslope_h |
---|
35 | use logic_mod, only: iflag_phys |
---|
36 | use mod_const_mpi, only: COMM_LMDZ |
---|
37 | use infotrac |
---|
38 | use geometry_mod, only: latitude_deg |
---|
39 | use conf_pem_mod, only: conf_pem |
---|
40 | use pemredem, only: pemdem0, pemdem1 |
---|
41 | use glaciers_mod, only: flow_co2glaciers, flow_h2oglaciers, co2ice_flow, h2oice_flow, inf_h2oice_threshold, & |
---|
42 | metam_h2oice_threshold, metam_co2ice_threshold, metam_h2oice, metam_co2ice |
---|
43 | use stopping_crit_mod, only: stopping_crit_h2o_ice, stopping_crit_co2 |
---|
44 | use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2 |
---|
45 | use evol_ice_mod, only: evol_co2_ice, evol_h2o_ice |
---|
46 | use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, & |
---|
47 | TI_PEM, & ! soil thermal inertia |
---|
48 | tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths |
---|
49 | fluxgeo ! Geothermal flux for the PEM and PCM |
---|
50 | use adsorption_mod, only: regolith_adsorption, adsorption_pem, & ! Bool to check if adsorption, main subroutine |
---|
51 | ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays |
---|
52 | co2_adsorbded_phys, h2o_adsorbded_phys ! Mass of co2 and h2O adsorbded |
---|
53 | use time_evol_mod, only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini |
---|
54 | use orbit_param_criterion_mod, only: orbit_param_criterion |
---|
55 | use recomp_orb_param_mod, only: recomp_orb_param |
---|
56 | use ice_table_mod, only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, & |
---|
57 | ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi |
---|
58 | use soil_thermalproperties_mod, only: update_soil_thermalproperties |
---|
59 | use time_phylmdz_mod, only: daysec, dtphys |
---|
60 | use abort_pem_mod, only: abort_pem |
---|
61 | use soil_settings_PEM_mod, only: soil_settings_PEM |
---|
62 | use compute_tend_mod, only: compute_tend |
---|
63 | use info_PEM_mod, only: info_PEM |
---|
64 | use nb_time_step_PCM_mod, only: nb_time_step_PCM |
---|
65 | use pemetat0_mod, only: pemetat0 |
---|
66 | use read_data_PCM_mod, only: read_data_PCM |
---|
67 | use recomp_tend_co2_slope_mod, only: recomp_tend_co2_slope |
---|
68 | use compute_soiltemp_mod, only: compute_tsoil_pem |
---|
69 | use writediagpem_mod, only: writediagpem, writediagsoilpem |
---|
70 | use co2condens_mod, only: CO2cond_ps |
---|
71 | use layering_mod, only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo |
---|
72 | use dyn_ss_ice_m_mod, only: dyn_ss_ice_m |
---|
73 | |
---|
74 | #ifndef CPP_STD |
---|
75 | use comsoil_h, only: tsoil, nsoilmx, ini_comsoil_h, inertiedat, mlayer, inertiesoil, flux_geo, nqsoil, qsoil |
---|
76 | use surfdat_h, only: tsurf, emis, qsurf, watercap, ini_surfdat_h, & |
---|
77 | albedodat, zmea, zstd, zsig, zgam, zthe, & |
---|
78 | albedo_h2o_frost,frost_albedo_threshold, & |
---|
79 | emissiv, watercaptag, perennial_co2ice, emisice, albedice |
---|
80 | use dimradmars_mod, only: totcloudfrac, albedo |
---|
81 | use dust_param_mod, only: tauscaling |
---|
82 | use tracer_mod, only: noms, igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses |
---|
83 | use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master |
---|
84 | use planete_h, only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit |
---|
85 | use comcstfi_h, only: pi, rad, g, mugaz, r |
---|
86 | use surfini_mod, only: surfini |
---|
87 | use comconst_mod, only: kappa, cpp |
---|
88 | #else |
---|
89 | use tracer_h, only: noms, igcm_h2o_ice, igcm_co2 ! Tracer names |
---|
90 | use phys_state_var_mod, only: cloudfrac, totcloudfrac, albedo_snow_SPECTV,HICE,RNAT, & |
---|
91 | PCTSRF_SIC, TSLAB, TSEA_ICE, SEA_ICE, ALBEDO_BAREGROUND, & |
---|
92 | ALBEDO_CO2_ICE_SPECTV, phys_state_var_init |
---|
93 | use aerosol_mod, only: iniaerosol |
---|
94 | use planete_mod, only: apoastr, periastr, year_day, peri_day, obliquit |
---|
95 | use comcstfi_mod, only: pi, rad, g, mugaz, r |
---|
96 | #endif |
---|
97 | |
---|
98 | #ifndef CPP_1D |
---|
99 | use iniphysiq_mod, only: iniphysiq |
---|
100 | use control_mod, only: iphysiq, day_step, nsplit_phys |
---|
101 | #else |
---|
102 | use time_phylmdz_mod, only: iphysiq, steps_per_sol |
---|
103 | use regular_lonlat_mod, only: init_regular_lonlat |
---|
104 | use physics_distribution_mod, only: init_physics_distribution |
---|
105 | use mod_grid_phy_lmdz, only: regular_lonlat |
---|
106 | use init_testphys1d_mod, only: init_testphys1d |
---|
107 | use comvert_mod, only: ap, bp |
---|
108 | use writerestart1D_mod, only: writerestart1D |
---|
109 | #endif |
---|
110 | |
---|
111 | implicit none |
---|
112 | |
---|
113 | include "dimensions.h" |
---|
114 | include "paramet.h" |
---|
115 | include "comgeom.h" |
---|
116 | include "iniprint.h" |
---|
117 | include "callkeys.h" |
---|
118 | |
---|
119 | integer ngridmx |
---|
120 | parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm) |
---|
121 | |
---|
122 | ! Same variable names as in the PCM |
---|
123 | integer, parameter :: nlayer = llm ! Number of vertical layer |
---|
124 | integer :: ngrid ! Number of physical grid points |
---|
125 | integer :: nq ! Number of tracer |
---|
126 | integer :: day_ini ! First day of the simulation |
---|
127 | real :: pday ! Physical day |
---|
128 | real :: time_phys ! Same as in PCM |
---|
129 | real :: ptimestep ! Same as in PCM |
---|
130 | real :: ztime_fin ! Same as in PCM |
---|
131 | |
---|
132 | ! Variables to read start.nc |
---|
133 | character(*), parameter :: start_name = "start.nc" ! Name of the file used to initialize the PEM |
---|
134 | |
---|
135 | ! Dynamic variables |
---|
136 | real, dimension(ip1jm,llm) :: vcov ! vents covariants |
---|
137 | real, dimension(ip1jmp1,llm) :: ucov ! vents covariants |
---|
138 | real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle |
---|
139 | real, dimension(:,:,:), allocatable :: q ! champs advectes |
---|
140 | real, dimension(ip1jmp1) :: ps ! pression au sol |
---|
141 | real, dimension(ip1jmp1) :: ps0 ! pression au sol initiale |
---|
142 | real, dimension(:), allocatable :: ps_start_PCM ! (ngrid) surface pressure |
---|
143 | real, dimension(:,:), allocatable :: ps_timeseries ! (ngrid x timelen) instantaneous surface pressure |
---|
144 | real, dimension(ip1jmp1,llm) :: masse ! masse d'air |
---|
145 | real, dimension(ip1jmp1) :: phis ! geopotentiel au sol |
---|
146 | real :: time_0 |
---|
147 | |
---|
148 | ! Variables to read starfi.nc |
---|
149 | character(*), parameter :: startfi_name = "startfi.nc" ! Name of the file used to initialize the PEM |
---|
150 | character(2) :: str2 |
---|
151 | integer :: ncid, status ! Variable for handling opening of files |
---|
152 | integer :: lonvarid, latvarid, areavarid, sdvarid ! Variable ID for Netcdf files |
---|
153 | integer :: apvarid, bpvarid ! Variable ID for Netcdf files |
---|
154 | |
---|
155 | ! Variables to read starfi.nc and write restartfi.nc |
---|
156 | real, dimension(:), allocatable :: longitude ! Longitude read in startfi_name and written in restartfi |
---|
157 | real, dimension(:), allocatable :: latitude ! Latitude read in startfi_name and written in restartfi |
---|
158 | real, dimension(:), allocatable :: cell_area ! Cell_area read in startfi_name and written in restartfi |
---|
159 | real :: Total_surface ! Total surface of the planet |
---|
160 | |
---|
161 | ! Variables for h2o_ice evolution |
---|
162 | real, dimension(:,:), allocatable :: h2o_ice ! h2o ice in the PEM |
---|
163 | real, dimension(:,:,:), allocatable :: min_h2o_ice ! Minima of h2o ice at each point for the PCM years [kg/m^2] |
---|
164 | real :: h2oice_ini_surf ! Initial surface of sublimating h2o ice |
---|
165 | logical, dimension(:,:), allocatable :: ini_h2oice_sublim ! Logical array to know if h2o ice is sublimating |
---|
166 | real :: global_avg_press_PCM ! constant: global average pressure retrieved in the PCM [Pa] |
---|
167 | real :: global_avg_press_old ! constant: Global average pressure of initial/previous time step |
---|
168 | real :: global_avg_press_new ! constant: Global average pressure of current time step |
---|
169 | real, dimension(:,:), allocatable :: zplev_new ! Physical x Atmospheric field: mass of the atmospheric layers in the pem at current time step [kg/m^2] |
---|
170 | real, dimension(:,:), allocatable :: zplev_PCM ! same but retrieved from the PCM [kg/m^2] |
---|
171 | real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2] |
---|
172 | real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step |
---|
173 | integer :: stopPEM ! which criterion is reached? 0 = no stopping; 1 = h2o ice surf; 2 = no h2o ice; 3 = co2 ice surf; 4 = ps; 5 = orb param; 6 = end of simu |
---|
174 | real, save :: A, B, mmean ! Molar mass: intermediate A, B for computations of the mean molar mass of the layer [mol/kg] |
---|
175 | real, dimension(:,:), allocatable :: q_h2o_PEM_phys ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg] |
---|
176 | integer :: timelen ! # time samples |
---|
177 | real :: ave ! intermediate varibale to compute average |
---|
178 | real, dimension(:,:), allocatable :: p ! Physics x Atmosphere: pressure to recompute and write in restart (ngrid,llmp1) |
---|
179 | real :: extra_mass ! Intermediate variables Extra mass of a tracer if it is greater than 1 |
---|
180 | |
---|
181 | ! Variables for co2_ice evolution |
---|
182 | real, dimension(:,:), allocatable :: co2_ice ! co2 ice in the PEM |
---|
183 | logical, dimension(:,:), allocatable :: is_co2ice_ini ! Was there co2 ice initially in the PEM? |
---|
184 | real, dimension(:,:,:), allocatable :: min_co2_ice ! Minimum of co2 ice at each point for the first year [kg/m^2] |
---|
185 | real :: co2ice_ini_surf ! Initial surface of sublimating co2 ice |
---|
186 | logical, dimension(:,:), allocatable :: ini_co2ice_sublim ! Logical array to know if co2 ice is sublimating |
---|
187 | real, dimension(:,:), allocatable :: vmr_co2_PCM ! Physics x Times co2 volume mixing ratio retrieve from the PCM [m^3/m^3] |
---|
188 | real, dimension(:,:), allocatable :: vmr_co2_PEM_phys ! Physics x Times co2 volume mixing ratio used in the PEM |
---|
189 | real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg] |
---|
190 | |
---|
191 | ! Variables for stratification (layering) evolution |
---|
192 | type(layering), dimension(:,:), allocatable :: stratif ! Layering (linked list of "stratum") for each grid point and slope |
---|
193 | type(ptrarray), dimension(:,:), allocatable :: current1, current2 ! Current active stratum in the layering |
---|
194 | logical, dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm |
---|
195 | |
---|
196 | ! Variables for slopes |
---|
197 | real, dimension(:,:), allocatable :: d_co2ice ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year |
---|
198 | real, dimension(:,:), allocatable :: d_co2ice_ini ! physical point x slope field: Tendency of evolution of perennial co2 ice over a year in the PCM |
---|
199 | real, dimension(:,:), allocatable :: d_h2oice ! physical point x slope field: Tendency of evolution of perennial h2o ice |
---|
200 | real, dimension(:,:), allocatable :: flag_co2flow ! (ngrid,nslope): Flag where there is a CO2 glacier flow |
---|
201 | real, dimension(:), allocatable :: flag_co2flow_mesh ! (ngrid) : Flag where there is a CO2 glacier flow |
---|
202 | real, dimension(:,:), allocatable :: flag_h2oflow ! (ngrid,nslope): Flag where there is a H2O glacier flow |
---|
203 | real, dimension(:), allocatable :: flag_h2oflow_mesh ! (ngrid) : Flag where there is a H2O glacier flow |
---|
204 | |
---|
205 | ! Variables for surface and soil |
---|
206 | real, dimension(:,:), allocatable :: tsurf_avg ! Physic x SLOPE field: Averaged Surface Temperature [K] |
---|
207 | real, dimension(:,:,:), allocatable :: tsoil_avg ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K] |
---|
208 | real, dimension(:,:,:), allocatable :: tsurf_PCM_timeseries ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K] |
---|
209 | real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K] |
---|
210 | real, dimension(:,:,:,:), allocatable :: tsoil_anom ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K] |
---|
211 | real, dimension(:,:), allocatable :: tsurf_avg_yr1 ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K] |
---|
212 | real, dimension(:,:), allocatable :: Tsoil_locslope ! Physic x Soil: Intermediate when computing Tsoil [K] |
---|
213 | real, dimension(:), allocatable :: Tsurf_locslope ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K] |
---|
214 | real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3] |
---|
215 | real, dimension(:,:), allocatable :: watersurf_density_avg ! Physic x Slope, water surface density, yearly averaged [kg/m^3] |
---|
216 | real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3] |
---|
217 | real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3] |
---|
218 | real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K] |
---|
219 | real, dimension(:), allocatable :: delta_co2_adsorbded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] |
---|
220 | real, dimension(:), allocatable :: delta_h2o_adsorbded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] |
---|
221 | real :: totmassco2_adsorbded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] |
---|
222 | real :: totmassh2o_adsorbded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] |
---|
223 | logical :: bool_sublim ! logical to check if there is sublimation or not |
---|
224 | logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep |
---|
225 | real, dimension(:,:), allocatable :: icetable_thickness_old ! ngrid x nslope: Thickness of the ice table at the previous iteration [m] |
---|
226 | real, dimension(:,:,:), allocatable :: ice_porefilling_old ! ngrid x nslope: Ice pore filling at the previous iteration [m] |
---|
227 | real, dimension(:), allocatable :: delta_h2o_icetablesublim ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg] |
---|
228 | real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table |
---|
229 | real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table |
---|
230 | |
---|
231 | ! Some variables for the PEM run |
---|
232 | real, parameter :: year_step = 1 ! Timestep for the pem |
---|
233 | real :: i_myear_leg ! Number of iteration |
---|
234 | real :: n_myear_leg ! Maximum number of iterations before stopping |
---|
235 | real :: i_myear ! Global number of Martian years of the chained simulations |
---|
236 | real :: n_myear ! Maximum number of Martian years of the chained simulations |
---|
237 | real :: timestep ! Timestep [s] |
---|
238 | character(20) :: job_id ! Job id provided as argument passed on the command line when the program was invoked |
---|
239 | logical :: timewall ! Flag to use the time limit stopping criterion in case of a PEM job |
---|
240 | integer(kind=8) :: cr ! Number of clock ticks per second (count rate) |
---|
241 | integer(kind=8) :: c1, c2 ! Counts of processor clock |
---|
242 | character(100) :: chtimelimit ! Time limit for the PEM job outputted by the SLURM command |
---|
243 | real :: timelimit ! Time limit for the PEM job in seconds |
---|
244 | real, parameter :: antetime = 1200 ! Anticipation time to prevent reaching the time limit: 1200 s = 20 min by default |
---|
245 | integer :: cstat, days, hours, minutes, seconds |
---|
246 | character(1) :: sep |
---|
247 | |
---|
248 | #ifdef CPP_STD |
---|
249 | real :: frost_albedo_threshold = 0.05 ! frost albedo threeshold to convert fresh frost to old ice |
---|
250 | real :: albedo_h2o_frost ! albedo of h2o frost |
---|
251 | real, dimension(:), allocatable :: tsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
252 | real, dimension(:,:), allocatable :: qsurf_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
253 | real, dimension(:,:), allocatable :: tsoil_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
254 | real, dimension(:), allocatable :: emis_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
255 | real, dimension(:,:), allocatable :: albedo_read_generic ! Temporary variable to do the subslope transfert dimension when reading form generic |
---|
256 | real, dimension(:,:), allocatable :: tsurf ! Subslope variable, only needed in the GENERIC case |
---|
257 | real, dimension(:,:,:), allocatable :: qsurf ! Subslope variable, only needed in the GENERIC case |
---|
258 | real, dimension(:,:,:), allocatable :: tsoil ! Subslope variable, only needed in the GENERIC case |
---|
259 | real, dimension(:,:), allocatable :: emis ! Subslope variable, only needed in the GENERIC case |
---|
260 | real, dimension(:,:), allocatable :: watercap ! Subslope variable, only needed in the GENERIC case =0 no watercap in generic model |
---|
261 | logical, dimension(:), allocatable :: watercaptag ! Subslope variable, only needed in the GENERIC case =false no watercaptag in generic model |
---|
262 | real, dimension(:,:,:), allocatable :: albedo ! Subslope variable, only needed in the GENERIC case |
---|
263 | real, dimension(:,:,:), allocatable :: inertiesoil ! Subslope variable, only needed in the GENERIC case |
---|
264 | #endif |
---|
265 | |
---|
266 | #ifdef CPP_1D |
---|
267 | integer :: nsplit_phys |
---|
268 | integer, parameter :: jjm_value = jjm - 1 |
---|
269 | integer :: day_step |
---|
270 | |
---|
271 | ! Dummy variables to use the subroutine 'init_testphys1d' |
---|
272 | logical :: therestart1D, therestartfi |
---|
273 | integer :: ndt, day0 |
---|
274 | real :: ptif, pks, day, gru, grv, atm_wat_profile, atm_wat_tau |
---|
275 | real, dimension(:), allocatable :: zqsat |
---|
276 | real, dimension(:,:,:), allocatable :: dq, dqdyn |
---|
277 | real, dimension(nlayer) :: play, w |
---|
278 | real, dimension(nlayer + 1) :: plev |
---|
279 | #else |
---|
280 | integer, parameter :: jjm_value = jjm |
---|
281 | real, dimension(:), allocatable :: ap ! Coefficient ap read in start_name and written in restart |
---|
282 | real, dimension(:), allocatable :: bp ! Coefficient bp read in start_name and written in restart |
---|
283 | #endif |
---|
284 | |
---|
285 | ! Loop variables |
---|
286 | integer :: i, l, ig, nnq, t, islope, ig_loop, islope_loop, isoil, icap |
---|
287 | |
---|
288 | ! Elapsed time with system clock |
---|
289 | call system_clock(count_rate = cr) |
---|
290 | call system_clock(c1) |
---|
291 | timewall = .true. |
---|
292 | timelimit = 86400 ! 86400 seconds = 24 h by default |
---|
293 | if (command_argument_count() > 0) then |
---|
294 | ! Read the job id passed as argument to the program |
---|
295 | call get_command_argument(1,job_id) |
---|
296 | ! Execute the system command |
---|
297 | call execute_command_line('squeue -j '//trim(job_id)//' -h --Format TimeLimit > tmp_cmdout.txt',cmdstat = cstat) |
---|
298 | if (cstat /= 0) then |
---|
299 | call execute_command_line('qstat -f '//trim(job_id)//' | grep "Walltime" | awk ''{print $3}'' > tmp_cmdout.txt',cmdstat = cstat) |
---|
300 | if (cstat > 0) then |
---|
301 | error stop 'pem: command execution failed!' |
---|
302 | else if (cstat < 0) then |
---|
303 | error stop 'pem: command execution not supported (neither SLURM nor PBS/TORQUE is installed)!' |
---|
304 | endif |
---|
305 | endif |
---|
306 | ! Read the output |
---|
307 | open(1,file = 'tmp_cmdout.txt',status = 'old') |
---|
308 | read(1,'(a)') chtimelimit |
---|
309 | close(1) |
---|
310 | chtimelimit = trim(chtimelimit) |
---|
311 | call execute_command_line('rm tmp_cmdout.txt',cmdstat = cstat) |
---|
312 | if (cstat > 0) then |
---|
313 | error stop 'pem: command execution failed!' |
---|
314 | else if (cstat < 0) then |
---|
315 | error stop 'pem: command execution not supported!' |
---|
316 | endif |
---|
317 | if (index(chtimelimit,'-') > 0) then ! 'chtimelimit' format is "D-HH:MM:SS" |
---|
318 | read(chtimelimit,'(i1,a1,i2,a1,i2,a1,i2)') days, sep, hours, sep, minutes, sep, seconds |
---|
319 | timelimit = days*86400 + hours*3600 + minutes*60 + seconds |
---|
320 | else if (index(chtimelimit,':') > 0 .and. len_trim(chtimelimit) > 5) then ! 'chtimelimit' format is "HH:MM:SS" |
---|
321 | read(chtimelimit,'(i2,a1,i2,a1,i2)') hours, sep, minutes, sep, seconds |
---|
322 | timelimit = hours*3600 + minutes*60 + seconds |
---|
323 | else ! 'chtimelimit' format is "MM:SS" |
---|
324 | read(chtimelimit,'(i2,a1,i2)') minutes, sep, seconds |
---|
325 | timelimit = minutes*60 + seconds |
---|
326 | endif |
---|
327 | else |
---|
328 | timewall = .false. |
---|
329 | endif |
---|
330 | |
---|
331 | ! Parallel variables |
---|
332 | #ifndef CPP_STD |
---|
333 | is_sequential = .true. |
---|
334 | is_parallel = .false. |
---|
335 | is_mpi_root = .true. |
---|
336 | is_omp_root = .true. |
---|
337 | is_master = .true. |
---|
338 | #endif |
---|
339 | |
---|
340 | ! Some constants |
---|
341 | day_ini = 0 ! test |
---|
342 | time_phys = 0. ! test |
---|
343 | ngrid = ngridmx |
---|
344 | A = (1/m_co2 - 1/m_noco2) |
---|
345 | B = 1/m_noco2 |
---|
346 | year_day = 669 |
---|
347 | daysec = 88775. |
---|
348 | timestep = year_day*daysec*year_step |
---|
349 | |
---|
350 | !----------------------------- INITIALIZATION -------------------------- |
---|
351 | !------------------------ |
---|
352 | ! I Initialization |
---|
353 | ! I_a Read the "run.def" |
---|
354 | !------------------------ |
---|
355 | #ifndef CPP_1D |
---|
356 | dtphys = daysec/48. ! Dummy value (overwritten in phyetat0) |
---|
357 | call conf_gcm(99,.true.) |
---|
358 | call infotrac_init |
---|
359 | nq = nqtot |
---|
360 | allocate(q(ip1jmp1,llm,nqtot)) |
---|
361 | allocate(longitude(ngrid),latitude(ngrid),cell_area(ngrid)) |
---|
362 | #else |
---|
363 | allocate(q(1,llm,nqtot)) |
---|
364 | allocate(longitude(1),latitude(1),cell_area(1)) |
---|
365 | |
---|
366 | therestart1D = .false. ! Default value |
---|
367 | inquire(file = 'start1D.txt',exist = therestart1D) |
---|
368 | if (.not. therestart1D) then |
---|
369 | write(*,*) 'There is no "start1D.txt" file!' |
---|
370 | error stop 'Initialization cannot be done for the 1D PEM.' |
---|
371 | endif |
---|
372 | therestartfi = .false. ! Default value |
---|
373 | inquire(file = 'startfi.nc',exist = therestartfi) |
---|
374 | if (.not. therestartfi) then |
---|
375 | write(*,*) 'There is no "startfi.nc" file!' |
---|
376 | error stop 'Initialization cannot be done for the 1D PEM.' |
---|
377 | endif |
---|
378 | |
---|
379 | call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,610.,nq,q, & |
---|
380 | time_0,ps(1),ucov,vcov,teta,ndt,ptif,pks,dtphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, & |
---|
381 | play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau) |
---|
382 | ps(2) = ps(1) |
---|
383 | nsplit_phys = 1 |
---|
384 | day_step = steps_per_sol |
---|
385 | #endif |
---|
386 | |
---|
387 | call conf_pem(i_myear,n_myear) |
---|
388 | |
---|
389 | !------------------------ |
---|
390 | ! I Initialization |
---|
391 | ! I_b Read of the "start.nc" and starfi_evol.nc |
---|
392 | !------------------------ |
---|
393 | ! I_b.1 Read "start.nc" |
---|
394 | allocate(ps_start_PCM(ngrid)) |
---|
395 | #ifndef CPP_1D |
---|
396 | call dynetat0(start_name,vcov,ucov,teta,q,masse,ps,phis,time_0) |
---|
397 | |
---|
398 | call gr_dyn_fi(1,iip1,jjp1,ngridmx,ps,ps_start_PCM) |
---|
399 | |
---|
400 | call iniconst !new |
---|
401 | call inigeom |
---|
402 | |
---|
403 | allocate(ap(nlayer + 1)) |
---|
404 | allocate(bp(nlayer + 1)) |
---|
405 | status = nf90_open(start_name,NF90_NOWRITE,ncid) |
---|
406 | status = nf90_inq_varid(ncid,"ap",apvarid) |
---|
407 | status = nf90_get_var(ncid,apvarid,ap) |
---|
408 | status = nf90_inq_varid(ncid,"bp",bpvarid) |
---|
409 | status = nf90_get_var(ncid,bpvarid,bp) |
---|
410 | status = nf90_close(ncid) |
---|
411 | |
---|
412 | call iniphysiq(iim,jjm,llm,(jjm-1)*iim+2,comm_lmdz,daysec,day_ini,dtphys/nsplit_phys,rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,iflag_phys) |
---|
413 | #else |
---|
414 | ps_start_PCM(1) = ps(1) |
---|
415 | #endif |
---|
416 | |
---|
417 | ! Save initial surface pressure |
---|
418 | ps0 = ps |
---|
419 | |
---|
420 | ! In the PCM, these values are given to the physic by the dynamic. |
---|
421 | ! Here we simply read them in the "startfi.nc" file |
---|
422 | status = nf90_open(startfi_name, NF90_NOWRITE, ncid) |
---|
423 | |
---|
424 | status = nf90_inq_varid(ncid,"longitude",lonvarid) |
---|
425 | status = nf90_get_var(ncid,lonvarid,longitude) |
---|
426 | |
---|
427 | status = nf90_inq_varid(ncid,"latitude",latvarid) |
---|
428 | status = nf90_get_var(ncid,latvarid,latitude) |
---|
429 | |
---|
430 | status = nf90_inq_varid(ncid,"area",areavarid) |
---|
431 | status = nf90_get_var(ncid,areavarid,cell_area) |
---|
432 | |
---|
433 | status = nf90_inq_varid(ncid,"soildepth",sdvarid) |
---|
434 | status = nf90_get_var(ncid,sdvarid,mlayer) |
---|
435 | |
---|
436 | status = nf90_close(ncid) |
---|
437 | |
---|
438 | ! I_b.2 Read the "startfi.nc" |
---|
439 | ! First we read the initial state (starfi.nc) |
---|
440 | #ifndef CPP_STD |
---|
441 | call phyetat0(startfi_name,0,0,nsoilmx,ngrid,nlayer,nq,nqsoil,day_ini,time_phys,tsurf, & |
---|
442 | tsoil,albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac,wstar, & |
---|
443 | watercap,perennial_co2ice,def_slope,def_slope_mean,subslope_dist) |
---|
444 | |
---|
445 | ! Remove unphysical values of surface tracer |
---|
446 | where (qsurf < 0.) qsurf = 0. |
---|
447 | |
---|
448 | call surfini(ngrid,nslope,qsurf) |
---|
449 | #else |
---|
450 | call phys_state_var_init(nq) |
---|
451 | if (.not. allocated(noms)) allocate(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on) |
---|
452 | call initracer(ngrid,nq) |
---|
453 | call iniaerosol() |
---|
454 | allocate(tsurf_read_generic(ngrid)) |
---|
455 | allocate(qsurf_read_generic(ngrid,nq)) |
---|
456 | allocate(tsoil_read_generic(ngrid,nsoilmx)) |
---|
457 | allocate(qsoil_read_generic(ngrid,nsoilmx,nqsoil,nslope)) |
---|
458 | allocate(emis_read_generic(ngrid)) |
---|
459 | allocate(tsurf(ngrid,1)) |
---|
460 | allocate(qsurf(ngrid,nq,1)) |
---|
461 | allocate(tsoil(ngrid,nsoilmx,1)) |
---|
462 | allocate(emis(ngrid,1)) |
---|
463 | allocate(watercap(ngrid,1)) |
---|
464 | allocate(watercaptag(ngrid)) |
---|
465 | allocate(albedo_read_generic(ngrid,2)) |
---|
466 | allocate(albedo(ngrid,2,1)) |
---|
467 | allocate(inertiesoil(ngrid,nsoilmx,1)) |
---|
468 | call phyetat0(.true.,ngrid,nlayer,startfi_name,0,0,nsoilmx,nq,nqsoil,day_ini,time_phys, & |
---|
469 | tsurf_read_generic,tsoil_read_generic,emis_read_generic,q2, & |
---|
470 | qsurf_read_generic,qsoil_read_generic,cloudfrac,totcloudfrac,hice, & |
---|
471 | rnat,pctsrf_sic,tslab,tsea_ice,sea_ice) |
---|
472 | call surfini(ngrid,nq,qsurf_read_generic,albedo_read_generic,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) |
---|
473 | |
---|
474 | nslope = 1 |
---|
475 | call ini_comslope_h(ngrid,1) |
---|
476 | |
---|
477 | qsurf(:,:,1) = qsurf_read_generic |
---|
478 | tsurf(:,1) = tsurf_read_generic |
---|
479 | tsoil(:,:,1) = tsoil_read_generic |
---|
480 | emis(:,1) = emis_read_generic |
---|
481 | watercap(:,1) = 0. |
---|
482 | watercaptag(:) = .false. |
---|
483 | albedo(:,1,1) = albedo_read_generic(:,1) |
---|
484 | albedo(:,2,1) = albedo_read_generic(:,2) |
---|
485 | inertiesoil(:,:,1) = inertiedat |
---|
486 | |
---|
487 | if (nslope == 1) then |
---|
488 | def_slope(1) = 0 |
---|
489 | def_slope(2) = 0 |
---|
490 | def_slope_mean = 0 |
---|
491 | subslope_dist(:,1) = 1. |
---|
492 | endif |
---|
493 | |
---|
494 | ! Remove unphysical values of surface tracer |
---|
495 | qsurf(:,:,1) = qsurf_read_generic |
---|
496 | where (qsurf < 0.) qsurf = 0. |
---|
497 | #endif |
---|
498 | |
---|
499 | do nnq = 1,nqtot ! Why not using ini_tracer ? |
---|
500 | if (noms(nnq) == "h2o_ice") igcm_h2o_ice = nnq |
---|
501 | if (noms(nnq) == "h2o_vap") then |
---|
502 | igcm_h2o_vap = nnq |
---|
503 | mmol(igcm_h2o_vap) = 18. |
---|
504 | endif |
---|
505 | if (noms(nnq) == "co2") igcm_co2 = nnq |
---|
506 | enddo |
---|
507 | r = 8.314511*1000./mugaz |
---|
508 | |
---|
509 | !------------------------ |
---|
510 | ! I Initialization |
---|
511 | ! I_c Subslope parametrisation |
---|
512 | !------------------------ |
---|
513 | ! Define some slope statistics |
---|
514 | iflat = 1 |
---|
515 | do islope = 2,nslope |
---|
516 | if (abs(def_slope_mean(islope)) < abs(def_slope_mean(iflat))) iflat = islope |
---|
517 | enddo |
---|
518 | |
---|
519 | write(*,*) 'Flat slope for islope = ',iflat |
---|
520 | write(*,*) 'corresponding criterium = ',def_slope_mean(iflat) |
---|
521 | |
---|
522 | allocate(flag_co2flow(ngrid,nslope)) |
---|
523 | allocate(flag_co2flow_mesh(ngrid)) |
---|
524 | allocate(flag_h2oflow(ngrid,nslope)) |
---|
525 | allocate(flag_h2oflow_mesh(ngrid)) |
---|
526 | |
---|
527 | flag_co2flow = 0 |
---|
528 | flag_co2flow_mesh = 0 |
---|
529 | flag_h2oflow = 0 |
---|
530 | flag_h2oflow_mesh = 0 |
---|
531 | |
---|
532 | !------------------------ |
---|
533 | ! I Initialization |
---|
534 | ! I_d Read the PCM data and convert them to the physical grid |
---|
535 | !------------------------ |
---|
536 | ! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value |
---|
537 | call nb_time_step_PCM("data_PCM_Y1.nc",timelen) |
---|
538 | |
---|
539 | allocate(tsoil_avg(ngrid,nsoilmx,nslope)) |
---|
540 | allocate(watersoil_density_PEM_avg(ngrid,nsoilmx_PEM,nslope)) |
---|
541 | allocate(vmr_co2_PCM(ngrid,timelen)) |
---|
542 | allocate(ps_timeseries(ngrid,timelen)) |
---|
543 | allocate(min_co2_ice(ngrid,nslope,2)) |
---|
544 | allocate(min_h2o_ice(ngrid,nslope,2)) |
---|
545 | allocate(tsurf_avg_yr1(ngrid,nslope)) |
---|
546 | allocate(tsurf_avg(ngrid,nslope)) |
---|
547 | allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen)) |
---|
548 | allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen)) |
---|
549 | allocate(q_co2_PEM_phys(ngrid,timelen)) |
---|
550 | allocate(q_h2o_PEM_phys(ngrid,timelen)) |
---|
551 | allocate(watersurf_density_avg(ngrid,nslope)) |
---|
552 | allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) |
---|
553 | allocate(Tsurfavg_before_saved(ngrid,nslope)) |
---|
554 | allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) |
---|
555 | allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) |
---|
556 | allocate(delta_co2_adsorbded(ngrid)) |
---|
557 | allocate(co2ice_disappeared(ngrid,nslope)) |
---|
558 | allocate(icetable_thickness_old(ngrid,nslope)) |
---|
559 | allocate(ice_porefilling_old(ngrid,nsoilmx_PEM,nslope)) |
---|
560 | allocate(delta_h2o_icetablesublim(ngrid)) |
---|
561 | allocate(delta_h2o_adsorbded(ngrid)) |
---|
562 | allocate(vmr_co2_PEM_phys(ngrid,timelen)) |
---|
563 | |
---|
564 | write(*,*) "Downloading data Y1..." |
---|
565 | call read_data_PCM("data_PCM_Y1.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,1),min_h2o_ice(:,:,1), & |
---|
566 | tsurf_avg_yr1,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, & |
---|
567 | watersurf_density_avg,watersoil_density_timeseries) |
---|
568 | write(*,*) "Downloading data Y1 done!" |
---|
569 | |
---|
570 | ! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value |
---|
571 | write(*,*) "Downloading data Y2..." |
---|
572 | call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_PCM,ps_timeseries,min_co2_ice(:,:,2),min_h2o_ice(:,:,2), & |
---|
573 | tsurf_avg,tsoil_avg,tsurf_PCM_timeseries,tsoil_anom,q_co2_PEM_phys,q_h2o_PEM_phys, & |
---|
574 | watersurf_density_avg,watersoil_density_timeseries) |
---|
575 | write(*,*) "Downloading data Y2 done!" |
---|
576 | |
---|
577 | !------------------------ |
---|
578 | ! I Initialization |
---|
579 | ! I_e Initialization of the PEM variables and soil |
---|
580 | !------------------------ |
---|
581 | call end_comsoil_h_PEM |
---|
582 | call ini_comsoil_h_PEM(ngrid,nslope) |
---|
583 | call end_adsorption_h_PEM |
---|
584 | call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM) |
---|
585 | call end_ice_table |
---|
586 | call ini_ice_table(ngrid,nslope,nsoilmx_PEM) |
---|
587 | |
---|
588 | if (soil_pem) then |
---|
589 | do t = 1,timelen |
---|
590 | tsoil_anom(:,:,:,t) = tsoil_anom(:,:,:,t) - tsoil_avg ! compute anomaly between Tsoil(t) in the startfi - <Tsoil> to recompute properly tsoil in the restart |
---|
591 | enddo |
---|
592 | call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM) |
---|
593 | tsoil_PEM(:,1:nsoilmx,:) = tsoil_avg |
---|
594 | tsoil_phys_PEM_timeseries(:,1:nsoilmx,:,:) = tsoil_anom |
---|
595 | watersoil_density_PEM_timeseries(:,1:nsoilmx,:,:) = watersoil_density_timeseries |
---|
596 | do l = nsoilmx + 1,nsoilmx_PEM |
---|
597 | tsoil_PEM(:,l,:) = tsoil_avg(:,nsoilmx,:) |
---|
598 | watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:) |
---|
599 | enddo |
---|
600 | watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen |
---|
601 | endif !soil_pem |
---|
602 | deallocate(tsoil_avg) |
---|
603 | |
---|
604 | !------------------------ |
---|
605 | ! I Initialization |
---|
606 | ! I_f Compute tendencies |
---|
607 | !------------------------ |
---|
608 | allocate(d_co2ice(ngrid,nslope),d_h2oice(ngrid,nslope)) |
---|
609 | allocate(d_co2ice_ini(ngrid,nslope)) |
---|
610 | |
---|
611 | ! Compute the tendencies of the evolution of ice over the years |
---|
612 | call compute_tend(ngrid,nslope,min_co2_ice,d_co2ice) |
---|
613 | call compute_tend(ngrid,nslope,min_h2o_ice,d_h2oice) |
---|
614 | d_co2ice_ini = d_co2ice |
---|
615 | deallocate(min_co2_ice,min_h2o_ice) |
---|
616 | |
---|
617 | !------------------------ |
---|
618 | ! I Initialization |
---|
619 | ! I_g Save the initial situation |
---|
620 | !------------------------ |
---|
621 | allocate(zplev_PCM(ngrid,nlayer + 1)) |
---|
622 | Total_surface = 0. |
---|
623 | do ig = 1,ngrid |
---|
624 | Total_surface = Total_surface + cell_area(ig) |
---|
625 | zplev_PCM(ig,:) = ap + bp*ps_start_PCM(ig) |
---|
626 | enddo |
---|
627 | global_avg_press_old = sum(cell_area*ps_start_PCM)/Total_surface |
---|
628 | global_avg_press_PCM = global_avg_press_old |
---|
629 | global_avg_press_new = global_avg_press_old |
---|
630 | write(*,*) "Total surface of the planet =", Total_surface |
---|
631 | write(*,*) "Initial global average pressure =", global_avg_press_PCM |
---|
632 | |
---|
633 | !------------------------ |
---|
634 | ! I Initialization |
---|
635 | ! I_h Read the "startpem.nc" |
---|
636 | !------------------------ |
---|
637 | allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope)) |
---|
638 | |
---|
639 | allocate(stratif(ngrid,nslope)) |
---|
640 | if (layering_algo) then |
---|
641 | do islope = 1,nslope |
---|
642 | do i = 1,ngrid |
---|
643 | call ini_layering(stratif(i,islope)) |
---|
644 | enddo |
---|
645 | enddo |
---|
646 | endif |
---|
647 | |
---|
648 | call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, & |
---|
649 | icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys, & |
---|
650 | ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice, & |
---|
651 | global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys, & |
---|
652 | delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif) |
---|
653 | |
---|
654 | ! We save the places where h2o ice is sublimating |
---|
655 | ! We compute the surface of h2o ice sublimating |
---|
656 | allocate(ini_co2ice_sublim(ngrid,nslope),ini_h2oice_sublim(ngrid,nslope),is_co2ice_ini(ngrid,nslope)) |
---|
657 | co2ice_ini_surf = 0. |
---|
658 | h2oice_ini_surf = 0. |
---|
659 | ini_co2ice_sublim = .false. |
---|
660 | ini_h2oice_sublim = .false. |
---|
661 | is_co2ice_ini = .false. |
---|
662 | co2ice_disappeared = .false. |
---|
663 | do i = 1,ngrid |
---|
664 | do islope = 1,nslope |
---|
665 | if (co2_ice(i,islope) > 0.) is_co2ice_ini(i,islope) = .true. |
---|
666 | if (d_co2ice(i,islope) < 0. .and. co2_ice(i,islope) > 0.) then |
---|
667 | ini_co2ice_sublim(i,islope) = .true. |
---|
668 | co2ice_ini_surf = co2ice_ini_surf + cell_area(i)*subslope_dist(i,islope) |
---|
669 | endif |
---|
670 | if (d_h2oice(i,islope) < 0. .and. h2o_ice(i,islope) > 0.) then |
---|
671 | ini_h2oice_sublim(i,islope) = .true. |
---|
672 | h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope) |
---|
673 | endif |
---|
674 | enddo |
---|
675 | enddo |
---|
676 | write(*,*) "Total initial surface of co2 ice sublimating (slope) =", co2ice_ini_surf |
---|
677 | write(*,*) "Total initial surface of h2o ice sublimating (slope) =", h2oice_ini_surf |
---|
678 | |
---|
679 | delta_h2o_icetablesublim = 0. |
---|
680 | |
---|
681 | if (adsorption_pem) then |
---|
682 | totmassco2_adsorbded = 0. |
---|
683 | totmassh2o_adsorbded = 0. |
---|
684 | do ig = 1,ngrid |
---|
685 | do islope = 1,nslope |
---|
686 | do l = 1,nsoilmx_PEM - 1 |
---|
687 | if (l==1) then |
---|
688 | totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & |
---|
689 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
690 | totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & |
---|
691 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
692 | else |
---|
693 | totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & |
---|
694 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
695 | totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & |
---|
696 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
697 | endif |
---|
698 | enddo |
---|
699 | enddo |
---|
700 | enddo |
---|
701 | write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded |
---|
702 | write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded |
---|
703 | endif ! adsorption |
---|
704 | deallocate(tsurf_avg_yr1) |
---|
705 | |
---|
706 | !------------------------ |
---|
707 | ! I Initialization |
---|
708 | ! I_i Compute orbit criterion |
---|
709 | !------------------------ |
---|
710 | #ifndef CPP_STD |
---|
711 | call iniorbit(aphelie,periheli,year_day,peri_day,obliquit) |
---|
712 | #else |
---|
713 | call iniorbit(apoastr,periastr,year_day,peri_day,obliquit) |
---|
714 | #endif |
---|
715 | |
---|
716 | n_myear_leg = Max_iter_pem |
---|
717 | if (evol_orbit_pem) call orbit_param_criterion(i_myear,n_myear_leg) |
---|
718 | |
---|
719 | !-------------------------- END INITIALIZATION ------------------------- |
---|
720 | |
---|
721 | !-------------------------------- RUN ---------------------------------- |
---|
722 | !------------------------ |
---|
723 | ! II Run |
---|
724 | ! II_a Update pressure, ice and tracers |
---|
725 | !------------------------ |
---|
726 | i_myear_leg = 0 |
---|
727 | stopPEM = 0 |
---|
728 | if (layering_algo) then |
---|
729 | allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope)) |
---|
730 | new_str = .true. |
---|
731 | new_lag1 = .true. |
---|
732 | new_lag2 = .true. |
---|
733 | do islope = 1,nslope |
---|
734 | do ig = 1,ngrid |
---|
735 | current1(ig,islope)%p => stratif(ig,islope)%top |
---|
736 | current2(ig,islope)%p => stratif(ig,islope)%top |
---|
737 | enddo |
---|
738 | enddo |
---|
739 | endif |
---|
740 | |
---|
741 | do while (i_myear_leg < n_myear_leg .and. i_myear < n_myear) |
---|
742 | ! II.a.1. Compute updated global pressure |
---|
743 | write(*,*) "Recomputing the new pressure..." |
---|
744 | do i = 1,ngrid |
---|
745 | do islope = 1,nslope |
---|
746 | global_avg_press_new = global_avg_press_new - CO2cond_ps*g*cell_area(i)*d_co2ice(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface |
---|
747 | enddo |
---|
748 | enddo |
---|
749 | |
---|
750 | if (adsorption_pem) then |
---|
751 | do i = 1,ngrid |
---|
752 | global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface |
---|
753 | enddo |
---|
754 | endif |
---|
755 | write(*,*) 'Global average pressure old time step',global_avg_press_old |
---|
756 | write(*,*) 'Global average pressure new time step',global_avg_press_new |
---|
757 | |
---|
758 | ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consumption) |
---|
759 | allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen)) |
---|
760 | write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..." |
---|
761 | do l = 1,nlayer + 1 |
---|
762 | do ig = 1,ngrid |
---|
763 | zplev_old_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) |
---|
764 | enddo |
---|
765 | enddo |
---|
766 | |
---|
767 | ! II.a.3. Surface pressure timeseries |
---|
768 | write(*,*) "Recomputing the surface pressure timeserie adapted to the new pressure..." |
---|
769 | do ig = 1,ngrid |
---|
770 | ps_timeseries(ig,:) = ps_timeseries(ig,:)*global_avg_press_new/global_avg_press_old |
---|
771 | enddo |
---|
772 | |
---|
773 | ! II.a.4. New pressure levels timeseries |
---|
774 | allocate(zplev_new_timeseries(ngrid,nlayer + 1,timelen)) |
---|
775 | write(*,*) "Recomputing the new pressure levels timeserie adapted to the new pressure..." |
---|
776 | do l = 1,nlayer + 1 |
---|
777 | do ig = 1,ngrid |
---|
778 | zplev_new_timeseries(ig,l,:) = ap(l) + bp(l)*ps_timeseries(ig,:) |
---|
779 | enddo |
---|
780 | enddo |
---|
781 | |
---|
782 | ! II.a.5. Tracers timeseries |
---|
783 | write(*,*) "Recomputing of tracer VMR timeseries for the new pressure..." |
---|
784 | |
---|
785 | l = 1 |
---|
786 | do ig = 1,ngrid |
---|
787 | do t = 1,timelen |
---|
788 | 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))/ & |
---|
789 | (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) |
---|
790 | if (q_h2o_PEM_phys(ig,t) < 0) then |
---|
791 | q_h2o_PEM_phys(ig,t) = 1.e-30 |
---|
792 | else if (q_h2o_PEM_phys(ig,t) > 1) then |
---|
793 | q_h2o_PEM_phys(ig,t) = 1. |
---|
794 | endif |
---|
795 | enddo |
---|
796 | enddo |
---|
797 | |
---|
798 | do ig = 1,ngrid |
---|
799 | do t = 1,timelen |
---|
800 | 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))/ & |
---|
801 | (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & |
---|
802 | + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) & |
---|
803 | - (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ & |
---|
804 | (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t)) |
---|
805 | if (q_co2_PEM_phys(ig,t) < 0) then |
---|
806 | q_co2_PEM_phys(ig,t) = 1.e-30 |
---|
807 | else if (q_co2_PEM_phys(ig,t) > 1) then |
---|
808 | q_co2_PEM_phys(ig,t) = 1. |
---|
809 | endif |
---|
810 | mmean=1/(A*q_co2_PEM_phys(ig,t) + B) |
---|
811 | vmr_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*mmean/m_co2 |
---|
812 | enddo |
---|
813 | enddo |
---|
814 | |
---|
815 | deallocate(zplev_new_timeseries,zplev_old_timeseries) |
---|
816 | |
---|
817 | !------------------------ |
---|
818 | ! II Run |
---|
819 | ! II_b Evolution of ice |
---|
820 | !------------------------ |
---|
821 | call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM) |
---|
822 | call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice) |
---|
823 | if (layering_algo) then |
---|
824 | do islope = 1,nslope |
---|
825 | do ig = 1,ngrid |
---|
826 | call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p) |
---|
827 | enddo |
---|
828 | enddo |
---|
829 | endif |
---|
830 | |
---|
831 | !------------------------ |
---|
832 | ! II Run |
---|
833 | ! II_c Flow of glaciers |
---|
834 | !------------------------ |
---|
835 | if (co2ice_flow .and. nslope > 1) call flow_co2glaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_PEM_phys,ps_timeseries, & |
---|
836 | global_avg_press_PCM,global_avg_press_new,co2_ice,flag_co2flow,flag_co2flow_mesh) |
---|
837 | if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow,flag_h2oflow_mesh) |
---|
838 | |
---|
839 | !------------------------ |
---|
840 | ! II Run |
---|
841 | ! II_d Update surface and soil temperatures |
---|
842 | !------------------------ |
---|
843 | ! II_d.1 Update Tsurf |
---|
844 | write(*,*) "Updating the new Tsurf" |
---|
845 | bool_sublim = .false. |
---|
846 | Tsurfavg_before_saved = tsurf_avg |
---|
847 | do ig = 1,ngrid |
---|
848 | do islope = 1,nslope |
---|
849 | if (is_co2ice_ini(ig,islope) .and. co2_ice(ig,islope) < 1.e-10 .and. .not. co2ice_disappeared(ig,islope)) then ! co2 ice disappeared, look for closest point without co2ice |
---|
850 | co2ice_disappeared(ig,islope) = .true. |
---|
851 | if (latitude_deg(ig) > 0) then |
---|
852 | do ig_loop = ig,ngrid |
---|
853 | do islope_loop = islope,iflat,-1 |
---|
854 | if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then |
---|
855 | tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop) |
---|
856 | bool_sublim = .true. |
---|
857 | exit |
---|
858 | endif |
---|
859 | enddo |
---|
860 | if (bool_sublim) exit |
---|
861 | enddo |
---|
862 | else |
---|
863 | do ig_loop = ig,1,-1 |
---|
864 | do islope_loop = islope,iflat |
---|
865 | if (.not. is_co2ice_ini(ig_loop,islope_loop) .and. co2_ice(ig_loop,islope_loop) < 1.e-10) then |
---|
866 | tsurf_avg(ig,islope) = tsurf_avg(ig_loop,islope_loop) |
---|
867 | bool_sublim = .true. |
---|
868 | exit |
---|
869 | endif |
---|
870 | enddo |
---|
871 | if (bool_sublim) exit |
---|
872 | enddo |
---|
873 | endif |
---|
874 | if ((co2_ice(ig,islope) < 1.e-10) .and. (h2o_ice(ig,islope) > frost_albedo_threshold)) then |
---|
875 | albedo(ig,1,islope) = albedo_h2o_frost |
---|
876 | albedo(ig,2,islope) = albedo_h2o_frost |
---|
877 | else |
---|
878 | albedo(ig,1,islope) = albedodat(ig) |
---|
879 | albedo(ig,2,islope) = albedodat(ig) |
---|
880 | emis(ig,islope) = emissiv |
---|
881 | endif |
---|
882 | else if ((co2_ice(ig,islope) > 1.e-3) .and. (d_co2ice(ig,islope) > 1.e-10)) then ! Put tsurf as tcond co2 |
---|
883 | ave = 0. |
---|
884 | do t = 1,timelen |
---|
885 | ave = ave + beta_clap_co2/(alpha_clap_co2 - log(vmr_co2_PEM_phys(ig,t)*ps_timeseries(ig,t)/100.)) |
---|
886 | enddo |
---|
887 | tsurf_avg(ig,islope) = ave/timelen |
---|
888 | endif |
---|
889 | enddo |
---|
890 | enddo |
---|
891 | |
---|
892 | do t = 1,timelen |
---|
893 | tsurf_PCM_timeseries(:,:,t) = tsurf_PCM_timeseries(:,:,t) + tsurf_avg - Tsurfavg_before_saved |
---|
894 | enddo |
---|
895 | ! for the start |
---|
896 | do ig = 1,ngrid |
---|
897 | do islope = 1,nslope |
---|
898 | tsurf(ig,islope) = tsurf(ig,islope) - (Tsurfavg_before_saved(ig,islope) - tsurf_avg(ig,islope)) |
---|
899 | enddo |
---|
900 | enddo |
---|
901 | |
---|
902 | if (soil_pem) then |
---|
903 | |
---|
904 | ! II_d.2 Update soil temperature |
---|
905 | write(*,*)"Updating soil temperature" |
---|
906 | allocate(Tsoil_locslope(ngrid,nsoilmx_PEM)) |
---|
907 | do islope = 1,nslope |
---|
908 | call compute_tsoil_pem(ngrid,nsoilmx_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) |
---|
909 | call compute_tsoil_pem(ngrid,nsoilmx_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_avg(:,islope),tsoil_PEM(:,:,islope)) |
---|
910 | |
---|
911 | do t = 1,timelen |
---|
912 | Tsoil_locslope(:,1:nsoilmx) = tsoil_PEM(:,1:nsoilmx,islope) + tsoil_anom(:,:,islope,t) |
---|
913 | Tsoil_locslope(:,nsoilmx + 1:) = tsoil_PEM(:,nsoilmx + 1:,islope) |
---|
914 | do ig = 1,ngrid |
---|
915 | do isoil = 1,nsoilmx_PEM |
---|
916 | watersoil_density_PEM_timeseries(ig,isoil,islope,t) = exp(beta_clap_h2o/Tsoil_locslope(ig,isoil) + alpha_clap_h2o)/Tsoil_locslope(ig,isoil)*mmol(igcm_h2o_vap)/(mugaz*r) |
---|
917 | if (isnan(tsoil_PEM(ig,isoil,islope))) call abort_pem("PEM - Update Tsoil","NaN detected in tsoil_PEM",1) |
---|
918 | enddo |
---|
919 | enddo |
---|
920 | enddo |
---|
921 | enddo |
---|
922 | watersoil_density_PEM_avg = sum(watersoil_density_PEM_timeseries,4)/timelen |
---|
923 | |
---|
924 | write(*,*) "Update of soil temperature done" |
---|
925 | |
---|
926 | deallocate(Tsoil_locslope) |
---|
927 | |
---|
928 | ! II_d.3 Update the ice table |
---|
929 | if (icetable_equilibrium) then |
---|
930 | write(*,*) "Compute ice table (equilibrium method)" |
---|
931 | icetable_thickness_old = icetable_thickness |
---|
932 | call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness) |
---|
933 | call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere |
---|
934 | else if (icetable_dynamic) then |
---|
935 | write(*,*) "Compute ice table (dynamic method)" |
---|
936 | ice_porefilling_old = ice_porefilling |
---|
937 | allocate(porefill(nsoilmx_PEM)) |
---|
938 | do ig = 1,ngrid |
---|
939 | do islope = 1,nslope |
---|
940 | call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),(/sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2)/),ice_porefilling(ig,:,islope),porefill,ssi_depth) |
---|
941 | icetable_depth(ig,islope) = ssi_depth |
---|
942 | ice_porefilling(ig,:,islope) = porefill |
---|
943 | enddo |
---|
944 | enddo |
---|
945 | deallocate(porefill) |
---|
946 | call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_old,ice_porefilling_old,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere |
---|
947 | endif |
---|
948 | |
---|
949 | ! II_d.4 Update the soil thermal properties |
---|
950 | call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) |
---|
951 | |
---|
952 | ! II_d.5 Update the mass of the regolith adsorbed |
---|
953 | if (adsorption_pem) then |
---|
954 | call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice, & |
---|
955 | h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & |
---|
956 | h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded) |
---|
957 | |
---|
958 | totmassco2_adsorbded = 0. |
---|
959 | totmassh2o_adsorbded = 0. |
---|
960 | do ig = 1,ngrid |
---|
961 | do islope = 1,nslope |
---|
962 | do l = 1,nsoilmx_PEM |
---|
963 | if (l == 1) then |
---|
964 | totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & |
---|
965 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
966 | totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* & |
---|
967 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
968 | else |
---|
969 | totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & |
---|
970 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
971 | totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & |
---|
972 | subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) |
---|
973 | endif |
---|
974 | enddo |
---|
975 | enddo |
---|
976 | enddo |
---|
977 | write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded |
---|
978 | write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded |
---|
979 | endif |
---|
980 | endif !soil_pem |
---|
981 | |
---|
982 | !------------------------ |
---|
983 | ! II Run |
---|
984 | ! II_e Outputs |
---|
985 | !------------------------ |
---|
986 | call writediagpem(ngrid,'ps_avg','Global average pressure','Pa',0,(/global_avg_press_new/)) |
---|
987 | do islope = 1,nslope |
---|
988 | write(str2(1:2),'(i2.2)') islope |
---|
989 | call writediagpem(ngrid,'h2o_ice_slope'//str2,'H2O ice','kg.m-2',2,h2o_ice(:,islope)) |
---|
990 | call writediagpem(ngrid,'co2_ice_slope'//str2,'CO2 ice','kg.m-2',2,co2_ice(:,islope)) |
---|
991 | call writediagpem(ngrid,'d_h2oice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,d_h2oice(:,islope)) |
---|
992 | call writediagpem(ngrid,'d_co2ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,d_co2ice(:,islope)) |
---|
993 | call writediagpem(ngrid,'Flow_co2ice_slope'//str2,'CO2 ice flow','Boolean',2,flag_co2flow(:,islope)) |
---|
994 | call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope)) |
---|
995 | if (icetable_equilibrium) then |
---|
996 | call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) |
---|
997 | call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table thickness','m',2,icetable_thickness(:,islope)) |
---|
998 | else if (icetable_dynamic) then |
---|
999 | call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope)) |
---|
1000 | endif |
---|
1001 | |
---|
1002 | if (soil_pem) then |
---|
1003 | call writediagsoilpem(ngrid,'tsoil_PEM_slope'//str2,'tsoil_PEM','K',3,tsoil_PEM(:,:,islope)) |
---|
1004 | call writediagsoilpem(ngrid,'inertiesoil_PEM_slope'//str2,'TI_PEM','K',3,TI_PEM(:,:,islope)) |
---|
1005 | if (icetable_dynamic) call writediagpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope)) |
---|
1006 | if (adsorption_pem) then |
---|
1007 | call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope)) |
---|
1008 | call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_phys(:,:,islope)) |
---|
1009 | endif |
---|
1010 | endif |
---|
1011 | enddo |
---|
1012 | |
---|
1013 | !------------------------ |
---|
1014 | ! II Run |
---|
1015 | ! II_f Update the tendencies |
---|
1016 | !------------------------ |
---|
1017 | call recomp_tend_co2_slope(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,global_avg_press_PCM,global_avg_press_new) |
---|
1018 | |
---|
1019 | !------------------------ |
---|
1020 | ! II Run |
---|
1021 | ! II_g Checking the stopping criterion |
---|
1022 | !------------------------ |
---|
1023 | |
---|
1024 | write(*,*) "Checking the stopping criteria..." |
---|
1025 | call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid) |
---|
1026 | call stopping_crit_co2(cell_area,co2ice_ini_surf,ini_co2ice_sublim,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope) |
---|
1027 | i_myear_leg = i_myear_leg + dt |
---|
1028 | i_myear = i_myear + dt |
---|
1029 | if (stopPEM <= 0 .and. i_myear_leg >= n_myear_leg) stopPEM = 5 |
---|
1030 | if (stopPEM <= 0 .and. i_myear >= n_myear) stopPEM = 6 |
---|
1031 | call system_clock(c2) |
---|
1032 | if (stopPEM <= 0 .and. timewall .and. real((c2 - c1)/cr) >= timelimit - antetime) stopPEM = 7 |
---|
1033 | if (stopPEM > 0) then |
---|
1034 | select case (stopPEM) |
---|
1035 | case(1) |
---|
1036 | write(*,*) "STOPPING because surface of h2o ice sublimating is too low:", stopPEM, "See message above." |
---|
1037 | case(2) |
---|
1038 | write(*,*) "STOPPING because tendencies on h2o ice = 0:", stopPEM, "See message above." |
---|
1039 | case(3) |
---|
1040 | write(*,*) "STOPPING because surface of co2 ice sublimating is too low:", stopPEM, "See message above." |
---|
1041 | case(4) |
---|
1042 | write(*,*) "STOPPING because surface global pressure changed too much:", stopPEM, "See message above." |
---|
1043 | case(5) |
---|
1044 | write(*,*) "STOPPING because maximum number of iterations is reached (possibly due to orbital parameters):", stopPEM |
---|
1045 | case(6) |
---|
1046 | write(*,*) "STOPPING because maximum number of Martian years to be simulated is reached:", stopPEM |
---|
1047 | case(7) |
---|
1048 | write(*,*) "STOPPING because the time limit for the PEM job will be reached soon:", stopPEM |
---|
1049 | case(8) |
---|
1050 | write(*,*) "STOPPING because the layering algorithm met an hasty end:", stopPEM |
---|
1051 | case default |
---|
1052 | write(*,*) "STOPPING with unknown stopping criterion code:", stopPEM |
---|
1053 | end select |
---|
1054 | exit |
---|
1055 | else |
---|
1056 | write(*,*) "The PEM can continue!" |
---|
1057 | write(*,*) "Number of iterations of the PEM: i_myear_leg =", i_myear_leg |
---|
1058 | write(*,*) "Number of simulated Martian years: i_myear =", i_myear |
---|
1059 | endif |
---|
1060 | |
---|
1061 | global_avg_press_old = global_avg_press_new |
---|
1062 | |
---|
1063 | enddo ! big time iteration loop of the pem |
---|
1064 | !------------------------------ END RUN -------------------------------- |
---|
1065 | |
---|
1066 | !------------------------------- OUTPUT -------------------------------- |
---|
1067 | !------------------------ |
---|
1068 | ! III Output |
---|
1069 | ! III_a Update surface value for the PCM start files |
---|
1070 | !------------------------ |
---|
1071 | ! III_a.1 Ice update (for startfi) |
---|
1072 | watercap = 0. |
---|
1073 | perennial_co2ice = co2_ice |
---|
1074 | do ig = 1,ngrid |
---|
1075 | ! H2O ice metamorphism |
---|
1076 | if (metam_h2oice .and. sum(qsurf(ig,igcm_h2o_ice,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_h2oice_threshold) then |
---|
1077 | h2o_ice(ig,:) = h2o_ice(ig,:) + qsurf(ig,igcm_h2o_ice,:) - metam_h2oice_threshold |
---|
1078 | qsurf(ig,igcm_h2o_ice,:) = metam_h2oice_threshold |
---|
1079 | endif |
---|
1080 | |
---|
1081 | ! Is H2O ice still considered as an infinite reservoir for the PCM? |
---|
1082 | if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then |
---|
1083 | ! There is enough ice to be considered as an infinite reservoir |
---|
1084 | watercaptag(ig) = .true. |
---|
1085 | else |
---|
1086 | ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost |
---|
1087 | watercaptag(ig) = .false. |
---|
1088 | qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:) |
---|
1089 | h2o_ice(ig,:) = 0. |
---|
1090 | endif |
---|
1091 | |
---|
1092 | ! CO2 ice metamorphism |
---|
1093 | if (metam_co2ice .and. sum(qsurf(ig,igcm_co2,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > metam_co2ice_threshold) then |
---|
1094 | perennial_co2ice(ig,:) = perennial_co2ice(ig,:) + qsurf(ig,igcm_co2,:) - metam_co2ice_threshold |
---|
1095 | qsurf(ig,igcm_co2,:) = metam_co2ice_threshold |
---|
1096 | endif |
---|
1097 | enddo |
---|
1098 | |
---|
1099 | ! III_a.2 Tsoil update (for startfi) |
---|
1100 | if (soil_pem) then |
---|
1101 | inertiesoil = TI_PEM(:,:nsoilmx,:) |
---|
1102 | tsoil = tsoil_PEM(:,1:nsoilmx,:) + tsoil_anom(:,:,:,timelen) |
---|
1103 | #ifndef CPP_STD |
---|
1104 | flux_geo = fluxgeo |
---|
1105 | #endif |
---|
1106 | endif |
---|
1107 | deallocate(tsoil_anom) |
---|
1108 | |
---|
1109 | ! III_a.4 Pressure (for start) |
---|
1110 | ps = ps*global_avg_press_new/global_avg_press_PCM |
---|
1111 | ps_start_PCM = ps_start_PCM*global_avg_press_new/global_avg_press_PCM |
---|
1112 | |
---|
1113 | ! III_a.5 Tracer (for start) |
---|
1114 | allocate(zplev_new(ngrid,nlayer + 1)) |
---|
1115 | |
---|
1116 | do l = 1,nlayer + 1 |
---|
1117 | zplev_new(:,l) = ap(l) + bp(l)*ps_start_PCM |
---|
1118 | enddo |
---|
1119 | |
---|
1120 | do nnq = 1,nqtot |
---|
1121 | if (noms(nnq) /= "co2") then |
---|
1122 | do l = 1,llm - 1 |
---|
1123 | do ig = 1,ngrid |
---|
1124 | q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) |
---|
1125 | enddo |
---|
1126 | q(:,llm,nnq) = q(:,llm - 1,nnq) |
---|
1127 | enddo |
---|
1128 | else |
---|
1129 | do l = 1,llm - 1 |
---|
1130 | do ig = 1,ngrid |
---|
1131 | q(ig,l,nnq) = q(ig,l,nnq)*(zplev_PCM(ig,l) - zplev_PCM(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) & |
---|
1132 | + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_PCM(ig,l) - zplev_PCM(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) |
---|
1133 | enddo |
---|
1134 | q(:,llm,nnq) = q(:,llm - 1,nnq) |
---|
1135 | enddo |
---|
1136 | endif |
---|
1137 | enddo |
---|
1138 | |
---|
1139 | ! Conserving the tracers mass for PCM start files |
---|
1140 | do nnq = 1,nqtot |
---|
1141 | do ig = 1,ngrid |
---|
1142 | do l = 1,llm - 1 |
---|
1143 | if (q(ig,l,nnq) > 1 .and. (noms(nnq) /= "dust_number") .and. (noms(nnq) /= "ccn_number") .and. (noms(nnq) /= "stormdust_number") .and. (noms(nnq) /= "topdust_number")) then |
---|
1144 | extra_mass = (q(ig,l,nnq) - 1)*(zplev_new(ig,l) - zplev_new(ig,l + 1)) |
---|
1145 | q(ig,l,nnq) = 1. |
---|
1146 | q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2)) |
---|
1147 | write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number" |
---|
1148 | endif |
---|
1149 | if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30 |
---|
1150 | enddo |
---|
1151 | enddo |
---|
1152 | enddo |
---|
1153 | |
---|
1154 | if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg) |
---|
1155 | |
---|
1156 | !------------------------ |
---|
1157 | ! III Output |
---|
1158 | ! III_b Write "restart.nc" and "restartfi.nc" |
---|
1159 | !------------------------ |
---|
1160 | ! III_b.1 Write "restart.nc" |
---|
1161 | ptimestep = iphysiq*daysec/real(day_step)/nsplit_phys ! dtphys/nsplit_phys |
---|
1162 | pday = day_ini |
---|
1163 | ztime_fin = time_phys |
---|
1164 | |
---|
1165 | allocate(p(ip1jmp1,nlayer + 1)) |
---|
1166 | #ifndef CPP_1D |
---|
1167 | ! Correction on teta due to surface pressure changes |
---|
1168 | do l = 1,nlayer |
---|
1169 | do i = 1,ip1jmp1 |
---|
1170 | teta(i,l)= teta(i,l)*(ps0(i)/ps(i))**kappa |
---|
1171 | enddo |
---|
1172 | enddo |
---|
1173 | ! Correction on atmospheric pressure |
---|
1174 | call pression(ip1jmp1,ap,bp,ps,p) |
---|
1175 | ! Correction on the mass of atmosphere |
---|
1176 | call massdair(p,masse) |
---|
1177 | call dynredem0("restart.nc",day_ini,phis) |
---|
1178 | call dynredem1("restart.nc",time_0,vcov,ucov,teta,q,masse,ps) |
---|
1179 | write(*,*) "restart.nc has been written" |
---|
1180 | #else |
---|
1181 | call writerestart1D('restart1D.txt',ps(1),tsurf(1,:),nlayer,size(tsurf,2),teta,ucov,vcov,nq,noms,qsurf(1,:,:),q) |
---|
1182 | write(*,*) "restart1D.txt has been written" |
---|
1183 | #endif |
---|
1184 | |
---|
1185 | ! III_b.2 Write the "restartfi.nc" |
---|
1186 | #ifndef CPP_STD |
---|
1187 | call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & |
---|
1188 | nlayer,nq,ptimestep,pday,0.,cell_area,albedodat, & |
---|
1189 | inertiedat,def_slope,subslope_dist) |
---|
1190 | call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & |
---|
1191 | ptimestep,ztime_fin,tsurf,tsoil,inertiesoil, & |
---|
1192 | albedo,emis,q2,qsurf,qsoil,tauscaling,totcloudfrac, & |
---|
1193 | wstar,watercap,perennial_co2ice) |
---|
1194 | #else |
---|
1195 | call physdem0("restartfi.nc",longitude,latitude,nsoilmx,ngrid, & |
---|
1196 | nlayer,nq,ptimestep,pday,time_phys,cell_area, & |
---|
1197 | albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe) |
---|
1198 | call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,nqsoil, & |
---|
1199 | ptimestep,ztime_fin,tsurf,tsoil,emis,q2,qsurf,qsoil, & |
---|
1200 | cloudfrac,totcloudfrac,hice,rnat,pctsrf_sic,tslab, & |
---|
1201 | tsea_ice,sea_ice) |
---|
1202 | #endif |
---|
1203 | write(*,*) "restartfi.nc has been written" |
---|
1204 | |
---|
1205 | !------------------------ |
---|
1206 | ! III Output |
---|
1207 | ! III_c Write the "restartpem.nc" |
---|
1208 | !------------------------ |
---|
1209 | if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings) |
---|
1210 | call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist) |
---|
1211 | call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, & |
---|
1212 | co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif) |
---|
1213 | write(*,*) "restartpem.nc has been written" |
---|
1214 | |
---|
1215 | call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear) |
---|
1216 | |
---|
1217 | write(*,*) "The PEM has run for", i_myear_leg, "Martian years." |
---|
1218 | write(*,*) "The chained simulation has run for", i_myear, "Martian years =", i_myear*convert_years, "Earth years." |
---|
1219 | write(*,*) "The reached date is now", (year_bp_ini + i_myear)*convert_years, "Earth years." |
---|
1220 | write(*,*) "PEM: so far, so good!" |
---|
1221 | |
---|
1222 | if (layering_algo) then |
---|
1223 | do islope = 1,nslope |
---|
1224 | do i = 1,ngrid |
---|
1225 | call del_layering(stratif(i,islope)) |
---|
1226 | enddo |
---|
1227 | enddo |
---|
1228 | deallocate(new_str,new_lag1,new_lag2,current1,current2) |
---|
1229 | endif |
---|
1230 | deallocate(vmr_co2_PCM,ps_timeseries,tsurf_PCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys) |
---|
1231 | deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved) |
---|
1232 | deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg) |
---|
1233 | deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim) |
---|
1234 | deallocate(icetable_thickness_old,ice_porefilling_old) |
---|
1235 | deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif) |
---|
1236 | !----------------------------- END OUTPUT ------------------------------ |
---|
1237 | |
---|
1238 | END PROGRAM pem |
---|