1 | module phyredem |
---|
2 | |
---|
3 | implicit none |
---|
4 | |
---|
5 | !======================================================================= |
---|
6 | contains |
---|
7 | !======================================================================= |
---|
8 | |
---|
9 | subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, & |
---|
10 | phystep,day_ini,time,airefi, & |
---|
11 | alb,ith,def_slope, & |
---|
12 | subslope_dist) |
---|
13 | |
---|
14 | ! create physics restart file and write time-independent variables |
---|
15 | use comsoil_h, only: inertiedat, volcapa, mlayer |
---|
16 | use geometry_mod, only: cell_area |
---|
17 | use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, & |
---|
18 | z0_default, albedice, emisice, emissiv, & |
---|
19 | iceradius, dtemisice, phisfi, z0, & |
---|
20 | hmons,summit,base,watercaptag |
---|
21 | use dimradmars_mod, only: tauvis |
---|
22 | use iostart, only : open_restartphy, close_restartphy, & |
---|
23 | put_var, put_field, length, ldscrpt, ndscrpt |
---|
24 | use mod_grid_phy_lmdz, only : klon_glo |
---|
25 | use planete_h, only: aphelie, emin_turb, lmixmin, obliquit, & |
---|
26 | peri_day, periheli, year_day |
---|
27 | use comcstfi_h, only: g, mugaz, omeg, rad, rcp |
---|
28 | use time_phylmdz_mod, only: daysec |
---|
29 | use comslope_mod, ONLY: nslope |
---|
30 | USE paleoclimate_mod, ONLY: paleoclimate, h2o_ice_depth, lag_co2_ice, d_coef |
---|
31 | |
---|
32 | implicit none |
---|
33 | |
---|
34 | character(len=*), intent(in) :: filename |
---|
35 | real,intent(in) :: lonfi(ngrid) |
---|
36 | real,intent(in) :: latfi(ngrid) |
---|
37 | integer,intent(in) :: nsoil |
---|
38 | integer,intent(in) :: ngrid |
---|
39 | integer,intent(in) :: nlay |
---|
40 | integer,intent(in) :: nq |
---|
41 | real,intent(in) :: phystep |
---|
42 | real,intent(in) :: day_ini |
---|
43 | real,intent(in) :: time |
---|
44 | real,intent(in) :: airefi(ngrid) |
---|
45 | real,intent(in) :: alb(ngrid) |
---|
46 | real,intent(in) :: ith(ngrid,nsoil) ! thermal inertia for present day |
---|
47 | real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes |
---|
48 | real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics |
---|
49 | |
---|
50 | character(ndscrpt), dimension(ldscrpt), parameter :: dscrpt_tab_cntrl = (/ & |
---|
51 | "(1) Total number of nodes on physics grid ", & |
---|
52 | "(2) Number of atmospheric layers ", & |
---|
53 | "(3) Initial day ", & |
---|
54 | "(4) Initial time of day ", & |
---|
55 | "(5) Radius of Mars (m) ~3397200 ", & |
---|
56 | "(6) Rotation rate (rad.s-1) ", & |
---|
57 | "(7) Gravity (m.s-2) ~3.72 ", & |
---|
58 | "(8) Molar mass of the atmosphere (g.mol-1) ~43.49", & |
---|
59 | "(9) = r/Cp ~0.256793 (=kappa in the dynamics) ", & |
---|
60 | "(10) Length of a sol (s) ~88775 ", & |
---|
61 | "(11) Time step in the physics ", & |
---|
62 | "(12) - ", & |
---|
63 | "(13) - ", & |
---|
64 | "(14) Length of year (sols) ~668.6 ", & |
---|
65 | "(15) Min. Sun-Mars distance (Mkm) ~206.66 ", & |
---|
66 | "(16) Max. Sun-Mars distance (Mkm) ~249.22 ", & |
---|
67 | "(17) Date of perihelion (sols since N. spring) ", & |
---|
68 | "(18) Obliquity of the planet (deg) ~23.98 ", & |
---|
69 | "(19) Default surface roughness (m) ~0.01 ", & |
---|
70 | "(20) Mixing length ~100 ", & |
---|
71 | "(21) Minimal energy ~1.e-8 ", & |
---|
72 | "(22) Albedo of northern cap ~0.5 ", & |
---|
73 | "(23) Albedo of southern cap ~0.5 ", & |
---|
74 | "(24) Emissivity of northern cap ~0.95 ", & |
---|
75 | "(25) Emissivity of southern cap ~0.95 ", & |
---|
76 | "(26) Emissivity of martian soil ~.95 ", & |
---|
77 | "(27) Mean visible optical depth ", & |
---|
78 | "(28) - ", & |
---|
79 | "(29) - ", & |
---|
80 | "(30) - ", & |
---|
81 | "(31) Mean scat radius of CO2 snow (north) ", & |
---|
82 | "(32) Mean scat radius of CO2 snow (south) ", & |
---|
83 | "(33) Time scale for snow metamorphism (north) ", & |
---|
84 | "(34) Time scale for snow metamorphism (south) ", & |
---|
85 | "(35) Soil volumetric heat capacity "/) |
---|
86 | real, dimension(length) :: tab_cntrl ! nb "length=100" defined in iostart module |
---|
87 | integer :: ig |
---|
88 | real, dimension(ngrid) :: watercaptag_tmp |
---|
89 | |
---|
90 | ! Create physics start file |
---|
91 | call open_restartphy(filename) |
---|
92 | |
---|
93 | ! Build tab_cntrl(:) array |
---|
94 | tab_cntrl(:)=0.0 |
---|
95 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
96 | ! Fill control array tab_cntrl(:) with parameters for this run |
---|
97 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
---|
98 | ! Informations on the physics grid |
---|
99 | tab_cntrl(1) = float(klon_glo) ! Total number of nodes on physics grid |
---|
100 | tab_cntrl(2) = float(nlay) ! Number of atmospheric layers |
---|
101 | tab_cntrl(3) = day_ini + int(time) ! Initial day |
---|
102 | tab_cntrl(4) = time -int(time) ! Initial time of day |
---|
103 | |
---|
104 | ! Informations about Mars, used by dynamics and physics |
---|
105 | tab_cntrl(5) = rad ! Radius of Mars (m) ~3397200 |
---|
106 | tab_cntrl(6) = omeg ! Rotation rate (rad.s-1) |
---|
107 | tab_cntrl(7) = g ! Gravity (m.s-2) ~3.72 |
---|
108 | tab_cntrl(8) = mugaz ! Molar mass of the atmosphere (g.mol-1) ~43.49 |
---|
109 | tab_cntrl(9) = rcp ! = r/cp ~0.256793 (=kappa dans dynamique) |
---|
110 | tab_cntrl(10) = daysec ! Length of a sol (s) ~88775 |
---|
111 | tab_cntrl(11) = phystep ! Time step in the physics |
---|
112 | |
---|
113 | ! Informations about Mars, only for physics |
---|
114 | tab_cntrl(14) = year_day ! Length of year (sols) ~668.6 |
---|
115 | tab_cntrl(15) = periheli ! Min. Sun-Mars distance (Mkm) ~206.66 |
---|
116 | tab_cntrl(16) = aphelie ! Max. SUn-Mars distance (Mkm) ~249.22 |
---|
117 | tab_cntrl(17) = peri_day ! Date of perihelion (sols since N. spring) |
---|
118 | tab_cntrl(18) = obliquit ! Obliquity of the planet (deg) ~23.98 |
---|
119 | |
---|
120 | ! Boundary layer and turbulence |
---|
121 | tab_cntrl(19) = z0_default ! Default surface roughness (m) ~0.01 |
---|
122 | tab_cntrl(20) = lmixmin ! Mixing length ~100 |
---|
123 | tab_cntrl(21) = emin_turb ! Minimal energy ~1.e-8 |
---|
124 | |
---|
125 | ! Optical properties of polar caps and ground emissivity |
---|
126 | tab_cntrl(22) = albedice(1) ! Albedo of northern cap ~0.5 |
---|
127 | tab_cntrl(23) = albedice(2) ! Albedo of southern cap ~0.5 |
---|
128 | tab_cntrl(24) = emisice(1) ! Emissivity of northern cap ~0.95 |
---|
129 | tab_cntrl(25) = emisice(2) ! Emissivity of southern cap ~0.95 |
---|
130 | tab_cntrl(26) = emissiv ! Emissivity of martian soil ~.95 |
---|
131 | tab_cntrl(31) = iceradius(1) ! Mean scat radius of CO2 snow (north) |
---|
132 | tab_cntrl(32) = iceradius(2) ! Mean scat radius of CO2 snow (south) |
---|
133 | tab_cntrl(33) = dtemisice(1) ! Time scale for snow metamorphism (north) |
---|
134 | tab_cntrl(34) = dtemisice(2) ! Time scale for snow metamorphism (south) |
---|
135 | |
---|
136 | ! dust aerosol properties |
---|
137 | tab_cntrl(27) = tauvis ! Mean visible optical depth |
---|
138 | |
---|
139 | ! Soil properties: |
---|
140 | tab_cntrl(35) = volcapa ! Soil volumetric heat capacity |
---|
141 | |
---|
142 | ! Write the controle array |
---|
143 | call put_var("controle","Control parameters",tab_cntrl) |
---|
144 | |
---|
145 | ! Write the controle array descriptor |
---|
146 | call put_var("controle_descriptor","Description of control parameters",dscrpt_tab_cntrl) |
---|
147 | |
---|
148 | ! Write the mid-layer depths |
---|
149 | call put_var("soildepth","Soil mid-layer depth",mlayer) |
---|
150 | |
---|
151 | ! Write longitudes |
---|
152 | call put_field("longitude","Longitudes of physics grid",lonfi) |
---|
153 | |
---|
154 | ! Write latitudes |
---|
155 | call put_field("latitude","Latitudes of physics grid",latfi) |
---|
156 | |
---|
157 | ! Write mesh areas |
---|
158 | call put_field("area","Mesh area",cell_area) |
---|
159 | |
---|
160 | ! Write surface geopotential |
---|
161 | call put_field("phisfi","Geopotential at the surface",phisfi) |
---|
162 | |
---|
163 | ! Write surface albedo |
---|
164 | call put_field("albedodat","Albedo of bare ground",alb) |
---|
165 | |
---|
166 | ! Subgrid topogaphy variables |
---|
167 | call put_field("ZMEA","Relief: mean relief",zmea) |
---|
168 | call put_field("ZSTD","Relief: standard deviation",zstd) |
---|
169 | call put_field("ZSIG","Relief: sigma parameter",zsig) |
---|
170 | call put_field("ZGAM","Relief: gamma parameter",zgam) |
---|
171 | call put_field("ZTHE","Relief: theta parameter",zthe) |
---|
172 | call put_field("hmons","Relief: hmons parameter (summit - base)",hmons) |
---|
173 | call put_field("summit","Relief: altitude from the aeroid of the summit of the highest subgrid topography",summit) |
---|
174 | call put_field("base","Relief: altitude from the aeroid of the base of the highest subgrid topography",base) |
---|
175 | |
---|
176 | ! Soil thermal inertia |
---|
177 | call put_field("inertiedat","Soil thermal inertia - present day",ith) |
---|
178 | |
---|
179 | ! Surface roughness length |
---|
180 | call put_field("z0","Surface roughness length",z0) |
---|
181 | |
---|
182 | ! Water cap |
---|
183 | watercaptag_tmp = 0 |
---|
184 | where (watercaptag) watercaptag_tmp = 1 |
---|
185 | call put_field("watercaptag","Infinite water reservoir",watercaptag_tmp) |
---|
186 | |
---|
187 | ! Sub grid scale slope parametrization |
---|
188 | call put_var("def_slope","slope criterium stages",def_slope) |
---|
189 | call put_field("subslope_dist","under mesh slope distribution",subslope_dist) |
---|
190 | |
---|
191 | ! Paleoclimate outputs |
---|
192 | if (paleoclimate) then |
---|
193 | call put_field("h2o_ice_depth","Depth to the fisrt H2O ice",h2o_ice_depth) |
---|
194 | call put_field("lag_co2_ice","Depth of the CO2 lags",lag_co2_ice) |
---|
195 | call put_field("d_coef","Diffusion coefficent",d_coef) |
---|
196 | endif |
---|
197 | |
---|
198 | ! Close file |
---|
199 | call close_restartphy |
---|
200 | |
---|
201 | end subroutine physdem0 |
---|
202 | |
---|
203 | !======================================================================= |
---|
204 | |
---|
205 | subroutine physdem1(filename,nsoil,ngrid,nlay,nq,nqsoil, & |
---|
206 | phystep,time,tsurf,tsoil,inertiesoil, & |
---|
207 | albedo,emis,q2,qsurf,qsoil,& |
---|
208 | tauscaling,totcloudfrac,wstar, & |
---|
209 | watercap,perennial_co2ice) |
---|
210 | ! write time-dependent variable to restart file |
---|
211 | use iostart, only : open_restartphy, close_restartphy, & |
---|
212 | put_var, put_field |
---|
213 | use tracer_mod, only: noms ! tracer names |
---|
214 | use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd |
---|
215 | use compute_dtau_mod, only: dtau |
---|
216 | use dust_rad_adjust_mod, only: dust_rad_adjust_prev,dust_rad_adjust_next |
---|
217 | use dust_param_mod, only: dustscaling_mode |
---|
218 | use comsoil_h,only: flux_geo,adsorption_soil,igcm_h2o_vap_soil, & |
---|
219 | igcm_h2o_ice_soil,igcm_h2o_vap_ads |
---|
220 | use comslope_mod, only: nslope |
---|
221 | use paleoclimate_mod, only: paleoclimate |
---|
222 | implicit none |
---|
223 | |
---|
224 | include "callkeys.h" |
---|
225 | |
---|
226 | character(len=*),intent(in) :: filename |
---|
227 | integer,intent(in) :: nsoil |
---|
228 | integer,intent(in) :: ngrid |
---|
229 | integer,intent(in) :: nlay |
---|
230 | integer,intent(in) :: nq |
---|
231 | integer,intent(in) :: nqsoil |
---|
232 | real,intent(in) :: phystep |
---|
233 | real,intent(in) :: time |
---|
234 | real,intent(in) :: tsurf(ngrid,nslope) |
---|
235 | real,intent(in) :: tsoil(ngrid,nsoil,nslope) |
---|
236 | real,intent(in) :: inertiesoil(ngrid,nsoil,nslope) |
---|
237 | real,intent(in) :: albedo(ngrid,2,nslope) |
---|
238 | real,intent(in) :: emis(ngrid,nslope) |
---|
239 | real,intent(in) :: q2(ngrid,nlay+1) |
---|
240 | real,intent(in) :: qsurf(ngrid,nq,nslope) |
---|
241 | real,intent(in) :: qsoil(ngrid,nsoil,nqsoil,nslope) |
---|
242 | real,intent(in) :: tauscaling(ngrid) |
---|
243 | real,intent(in) :: totcloudfrac(ngrid) |
---|
244 | real,intent(in) :: wstar(ngrid) |
---|
245 | real,intent(in) :: watercap(ngrid,nslope) |
---|
246 | real,intent(in) :: perennial_co2ice(ngrid,nslope) |
---|
247 | |
---|
248 | integer :: iq |
---|
249 | character(len=30) :: txt ! to store some text |
---|
250 | ! indexes of water vapour & water ice tracers (if any): |
---|
251 | integer :: i_h2o_vap=0 |
---|
252 | integer :: i_h2o_ice=0 |
---|
253 | integer :: i_hdo_vap=0 |
---|
254 | integer :: i_hdo_ice=0 |
---|
255 | |
---|
256 | |
---|
257 | ! Open file |
---|
258 | call open_restartphy(filename) |
---|
259 | |
---|
260 | ! First variable to write must be "Time", in order to correctly |
---|
261 | ! set time counter in file |
---|
262 | call put_var("Time","Temps de simulation",time) |
---|
263 | |
---|
264 | ! Water ice layer |
---|
265 | call put_field("watercap","H2O ice cover",watercap,time) |
---|
266 | |
---|
267 | ! Perennial CO2 ice layer |
---|
268 | if (paleoclimate) call put_field("perennial_co2ice","CO2 ice cover",perennial_co2ice,time) |
---|
269 | |
---|
270 | ! Surface temperature |
---|
271 | call put_field("tsurf","Surface temperature",tsurf,time) |
---|
272 | |
---|
273 | ! Soil temperature |
---|
274 | call put_field("inertiesoil","Soil thermal inertia",inertiesoil,time) |
---|
275 | |
---|
276 | ! Soil temperature |
---|
277 | call put_field("tsoil","Soil temperature",tsoil,time) |
---|
278 | |
---|
279 | ! Albedo of the surface |
---|
280 | call put_field("albedo","Surface albedo",albedo(:,1,:),time) |
---|
281 | |
---|
282 | ! Emissivity of the surface |
---|
283 | call put_field("emis","Surface emissivity",emis,time) |
---|
284 | |
---|
285 | ! Planetary Boundary Layer |
---|
286 | call put_field("q2","pbl wind variance",q2,time) |
---|
287 | |
---|
288 | ! Sub-grid cloud fraction |
---|
289 | call put_field("totcloudfrac","Total cloud fraction",totcloudfrac,time) |
---|
290 | |
---|
291 | ! Dust conversion factor |
---|
292 | ! Only to be read by newstart to convert to actual dust values |
---|
293 | ! Or by any user who wants to reconstruct dust, opacity from the start files. |
---|
294 | call put_field("tauscaling","dust conversion factor",tauscaling,time) |
---|
295 | |
---|
296 | ! Radiative scaling coefficients |
---|
297 | if (dustscaling_mode==2) then |
---|
298 | call put_field("dust_rad_adjust_prev", & |
---|
299 | "radiative dust adjustement factor prev. sol", & |
---|
300 | dust_rad_adjust_prev,time) |
---|
301 | call put_field("dust_rad_adjust_next", & |
---|
302 | "radiative dust adjustement factor next sol", & |
---|
303 | dust_rad_adjust_next,time) |
---|
304 | endif |
---|
305 | |
---|
306 | if (dustinjection.gt.0) then |
---|
307 | call put_field("dtau","dust opacity difference between GCM and scenario",& |
---|
308 | dtau,time) |
---|
309 | endif |
---|
310 | |
---|
311 | if (calltherm) then |
---|
312 | call put_field("wstar","Max vertical velocity in thermals",wstar,time) |
---|
313 | endif |
---|
314 | |
---|
315 | ! Tracers on the surface |
---|
316 | ! preliminary stuff: look for water vapour & water ice tracers (if any) |
---|
317 | do iq=1,nq |
---|
318 | if (noms(iq).eq."h2o_vap") then |
---|
319 | i_h2o_vap=iq |
---|
320 | endif |
---|
321 | if (noms(iq).eq."h2o_ice") then |
---|
322 | i_h2o_ice=iq |
---|
323 | endif |
---|
324 | |
---|
325 | ! look for HDO vapour & HDO ice tracers (if any) |
---|
326 | if (noms(iq).eq."hdo_vap") then |
---|
327 | i_hdo_vap=iq |
---|
328 | endif |
---|
329 | if (noms(iq).eq."hdo_ice") then |
---|
330 | i_hdo_ice=iq |
---|
331 | endif |
---|
332 | enddo |
---|
333 | |
---|
334 | |
---|
335 | if (nq.gt.0) then |
---|
336 | do iq=1,nq |
---|
337 | txt=noms(iq) |
---|
338 | ! Exception: there is no water vapour surface tracer |
---|
339 | if (txt.eq."h2o_vap") then |
---|
340 | write(*,*)"physdem1: skipping water vapour tracer" |
---|
341 | if (i_h2o_ice.eq.i_h2o_vap) then |
---|
342 | ! then there is no "water ice" tracer; but still |
---|
343 | ! there is some water ice on the surface |
---|
344 | write(*,*)" writing water ice instead" |
---|
345 | txt="h2o_ice" |
---|
346 | else |
---|
347 | ! there is a "water ice" tracer which has been / will be |
---|
348 | ! delt with in due time |
---|
349 | cycle |
---|
350 | endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap) |
---|
351 | endif ! of if (txt.eq."h2o_vap") |
---|
352 | |
---|
353 | if (txt.eq."hdo_vap") then |
---|
354 | write(*,*)"physdem1: skipping HDO vapour tracer" |
---|
355 | if (i_hdo_ice.eq.i_hdo_vap) then |
---|
356 | ! then there is no "water ice" tracer; but still |
---|
357 | ! there is some water ice on the surface |
---|
358 | write(*,*)" writing HDO ice instead" |
---|
359 | txt="hdo_ice" |
---|
360 | else |
---|
361 | ! there is a "water ice" tracer which has been / will be |
---|
362 | ! delt with in due time |
---|
363 | cycle |
---|
364 | endif ! of if (igcm_hdo_ice.eq.igcm_hdo_vap) |
---|
365 | endif ! of if (txt.eq."hdo_vap") |
---|
366 | |
---|
367 | ! co2_ice has been added to qsurf(:,igcm_co2) in co2condens4micro |
---|
368 | if (txt.eq."co2_ice") then |
---|
369 | write(*,*)"physdem1: skipping co2_ice tracer" |
---|
370 | cycle |
---|
371 | end if |
---|
372 | |
---|
373 | call put_field(trim(txt),"tracer on surface",qsurf(:,iq,:),time) |
---|
374 | enddo |
---|
375 | endif |
---|
376 | |
---|
377 | ! Non-orographic gavity waves |
---|
378 | if (calllott_nonoro) then |
---|
379 | call put_field("du_nonoro_gwd","Zonal wind tendency due to GW",du_nonoro_gwd,time) |
---|
380 | call put_field("dv_nonoro_gwd","Meridional wind tendency due to GW",dv_nonoro_gwd,time) |
---|
381 | endif |
---|
382 | |
---|
383 | ! Geothermal Flux |
---|
384 | call put_field('flux_geo','Geothermal flux',flux_geo,time) |
---|
385 | |
---|
386 | ! Adsorption |
---|
387 | if (adsorption_soil) then |
---|
388 | call put_field("h2o_vap_soil","subsurface water vapour", & |
---|
389 | qsoil(:,:,igcm_h2o_vap_soil,:), time) |
---|
390 | call put_field("h2o_ice_soil","subsurface water ice", & |
---|
391 | qsoil(:,:,igcm_h2o_ice_soil,:), time) |
---|
392 | call put_field("h2o_vap_ads", "adsorbed water", & |
---|
393 | qsoil(:,:,igcm_h2o_vap_ads,:), time) |
---|
394 | endif |
---|
395 | |
---|
396 | ! Close file |
---|
397 | call close_restartphy |
---|
398 | |
---|
399 | end subroutine physdem1 |
---|
400 | |
---|
401 | end module phyredem |
---|