1 | MODULE physiq_mod |
---|
2 | |
---|
3 | IMPLICIT NONE |
---|
4 | |
---|
5 | CONTAINS |
---|
6 | |
---|
7 | SUBROUTINE physiq( |
---|
8 | $ ngrid,nlayer,nq |
---|
9 | $ ,firstcall,lastcall |
---|
10 | $ ,pday,ptime,ptimestep |
---|
11 | $ ,pplev,pplay,pphi |
---|
12 | $ ,pu,pv,pt,pq |
---|
13 | $ ,flxw |
---|
14 | $ ,pdu,pdv,pdt,pdq,pdpsrf) |
---|
15 | |
---|
16 | use watercloud_mod, only: watercloud, zdqcloud, zdqscloud |
---|
17 | use calchim_mod, only: calchim, ichemistry, zdqchim, zdqschim |
---|
18 | use watersat_mod, only: watersat |
---|
19 | use co2condens_mod, only: co2condens |
---|
20 | use co2cloud_mod, only: co2cloud |
---|
21 | use callradite_mod, only: callradite |
---|
22 | use callsedim_mod, only: callsedim |
---|
23 | use rocketduststorm_mod, only: rocketduststorm, dustliftday |
---|
24 | use calcstormfract_mod, only: calcstormfract |
---|
25 | use topmons_mod, only: topmons,topmons_setup |
---|
26 | use nltecool_mod, only: nltecool |
---|
27 | use nlte_tcool_mod, only: nlte_tcool |
---|
28 | use blendrad_mod, only: blendrad |
---|
29 | use nlthermeq_mod, only: nlthermeq |
---|
30 | use thermosphere_mod, only: thermosphere |
---|
31 | use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice, |
---|
32 | & igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice, |
---|
33 | & igcm_hdo_vap, igcm_hdo_ice, |
---|
34 | & igcm_ccn_mass, igcm_ccn_number, |
---|
35 | & igcm_ccnco2_mass, igcm_ccnco2_number, |
---|
36 | & igcm_ccnco2_h2o_mass_ice, |
---|
37 | & igcm_ccnco2_h2o_mass_ccn, |
---|
38 | & igcm_ccnco2_h2o_number, |
---|
39 | & igcm_ccnco2_meteor_mass, |
---|
40 | & igcm_ccnco2_meteor_number, |
---|
41 | & rho_ice_co2,nuiceco2_sed,nuiceco2_ref, |
---|
42 | & igcm_dust_mass, igcm_dust_number, igcm_h2o2, |
---|
43 | & nuice_ref, rho_ice, rho_dust, ref_r0, |
---|
44 | & igcm_he, igcm_stormdust_mass, |
---|
45 | & igcm_stormdust_number, igcm_topdust_mass, |
---|
46 | & igcm_topdust_number, |
---|
47 | & qperemin |
---|
48 | use comsoil_h, only: inertiedat, inertiesoil,! dat: soil thermal inertia for present climate, inertiesoil is the TI read in the start |
---|
49 | & tsoil, nsoilmx,!number of subsurface layers |
---|
50 | & mlayer,layer ! soil mid layer depths |
---|
51 | use geometry_mod, only: longitude, latitude, cell_area, |
---|
52 | & longitude_deg |
---|
53 | use comgeomfi_h, only: sinlon, coslon, sinlat, coslat |
---|
54 | use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, |
---|
55 | & zthe, z0, albedo_h2o_cap,albedo_h2o_frost, |
---|
56 | & frost_albedo_threshold,frost_metam_threshold, |
---|
57 | & tsurf, emis, |
---|
58 | & capcal, fluxgrd, qsurf, |
---|
59 | & hmons,summit,base,watercap,watercaptag, |
---|
60 | & perenial_co2ice |
---|
61 | use comsaison_h, only: dist_sol, declin, zls, |
---|
62 | & mu0, fract, local_time |
---|
63 | use nirdata_mod, only: NIR_leedat |
---|
64 | use nirco2abs_mod, only: nirco2abs |
---|
65 | use slope_mod, only: theta_sl, psi_sl |
---|
66 | use conc_mod, only: rnew, cpnew, mmean |
---|
67 | use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec |
---|
68 | use dimradmars_mod, only: aerosol, totcloudfrac, |
---|
69 | & dtrad, fluxrad_sky, fluxrad, albedo, |
---|
70 | & naerkind, iaer_dust_doubleq, |
---|
71 | & iaer_stormdust_doubleq, iaer_h2o_ice, |
---|
72 | & flux_1AU |
---|
73 | use dust_param_mod, only: doubleq, lifting, callddevil, |
---|
74 | & tauscaling, odpref, dustbin, |
---|
75 | & dustscaling_mode, dust_rad_adjust, |
---|
76 | & freedust, reff_driven_IRtoVIS_scenario |
---|
77 | use turb_mod, only: q2, wstar, ustar, sensibFlux, |
---|
78 | & zmax_th, hfmax_th, turb_resolved |
---|
79 | use planete_h, only: aphelie, periheli, year_day, peri_day, |
---|
80 | & obliquit |
---|
81 | USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad |
---|
82 | USE calldrag_noro_mod, ONLY: calldrag_noro |
---|
83 | USE vdifc_mod, ONLY: vdifc |
---|
84 | use param_v4_h, only: nreact,n_avog, |
---|
85 | & fill_data_thermos, allocate_param_thermos |
---|
86 | use iono_h, only: allocate_param_iono |
---|
87 | use compute_dtau_mod, only: compute_dtau |
---|
88 | use nonoro_gwd_ran_mod, only: nonoro_gwd_ran |
---|
89 | use check_fields_mod, only: check_physics_fields |
---|
90 | #ifdef MESOSCALE |
---|
91 | use comsoil_h, only: mlayer,layer |
---|
92 | use surfdat_h, only: z0_default |
---|
93 | use comm_wrf |
---|
94 | #else |
---|
95 | USE planetwide_mod, ONLY: planetwide_maxval, planetwide_minval, |
---|
96 | & planetwide_sumval |
---|
97 | use phyredem, only: physdem0, physdem1 |
---|
98 | use phyetat0_mod, only: phyetat0, tab_cntrl_mod |
---|
99 | use wstats_mod, only: callstats, wstats, mkstats |
---|
100 | use eofdump_mod, only: eofdump |
---|
101 | USE vertical_layers_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt |
---|
102 | USE mod_phys_lmdz_omp_data, ONLY: is_omp_master |
---|
103 | USE time_phylmdz_mod, ONLY: day_end |
---|
104 | #endif |
---|
105 | |
---|
106 | #ifdef CPP_XIOS |
---|
107 | use xios_output_mod, only: initialize_xios_output, |
---|
108 | & update_xios_timestep, |
---|
109 | & send_xios_field |
---|
110 | use wxios, only: wxios_context_init, xios_context_finalize |
---|
111 | #endif |
---|
112 | USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured |
---|
113 | use ioipsl_getin_p_mod, only: getin_p |
---|
114 | use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, |
---|
115 | & subslope_dist,iflat,sky_slope, |
---|
116 | & major_slope,compute_meshgridavg, |
---|
117 | & ini_comslope_h |
---|
118 | use write_output_mod, only: write_output |
---|
119 | IMPLICIT NONE |
---|
120 | c======================================================================= |
---|
121 | c |
---|
122 | c subject: |
---|
123 | c -------- |
---|
124 | c |
---|
125 | c Organisation of the physical parametrisations of the LMD |
---|
126 | c martian atmospheric general circulation model. |
---|
127 | c |
---|
128 | c The GCM can be run without or with tracer transport |
---|
129 | c depending on the value of Logical "tracer" in file "callphys.def" |
---|
130 | c Tracers may be water vapor, ice OR chemical species OR dust particles |
---|
131 | c |
---|
132 | c SEE comments in initracer.F about numbering of tracer species... |
---|
133 | c |
---|
134 | c It includes: |
---|
135 | c |
---|
136 | c 1. Initialization: |
---|
137 | c 1.1 First call initializations |
---|
138 | c 1.2 Initialization for every call to physiq |
---|
139 | c 1.2.5 Compute mean mass and cp, R and thermal conduction coeff. |
---|
140 | c 2. Compute radiative transfer tendencies |
---|
141 | c (longwave and shortwave) for CO2 and aerosols. |
---|
142 | c 3. Gravity wave and subgrid scale topography drag : |
---|
143 | c 4. Vertical diffusion (turbulent mixing): |
---|
144 | c 5. Convective adjustment |
---|
145 | c 6. Condensation and sublimation of carbon dioxide. |
---|
146 | c 7. TRACERS : |
---|
147 | c 7a. water, water ice, co2 ice (clouds) |
---|
148 | c 7b. call for photochemistry when tracers are chemical species |
---|
149 | c 7c. other scheme for tracer (dust) transport (lifting, sedimentation) |
---|
150 | c 7d. updates (CO2 pressure variations, surface budget) |
---|
151 | c 8. Contribution to tendencies due to thermosphere |
---|
152 | c 9. Surface and sub-surface temperature calculations |
---|
153 | c 10. Write outputs : |
---|
154 | c - "startfi", "histfi" (if it's time) |
---|
155 | c - Saving statistics (if "callstats = .true.") |
---|
156 | c - Dumping eof (if "calleofdump = .true.") |
---|
157 | c - Output any needed variables in "diagfi" |
---|
158 | c 11. Diagnostic: mass conservation of tracers |
---|
159 | c |
---|
160 | c author: |
---|
161 | c ------- |
---|
162 | c Frederic Hourdin 15/10/93 |
---|
163 | c Francois Forget 1994 |
---|
164 | c Christophe Hourdin 02/1997 |
---|
165 | c Subroutine completly rewritten by F.Forget (01/2000) |
---|
166 | c Introduction of the photochemical module: S. Lebonnois (11/2002) |
---|
167 | c Introduction of the thermosphere module: M. Angelats i Coll (2002) |
---|
168 | c Water ice clouds: Franck Montmessin (update 06/2003) |
---|
169 | c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009) |
---|
170 | c Nb: See callradite.F for more information. |
---|
171 | c Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags |
---|
172 | c jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization |
---|
173 | c |
---|
174 | c 10/16 J. Audouard: modifications for CO2 clouds scheme |
---|
175 | |
---|
176 | c arguments: |
---|
177 | c ---------- |
---|
178 | c |
---|
179 | c input: |
---|
180 | c ------ |
---|
181 | c ecri period (in dynamical timestep) to write output |
---|
182 | c ngrid Size of the horizontal grid. |
---|
183 | c All internal loops are performed on that grid. |
---|
184 | c nlayer Number of vertical layers. |
---|
185 | c nq Number of advected fields |
---|
186 | c firstcall True at the first call |
---|
187 | c lastcall True at the last call |
---|
188 | c pday Number of days counted from the North. Spring |
---|
189 | c equinoxe. |
---|
190 | c ptime Universal time (0<ptime<1): ptime=0.5 at 12:00 UT |
---|
191 | c ptimestep timestep (s) |
---|
192 | c pplay(ngrid,nlayer) Pressure at the middle of the layers (Pa) |
---|
193 | c pplev(ngrid,nlayer+1) intermediate pressure levels (pa) |
---|
194 | c pphi(ngrid,nlayer) Geopotential at the middle of the layers (m2s-2) |
---|
195 | c pu(ngrid,nlayer) u component of the wind (ms-1) |
---|
196 | c pv(ngrid,nlayer) v component of the wind (ms-1) |
---|
197 | c pt(ngrid,nlayer) Temperature (K) |
---|
198 | c pq(ngrid,nlayer,nq) Advected fields |
---|
199 | c pudyn(ngrid,nlayer) | |
---|
200 | c pvdyn(ngrid,nlayer) | Dynamical temporal derivative for the |
---|
201 | c ptdyn(ngrid,nlayer) | corresponding variables |
---|
202 | c pqdyn(ngrid,nlayer,nq) | |
---|
203 | c flxw(ngrid,nlayer) vertical mass flux (kg/s) at layer lower boundary |
---|
204 | c |
---|
205 | c output: |
---|
206 | c ------- |
---|
207 | c |
---|
208 | c pdu(ngrid,nlayer) | |
---|
209 | c pdv(ngrid,nlayer) | Temporal derivative of the corresponding |
---|
210 | c pdt(ngrid,nlayer) | variables due to physical processes. |
---|
211 | c pdq(ngrid,nlayer,nq) | |
---|
212 | c pdpsrf(ngrid) | |
---|
213 | |
---|
214 | c |
---|
215 | c======================================================================= |
---|
216 | c |
---|
217 | c 0. Declarations : |
---|
218 | c ------------------ |
---|
219 | |
---|
220 | include "callkeys.h" |
---|
221 | include "netcdf.inc" |
---|
222 | |
---|
223 | c Arguments : |
---|
224 | c ----------- |
---|
225 | |
---|
226 | c inputs: |
---|
227 | c ------- |
---|
228 | INTEGER,INTENT(in) :: ngrid ! number of atmospheric columns |
---|
229 | INTEGER,INTENT(in) :: nlayer ! number of atmospheric layers |
---|
230 | INTEGER,INTENT(in) :: nq ! number of tracers |
---|
231 | LOGICAL,INTENT(in) :: firstcall ! signals first call to physics |
---|
232 | LOGICAL,INTENT(in) :: lastcall ! signals last call to physics |
---|
233 | REAL,INTENT(in) :: pday ! number of elapsed sols since reference Ls=0 |
---|
234 | REAL,INTENT(in) :: ptime ! "universal time", given as fraction of sol (e.g.: 0.5 for noon) |
---|
235 | REAL,INTENT(in) :: ptimestep ! physics timestep (s) |
---|
236 | REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) |
---|
237 | REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa) |
---|
238 | REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential at mid-layer (m2s-2) |
---|
239 | REAL,INTENT(in) :: pu(ngrid,nlayer) ! zonal wind component (m/s) |
---|
240 | REAL,INTENT(in) :: pv(ngrid,nlayer) ! meridional wind component (m/s) |
---|
241 | REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K) |
---|
242 | REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air) |
---|
243 | REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s) |
---|
244 | ! at lower boundary of layer |
---|
245 | |
---|
246 | c outputs: |
---|
247 | c -------- |
---|
248 | c physical tendencies |
---|
249 | REAL,INTENT(out) :: pdu(ngrid,nlayer) ! zonal wind tendency (m/s/s) |
---|
250 | REAL,INTENT(out) :: pdv(ngrid,nlayer) ! meridional wind tendency (m/s/s) |
---|
251 | REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s) |
---|
252 | REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s) |
---|
253 | REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s) |
---|
254 | |
---|
255 | |
---|
256 | |
---|
257 | c Local saved variables: |
---|
258 | c ---------------------- |
---|
259 | INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0) |
---|
260 | INTEGER,SAVE :: icount ! counter of calls to physiq during the run. |
---|
261 | REAL,SAVE :: time_phys |
---|
262 | |
---|
263 | !$OMP THREADPRIVATE(day_ini,icount,time_phys) |
---|
264 | |
---|
265 | #ifdef DUSTSTORM |
---|
266 | REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb |
---|
267 | #endif |
---|
268 | |
---|
269 | c Variables used by the water ice microphysical scheme: |
---|
270 | REAL rice(ngrid,nlayer) ! Water ice geometric mean radius (m) |
---|
271 | REAL nuice(ngrid,nlayer) ! Estimated effective variance |
---|
272 | ! of the size distribution |
---|
273 | real rsedcloud(ngrid,nlayer) ! Cloud sedimentation radius |
---|
274 | real rhocloud(ngrid,nlayer) ! Cloud density (kg.m-3) |
---|
275 | real rsedcloudco2(ngrid,nlayer) ! CO2 Cloud sedimentation radius |
---|
276 | real rhocloudco2(ngrid,nlayer) ! CO2 Cloud density (kg.m-3) |
---|
277 | real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the |
---|
278 | ! size distribution |
---|
279 | REAL inertiesoil_tifeedback(ngrid,nsoilmx,nslope) ! Time varying subsurface |
---|
280 | ! thermal inertia (J.s-1/2.m-2.K-1) |
---|
281 | ! (used only when tifeedback=.true.) |
---|
282 | c Variables used by the CO2 clouds microphysical scheme: |
---|
283 | DOUBLE PRECISION riceco2(ngrid,nlayer) ! co2 ice geometric mean radius (m) |
---|
284 | real zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1) |
---|
285 | real zdqssed_ccn(ngrid,nq) ! CCN flux at the surface (kg.m-2.s-1) |
---|
286 | real zcondicea_co2microp(ngrid,nlayer) |
---|
287 | c Variables used by the photochemistry |
---|
288 | REAL surfdust(ngrid,nlayer) ! dust surface area (m2/m3, if photochemistry) |
---|
289 | REAL surfice(ngrid,nlayer) ! ice surface area (m2/m3, if photochemistry) |
---|
290 | c Variables used by the slope model |
---|
291 | REAL sl_ls, sl_lct, sl_lat |
---|
292 | REAL sl_tau, sl_alb, sl_the, sl_psi |
---|
293 | REAL sl_fl0, sl_flu |
---|
294 | REAL sl_ra, sl_di0 |
---|
295 | REAL sky |
---|
296 | REAL fluxsurf_dir_dn_sw(ngrid) ! Incident direct solar flux on Mars at surface (W.m-2) |
---|
297 | |
---|
298 | REAL,PARAMETER :: stephan = 5.67e-08 ! Stephan Boltzman constant |
---|
299 | |
---|
300 | c Local variables : |
---|
301 | c ----------------- |
---|
302 | |
---|
303 | REAL CBRT |
---|
304 | EXTERNAL CBRT |
---|
305 | |
---|
306 | ! CHARACTER*80 fichier |
---|
307 | INTEGER l,ig,ierr,igout,iq,tapphys,isoil |
---|
308 | |
---|
309 | REAL fluxsurf_lw(ngrid,nslope) !incident LW (IR) surface flux (W.m-2) |
---|
310 | REAL fluxsurf_dn_sw(ngrid,2,nslope) ! Incident SW (solar) surface flux (W.m-2) |
---|
311 | REAL fluxsurf_up_sw(ngrid,2) ! Reflected SW (solar) surface flux (W.m-2) |
---|
312 | REAL fluxtop_lw(ngrid) !Outgoing LW (IR) flux to space (W.m-2) |
---|
313 | REAL fluxtop_dn_sw(ngrid,2) ! Incoming SW (solar) flux from space (W.m-2) |
---|
314 | REAL fluxtop_up_sw(ngrid,2) ! Outgoing SW (solar) flux to space (W.m-2) |
---|
315 | REAL tau_pref_scenario(ngrid) ! prescribed dust column visible opacity |
---|
316 | ! at odpref |
---|
317 | REAL IRtoVIScoef(ngrid) ! conversion coefficient to apply on |
---|
318 | ! scenario absorption IR (9.3um) CDOD |
---|
319 | ! = tau_pref_gcm_VIS / tau_pref_gcm_IR |
---|
320 | REAL tau_pref_gcm(ngrid) ! dust column visible opacity at odpref in the GCM |
---|
321 | c rocket dust storm |
---|
322 | REAL totstormfract(ngrid) ! fraction of the mesh where the dust storm is contained |
---|
323 | logical clearatm ! clearatm used to calculate twice the radiative |
---|
324 | ! transfer when rdstorm is active : |
---|
325 | ! - in a mesh with stormdust and background dust (false) |
---|
326 | ! - in a mesh with background dust only (true) |
---|
327 | c entrainment by mountain top dust flows |
---|
328 | logical nohmons ! nohmons used to calculate twice the radiative |
---|
329 | ! transfer when topflows is active : |
---|
330 | ! - in a mesh with topdust and background dust (false) |
---|
331 | ! - in a mesh with background dust only (true) |
---|
332 | |
---|
333 | REAL tau(ngrid,naerkind) ! Column dust optical depth at each point |
---|
334 | ! AS: TBD: this one should be in a module ! |
---|
335 | REAL zday ! date (time since Ls=0, in martian days) |
---|
336 | REAL zzlay(ngrid,nlayer) ! altitude at the middle of the layers |
---|
337 | REAL zzlev(ngrid,nlayer+1) ! altitude at layer boundaries |
---|
338 | ! REAL latvl1,lonvl1 ! Viking Lander 1 point (for diagnostic) |
---|
339 | |
---|
340 | c Tendancies due to various processes: |
---|
341 | REAL dqsurf(ngrid,nq,nslope) ! tendency for tracers on surface (Kg/m2/s) |
---|
342 | REAL zdtlw(ngrid,nlayer) ! (K/s) |
---|
343 | REAL zdtsw(ngrid,nlayer) ! (K/s) |
---|
344 | REAL pdqrds(ngrid,nlayer,nq) ! tendency for dust after rocketduststorm |
---|
345 | |
---|
346 | REAL zdtnirco2(ngrid,nlayer) ! (K/s) |
---|
347 | REAL zdtnlte(ngrid,nlayer) ! (K/s) |
---|
348 | REAL zdtsurf(ngrid,nslope) ! (K/s) |
---|
349 | REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer) |
---|
350 | REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer) ! (m.s-2) |
---|
351 | REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid,nslope) ! (K/s) |
---|
352 | REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer) ! (m.s-2) |
---|
353 | REAL zdhadj(ngrid,nlayer) ! (K/s) |
---|
354 | REAL zdtgw(ngrid,nlayer) ! (K/s) |
---|
355 | REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer) ! (m.s-2) |
---|
356 | REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid,nslope) |
---|
357 | REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer) |
---|
358 | |
---|
359 | REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq,nslope) |
---|
360 | REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq) |
---|
361 | REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq) |
---|
362 | REAL zdqadj(ngrid,nlayer,nq) |
---|
363 | REAL zdqc(ngrid,nlayer,nq) |
---|
364 | REAL zdqcloudco2(ngrid,nlayer,nq) |
---|
365 | REAL zdqsc(ngrid,nq,nslope) |
---|
366 | |
---|
367 | REAL zdteuv(ngrid,nlayer) ! (K/s) |
---|
368 | REAL zdtconduc(ngrid,nlayer) ! (K/s) |
---|
369 | REAL zdumolvis(ngrid,nlayer) |
---|
370 | REAL zdvmolvis(ngrid,nlayer) |
---|
371 | real zdqmoldiff(ngrid,nlayer,nq) |
---|
372 | real*8 PhiEscH,PhiEscH2,PhiEscD |
---|
373 | |
---|
374 | REAL dwatercap(ngrid,nslope), dwatercap_dif(ngrid,nslope) ! (kg/m-2) |
---|
375 | |
---|
376 | c Local variable for local intermediate calcul: |
---|
377 | REAL zflubid(ngrid,nslope) |
---|
378 | REAL zplanck(ngrid),zpopsk(ngrid,nlayer) |
---|
379 | REAL zdum1(ngrid,nlayer) |
---|
380 | REAL zdum2(ngrid,nlayer) |
---|
381 | REAL ztim1,ztim2,ztim3, z1,z2 |
---|
382 | REAL ztime_fin |
---|
383 | REAL zdh(ngrid,nlayer) |
---|
384 | REAL zh(ngrid,nlayer) ! potential temperature (K) |
---|
385 | REAL pw(ngrid,nlayer) ! vertical velocity (m/s) (>0 when downwards) |
---|
386 | INTEGER length |
---|
387 | PARAMETER (length=100) |
---|
388 | |
---|
389 | c Variables for the total dust for diagnostics |
---|
390 | REAL qdusttotal(ngrid,nlayer) !it equals to dust + stormdust |
---|
391 | |
---|
392 | INTEGER iaer |
---|
393 | |
---|
394 | c local variables only used for diagnostic (output in file "diagfi" or "stats") |
---|
395 | c ----------------------------------------------------------------------------- |
---|
396 | REAL ps(ngrid), zt(ngrid,nlayer) |
---|
397 | REAL zu(ngrid,nlayer),zv(ngrid,nlayer) |
---|
398 | REAL zq(ngrid,nlayer,nq) |
---|
399 | |
---|
400 | REAL fluxtop_dn_sw_tot(ngrid), fluxtop_up_sw_tot(ngrid) |
---|
401 | REAL fluxsurf_dn_sw_tot(ngrid,nslope), fluxsurf_up_sw_tot(ngrid) |
---|
402 | character*2 str2 |
---|
403 | ! character*5 str5 |
---|
404 | real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer) |
---|
405 | real rdust(ngrid,nlayer) ! dust geometric mean radius (m) |
---|
406 | real rstormdust(ngrid,nlayer) ! stormdust geometric mean radius (m) |
---|
407 | real rtopdust(ngrid,nlayer) ! topdust geometric mean radius (m) |
---|
408 | integer igmin, lmin |
---|
409 | logical tdiag |
---|
410 | |
---|
411 | real co2col(ngrid) ! CO2 column |
---|
412 | ! pplev and pplay are dynamical inputs and must not be modified in the physics. |
---|
413 | ! instead, use zplay and zplev : |
---|
414 | REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer) |
---|
415 | ! REAL zstress(ngrid),cd |
---|
416 | real tmean, zlocal(nlayer) |
---|
417 | real rho(ngrid,nlayer) ! density |
---|
418 | real vmr(ngrid,nlayer) ! volume mixing ratio |
---|
419 | real rhopart(ngrid,nlayer) ! number density of a given species |
---|
420 | real colden(ngrid,nq) ! vertical column of tracers |
---|
421 | real mass(nq) ! global mass of tracers (g) |
---|
422 | REAL mtot(ngrid) ! Total mass of water vapor (kg/m2) |
---|
423 | REAL mstormdtot(ngrid) ! Total mass of stormdust tracer (kg/m2) |
---|
424 | REAL mdusttot(ngrid) ! Total mass of dust tracer (kg/m2) |
---|
425 | REAL icetot(ngrid) ! Total mass of water ice (kg/m2) |
---|
426 | REAL mtotco2(ngrid) ! Total mass of co2, including ice at the surface (kg/m2) |
---|
427 | REAL vaptotco2(ngrid) ! Total mass of co2 vapor (kg/m2) |
---|
428 | REAL icetotco2(ngrid) ! Total mass of co2 ice (kg/m2) |
---|
429 | REAL Nccntot(ngrid) ! Total number of ccn (nbr/m2) |
---|
430 | REAL NccnCO2tot(ngrid) ! Total number of ccnCO2 (nbr/m2) |
---|
431 | REAL Mccntot(ngrid) ! Total mass of ccn (kg/m2) |
---|
432 | REAL rave(ngrid) ! Mean water ice effective radius (m) |
---|
433 | REAL opTES(ngrid,nlayer) ! abs optical depth at 825 cm-1 |
---|
434 | REAL tauTES(ngrid) ! column optical depth at 825 cm-1 |
---|
435 | REAL Qabsice ! Water ice absorption coefficient |
---|
436 | REAL taucloudtes(ngrid) ! Cloud opacity at infrared |
---|
437 | ! reference wavelength using |
---|
438 | ! Qabs instead of Qext |
---|
439 | ! (direct comparison with TES) |
---|
440 | REAL mtotD(ngrid) ! Total mass of HDO vapor (kg/m2) |
---|
441 | REAL icetotD(ngrid) ! Total mass of HDO ice (kg/m2) |
---|
442 | REAL DoH_vap(ngrid,nlayer) !D/H ratio |
---|
443 | REAL DoH_ice(ngrid,nlayer) !D/H ratio |
---|
444 | REAL DoH_surf(ngrid) !D/H ratio surface |
---|
445 | |
---|
446 | REAL dqdustsurf(ngrid) ! surface q dust flux (kg/m2/s) |
---|
447 | REAL dndustsurf(ngrid) ! surface n dust flux (number/m2/s) |
---|
448 | REAL ndust(ngrid,nlayer) ! true n dust (kg/kg) |
---|
449 | REAL qdust(ngrid,nlayer) ! true q dust (kg/kg) |
---|
450 | REAL nccn(ngrid,nlayer) ! true n ccn (kg/kg) |
---|
451 | REAL qccn(ngrid,nlayer) ! true q ccn (kg/kg) |
---|
452 | c definition tendancies of stormdust tracers |
---|
453 | REAL rdsdqdustsurf(ngrid) ! surface q stormdust flux (kg/m2/s) |
---|
454 | REAL rdsdndustsurf(ngrid) ! surface n stormdust flux (number/m2/s) |
---|
455 | REAL rdsndust(ngrid,nlayer) ! true n stormdust (kg/kg) |
---|
456 | REAL rdsqdust(ngrid,nlayer) ! true q stormdust (kg/kg) |
---|
457 | REAL wspeed(ngrid,nlayer+1) ! vertical velocity stormdust tracer |
---|
458 | REAL wtop(ngrid,nlayer+1) ! vertical velocity topdust tracer |
---|
459 | |
---|
460 | REAL dsodust(ngrid,nlayer) ! density scaled opacity for background dust |
---|
461 | REAL dsords(ngrid,nlayer) ! density scaled opacity for stormdust |
---|
462 | REAL dsotop(ngrid,nlayer) ! density scaled opacity for topdust |
---|
463 | |
---|
464 | c Test 1d/3d scavenging |
---|
465 | real h2otot(ngrid) |
---|
466 | real hdotot(ngrid) |
---|
467 | REAL satu(ngrid,nlayer) ! satu ratio for output |
---|
468 | REAL zqsat(ngrid,nlayer) ! saturation |
---|
469 | REAL satuco2(ngrid,nlayer) ! co2 satu ratio for output |
---|
470 | REAL zqsatco2(ngrid,nlayer) ! saturation co2 |
---|
471 | |
---|
472 | |
---|
473 | ! Added for new NLTE scheme |
---|
474 | |
---|
475 | real co2vmr_gcm(ngrid,nlayer) |
---|
476 | real n2vmr_gcm(ngrid,nlayer) |
---|
477 | real ovmr_gcm(ngrid,nlayer) |
---|
478 | real covmr_gcm(ngrid,nlayer) |
---|
479 | integer ierr_nlte |
---|
480 | real*8 varerr |
---|
481 | |
---|
482 | C Non-oro GW drag & Calcul of Brunt-Vaisala freq. (BV2) |
---|
483 | REAL ztetalev(ngrid,nlayer) |
---|
484 | real zdtetalev(ngrid,nlayer), zdzlev(ngrid,nlayer) |
---|
485 | REAL bv2(ngrid,nlayer) ! BV2 at zlev |
---|
486 | c Non-oro GW tendencies |
---|
487 | REAL d_u_hin(ngrid,nlayer), d_v_hin(ngrid,nlayer) |
---|
488 | REAL d_t_hin(ngrid,nlayer) |
---|
489 | c Diagnostics 2D of gw_nonoro |
---|
490 | REAL zustrhi(ngrid), zvstrhi(ngrid) |
---|
491 | c Variables for PBL |
---|
492 | REAL zz1(ngrid) |
---|
493 | REAL lmax_th_out(ngrid) |
---|
494 | REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer) |
---|
495 | REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq) |
---|
496 | INTEGER lmax_th(ngrid),dimout,n_out,n |
---|
497 | CHARACTER(50) zstring |
---|
498 | REAL dtke_th(ngrid,nlayer+1) |
---|
499 | REAL zcdv(ngrid), zcdh(ngrid) |
---|
500 | REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out |
---|
501 | REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out |
---|
502 | REAL u_out1(ngrid) |
---|
503 | REAL T_out1(ngrid) |
---|
504 | REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters] |
---|
505 | REAL tstar(ngrid) ! friction velocity and friction potential temp |
---|
506 | REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid) |
---|
507 | real qdustrds0(ngrid,nlayer),qdustrds1(ngrid,nlayer) |
---|
508 | real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer) |
---|
509 | real qdusttotal0(ngrid),qdusttotal1(ngrid) |
---|
510 | |
---|
511 | c sub-grid scale water ice clouds (A. Pottier 2013) |
---|
512 | logical clearsky |
---|
513 | ! flux for the part without clouds |
---|
514 | real zdtswclf(ngrid,nlayer) |
---|
515 | real zdtlwclf(ngrid,nlayer) |
---|
516 | real fluxsurf_lwclf(ngrid) |
---|
517 | real fluxsurf_dn_swclf(ngrid,2),fluxsurf_up_swclf(ngrid,2) |
---|
518 | real fluxtop_lwclf(ngrid) |
---|
519 | real fluxtop_dn_swclf(ngrid,2),fluxtop_up_swclf(ngrid,2) |
---|
520 | real taucloudtesclf(ngrid) |
---|
521 | real tf_clf, ntf_clf ! tf: fraction of clouds, ntf: fraction without clouds |
---|
522 | real rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m) |
---|
523 | C test de conservation de la masse de CO2 |
---|
524 | REAL co2totA |
---|
525 | REAL co2totB |
---|
526 | REAL co2conservation |
---|
527 | |
---|
528 | c entrainment by mountain top dust flows above sub-grid scale topography |
---|
529 | REAL pdqtop(ngrid,nlayer,nq) ! tendency for dust after topmons |
---|
530 | |
---|
531 | c when no startfi file is asked for init |
---|
532 | real alpha,lay1 ! coefficients for building layers |
---|
533 | integer iloop |
---|
534 | |
---|
535 | ! flags to trigger extra sanity checks |
---|
536 | logical,save :: check_physics_inputs=.false. |
---|
537 | logical,save :: check_physics_outputs=.false. |
---|
538 | |
---|
539 | !$OMP THREADPRIVATE(check_physics_inputs,check_physics_outputs) |
---|
540 | |
---|
541 | c Sub-grid scale slopes |
---|
542 | real :: tsurf_meshavg(ngrid) ! Surface temperature grid box averaged [K] |
---|
543 | real :: albedo_meshavg(ngrid,2) ! albedo temperature grid box averaged [1] |
---|
544 | real :: emis_meshavg(ngrid,2) ! emis temperature grid box averaged [1] |
---|
545 | real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2] |
---|
546 | real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie |
---|
547 | integer :: islope |
---|
548 | real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting |
---|
549 | |
---|
550 | logical :: write_restart |
---|
551 | |
---|
552 | ! Variable for ice table |
---|
553 | REAL :: rhowater_surf(ngrid,nslope) ! Water density at the surface [kg/m^3] |
---|
554 | REAL :: rhowater_surf_sat(ngrid,nslope) ! Water density at the surface at saturation [kg/m^3] |
---|
555 | REAL :: rhowater_soil(ngrid,nsoilmx,nslope) ! Water density in soil layers [kg/m^3] |
---|
556 | REAL,PARAMETER :: alpha_clap_h2o = 28.9074 ! Coeff for Clapeyron law [/] |
---|
557 | REAL,PARAMETER :: beta_clap_h2o = -6143.7 ! Coeff for Clapeyron law [K] |
---|
558 | REAL :: pvap_surf(ngrid) ! Water vapor partial pressure in first layer [Pa] |
---|
559 | REAL,PARAMETER :: m_co2 = 44.01E-3 ! CO2 molecular mass [kg/mol] |
---|
560 | REAL,PARAMETER :: m_noco2 = 33.37E-3 ! Non condensible mol mass [kg/mol] |
---|
561 | REAL :: ztmp1,ztmp2 ! intermediate variables to compute the mean molar mass of the layer |
---|
562 | |
---|
563 | c======================================================================= |
---|
564 | pdq(:,:,:) = 0. |
---|
565 | |
---|
566 | c 1. Initialisation: |
---|
567 | c ----------------- |
---|
568 | c 1.1 Initialisation only at first call |
---|
569 | c --------------------------------------- |
---|
570 | |
---|
571 | IF (firstcall) THEN |
---|
572 | |
---|
573 | call getin_p("check_physics_inputs",check_physics_inputs) |
---|
574 | call getin_p("check_physics_outputs",check_physics_outputs) |
---|
575 | |
---|
576 | c variables set to 0 |
---|
577 | c ~~~~~~~~~~~~~~~~~~ |
---|
578 | aerosol(:,:,:)=0 |
---|
579 | dtrad(:,:)=0 |
---|
580 | |
---|
581 | #ifndef MESOSCALE |
---|
582 | fluxrad(:,:)=0 |
---|
583 | wstar(:)=0. |
---|
584 | #endif |
---|
585 | |
---|
586 | #ifdef CPP_XIOS |
---|
587 | ! Initialize XIOS context |
---|
588 | write(*,*) "physiq: call wxios_context_init" |
---|
589 | CALL wxios_context_init |
---|
590 | #endif |
---|
591 | |
---|
592 | c read startfi |
---|
593 | c ~~~~~~~~~~~~ |
---|
594 | #ifndef MESOSCALE |
---|
595 | |
---|
596 | ! GCM. Read netcdf initial physical parameters. |
---|
597 | CALL phyetat0 ("startfi.nc",0,0, |
---|
598 | & nsoilmx,ngrid,nlayer,nq, |
---|
599 | & day_ini,time_phys, |
---|
600 | & tsurf,tsoil,albedo,emis, |
---|
601 | & q2,qsurf,tauscaling,totcloudfrac,wstar, |
---|
602 | & watercap,perenial_co2ice, |
---|
603 | & def_slope,def_slope_mean,subslope_dist) |
---|
604 | |
---|
605 | ! Sky view: |
---|
606 | DO islope=1,nslope |
---|
607 | sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2. |
---|
608 | END DO |
---|
609 | ! Determine the 'flatest' slopes |
---|
610 | iflat = 1 |
---|
611 | DO islope=2,nslope |
---|
612 | IF(abs(def_slope_mean(islope)).lt. |
---|
613 | & abs(def_slope_mean(iflat)))THEN |
---|
614 | iflat = islope |
---|
615 | ENDIF |
---|
616 | ENDDO |
---|
617 | PRINT*,'Flat slope for islope = ',iflat |
---|
618 | PRINT*,'corresponding criterium = ',def_slope_mean(iflat) |
---|
619 | |
---|
620 | #else |
---|
621 | ! MESOSCALE. Supposedly everything is already set in modules. |
---|
622 | ! So we just check. And we fill day_ini |
---|
623 | print*,"check: --- in physiq.F" |
---|
624 | print*,"check: rad,cpp,g,r,rcp,daysec" |
---|
625 | print*,rad,cpp,g,r,rcp,daysec |
---|
626 | PRINT*,'check: tsurf ',tsurf(1,:),tsurf(ngrid,:) |
---|
627 | PRINT*,'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:) |
---|
628 | PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx) |
---|
629 | PRINT*,'check: midlayer,layer ', mlayer(:),layer(:) |
---|
630 | PRINT*,'check: tracernames ', noms |
---|
631 | PRINT*,'check: emis ',emis(1,:),emis(ngrid,:) |
---|
632 | PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1) |
---|
633 | PRINT*,'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:) |
---|
634 | PRINT*,'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:) |
---|
635 | !!! |
---|
636 | day_ini = pday |
---|
637 | !!! a couple initializations (dummy for mesoscale) done in phyetat0 |
---|
638 | !!! --- maybe this should be done in update_inputs_physiq_mod |
---|
639 | |
---|
640 | tauscaling(:)=1.0 !! probably important |
---|
641 | totcloudfrac(:)=1.0 |
---|
642 | DO islope = 1,nslope |
---|
643 | albedo(:,1,islope)=albedodat(:) |
---|
644 | albedo(:,2,islope)=albedo(:,1,islope) |
---|
645 | inertiesoil(:,:,islope) = inertiedat(:,:) |
---|
646 | watercap(:,:)=0.0 |
---|
647 | ENDDO |
---|
648 | #endif |
---|
649 | #ifndef MESOSCALE |
---|
650 | if (.not.startphy_file) then |
---|
651 | ! starting without startfi.nc and with callsoil |
---|
652 | ! is not yet possible as soildepth default is not defined |
---|
653 | if (callsoil) then |
---|
654 | ! default mlayer distribution, following a power law: |
---|
655 | ! mlayer(k)=lay1*alpha**(k-1/2) |
---|
656 | lay1=2.e-4 |
---|
657 | alpha=2 |
---|
658 | do iloop=0,nsoilmx-1 |
---|
659 | mlayer(iloop)=lay1*(alpha**(iloop-0.5)) |
---|
660 | enddo |
---|
661 | lay1=sqrt(mlayer(0)*mlayer(1)) |
---|
662 | alpha=mlayer(1)/mlayer(0) |
---|
663 | do iloop=1,nsoilmx |
---|
664 | layer(iloop)=lay1*(alpha**(iloop-1)) |
---|
665 | enddo |
---|
666 | endif |
---|
667 | ! additionnal "academic" initialization of physics |
---|
668 | do islope = 1,nslope |
---|
669 | tsurf(:,islope)=pt(:,1) |
---|
670 | enddo |
---|
671 | write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!" |
---|
672 | do isoil=1,nsoilmx |
---|
673 | tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:) |
---|
674 | enddo |
---|
675 | write(*,*) "Physiq: initializing inertiedat !!" |
---|
676 | inertiedat(:,:)=400. |
---|
677 | inertiesoil(:,:,:)=400. |
---|
678 | write(*,*) "Physiq: initializing day_ini to pdat !" |
---|
679 | day_ini=pday |
---|
680 | endif |
---|
681 | #endif |
---|
682 | if (pday.ne.day_ini) then |
---|
683 | write(*,*) "PHYSIQ: ERROR: bad synchronization between ", |
---|
684 | & "physics and dynamics" |
---|
685 | write(*,*) "dynamics day [pday]: ",pday |
---|
686 | write(*,*) "physics day [day_ini]: ",day_ini |
---|
687 | call abort_physic("physiq","dynamics day /= physics day",1) |
---|
688 | endif |
---|
689 | |
---|
690 | write (*,*) 'In physiq day_ini =', day_ini |
---|
691 | |
---|
692 | c initialize tracers |
---|
693 | c ~~~~~~~~~~~~~~~~~~ |
---|
694 | CALL initracer(ngrid,nq,qsurf) |
---|
695 | |
---|
696 | c Initialize albedo and orbital calculation |
---|
697 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
698 | CALL surfini(ngrid,nslope,qsurf) |
---|
699 | CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit) |
---|
700 | c initialize soil |
---|
701 | c ~~~~~~~~~~~~~~~ |
---|
702 | IF (callsoil) THEN |
---|
703 | c Thermal inertia feedback: |
---|
704 | IF (tifeedback) THEN |
---|
705 | DO islope = 1,nslope |
---|
706 | CALL soil_tifeedback(ngrid,nsoilmx, |
---|
707 | s qsurf(:,:,islope), |
---|
708 | s inertiesoil_tifeedback(:,:,islope)) |
---|
709 | ENDDO |
---|
710 | CALL soil(ngrid,nsoilmx,firstcall, |
---|
711 | s inertiesoil_tifeedback, |
---|
712 | s ptimestep,tsurf,tsoil,capcal,fluxgrd) |
---|
713 | ELSE |
---|
714 | CALL soil(ngrid,nsoilmx,firstcall,inertiesoil, |
---|
715 | s ptimestep,tsurf,tsoil,capcal,fluxgrd) |
---|
716 | ENDIF ! of IF (tifeedback) |
---|
717 | ELSE |
---|
718 | PRINT*, |
---|
719 | & 'PHYSIQ WARNING! Thermal conduction in the soil turned off' |
---|
720 | DO ig=1,ngrid |
---|
721 | capcal(ig,:)=1.e5 |
---|
722 | fluxgrd(ig,:)=0. |
---|
723 | ENDDO |
---|
724 | ENDIF |
---|
725 | icount=1 |
---|
726 | |
---|
727 | #ifndef MESOSCALE |
---|
728 | c Initialize thermospheric parameters |
---|
729 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
730 | |
---|
731 | if (callthermos) then |
---|
732 | call fill_data_thermos |
---|
733 | call allocate_param_thermos(nlayer) |
---|
734 | call allocate_param_iono(nlayer,nreact) |
---|
735 | call param_read_e107 |
---|
736 | endif |
---|
737 | #endif |
---|
738 | c Initialize R and Cp as constant |
---|
739 | |
---|
740 | if (.not.callthermos .and. .not.photochem) then |
---|
741 | do l=1,nlayer |
---|
742 | do ig=1,ngrid |
---|
743 | rnew(ig,l)=r |
---|
744 | cpnew(ig,l)=cpp |
---|
745 | mmean(ig,l)=mugaz |
---|
746 | enddo |
---|
747 | enddo |
---|
748 | endif |
---|
749 | |
---|
750 | if(callnlte.and.nltemodel.eq.2) call nlte_setup |
---|
751 | if(callnirco2.and.nircorr.eq.1) call NIR_leedat |
---|
752 | |
---|
753 | |
---|
754 | IF (water.AND.(ngrid.NE.1)) THEN |
---|
755 | write(*,*)"physiq: water_param Surface water frost albedo:", |
---|
756 | . albedo_h2o_frost |
---|
757 | write(*,*)"physiq: water_param Surface watercap albedo:", |
---|
758 | . albedo_h2o_cap |
---|
759 | ENDIF |
---|
760 | |
---|
761 | #ifndef MESOSCALE |
---|
762 | if(ngrid.ne.1) then |
---|
763 | ! no need to compute slopes when in 1D; it is an input |
---|
764 | if (callslope) call getslopes(ngrid,phisfi) |
---|
765 | endif !1D |
---|
766 | if (ecritstart.GT.0) then |
---|
767 | call physdem0("restartfi.nc",longitude,latitude, |
---|
768 | & nsoilmx,ngrid,nlayer,nq, |
---|
769 | & ptimestep,pday,0.,cell_area, |
---|
770 | & albedodat,inertiedat,def_slope, |
---|
771 | & subslope_dist) |
---|
772 | else |
---|
773 | call physdem0("restartfi.nc",longitude,latitude, |
---|
774 | & nsoilmx,ngrid,nlayer,nq, |
---|
775 | & ptimestep,float(day_end),0.,cell_area, |
---|
776 | & albedodat,inertiedat,def_slope, |
---|
777 | & subslope_dist) |
---|
778 | endif |
---|
779 | |
---|
780 | c Initialize mountain mesh fraction for the entrainment by top flows param. |
---|
781 | c ~~~~~~~~~~~~~~~ |
---|
782 | if (topflows) call topmons_setup(ngrid) |
---|
783 | |
---|
784 | #endif |
---|
785 | |
---|
786 | #ifdef CPP_XIOS |
---|
787 | ! XIOS outputs |
---|
788 | write(*,*) "physiq firstcall: call initialize_xios_output" |
---|
789 | call initialize_xios_output(pday,ptime,ptimestep,daysec, |
---|
790 | & presnivs,pseudoalt,mlayer) |
---|
791 | #endif |
---|
792 | ENDIF ! (end of "if firstcall") |
---|
793 | |
---|
794 | if (check_physics_inputs) then |
---|
795 | ! Check the validity of input fields coming from the dynamics |
---|
796 | call check_physics_fields("begin physiq:",pt,pu,pv,pplev,pq) |
---|
797 | endif |
---|
798 | |
---|
799 | c --------------------------------------------------- |
---|
800 | c 1.2 Initializations done at every physical timestep: |
---|
801 | c --------------------------------------------------- |
---|
802 | c |
---|
803 | |
---|
804 | #ifdef CPP_XIOS |
---|
805 | ! update XIOS time/calendar |
---|
806 | call update_xios_timestep |
---|
807 | #endif |
---|
808 | |
---|
809 | |
---|
810 | |
---|
811 | |
---|
812 | |
---|
813 | c Initialize various variables |
---|
814 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
815 | pdv(:,:)=0 |
---|
816 | pdu(:,:)=0 |
---|
817 | pdt(:,:)=0 |
---|
818 | pdq(:,:,:)=0 |
---|
819 | pdpsrf(:)=0 |
---|
820 | zflubid(:,:)=0 |
---|
821 | zdtsurf(:,:)=0 |
---|
822 | dqsurf(:,:,:)=0 |
---|
823 | dsodust(:,:)=0. |
---|
824 | dsords(:,:)=0. |
---|
825 | dsotop(:,:)=0. |
---|
826 | dwatercap(:,:)=0 |
---|
827 | |
---|
828 | call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf, |
---|
829 | & albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg) |
---|
830 | |
---|
831 | ! Dust scenario conversion coefficient from IRabs to VISext |
---|
832 | IRtoVIScoef(1:ngrid)=2.6 ! initialized with former value from Montabone et al 2015 |
---|
833 | ! recomputed in aeropacity if reff_driven_IRtoVIS_scenario=.true. |
---|
834 | |
---|
835 | #ifdef DUSTSTORM |
---|
836 | pq_tmp(:,:,:)=0 |
---|
837 | #endif |
---|
838 | igout=ngrid/2+1 |
---|
839 | |
---|
840 | |
---|
841 | zday=pday+ptime ! compute time, in sols (and fraction thereof) |
---|
842 | ! Compute local time at each grid point |
---|
843 | DO ig=1,ngrid |
---|
844 | local_time(ig)=modulo(1.+(zday-INT(zday)) |
---|
845 | & +(longitude_deg(ig)/15)/24,1.) |
---|
846 | ENDDO |
---|
847 | c Compute Solar Longitude (Ls) : |
---|
848 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
849 | if (season) then |
---|
850 | call solarlong(zday,zls) |
---|
851 | else |
---|
852 | call solarlong(float(day_ini),zls) |
---|
853 | end if |
---|
854 | |
---|
855 | c Initialize pressure levels |
---|
856 | c ~~~~~~~~~~~~~~~~~~ |
---|
857 | zplev(:,:) = pplev(:,:) |
---|
858 | zplay(:,:) = pplay(:,:) |
---|
859 | ps(:) = pplev(:,1) |
---|
860 | |
---|
861 | c Compute geopotential at interlayers |
---|
862 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
863 | c ponderation des altitudes au niveau des couches en dp/p |
---|
864 | |
---|
865 | DO l=1,nlayer |
---|
866 | DO ig=1,ngrid |
---|
867 | zzlay(ig,l)=pphi(ig,l)/g |
---|
868 | ENDDO |
---|
869 | ENDDO |
---|
870 | DO ig=1,ngrid |
---|
871 | zzlev(ig,1)=0. |
---|
872 | zzlev(ig,nlayer+1)=1.e7 ! dummy top of last layer above 10000 km... |
---|
873 | ENDDO |
---|
874 | DO l=2,nlayer |
---|
875 | DO ig=1,ngrid |
---|
876 | z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l)) |
---|
877 | z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l)) |
---|
878 | zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) |
---|
879 | ENDDO |
---|
880 | ENDDO |
---|
881 | |
---|
882 | |
---|
883 | ! Potential temperature calculation not the same in physiq and dynamic |
---|
884 | |
---|
885 | c Compute potential temperature |
---|
886 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
887 | DO l=1,nlayer |
---|
888 | DO ig=1,ngrid |
---|
889 | zpopsk(ig,l)=(zplay(ig,l)/zplev(ig,1))**rcp |
---|
890 | zh(ig,l)=pt(ig,l)/zpopsk(ig,l) |
---|
891 | ENDDO |
---|
892 | ENDDO |
---|
893 | |
---|
894 | #ifndef MESOSCALE |
---|
895 | c----------------------------------------------------------------------- |
---|
896 | c 1.2.5 Compute mean mass, cp, and R |
---|
897 | c -------------------------------- |
---|
898 | |
---|
899 | if(photochem.or.callthermos) then |
---|
900 | call concentrations(ngrid,nlayer,nq, |
---|
901 | & zplay,pt,pdt,pq,pdq,ptimestep) |
---|
902 | endif |
---|
903 | #endif |
---|
904 | |
---|
905 | ! Compute vertical velocity (m/s) from vertical mass flux |
---|
906 | ! w = F / (rho*area) and rho = P/(r*T) |
---|
907 | ! but first linearly interpolate mass flux to mid-layers |
---|
908 | do l=1,nlayer-1 |
---|
909 | pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1)) |
---|
910 | enddo |
---|
911 | pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0 |
---|
912 | do l=1,nlayer |
---|
913 | pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) / |
---|
914 | & (pplay(1:ngrid,l)*cell_area(1:ngrid)) |
---|
915 | ! NB: here we use r and not rnew since this diagnostic comes |
---|
916 | ! from the dynamics |
---|
917 | enddo |
---|
918 | |
---|
919 | ! test for co2 conservation with co2 microphysics |
---|
920 | if (igcm_co2_ice.ne.0) then |
---|
921 | ! calculates the amount of co2 at the beginning of physics |
---|
922 | co2totA = 0. |
---|
923 | do ig=1,ngrid |
---|
924 | do l=1,nlayer |
---|
925 | co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g* |
---|
926 | & (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice) |
---|
927 | & +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep) |
---|
928 | end do |
---|
929 | do islope = 1,nslope |
---|
930 | co2totA = co2totA + qsurf(ig,igcm_co2,islope)* |
---|
931 | & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) |
---|
932 | enddo |
---|
933 | end do |
---|
934 | else |
---|
935 | co2totA = 0. |
---|
936 | do ig=1,ngrid |
---|
937 | do l=1,nlayer |
---|
938 | co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g* |
---|
939 | & (pq(ig,l,igcm_co2) |
---|
940 | & +pdq(ig,l,igcm_co2)*ptimestep) |
---|
941 | end do |
---|
942 | do islope = 1,nslope |
---|
943 | co2totA = co2totA + qsurf(ig,igcm_co2,islope)* |
---|
944 | & subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) |
---|
945 | enddo |
---|
946 | end do |
---|
947 | endif ! of if (igcm_co2_ice.ne.0) |
---|
948 | c----------------------------------------------------------------------- |
---|
949 | c 2. Compute radiative tendencies : |
---|
950 | c------------------------------------ |
---|
951 | |
---|
952 | IF (callrad) THEN |
---|
953 | |
---|
954 | c Local Solar zenith angle |
---|
955 | c ~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
956 | |
---|
957 | CALL orbite(zls,dist_sol,declin) |
---|
958 | |
---|
959 | IF (diurnal) THEN |
---|
960 | ztim1=SIN(declin) |
---|
961 | ztim2=COS(declin)*COS(2.*pi*(zday-.5)) |
---|
962 | ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) |
---|
963 | |
---|
964 | CALL solang(ngrid,sinlon,coslon,sinlat,coslat, |
---|
965 | & ztim1,ztim2,ztim3, mu0,fract) |
---|
966 | |
---|
967 | ELSE |
---|
968 | CALL mucorr(ngrid,declin,latitude,mu0,fract,10000.,rad) |
---|
969 | ENDIF ! of IF (diurnal) |
---|
970 | |
---|
971 | IF( MOD(icount-1,iradia).EQ.0) THEN |
---|
972 | |
---|
973 | c NLTE cooling from CO2 emission |
---|
974 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
975 | IF(callnlte) then |
---|
976 | if(nltemodel.eq.0.or.nltemodel.eq.1) then |
---|
977 | CALL nltecool(ngrid,nlayer,nq,zplay,pt,pq,zdtnlte) |
---|
978 | else if(nltemodel.eq.2) then |
---|
979 | co2vmr_gcm(1:ngrid,1:nlayer)= |
---|
980 | & pq(1:ngrid,1:nlayer,igcm_co2)* |
---|
981 | & mmean(1:ngrid,1:nlayer)/mmol(igcm_co2) |
---|
982 | n2vmr_gcm(1:ngrid,1:nlayer)= |
---|
983 | & pq(1:ngrid,1:nlayer,igcm_n2)* |
---|
984 | & mmean(1:ngrid,1:nlayer)/mmol(igcm_n2) |
---|
985 | covmr_gcm(1:ngrid,1:nlayer)= |
---|
986 | & pq(1:ngrid,1:nlayer,igcm_co)* |
---|
987 | & mmean(1:ngrid,1:nlayer)/mmol(igcm_co) |
---|
988 | ovmr_gcm(1:ngrid,1:nlayer)= |
---|
989 | & pq(1:ngrid,1:nlayer,igcm_o)* |
---|
990 | & mmean(1:ngrid,1:nlayer)/mmol(igcm_o) |
---|
991 | |
---|
992 | CALL nlte_tcool(ngrid,nlayer,zplay*9.869e-6, |
---|
993 | $ pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, |
---|
994 | $ ovmr_gcm, zdtnlte,ierr_nlte,varerr ) |
---|
995 | if(ierr_nlte.gt.0) then |
---|
996 | write(*,*) |
---|
997 | $ 'WARNING: nlte_tcool output with error message', |
---|
998 | $ 'ierr_nlte=',ierr_nlte,'varerr=',varerr |
---|
999 | write(*,*)'I will continue anyway' |
---|
1000 | endif |
---|
1001 | |
---|
1002 | zdtnlte(1:ngrid,1:nlayer)= |
---|
1003 | & zdtnlte(1:ngrid,1:nlayer)/86400. |
---|
1004 | endif |
---|
1005 | ELSE |
---|
1006 | zdtnlte(:,:)=0. |
---|
1007 | ENDIF !end callnlte |
---|
1008 | |
---|
1009 | c Find number of layers for LTE radiation calculations |
---|
1010 | IF(MOD(iphysiq*(icount-1),day_step).EQ.0) |
---|
1011 | & CALL nlthermeq(ngrid,nlayer,zplev,zplay) |
---|
1012 | |
---|
1013 | c rocketstorm : compute dust storm mesh fraction |
---|
1014 | IF (rdstorm) THEN |
---|
1015 | CALL calcstormfract(ngrid,nlayer,nq,pq, |
---|
1016 | & totstormfract) |
---|
1017 | ENDIF |
---|
1018 | |
---|
1019 | c Note: Dustopacity.F has been transferred to callradite.F |
---|
1020 | |
---|
1021 | #ifdef DUSTSTORM |
---|
1022 | !! specific case: save the quantity of dust before adding perturbation |
---|
1023 | |
---|
1024 | if (firstcall) then |
---|
1025 | pq_tmp(1:ngrid,1:nlayer,1)=pq(1:ngrid,1:nlayer,igcm_dust_mass) |
---|
1026 | pq_tmp(1:ngrid,1:nlayer,2)=pq(1:ngrid,1:nlayer,igcm_dust_number) |
---|
1027 | endif |
---|
1028 | #endif |
---|
1029 | |
---|
1030 | c Call main radiative transfer scheme |
---|
1031 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
1032 | c Transfer through CO2 (except NIR CO2 absorption) |
---|
1033 | c and aerosols (dust and water ice) |
---|
1034 | ! callradite for background dust (out of the rdstorm fraction) |
---|
1035 | clearatm=.true. |
---|
1036 | !! callradite for background dust (out of the topflows fraction) |
---|
1037 | nohmons=.true. |
---|
1038 | |
---|
1039 | c Radiative transfer |
---|
1040 | c ------------------ |
---|
1041 | ! callradite for the part with clouds |
---|
1042 | clearsky=.false. ! part with clouds for both cases CLFvarying true/false |
---|
1043 | CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq, |
---|
1044 | & albedo_meshavg,emis_meshavg, |
---|
1045 | & mu0,zplev,zplay,pt,tsurf_meshavg,fract,dist_sol,igout, |
---|
1046 | & zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat), |
---|
1047 | & fluxsurf_up_sw, |
---|
1048 | & fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw, |
---|
1049 | & tau_pref_scenario,tau_pref_gcm, |
---|
1050 | & tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, |
---|
1051 | & taucloudtes,rdust,rice,nuice,riceco2,nuiceco2, |
---|
1052 | & qsurf_meshavg(:,igcm_co2),rstormdust,rtopdust,totstormfract, |
---|
1053 | & clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac) |
---|
1054 | |
---|
1055 | DO islope=1,nslope |
---|
1056 | fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat) |
---|
1057 | fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat) |
---|
1058 | ENDDO |
---|
1059 | |
---|
1060 | ! case of sub-grid water ice clouds: callradite for the clear case |
---|
1061 | IF (CLFvarying) THEN |
---|
1062 | ! ---> PROBLEMS WITH ALLOCATED ARRAYS |
---|
1063 | ! (temporary solution in callcorrk: do not deallocate |
---|
1064 | ! if |
---|
1065 | ! CLFvarying ...) ?? AP ?? |
---|
1066 | clearsky=.true. |
---|
1067 | CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq, |
---|
1068 | & albedo_meshavg,emis_meshavg,mu0,zplev,zplay,pt, |
---|
1069 | & tsurf_meshavg,fract, |
---|
1070 | & dist_sol,igout,zdtlwclf,zdtswclf, |
---|
1071 | & fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf, |
---|
1072 | & fluxtop_lwclf,fluxtop_dn_swclf,fluxtop_up_swclf, |
---|
1073 | & tau_pref_scenario,tau_pref_gcm,tau,aerosol, |
---|
1074 | & dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, |
---|
1075 | & taucloudtesclf,rdust, |
---|
1076 | & rice,nuice,riceco2, nuiceco2, |
---|
1077 | & qsurf_meshavg(:,igcm_co2), |
---|
1078 | & rstormdust,rtopdust,totstormfract, |
---|
1079 | & clearatm,dsords,dsotop, |
---|
1080 | & nohmons,clearsky,totcloudfrac) |
---|
1081 | clearsky = .false. ! just in case. |
---|
1082 | ! Sum the fluxes and heating rates from cloudy/clear |
---|
1083 | ! cases |
---|
1084 | DO ig=1,ngrid |
---|
1085 | tf_clf=totcloudfrac(ig) |
---|
1086 | ntf_clf=1.-tf_clf |
---|
1087 | DO islope=1,nslope |
---|
1088 | fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig) |
---|
1089 | & + tf_clf*fluxsurf_lw(ig,islope) |
---|
1090 | fluxsurf_dn_sw(ig,1:2,islope) = |
---|
1091 | & ntf_clf*fluxsurf_dn_swclf(ig,1:2) |
---|
1092 | & + tf_clf*fluxsurf_dn_sw(ig,1:2,islope) |
---|
1093 | ENDDO |
---|
1094 | fluxsurf_up_sw(ig,1:2) = |
---|
1095 | & ntf_clf*fluxsurf_up_swclf(ig,1:2) |
---|
1096 | & + tf_clf*fluxsurf_up_sw(ig,1:2) |
---|
1097 | fluxtop_lw(ig) = ntf_clf*fluxtop_lwclf(ig) |
---|
1098 | & + tf_clf*fluxtop_lw(ig) |
---|
1099 | fluxtop_dn_sw(ig,1:2)=ntf_clf*fluxtop_dn_swclf(ig,1:2) |
---|
1100 | & + tf_clf*fluxtop_dn_sw(ig,1:2) |
---|
1101 | fluxtop_up_sw(ig,1:2)=ntf_clf*fluxtop_up_swclf(ig,1:2) |
---|
1102 | & + tf_clf*fluxtop_up_sw(ig,1:2) |
---|
1103 | taucloudtes(ig) = ntf_clf*taucloudtesclf(ig) |
---|
1104 | & + tf_clf*taucloudtes(ig) |
---|
1105 | zdtlw(ig,1:nlayer) = ntf_clf*zdtlwclf(ig,1:nlayer) |
---|
1106 | & + tf_clf*zdtlw(ig,1:nlayer) |
---|
1107 | zdtsw(ig,1:nlayer) = ntf_clf*zdtswclf(ig,1:nlayer) |
---|
1108 | & + tf_clf*zdtsw(ig,1:nlayer) |
---|
1109 | ENDDO |
---|
1110 | |
---|
1111 | ENDIF ! (CLFvarying) |
---|
1112 | |
---|
1113 | !============================================================================ |
---|
1114 | |
---|
1115 | #ifdef DUSTSTORM |
---|
1116 | !! specific case: compute the added quantity of dust for perturbation |
---|
1117 | if (firstcall) then |
---|
1118 | pdq(1:ngrid,1:nlayer,igcm_dust_mass)= |
---|
1119 | & pdq(1:ngrid,1:nlayer,igcm_dust_mass) |
---|
1120 | & - pq_tmp(1:ngrid,1:nlayer,1) |
---|
1121 | & + pq(1:ngrid,1:nlayer,igcm_dust_mass) |
---|
1122 | pdq(1:ngrid,1:nlayer,igcm_dust_number)= |
---|
1123 | & pdq(1:ngrid,1:nlayer,igcm_dust_number) |
---|
1124 | & - pq_tmp(1:ngrid,1:nlayer,2) |
---|
1125 | & + pq(1:ngrid,1:nlayer,igcm_dust_number) |
---|
1126 | endif |
---|
1127 | #endif |
---|
1128 | |
---|
1129 | c Outputs for basic check (middle of domain) |
---|
1130 | c ------------------------------------------ |
---|
1131 | write(*,'("Ls =",f11.6," check lat =",f10.6, |
---|
1132 | & " lon =",f11.6)') |
---|
1133 | & zls*180./pi,latitude(igout)*180/pi, |
---|
1134 | & longitude(igout)*180/pi |
---|
1135 | |
---|
1136 | write(*,'(" tau_pref_gcm(",f4.0," Pa) =",f9.6, |
---|
1137 | & " tau(",f4.0," Pa) =",f9.6)') |
---|
1138 | & odpref,tau_pref_gcm(igout), |
---|
1139 | & odpref,tau(igout,1)*odpref/zplev(igout,1) |
---|
1140 | |
---|
1141 | |
---|
1142 | c --------------------------------------------------------- |
---|
1143 | c Call slope parameterization for direct and scattered flux |
---|
1144 | c --------------------------------------------------------- |
---|
1145 | IF(callslope) THEN |
---|
1146 | ! assume that in this case, nslope = 1 |
---|
1147 | if(nslope.ne.1) then |
---|
1148 | call abort_physic( |
---|
1149 | & "physiq","callslope=true but nslope.ne.1",1) |
---|
1150 | endif |
---|
1151 | print *, 'Slope scheme is on and computing...' |
---|
1152 | DO ig=1,ngrid |
---|
1153 | sl_the = theta_sl(ig) |
---|
1154 | IF (sl_the .ne. 0.) THEN |
---|
1155 | ztim1=fluxsurf_dn_sw(ig,1,iflat) |
---|
1156 | & +fluxsurf_dn_sw(ig,2,iflat) |
---|
1157 | DO l=1,2 |
---|
1158 | sl_lct = ptime*24. + 180.*longitude(ig)/pi/15. |
---|
1159 | sl_ra = pi*(1.0-sl_lct/12.) |
---|
1160 | sl_lat = 180.*latitude(ig)/pi |
---|
1161 | sl_tau = tau(ig,1) !il faudrait iaerdust(iaer) |
---|
1162 | sl_alb = albedo(ig,l,iflat) |
---|
1163 | sl_psi = psi_sl(ig) |
---|
1164 | sl_fl0 = fluxsurf_dn_sw(ig,l,iflat) |
---|
1165 | sl_di0 = 0. |
---|
1166 | if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then |
---|
1167 | sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) |
---|
1168 | sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol |
---|
1169 | sl_di0 = sl_di0/ztim1 |
---|
1170 | sl_di0 = fluxsurf_dn_sw(ig,l,iflat)*sl_di0 |
---|
1171 | endif |
---|
1172 | ! you never know (roundup concern...) |
---|
1173 | if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 |
---|
1174 | !!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
1175 | CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, |
---|
1176 | & sl_tau, sl_alb, sl_the, sl_psi, |
---|
1177 | & sl_di0, sl_fl0, sl_flu ) |
---|
1178 | !!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
1179 | fluxsurf_dn_sw(ig,l,1) = sl_flu |
---|
1180 | ENDDO |
---|
1181 | !!! compute correction on IR flux as well |
---|
1182 | sky= (1.+cos(pi*theta_sl(ig)/180.))/2. |
---|
1183 | fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky |
---|
1184 | ENDIF |
---|
1185 | ENDDO |
---|
1186 | ELSE ! not calling subslope, nslope might be > 1 |
---|
1187 | DO islope = 1,nslope |
---|
1188 | sl_the=abs(def_slope_mean(islope)) |
---|
1189 | IF (sl_the .gt. 1e-6) THEN |
---|
1190 | IF(def_slope_mean(islope).ge.0.) THEN |
---|
1191 | psi_sl(:) = 0. !Northward slope |
---|
1192 | ELSE |
---|
1193 | psi_sl(:) = 180. !Southward slope |
---|
1194 | ENDIF |
---|
1195 | DO ig=1,ngrid |
---|
1196 | ztim1=fluxsurf_dn_sw(ig,1,islope) |
---|
1197 | s +fluxsurf_dn_sw(ig,2,islope) |
---|
1198 | DO l=1,2 |
---|
1199 | sl_lct = ptime*24. + 180.*longitude(ig)/pi/15. |
---|
1200 | sl_ra = pi*(1.0-sl_lct/12.) |
---|
1201 | sl_lat = 180.*latitude(ig)/pi |
---|
1202 | sl_tau = tau(ig,1) !il faudrait iaerdust(iaer) |
---|
1203 | sl_alb = albedo(ig,l,islope) |
---|
1204 | sl_psi = psi_sl(ig) |
---|
1205 | sl_fl0 = fluxsurf_dn_sw(ig,l,islope) |
---|
1206 | sl_di0 = 0. |
---|
1207 | if (mu0(ig) .gt. 0.) then |
---|
1208 | sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) |
---|
1209 | sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol |
---|
1210 | sl_di0 = sl_di0/ztim1 |
---|
1211 | sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0 |
---|
1212 | endif |
---|
1213 | ! you never know (roundup concern...) |
---|
1214 | if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 |
---|
1215 | !!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
1216 | CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, |
---|
1217 | & sl_tau, sl_alb, sl_the, sl_psi, |
---|
1218 | & sl_di0, sl_fl0, sl_flu ) |
---|
1219 | !!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
1220 | fluxsurf_dn_sw(ig,l,islope) = sl_flu |
---|
1221 | ENDDO |
---|
1222 | !!! compute correction on IR flux as well |
---|
1223 | |
---|
1224 | fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope) |
---|
1225 | & *sky_slope(islope) |
---|
1226 | ENDDO |
---|
1227 | ENDIF ! sub grid is not flat |
---|
1228 | ENDDO ! islope = 1,nslope |
---|
1229 | ENDIF ! callslope |
---|
1230 | |
---|
1231 | c CO2 near infrared absorption |
---|
1232 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
1233 | zdtnirco2(:,:)=0 |
---|
1234 | if (callnirco2) then |
---|
1235 | call nirco2abs (ngrid,nlayer,zplay,dist_sol,nq,pq, |
---|
1236 | . mu0,fract,declin, zdtnirco2) |
---|
1237 | endif |
---|
1238 | |
---|
1239 | c Radiative flux from the sky absorbed by the surface (W.m-2) |
---|
1240 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
1241 | DO ig=1,ngrid |
---|
1242 | DO islope = 1,nslope |
---|
1243 | fluxrad_sky(ig,islope) = |
---|
1244 | $ emis(ig,islope)*fluxsurf_lw(ig,islope) |
---|
1245 | $ +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope)) |
---|
1246 | $ +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope)) |
---|
1247 | ENDDO |
---|
1248 | ENDDO |
---|
1249 | |
---|
1250 | c Net atmospheric radiative heating rate (K.s-1) |
---|
1251 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
1252 | IF(callnlte) THEN |
---|
1253 | CALL blendrad(ngrid, nlayer, zplay, |
---|
1254 | & zdtsw, zdtlw, zdtnirco2, zdtnlte, dtrad) |
---|
1255 | ELSE |
---|
1256 | DO l=1,nlayer |
---|
1257 | DO ig=1,ngrid |
---|
1258 | dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) |
---|
1259 | & +zdtnirco2(ig,l) |
---|
1260 | ENDDO |
---|
1261 | ENDDO |
---|
1262 | ENDIF |
---|
1263 | |
---|
1264 | ENDIF ! of if(mod(icount-1,iradia).eq.0) |
---|
1265 | |
---|
1266 | c Transformation of the radiative tendencies: |
---|
1267 | c ------------------------------------------- |
---|
1268 | |
---|
1269 | c Net radiative surface flux (W.m-2) |
---|
1270 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
1271 | |
---|
1272 | c |
---|
1273 | DO ig=1,ngrid |
---|
1274 | DO islope = 1,nslope |
---|
1275 | zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope) |
---|
1276 | zplanck(ig)=emis(ig,islope)* |
---|
1277 | $ stephan*zplanck(ig)*zplanck(ig) |
---|
1278 | fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig) |
---|
1279 | IF(callslope) THEN |
---|
1280 | sky= (1.+cos(pi*theta_sl(ig)/180.))/2. |
---|
1281 | fluxrad(ig,nslope)=fluxrad(ig,nslope)+ |
---|
1282 | $ (1.-sky)*zplanck(ig) |
---|
1283 | ELSE |
---|
1284 | fluxrad(ig,islope)=fluxrad(ig,islope) + |
---|
1285 | $ (1.-sky_slope(iflat))*emis(ig,iflat)* |
---|
1286 | $ stephan*tsurf(ig,iflat)**4 |
---|
1287 | ENDIF |
---|
1288 | ENDDO |
---|
1289 | ENDDO |
---|
1290 | |
---|
1291 | DO l=1,nlayer |
---|
1292 | DO ig=1,ngrid |
---|
1293 | pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) |
---|
1294 | ENDDO |
---|
1295 | ENDDO |
---|
1296 | |
---|
1297 | ENDIF ! of IF (callrad) |
---|
1298 | |
---|
1299 | c 3.1 Rocket dust storm |
---|
1300 | c ------------------------------------------- |
---|
1301 | IF (rdstorm) THEN |
---|
1302 | clearatm=.false. |
---|
1303 | pdqrds(:,:,:)=0. |
---|
1304 | qdusttotal0(:)=0. |
---|
1305 | qdusttotal1(:)=0. |
---|
1306 | do ig=1,ngrid |
---|
1307 | do l=1,nlayer |
---|
1308 | zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ! updated potential |
---|
1309 | ! temperature tendency |
---|
1310 | ! for diagnostics |
---|
1311 | ! qdustrds0(ig,l)=pq(ig,l,igcm_dust_mass)+ |
---|
1312 | ! & pdq(ig,l,igcm_dust_mass)*ptimestep |
---|
1313 | ! qstormrds0(ig,l)=pq(ig,l,igcm_stormdust_mass)+ |
---|
1314 | ! & pdq(ig,l,igcm_stormdust_mass)*ptimestep |
---|
1315 | ! qdusttotal0(ig)=qdusttotal0(ig)+(qdustrds0(ig,l)+ |
---|
1316 | ! & qstormrds0(ig,l))*(zplev(ig,l)- |
---|
1317 | ! & zplev(ig,l+1))/g |
---|
1318 | enddo |
---|
1319 | enddo |
---|
1320 | ! call write_output('qdustrds0','qdust before rds', |
---|
1321 | ! & 'kg/kg ',qdustrds0(:,:)) |
---|
1322 | ! call write_output('qstormrds0','qstorm before rds', |
---|
1323 | ! & 'kg/kg ',qstormrds0(:,:)) |
---|
1324 | |
---|
1325 | CALL rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep, |
---|
1326 | & pq,pdq,pt,pdt,zplev,zplay,zzlev, |
---|
1327 | & zzlay,zdtsw,zdtlw, |
---|
1328 | c for radiative transfer |
---|
1329 | & clearatm,icount,zday,zls, |
---|
1330 | & tsurf_meshavg,qsurf_meshavg(:,igcm_co2), |
---|
1331 | & igout,totstormfract,tauscaling, |
---|
1332 | & dust_rad_adjust,IRtoVIScoef, |
---|
1333 | & albedo_meshavg,emis_meshavg, |
---|
1334 | c input sub-grid scale cloud |
---|
1335 | & clearsky,totcloudfrac, |
---|
1336 | c input sub-grid scale topography |
---|
1337 | & nohmons, |
---|
1338 | c output |
---|
1339 | & pdqrds,wspeed,dsodust,dsords,dsotop, |
---|
1340 | & tau_pref_scenario,tau_pref_gcm) |
---|
1341 | |
---|
1342 | c update the tendencies of both dust after vertical transport |
---|
1343 | DO l=1,nlayer |
---|
1344 | DO ig=1,ngrid |
---|
1345 | pdq(ig,l,igcm_stormdust_mass)= |
---|
1346 | & pdq(ig,l,igcm_stormdust_mass)+ |
---|
1347 | & pdqrds(ig,l,igcm_stormdust_mass) |
---|
1348 | pdq(ig,l,igcm_stormdust_number)= |
---|
1349 | & pdq(ig,l,igcm_stormdust_number)+ |
---|
1350 | & pdqrds(ig,l,igcm_stormdust_number) |
---|
1351 | |
---|
1352 | pdq(ig,l,igcm_dust_mass)= |
---|
1353 | & pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass) |
---|
1354 | pdq(ig,l,igcm_dust_number)= |
---|
1355 | & pdq(ig,l,igcm_dust_number)+ |
---|
1356 | & pdqrds(ig,l,igcm_dust_number) |
---|
1357 | |
---|
1358 | ENDDO |
---|
1359 | ENDDO |
---|
1360 | do l=1,nlayer |
---|
1361 | do ig=1,ngrid |
---|
1362 | qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+ |
---|
1363 | & pdq(ig,l,igcm_dust_mass)*ptimestep |
---|
1364 | qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+ |
---|
1365 | & pdq(ig,l,igcm_stormdust_mass)*ptimestep |
---|
1366 | qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+ |
---|
1367 | & qstormrds1(ig,l))*(zplev(ig,l)- |
---|
1368 | & zplev(ig,l+1))/g |
---|
1369 | enddo |
---|
1370 | enddo |
---|
1371 | |
---|
1372 | c for diagnostics |
---|
1373 | ! call write_output('qdustrds1','qdust after rds', |
---|
1374 | ! & 'kg/kg ',qdustrds1(:,:)) |
---|
1375 | ! call write_output('qstormrds1','qstorm after rds', |
---|
1376 | ! & 'kg/kg ',qstormrds1(:,:)) |
---|
1377 | ! |
---|
1378 | ! call write_output('qdusttotal0','q sum before rds', |
---|
1379 | ! & 'kg/m2 ',qdusttotal0(:)) |
---|
1380 | ! call write_output('qdusttotal1','q sum after rds', |
---|
1381 | ! & 'kg/m2 ',qdusttotal1(:)) |
---|
1382 | |
---|
1383 | ENDIF ! end of if(rdstorm) |
---|
1384 | |
---|
1385 | c 3.2 Dust entrained from the PBL up to the top of sub-grid scale topography |
---|
1386 | c ------------------------------------------- |
---|
1387 | IF (topflows) THEN |
---|
1388 | clearatm=.true. ! stormdust is not accounted in the extra heating on top of the mountains |
---|
1389 | nohmons=.false. |
---|
1390 | pdqtop(:,:,:)=0. |
---|
1391 | CALL topmons(ngrid,nlayer,nq,ptime,ptimestep, |
---|
1392 | & pq,pdq,pt,pdt,zplev,zplay,zzlev, |
---|
1393 | & zzlay,zdtsw,zdtlw, |
---|
1394 | & icount,zday,zls,tsurf(:,iflat), |
---|
1395 | & qsurf_meshavg(:,igcm_co2), |
---|
1396 | & igout,aerosol,tauscaling,dust_rad_adjust, |
---|
1397 | & IRtoVIScoef,albedo_meshavg,emis_meshavg, |
---|
1398 | & totstormfract,clearatm, |
---|
1399 | & clearsky,totcloudfrac, |
---|
1400 | & nohmons, |
---|
1401 | & pdqtop,wtop,dsodust,dsords,dsotop, |
---|
1402 | & tau_pref_scenario,tau_pref_gcm) |
---|
1403 | |
---|
1404 | c update the tendencies of both dust after vertical transport |
---|
1405 | DO l=1,nlayer |
---|
1406 | DO ig=1,ngrid |
---|
1407 | pdq(ig,l,igcm_topdust_mass)= |
---|
1408 | & pdq(ig,l,igcm_topdust_mass)+ |
---|
1409 | & pdqtop(ig,l,igcm_topdust_mass) |
---|
1410 | pdq(ig,l,igcm_topdust_number)= |
---|
1411 | & pdq(ig,l,igcm_topdust_number)+ |
---|
1412 | & pdqtop(ig,l,igcm_topdust_number) |
---|
1413 | pdq(ig,l,igcm_dust_mass)= |
---|
1414 | & pdq(ig,l,igcm_dust_mass)+ pdqtop(ig,l,igcm_dust_mass) |
---|
1415 | pdq(ig,l,igcm_dust_number)= |
---|
1416 | & pdq(ig,l,igcm_dust_number)+pdqtop(ig,l,igcm_dust_number) |
---|
1417 | |
---|
1418 | ENDDO |
---|
1419 | ENDDO |
---|
1420 | |
---|
1421 | ENDIF ! end of if (topflows) |
---|
1422 | |
---|
1423 | c 3.3 Dust injection from the surface |
---|
1424 | c ------------------------------------------- |
---|
1425 | if (dustinjection.gt.0) then |
---|
1426 | |
---|
1427 | CALL compute_dtau(ngrid,nlayer, |
---|
1428 | & zday,pplev,tau_pref_gcm, |
---|
1429 | & ptimestep,local_time,IRtoVIScoef, |
---|
1430 | & dustliftday) |
---|
1431 | endif ! end of if (dustinjection.gt.0) |
---|
1432 | |
---|
1433 | c----------------------------------------------------------------------- |
---|
1434 | c 4. Gravity wave and subgrid scale topography drag : |
---|
1435 | c ------------------------------------------------- |
---|
1436 | |
---|
1437 | |
---|
1438 | IF(calllott)THEN |
---|
1439 | CALL calldrag_noro(ngrid,nlayer,ptimestep, |
---|
1440 | & zplay,zplev,pt,pu,pv,zdtgw,zdugw,zdvgw) |
---|
1441 | |
---|
1442 | DO l=1,nlayer |
---|
1443 | DO ig=1,ngrid |
---|
1444 | pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l) |
---|
1445 | pdu(ig,l)=pdu(ig,l)+zdugw(ig,l) |
---|
1446 | pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l) |
---|
1447 | ENDDO |
---|
1448 | ENDDO |
---|
1449 | ENDIF |
---|
1450 | |
---|
1451 | c----------------------------------------------------------------------- |
---|
1452 | c 5. Vertical diffusion (turbulent mixing): |
---|
1453 | c ----------------------------------------- |
---|
1454 | |
---|
1455 | IF (calldifv) THEN |
---|
1456 | DO ig=1,ngrid |
---|
1457 | DO islope = 1,nslope |
---|
1458 | zflubid(ig,islope)=fluxrad(ig,islope) |
---|
1459 | & +fluxgrd(ig,islope) |
---|
1460 | ENDDO |
---|
1461 | ENDDO |
---|
1462 | zdum1(:,:)=0 |
---|
1463 | zdum2(:,:)=0 |
---|
1464 | do l=1,nlayer |
---|
1465 | do ig=1,ngrid |
---|
1466 | zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) |
---|
1467 | enddo |
---|
1468 | enddo |
---|
1469 | |
---|
1470 | c ---------------------- |
---|
1471 | c Treatment of a special case : using new surface layer (Richardson based) |
---|
1472 | c without using the thermals in gcm and mesoscale can yield problems in |
---|
1473 | c weakly unstable situations when winds are near to 0. For those cases, we add |
---|
1474 | c a unit subgrid gustiness. Remember that thermals should be used we using the |
---|
1475 | c Richardson based surface layer model. |
---|
1476 | IF ( .not.calltherm |
---|
1477 | . .and. callrichsl |
---|
1478 | . .and. .not.turb_resolved) THEN |
---|
1479 | |
---|
1480 | DO ig=1, ngrid |
---|
1481 | IF (zh(ig,1) .lt. tsurf_meshavg(ig)) THEN |
---|
1482 | wstar(ig)=1. |
---|
1483 | hfmax_th(ig)=0.2 |
---|
1484 | ELSE |
---|
1485 | wstar(ig)=0. |
---|
1486 | hfmax_th(ig)=0. |
---|
1487 | ENDIF |
---|
1488 | ENDDO |
---|
1489 | ENDIF |
---|
1490 | |
---|
1491 | c ---------------------- |
---|
1492 | |
---|
1493 | IF (tke_heat_flux .ne. 0.) THEN |
---|
1494 | |
---|
1495 | zz1(:)=(pt(:,1)+pdt(:,1)*ptimestep)*(r/g)* |
---|
1496 | & (-alog(zplay(:,1)/zplev(:,1))) |
---|
1497 | pdt(:,1)=pdt(:,1) + (tke_heat_flux/zz1(:))*zpopsk(:,1) |
---|
1498 | ENDIF |
---|
1499 | |
---|
1500 | c Calling vdif (Martian version WITH CO2 condensation) |
---|
1501 | dwatercap_dif(:,:) = 0. |
---|
1502 | zcdh(:) = 0. |
---|
1503 | zcdv(:) = 0. |
---|
1504 | |
---|
1505 | CALL vdifc(ngrid,nlayer,nq,zpopsk, |
---|
1506 | $ ptimestep,capcal,lwrite, |
---|
1507 | $ zplay,zplev,zzlay,zzlev,z0, |
---|
1508 | $ pu,pv,zh,pq,tsurf,emis,qsurf, |
---|
1509 | $ zdum1,zdum2,zdh,pdq,zflubid, |
---|
1510 | $ zdudif,zdvdif,zdhdif,zdtsdif,q2, |
---|
1511 | & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th, |
---|
1512 | & zcondicea_co2microp,sensibFlux, |
---|
1513 | & dustliftday,local_time,watercap,dwatercap_dif) |
---|
1514 | |
---|
1515 | DO ig=1,ngrid |
---|
1516 | zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:) |
---|
1517 | dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:) |
---|
1518 | ENDDO |
---|
1519 | |
---|
1520 | call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif, |
---|
1521 | & albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp) |
---|
1522 | IF (.not.turb_resolved) THEN |
---|
1523 | DO l=1,nlayer |
---|
1524 | DO ig=1,ngrid |
---|
1525 | pdv(ig,l)=pdv(ig,l)+zdvdif(ig,l) |
---|
1526 | pdu(ig,l)=pdu(ig,l)+zdudif(ig,l) |
---|
1527 | pdt(ig,l)=pdt(ig,l)+zdhdif(ig,l)*zpopsk(ig,l) |
---|
1528 | |
---|
1529 | zdtdif(ig,l)=zdhdif(ig,l)*zpopsk(ig,l) ! for diagnostic only |
---|
1530 | ENDDO |
---|
1531 | ENDDO |
---|
1532 | |
---|
1533 | DO iq=1, nq |
---|
1534 | DO l=1,nlayer |
---|
1535 | DO ig=1,ngrid |
---|
1536 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) |
---|
1537 | ENDDO |
---|
1538 | ENDDO |
---|
1539 | ENDDO |
---|
1540 | DO iq=1, nq |
---|
1541 | DO ig=1,ngrid |
---|
1542 | dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:) |
---|
1543 | ENDDO |
---|
1544 | ENDDO |
---|
1545 | |
---|
1546 | ELSE |
---|
1547 | write (*,*) '******************************************' |
---|
1548 | write (*,*) '** LES mode: the difv part is only used to' |
---|
1549 | write (*,*) '** - provide HFX and UST to the dynamics' |
---|
1550 | write (*,*) '** - update TSURF' |
---|
1551 | write (*,*) '******************************************' |
---|
1552 | !! Specific treatment for lifting in turbulent-resolving mode (AC) |
---|
1553 | IF (lifting .and. doubleq) THEN |
---|
1554 | !! lifted dust is injected in the first layer. |
---|
1555 | !! Sedimentation must be called after turbulent mixing, i.e. on next step, after WRF. |
---|
1556 | !! => lifted dust is not incremented before the sedimentation step. |
---|
1557 | zdqdif(1:ngrid,1,1:nq)=0. |
---|
1558 | zdqdif(1:ngrid,1,igcm_dust_number) = |
---|
1559 | . -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_number) |
---|
1560 | zdqdif(1:ngrid,1,igcm_dust_mass) = |
---|
1561 | . -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_mass) |
---|
1562 | zdqdif(1:ngrid,2:nlayer,1:nq) = 0. |
---|
1563 | DO iq=1, nq |
---|
1564 | IF ((iq .ne. igcm_dust_mass) |
---|
1565 | & .and. (iq .ne. igcm_dust_number)) THEN |
---|
1566 | zdqsdif(:,iq,:)=0. |
---|
1567 | ENDIF |
---|
1568 | ENDDO |
---|
1569 | ELSE |
---|
1570 | zdqdif(1:ngrid,1:nlayer,1:nq) = 0. |
---|
1571 | zdqsdif(1:ngrid,1:nq,1:nslope) = 0. |
---|
1572 | ENDIF |
---|
1573 | ENDIF |
---|
1574 | ELSE |
---|
1575 | DO ig=1,ngrid |
---|
1576 | DO islope=1,nslope |
---|
1577 | zdtsurf(ig,islope)=zdtsurf(ig,islope)+ |
---|
1578 | s (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope) |
---|
1579 | ENDDO |
---|
1580 | ENDDO |
---|
1581 | |
---|
1582 | IF (turb_resolved) THEN |
---|
1583 | write(*,*) 'Turbulent-resolving mode !' |
---|
1584 | write(*,*) 'Please set calldifv to T in callphys.def' |
---|
1585 | call abort_physic("physiq","turbulent-resolving mode",1) |
---|
1586 | ENDIF |
---|
1587 | ENDIF ! of IF (calldifv) |
---|
1588 | |
---|
1589 | c----------------------------------------------------------------------- |
---|
1590 | c 6. Thermals : |
---|
1591 | c ----------------------------- |
---|
1592 | |
---|
1593 | if(calltherm .and. .not.turb_resolved) then |
---|
1594 | |
---|
1595 | call calltherm_interface(ngrid,nlayer,nq,igcm_co2, |
---|
1596 | $ zzlev,zzlay, |
---|
1597 | $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, |
---|
1598 | $ zplay,zplev,pphi,zpopsk, |
---|
1599 | $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th, |
---|
1600 | $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux) |
---|
1601 | |
---|
1602 | DO l=1,nlayer |
---|
1603 | DO ig=1,ngrid |
---|
1604 | pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l) |
---|
1605 | pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l) |
---|
1606 | pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l) |
---|
1607 | q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep |
---|
1608 | ENDDO |
---|
1609 | ENDDO |
---|
1610 | |
---|
1611 | DO ig=1,ngrid |
---|
1612 | q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep |
---|
1613 | ENDDO |
---|
1614 | |
---|
1615 | DO iq=1,nq |
---|
1616 | DO l=1,nlayer |
---|
1617 | DO ig=1,ngrid |
---|
1618 | pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq) |
---|
1619 | ENDDO |
---|
1620 | ENDDO |
---|
1621 | ENDDO |
---|
1622 | |
---|
1623 | lmax_th_out(:)=real(lmax_th(:)) |
---|
1624 | |
---|
1625 | else !of if calltherm |
---|
1626 | lmax_th(:)=0 |
---|
1627 | wstar(:)=0. |
---|
1628 | hfmax_th(:)=0. |
---|
1629 | lmax_th_out(:)=0. |
---|
1630 | end if |
---|
1631 | |
---|
1632 | c----------------------------------------------------------------------- |
---|
1633 | c 7. Dry convective adjustment: |
---|
1634 | c ----------------------------- |
---|
1635 | |
---|
1636 | IF(calladj) THEN |
---|
1637 | |
---|
1638 | DO l=1,nlayer |
---|
1639 | DO ig=1,ngrid |
---|
1640 | zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) |
---|
1641 | ENDDO |
---|
1642 | ENDDO |
---|
1643 | zduadj(:,:)=0 |
---|
1644 | zdvadj(:,:)=0 |
---|
1645 | zdhadj(:,:)=0 |
---|
1646 | |
---|
1647 | CALL convadj(ngrid,nlayer,nq,ptimestep, |
---|
1648 | $ zplay,zplev,zpopsk,lmax_th, |
---|
1649 | $ pu,pv,zh,pq, |
---|
1650 | $ pdu,pdv,zdh,pdq, |
---|
1651 | $ zduadj,zdvadj,zdhadj, |
---|
1652 | $ zdqadj) |
---|
1653 | |
---|
1654 | DO l=1,nlayer |
---|
1655 | DO ig=1,ngrid |
---|
1656 | pdu(ig,l)=pdu(ig,l)+zduadj(ig,l) |
---|
1657 | pdv(ig,l)=pdv(ig,l)+zdvadj(ig,l) |
---|
1658 | pdt(ig,l)=pdt(ig,l)+zdhadj(ig,l)*zpopsk(ig,l) |
---|
1659 | |
---|
1660 | zdtadj(ig,l)=zdhadj(ig,l)*zpopsk(ig,l) ! for diagnostic only |
---|
1661 | ENDDO |
---|
1662 | ENDDO |
---|
1663 | |
---|
1664 | DO iq=1, nq |
---|
1665 | DO l=1,nlayer |
---|
1666 | DO ig=1,ngrid |
---|
1667 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqadj(ig,l,iq) |
---|
1668 | ENDDO |
---|
1669 | ENDDO |
---|
1670 | ENDDO |
---|
1671 | ENDIF ! of IF(calladj) |
---|
1672 | |
---|
1673 | c----------------------------------------------------- |
---|
1674 | c 8. Non orographic Gravity waves : |
---|
1675 | c ------------------------------------------------- |
---|
1676 | |
---|
1677 | IF (calllott_nonoro) THEN |
---|
1678 | |
---|
1679 | CALL nonoro_gwd_ran(ngrid,nlayer,ptimestep, |
---|
1680 | & cpnew,rnew, |
---|
1681 | & zplay, |
---|
1682 | & zmax_th, ! max altitude reached by thermals (m) |
---|
1683 | & pt, pu, pv, |
---|
1684 | & pdt, pdu, pdv, |
---|
1685 | & zustrhi,zvstrhi, |
---|
1686 | & d_t_hin, d_u_hin, d_v_hin) |
---|
1687 | |
---|
1688 | ! Update tendencies |
---|
1689 | pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer) |
---|
1690 | & +d_t_hin(1:ngrid,1:nlayer) |
---|
1691 | pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer) |
---|
1692 | & +d_u_hin(1:ngrid,1:nlayer) |
---|
1693 | pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer) |
---|
1694 | & +d_v_hin(1:ngrid,1:nlayer) |
---|
1695 | |
---|
1696 | ENDIF ! of IF (calllott_nonoro) |
---|
1697 | |
---|
1698 | c----------------------------------------------------------------------- |
---|
1699 | c 9. Specific parameterizations for tracers |
---|
1700 | c: ----------------------------------------- |
---|
1701 | |
---|
1702 | |
---|
1703 | c 9a. Water and ice |
---|
1704 | c --------------- |
---|
1705 | |
---|
1706 | c --------------------------------------- |
---|
1707 | c Water ice condensation in the atmosphere |
---|
1708 | c ---------------------------------------- |
---|
1709 | IF (water) THEN |
---|
1710 | |
---|
1711 | call watercloud(ngrid,nlayer,ptimestep, |
---|
1712 | & zplev,zplay,pdpsrf,zzlay, pt,pdt, |
---|
1713 | & pq,pdq,zdqcloud,zdtcloud, |
---|
1714 | & nq,tau,tauscaling,rdust,rice,nuice, |
---|
1715 | & rsedcloud,rhocloud,totcloudfrac) |
---|
1716 | c Temperature variation due to latent heat release |
---|
1717 | if (activice) then |
---|
1718 | pdt(1:ngrid,1:nlayer) = |
---|
1719 | & pdt(1:ngrid,1:nlayer) + |
---|
1720 | & zdtcloud(1:ngrid,1:nlayer) |
---|
1721 | endif |
---|
1722 | |
---|
1723 | ! increment water vapour and ice atmospheric tracers tendencies |
---|
1724 | pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = |
---|
1725 | & pdq(1:ngrid,1:nlayer,igcm_h2o_vap) + |
---|
1726 | & zdqcloud(1:ngrid,1:nlayer,igcm_h2o_vap) |
---|
1727 | pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = |
---|
1728 | & pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + |
---|
1729 | & zdqcloud(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
1730 | |
---|
1731 | if (hdo) then |
---|
1732 | ! increment HDO vapour and ice atmospheric tracers tendencies |
---|
1733 | pdq(1:ngrid,1:nlayer,igcm_hdo_vap) = |
---|
1734 | & pdq(1:ngrid,1:nlayer,igcm_hdo_vap) + |
---|
1735 | & zdqcloud(1:ngrid,1:nlayer,igcm_hdo_vap) |
---|
1736 | pdq(1:ngrid,1:nlayer,igcm_hdo_ice) = |
---|
1737 | & pdq(1:ngrid,1:nlayer,igcm_hdo_ice) + |
---|
1738 | & zdqcloud(1:ngrid,1:nlayer,igcm_hdo_ice) |
---|
1739 | endif !hdo |
---|
1740 | |
---|
1741 | ! increment dust and ccn masses and numbers |
---|
1742 | ! We need to check that we have Nccn & Ndust > 0 |
---|
1743 | ! This is due to single precision rounding problems |
---|
1744 | if (microphys) then |
---|
1745 | pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = |
---|
1746 | & pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + |
---|
1747 | & zdqcloud(1:ngrid,1:nlayer,igcm_ccn_mass) |
---|
1748 | pdq(1:ngrid,1:nlayer,igcm_ccn_number) = |
---|
1749 | & pdq(1:ngrid,1:nlayer,igcm_ccn_number) + |
---|
1750 | & zdqcloud(1:ngrid,1:nlayer,igcm_ccn_number) |
---|
1751 | where (pq(:,:,igcm_ccn_mass) + |
---|
1752 | & ptimestep*pdq(:,:,igcm_ccn_mass) < 0.) |
---|
1753 | pdq(:,:,igcm_ccn_mass) = |
---|
1754 | & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 |
---|
1755 | pdq(:,:,igcm_ccn_number) = |
---|
1756 | & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 |
---|
1757 | end where |
---|
1758 | where (pq(:,:,igcm_ccn_number) + |
---|
1759 | & ptimestep*pdq(:,:,igcm_ccn_number) < 0.) |
---|
1760 | pdq(:,:,igcm_ccn_mass) = |
---|
1761 | & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 |
---|
1762 | pdq(:,:,igcm_ccn_number) = |
---|
1763 | & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 |
---|
1764 | end where |
---|
1765 | endif |
---|
1766 | |
---|
1767 | if (scavenging) then |
---|
1768 | pdq(1:ngrid,1:nlayer,igcm_dust_mass) = |
---|
1769 | & pdq(1:ngrid,1:nlayer,igcm_dust_mass) + |
---|
1770 | & zdqcloud(1:ngrid,1:nlayer,igcm_dust_mass) |
---|
1771 | pdq(1:ngrid,1:nlayer,igcm_dust_number) = |
---|
1772 | & pdq(1:ngrid,1:nlayer,igcm_dust_number) + |
---|
1773 | & zdqcloud(1:ngrid,1:nlayer,igcm_dust_number) |
---|
1774 | where (pq(:,:,igcm_dust_mass) + |
---|
1775 | & ptimestep*pdq(:,:,igcm_dust_mass) < 0.) |
---|
1776 | pdq(:,:,igcm_dust_mass) = |
---|
1777 | & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 |
---|
1778 | pdq(:,:,igcm_dust_number) = |
---|
1779 | & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 |
---|
1780 | end where |
---|
1781 | where (pq(:,:,igcm_dust_number) + |
---|
1782 | & ptimestep*pdq(:,:,igcm_dust_number) < 0.) |
---|
1783 | pdq(:,:,igcm_dust_mass) = |
---|
1784 | & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 |
---|
1785 | pdq(:,:,igcm_dust_number) = |
---|
1786 | & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 |
---|
1787 | end where |
---|
1788 | endif ! of if scavenging |
---|
1789 | |
---|
1790 | END IF ! of IF (water) |
---|
1791 | |
---|
1792 | c 9a bis. CO2 clouds (CL & JA) |
---|
1793 | c --------------------------------------- |
---|
1794 | c CO2 ice cloud condensation in the atmosphere |
---|
1795 | c ---------------------------------------- |
---|
1796 | c flag needed in callphys.def: |
---|
1797 | c co2clouds=.true. is mandatory (default is .false.) |
---|
1798 | c co2useh2o=.true. if you want to allow co2 condensation |
---|
1799 | c on water ice particles |
---|
1800 | c meteo_flux=.true. if you want to add a meteoritic |
---|
1801 | c supply of CCN |
---|
1802 | c CLFvaryingCO2=.true. if you want to have a sub-grid |
---|
1803 | c temperature distribution |
---|
1804 | c spantCO2=integer (i.e. 3) amplitude of the sub-grid T disti |
---|
1805 | c nuiceco2_sed=0.2 variance of the size distribution for the |
---|
1806 | c sedimentation |
---|
1807 | c nuiceco2_ref=0.2 variance of the size distribution for the |
---|
1808 | c nucleation |
---|
1809 | c imicroco2=50 micro-timestep is 1/50 of physical timestep |
---|
1810 | zdqssed_co2(:) = 0. |
---|
1811 | zdqssed_ccn(:,:) = 0. |
---|
1812 | |
---|
1813 | IF (co2clouds) THEN |
---|
1814 | call co2cloud(ngrid,nlayer,ptimestep, |
---|
1815 | & zplev,zplay,pdpsrf,zzlay,pt,pdt, |
---|
1816 | & pq,pdq,zdqcloudco2,zdtcloudco2, |
---|
1817 | & nq,tau,tauscaling,rdust,rice,riceco2,nuice, |
---|
1818 | & rhocloud, rsedcloudco2,rhocloudco2,zzlev,zdqssed_co2, |
---|
1819 | & zdqssed_ccn,pdu,pu,zcondicea_co2microp) |
---|
1820 | |
---|
1821 | DO iq=1, nq |
---|
1822 | DO ig=1,ngrid |
---|
1823 | DO islope = 1,nslope |
---|
1824 | dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+ |
---|
1825 | & zdqssed_ccn(ig,iq)*cos(pi*def_slope_mean(islope)/180.) |
---|
1826 | ENDDO !(islope) |
---|
1827 | ENDDO ! (ig) |
---|
1828 | ENDDO ! (iq)q) |
---|
1829 | c Temperature variation due to latent heat release |
---|
1830 | pdt(1:ngrid,1:nlayer) = |
---|
1831 | & pdt(1:ngrid,1:nlayer) + |
---|
1832 | & zdtcloudco2(1:ngrid,1:nlayer) |
---|
1833 | |
---|
1834 | ! increment dust and ccn masses and numbers |
---|
1835 | ! We need to check that we have Nccn & Ndust > 0 |
---|
1836 | ! This is due to single precision rounding problems |
---|
1837 | ! increment dust tracers tendancies |
---|
1838 | pdq(:,:,igcm_dust_mass) = pdq(:,:,igcm_dust_mass) |
---|
1839 | & + zdqcloudco2(:,:,igcm_dust_mass) |
---|
1840 | |
---|
1841 | pdq(:,:,igcm_dust_number) = pdq(:,:,igcm_dust_number) |
---|
1842 | & + zdqcloudco2(:,:,igcm_dust_number) |
---|
1843 | |
---|
1844 | pdq(:,:,igcm_co2) = pdq(:,:,igcm_co2) |
---|
1845 | & + zdqcloudco2(:,:,igcm_co2) |
---|
1846 | |
---|
1847 | pdq(:,:,igcm_co2_ice) = pdq(:,:,igcm_co2_ice) |
---|
1848 | & + zdqcloudco2(:,:,igcm_co2_ice) |
---|
1849 | |
---|
1850 | pdq(:,:,igcm_ccnco2_mass) = pdq(:,:,igcm_ccnco2_mass) |
---|
1851 | & + zdqcloudco2(:,:,igcm_ccnco2_mass) |
---|
1852 | |
---|
1853 | pdq(:,:,igcm_ccnco2_number) = pdq(:,:,igcm_ccnco2_number) |
---|
1854 | & + zdqcloudco2(:,:,igcm_ccnco2_number) |
---|
1855 | |
---|
1856 | if (meteo_flux) then |
---|
1857 | pdq(:,:,igcm_ccnco2_meteor_mass) = |
---|
1858 | & pdq(:,:,igcm_ccnco2_meteor_mass) + |
---|
1859 | & zdqcloudco2(:,:,igcm_ccnco2_meteor_mass) |
---|
1860 | |
---|
1861 | pdq(:,:,igcm_ccnco2_meteor_number) = |
---|
1862 | & pdq(:,:,igcm_ccnco2_meteor_number) |
---|
1863 | & + zdqcloudco2(:,:,igcm_ccnco2_meteor_number) |
---|
1864 | end if |
---|
1865 | !Update water ice clouds values as well |
---|
1866 | if (co2useh2o) then |
---|
1867 | pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = |
---|
1868 | & pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + |
---|
1869 | & zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
1870 | pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = |
---|
1871 | & pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + |
---|
1872 | & zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass) |
---|
1873 | pdq(1:ngrid,1:nlayer,igcm_ccn_number) = |
---|
1874 | & pdq(1:ngrid,1:nlayer,igcm_ccn_number) + |
---|
1875 | & zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number) |
---|
1876 | |
---|
1877 | pdq(:,:,igcm_ccnco2_h2o_mass_ice) = |
---|
1878 | & pdq(:,:,igcm_ccnco2_h2o_mass_ice) + |
---|
1879 | & zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice) |
---|
1880 | |
---|
1881 | pdq(:,:,igcm_ccnco2_h2o_mass_ccn) = |
---|
1882 | & pdq(:,:,igcm_ccnco2_h2o_mass_ccn) + |
---|
1883 | & zdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn) |
---|
1884 | |
---|
1885 | pdq(:,:,igcm_ccnco2_h2o_number) = |
---|
1886 | & pdq(:,:,igcm_ccnco2_h2o_number) + |
---|
1887 | & zdqcloudco2(:,:,igcm_ccnco2_h2o_number) |
---|
1888 | |
---|
1889 | c Negative values? |
---|
1890 | where (pq(:,:,igcm_ccn_mass) + |
---|
1891 | & ptimestep*pdq(:,:,igcm_ccn_mass) < 0.) |
---|
1892 | pdq(:,:,igcm_ccn_mass) = |
---|
1893 | & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 |
---|
1894 | pdq(:,:,igcm_ccn_number) = |
---|
1895 | & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 |
---|
1896 | end where |
---|
1897 | c Negative values? |
---|
1898 | where (pq(:,:,igcm_ccn_number) + |
---|
1899 | & ptimestep*pdq(:,:,igcm_ccn_number) < 0.) |
---|
1900 | pdq(:,:,igcm_ccn_mass) = |
---|
1901 | & - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30 |
---|
1902 | pdq(:,:,igcm_ccn_number) = |
---|
1903 | & - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30 |
---|
1904 | end where |
---|
1905 | where (pq(:,:,igcm_ccnco2_h2o_mass_ice) + |
---|
1906 | & pq(:,:,igcm_ccnco2_h2o_mass_ccn) + |
---|
1907 | & (pdq(:,:,igcm_ccnco2_h2o_mass_ice) + |
---|
1908 | & pdq(:,:,igcm_ccnco2_h2o_mass_ccn) |
---|
1909 | & )*ptimestep < 0.) |
---|
1910 | pdq(:,:,igcm_ccnco2_h2o_mass_ice) = |
---|
1911 | & - pq(:,:,igcm_ccnco2_h2o_mass_ice) |
---|
1912 | & /ptimestep + 1.e-30 |
---|
1913 | pdq(:,:,igcm_ccnco2_h2o_mass_ccn) = |
---|
1914 | & - pq(:,:,igcm_ccnco2_h2o_mass_ccn) |
---|
1915 | & /ptimestep + 1.e-30 |
---|
1916 | pdq(:,:,igcm_ccnco2_h2o_number) = |
---|
1917 | & - pq(:,:,igcm_ccnco2_h2o_number) |
---|
1918 | & /ptimestep + 1.e-30 |
---|
1919 | end where |
---|
1920 | |
---|
1921 | where (pq(:,:,igcm_ccnco2_h2o_number) + |
---|
1922 | & (pdq(:,:,igcm_ccnco2_h2o_number) |
---|
1923 | & )*ptimestep < 0.) |
---|
1924 | pdq(:,:,igcm_ccnco2_h2o_mass_ice) = |
---|
1925 | & - pq(:,:,igcm_ccnco2_h2o_mass_ice) |
---|
1926 | & /ptimestep + 1.e-30 |
---|
1927 | pdq(:,:,igcm_ccnco2_h2o_mass_ccn) = |
---|
1928 | & - pq(:,:,igcm_ccnco2_h2o_mass_ccn) |
---|
1929 | & /ptimestep + 1.e-30 |
---|
1930 | pdq(:,:,igcm_ccnco2_h2o_number) = |
---|
1931 | & - pq(:,:,igcm_ccnco2_h2o_number) |
---|
1932 | & /ptimestep + 1.e-30 |
---|
1933 | end where |
---|
1934 | endif ! of if (co2useh2o) |
---|
1935 | c Negative values? |
---|
1936 | where (pq(:,:,igcm_ccnco2_mass) + |
---|
1937 | & ptimestep*pdq(:,:,igcm_ccnco2_mass) < 0.) |
---|
1938 | pdq(:,:,igcm_ccnco2_mass) = |
---|
1939 | & - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30 |
---|
1940 | pdq(:,:,igcm_ccnco2_number) = |
---|
1941 | & - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30 |
---|
1942 | end where |
---|
1943 | where (pq(:,:,igcm_ccnco2_number) + |
---|
1944 | & ptimestep*pdq(:,:,igcm_ccnco2_number) < 0.) |
---|
1945 | pdq(:,:,igcm_ccnco2_mass) = |
---|
1946 | & - pq(:,:,igcm_ccnco2_mass)/ptimestep + 1.e-30 |
---|
1947 | pdq(:,:,igcm_ccnco2_number) = |
---|
1948 | & - pq(:,:,igcm_ccnco2_number)/ptimestep + 1.e-30 |
---|
1949 | end where |
---|
1950 | |
---|
1951 | c Negative values? |
---|
1952 | where (pq(:,:,igcm_dust_mass) + |
---|
1953 | & ptimestep*pdq(:,:,igcm_dust_mass) < 0.) |
---|
1954 | pdq(:,:,igcm_dust_mass) = |
---|
1955 | & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 |
---|
1956 | pdq(:,:,igcm_dust_number) = |
---|
1957 | & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 |
---|
1958 | end where |
---|
1959 | where (pq(:,:,igcm_dust_number) + |
---|
1960 | & ptimestep*pdq(:,:,igcm_dust_number) < 0.) |
---|
1961 | pdq(:,:,igcm_dust_mass) = |
---|
1962 | & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 |
---|
1963 | pdq(:,:,igcm_dust_number) = |
---|
1964 | & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 |
---|
1965 | end where |
---|
1966 | if (meteo_flux) then |
---|
1967 | where (pq(:,:,igcm_ccnco2_meteor_mass) + |
---|
1968 | & ptimestep*pdq(:,:,igcm_ccnco2_meteor_mass) < 0.) |
---|
1969 | pdq(:,:,igcm_ccnco2_meteor_mass) = |
---|
1970 | & - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30 |
---|
1971 | pdq(:,:,igcm_ccnco2_meteor_number) = |
---|
1972 | & - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30 |
---|
1973 | end where |
---|
1974 | where (pq(:,:,igcm_ccnco2_meteor_number) + |
---|
1975 | & ptimestep*pdq(:,:,igcm_ccnco2_meteor_number) < 0.) |
---|
1976 | pdq(:,:,igcm_ccnco2_meteor_mass) = |
---|
1977 | & - pq(:,:,igcm_ccnco2_meteor_mass)/ptimestep + 1.e-30 |
---|
1978 | pdq(:,:,igcm_ccnco2_meteor_number) = |
---|
1979 | & - pq(:,:,igcm_ccnco2_meteor_number)/ptimestep + 1.e-30 |
---|
1980 | end where |
---|
1981 | end if |
---|
1982 | END IF ! of IF (co2clouds) |
---|
1983 | |
---|
1984 | c 9b. Aerosol particles |
---|
1985 | c ------------------- |
---|
1986 | c ---------- |
---|
1987 | c Dust devil : |
---|
1988 | c ---------- |
---|
1989 | IF(callddevil) then |
---|
1990 | call dustdevil(ngrid,nlayer,nq, zplev,pu,pv,pt, tsurf,q2, |
---|
1991 | & zdqdev,zdqsdev) |
---|
1992 | |
---|
1993 | if (dustbin.ge.1) then |
---|
1994 | do iq=1,nq |
---|
1995 | DO l=1,nlayer |
---|
1996 | DO ig=1,ngrid |
---|
1997 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdev(ig,l,iq) |
---|
1998 | ENDDO |
---|
1999 | ENDDO |
---|
2000 | enddo |
---|
2001 | do iq=1,nq |
---|
2002 | DO ig=1,ngrid |
---|
2003 | DO islope = 1,nslope |
---|
2004 | dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) + |
---|
2005 | & zdqsdev(ig,iq)*cos(pi*def_slope_mean(islope)/180.) |
---|
2006 | ENDDO |
---|
2007 | ENDDO |
---|
2008 | enddo |
---|
2009 | endif ! of if (dustbin.ge.1) |
---|
2010 | |
---|
2011 | END IF ! of IF (callddevil) |
---|
2012 | |
---|
2013 | c ------------- |
---|
2014 | c Sedimentation : acts also on water ice |
---|
2015 | c ------------- |
---|
2016 | IF (sedimentation) THEN |
---|
2017 | zdqsed(1:ngrid,1:nlayer,1:nq)=0 |
---|
2018 | zdqssed(1:ngrid,1:nq)=0 |
---|
2019 | |
---|
2020 | c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep |
---|
2021 | c Zdqssed isn't |
---|
2022 | |
---|
2023 | call callsedim(ngrid,nlayer,ptimestep, |
---|
2024 | & zplev,zzlev,zzlay,pt,pdt, |
---|
2025 | & rdust,rstormdust,rtopdust, |
---|
2026 | & rice,rsedcloud,rhocloud, |
---|
2027 | & pq,pdq,zdqsed,zdqssed,nq, |
---|
2028 | & tau,tauscaling) |
---|
2029 | |
---|
2030 | c Flux at the surface of co2 ice computed in co2cloud microtimestep |
---|
2031 | IF (rdstorm) THEN |
---|
2032 | c Storm dust cannot sediment to the surface |
---|
2033 | DO ig=1,ngrid |
---|
2034 | zdqsed(ig,1,igcm_stormdust_mass)= |
---|
2035 | & zdqsed(ig,1,igcm_stormdust_mass)+ |
---|
2036 | & zdqssed(ig,igcm_stormdust_mass) / |
---|
2037 | & ((pplev(ig,1)-pplev(ig,2))/g) |
---|
2038 | zdqsed(ig,1,igcm_stormdust_number)= |
---|
2039 | & zdqsed(ig,1,igcm_stormdust_number)+ |
---|
2040 | & zdqssed(ig,igcm_stormdust_number) / |
---|
2041 | & ((pplev(ig,1)-pplev(ig,2))/g) |
---|
2042 | zdqssed(ig,igcm_stormdust_mass)=0. |
---|
2043 | zdqssed(ig,igcm_stormdust_number)=0. |
---|
2044 | ENDDO |
---|
2045 | ENDIF !rdstorm |
---|
2046 | |
---|
2047 | DO iq=1, nq |
---|
2048 | DO l=1,nlayer |
---|
2049 | DO ig=1,ngrid |
---|
2050 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqsed(ig,l,iq) |
---|
2051 | ENDDO |
---|
2052 | ENDDO |
---|
2053 | ENDDO |
---|
2054 | DO iq=1, nq |
---|
2055 | DO ig=1,ngrid |
---|
2056 | DO islope = 1,nslope |
---|
2057 | dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) + |
---|
2058 | & zdqssed(ig,iq)*cos(pi*def_slope_mean(islope)/180.) |
---|
2059 | ENDDO |
---|
2060 | ENDDO |
---|
2061 | ENDDO |
---|
2062 | |
---|
2063 | END IF ! of IF (sedimentation) |
---|
2064 | |
---|
2065 | c Add lifted dust to tendancies after sedimentation in the LES (AC) |
---|
2066 | IF (turb_resolved) THEN |
---|
2067 | DO iq=1, nq |
---|
2068 | DO l=1,nlayer |
---|
2069 | DO ig=1,ngrid |
---|
2070 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqdif(ig,l,iq) |
---|
2071 | ENDDO |
---|
2072 | ENDDO |
---|
2073 | ENDDO |
---|
2074 | DO iq=1, nq |
---|
2075 | DO ig=1,ngrid |
---|
2076 | dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:) |
---|
2077 | ENDDO |
---|
2078 | ENDDO |
---|
2079 | ENDIF |
---|
2080 | c |
---|
2081 | c 9c. Chemical species |
---|
2082 | c ------------------ |
---|
2083 | |
---|
2084 | #ifndef MESOSCALE |
---|
2085 | c -------------- |
---|
2086 | c photochemistry : |
---|
2087 | c -------------- |
---|
2088 | IF (photochem) then |
---|
2089 | |
---|
2090 | if (modulo(icount-1,ichemistry).eq.0) then |
---|
2091 | ! compute chemistry every ichemistry physics step |
---|
2092 | |
---|
2093 | ! dust and ice surface area |
---|
2094 | call surfacearea(ngrid, nlayer, naerkind, |
---|
2095 | $ ptimestep, zplay, zzlay, |
---|
2096 | $ pt, pq, pdq, nq, |
---|
2097 | $ rdust, rice, tau, tauscaling, |
---|
2098 | $ surfdust, surfice) |
---|
2099 | ! call photochemistry |
---|
2100 | DO ig = 1,ngrid |
---|
2101 | qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig)) |
---|
2102 | ENDDO |
---|
2103 | call calchim(ngrid,nlayer,nq, |
---|
2104 | & ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0, |
---|
2105 | $ zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim, |
---|
2106 | $ zdqcloud,zdqscloud,tau(:,1), |
---|
2107 | $ qsurf_tmp(:,igcm_co2), |
---|
2108 | $ pu,pdu,pv,pdv,surfdust,surfice) |
---|
2109 | |
---|
2110 | endif ! of if (modulo(icount-1,ichemistry).eq.0) |
---|
2111 | |
---|
2112 | ! increment values of tracers: |
---|
2113 | DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry |
---|
2114 | ! tracers is zero anyways |
---|
2115 | DO l=1,nlayer |
---|
2116 | DO ig=1,ngrid |
---|
2117 | pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) |
---|
2118 | ENDDO |
---|
2119 | ENDDO |
---|
2120 | ENDDO ! of DO iq=1,nq |
---|
2121 | |
---|
2122 | ! add condensation tendency for H2O2 |
---|
2123 | if (igcm_h2o2.ne.0) then |
---|
2124 | DO l=1,nlayer |
---|
2125 | DO ig=1,ngrid |
---|
2126 | pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2) |
---|
2127 | & +zdqcloud(ig,l,igcm_h2o2) |
---|
2128 | ENDDO |
---|
2129 | ENDDO |
---|
2130 | endif |
---|
2131 | |
---|
2132 | ! increment surface values of tracers: |
---|
2133 | DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry |
---|
2134 | ! tracers is zero anyways |
---|
2135 | DO ig=1,ngrid |
---|
2136 | DO islope = 1,nslope |
---|
2137 | dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) + |
---|
2138 | & zdqschim(ig,iq)*cos(pi*def_slope_mean(islope)/180.) |
---|
2139 | ENDDO |
---|
2140 | ENDDO |
---|
2141 | ENDDO ! of DO iq=1,nq |
---|
2142 | |
---|
2143 | ! add condensation tendency for H2O2 |
---|
2144 | if (igcm_h2o2.ne.0) then |
---|
2145 | DO ig=1,ngrid |
---|
2146 | DO islope = 1,nslope |
---|
2147 | dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+ |
---|
2148 | & zdqscloud(ig,igcm_h2o2)*cos(pi*def_slope_mean(islope)/180.) |
---|
2149 | ENDDO |
---|
2150 | ENDDO |
---|
2151 | endif |
---|
2152 | |
---|
2153 | END IF ! of IF (photochem) |
---|
2154 | #endif |
---|
2155 | |
---|
2156 | |
---|
2157 | #ifndef MESOSCALE |
---|
2158 | c----------------------------------------------------------------------- |
---|
2159 | c 10. THERMOSPHERE CALCULATION |
---|
2160 | c----------------------------------------------------------------------- |
---|
2161 | |
---|
2162 | if (callthermos) then |
---|
2163 | call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol, |
---|
2164 | $ mu0,ptimestep,ptime,zday,tsurf_meshavg,zzlev,zzlay, |
---|
2165 | & pt,pq,pu,pv,pdt,pdq, |
---|
2166 | $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff, |
---|
2167 | $ PhiEscH,PhiEscH2,PhiEscD) |
---|
2168 | |
---|
2169 | DO l=1,nlayer |
---|
2170 | DO ig=1,ngrid |
---|
2171 | dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l) |
---|
2172 | pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l)+zdteuv(ig,l) |
---|
2173 | pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l) |
---|
2174 | pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l) |
---|
2175 | DO iq=1, nq |
---|
2176 | pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq) |
---|
2177 | ENDDO |
---|
2178 | ENDDO |
---|
2179 | ENDDO |
---|
2180 | |
---|
2181 | endif ! of if (callthermos) |
---|
2182 | #endif |
---|
2183 | |
---|
2184 | c----------------------------------------------------------------------- |
---|
2185 | c 11. Carbon dioxide condensation-sublimation: |
---|
2186 | c (should be the last atmospherical physical process to be computed) |
---|
2187 | c ------------------------------------------- |
---|
2188 | IF (tituscap) THEN |
---|
2189 | !!! get the actual co2 seasonal cap from Titus observations |
---|
2190 | CALL geticecover(ngrid, 180.*zls/pi, |
---|
2191 | . 180.*longitude/pi, 180.*latitude/pi, |
---|
2192 | . qsurf_tmp(:,igcm_co2) ) |
---|
2193 | qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000. |
---|
2194 | ENDIF |
---|
2195 | |
---|
2196 | |
---|
2197 | IF (callcond) THEN |
---|
2198 | zdtc(:,:) = 0. |
---|
2199 | zdtsurfc(:,:) = 0. |
---|
2200 | zduc(:,:) = 0. |
---|
2201 | zdvc(:,:) = 0. |
---|
2202 | zdqc(:,:,:) = 0. |
---|
2203 | zdqsc(:,:,:) = 0. |
---|
2204 | CALL co2condens(ngrid,nlayer,nq,nslope,ptimestep, |
---|
2205 | $ capcal,zplay,zplev,tsurf,pt, |
---|
2206 | $ pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq, |
---|
2207 | $ qsurf(:,igcm_co2,:),perenial_co2ice, |
---|
2208 | $ albedo,emis,rdust, |
---|
2209 | $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, |
---|
2210 | $ fluxsurf_dn_sw,zls, |
---|
2211 | $ zdqssed_co2,zcondicea_co2microp, |
---|
2212 | & zdqsc) |
---|
2213 | |
---|
2214 | DO iq=1, nq |
---|
2215 | DO ig=1,ngrid |
---|
2216 | dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:) |
---|
2217 | ENDDO ! (ig) |
---|
2218 | ENDDO ! (iq) |
---|
2219 | DO l=1,nlayer |
---|
2220 | DO ig=1,ngrid |
---|
2221 | pdt(ig,l)=pdt(ig,l)+zdtc(ig,l) |
---|
2222 | pdv(ig,l)=pdv(ig,l)+zdvc(ig,l) |
---|
2223 | pdu(ig,l)=pdu(ig,l)+zduc(ig,l) |
---|
2224 | ENDDO |
---|
2225 | ENDDO |
---|
2226 | DO ig=1,ngrid |
---|
2227 | zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:) |
---|
2228 | ENDDO |
---|
2229 | |
---|
2230 | DO iq=1, nq |
---|
2231 | DO l=1,nlayer |
---|
2232 | DO ig=1,ngrid |
---|
2233 | pdq(ig,l,iq)=pdq(ig,l,iq)+ zdqc(ig,l,iq) |
---|
2234 | ENDDO |
---|
2235 | ENDDO |
---|
2236 | ENDDO |
---|
2237 | |
---|
2238 | #ifndef MESOSCALE |
---|
2239 | ! update surface pressure |
---|
2240 | DO ig=1,ngrid |
---|
2241 | ps(ig) = zplev(ig,1) + pdpsrf(ig)*ptimestep |
---|
2242 | ENDDO |
---|
2243 | ! update pressure levels |
---|
2244 | DO l=1,nlayer |
---|
2245 | DO ig=1,ngrid |
---|
2246 | zplay(ig,l) = aps(l) + bps(l)*ps(ig) |
---|
2247 | zplev(ig,l) = ap(l) + bp(l)*ps(ig) |
---|
2248 | ENDDO |
---|
2249 | ENDDO |
---|
2250 | zplev(:,nlayer+1) = 0. |
---|
2251 | ! update layers altitude |
---|
2252 | DO l=2,nlayer |
---|
2253 | DO ig=1,ngrid |
---|
2254 | z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l)) |
---|
2255 | z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l)) |
---|
2256 | zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) |
---|
2257 | ENDDO |
---|
2258 | ENDDO |
---|
2259 | #endif |
---|
2260 | ENDIF ! of IF (callcond) |
---|
2261 | |
---|
2262 | c----------------------------------------------------------------------- |
---|
2263 | c Updating tracer budget on surface |
---|
2264 | c----------------------------------------------------------------------- |
---|
2265 | DO iq=1, nq |
---|
2266 | DO ig=1,ngrid |
---|
2267 | DO islope = 1,nslope |
---|
2268 | qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+ |
---|
2269 | & ptimestep*dqsurf(ig,iq,islope) |
---|
2270 | ENDDO |
---|
2271 | ENDDO ! (ig) |
---|
2272 | ENDDO ! (iq) |
---|
2273 | c----------------------------------------------------------------------- |
---|
2274 | c 12. Surface and sub-surface soil temperature |
---|
2275 | c----------------------------------------------------------------------- |
---|
2276 | c |
---|
2277 | c |
---|
2278 | c 12.1 Increment Surface temperature: |
---|
2279 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
2280 | |
---|
2281 | DO ig=1,ngrid |
---|
2282 | DO islope = 1,nslope |
---|
2283 | tsurf(ig,islope)=tsurf(ig,islope)+ |
---|
2284 | & ptimestep*zdtsurf(ig,islope) |
---|
2285 | ENDDO |
---|
2286 | ENDDO |
---|
2287 | |
---|
2288 | c Prescribe a cold trap at south pole (except at high obliquity !!) |
---|
2289 | c Temperature at the surface is set there to be the temperature |
---|
2290 | c corresponding to equilibrium temperature between phases of CO2 |
---|
2291 | |
---|
2292 | |
---|
2293 | IF (water) THEN |
---|
2294 | !#ifndef MESOSCALE |
---|
2295 | ! if (caps.and.(obliquit.lt.27.)) then => now done in co2condens |
---|
2296 | ! NB: Updated surface pressure, at grid point 'ngrid', is |
---|
2297 | ! ps(ngrid)=zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep |
---|
2298 | ! tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095* |
---|
2299 | ! & (zplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) |
---|
2300 | ! tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095*ps(ngrid))) |
---|
2301 | ! endif |
---|
2302 | !#endif |
---|
2303 | c ------------------------------------------------------------- |
---|
2304 | c Change of surface albedo in case of ground frost |
---|
2305 | c everywhere except on the north permanent cap and in regions |
---|
2306 | c covered by dry ice. |
---|
2307 | c ALWAYS PLACE these lines after co2condens !!! |
---|
2308 | c ------------------------------------------------------------- |
---|
2309 | do ig=1,ngrid |
---|
2310 | do islope = 1,nslope |
---|
2311 | if ((qsurf(ig,igcm_co2,islope).eq.0).and. |
---|
2312 | & (qsurf(ig,igcm_h2o_ice,islope) |
---|
2313 | & .gt.frost_albedo_threshold)) then |
---|
2314 | if ((watercaptag(ig)).and.(cst_cap_albedo)) then |
---|
2315 | albedo(ig,1,islope) = albedo_h2o_cap |
---|
2316 | albedo(ig,2,islope) = albedo_h2o_cap |
---|
2317 | else |
---|
2318 | albedo(ig,1,islope) = albedo_h2o_frost |
---|
2319 | albedo(ig,2,islope) = albedo_h2o_frost |
---|
2320 | endif !((watercaptag(ig)).and.(cst_cap_albedo)) then |
---|
2321 | c write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice) |
---|
2322 | c write(*,*) "physiq.F frost :" |
---|
2323 | c & ,latitude(ig)*180./pi, longitude(ig)*180./pi |
---|
2324 | endif |
---|
2325 | enddo ! islope |
---|
2326 | enddo ! of do ig=1,ngrid |
---|
2327 | ENDIF ! of IF (water) |
---|
2328 | |
---|
2329 | c |
---|
2330 | c 12.2 Compute soil temperatures and subsurface heat flux: |
---|
2331 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
2332 | IF (callsoil) THEN |
---|
2333 | c Thermal inertia feedback |
---|
2334 | IF (tifeedback) THEN |
---|
2335 | DO islope = 1,nslope |
---|
2336 | CALL soil_tifeedback(ngrid,nsoilmx, |
---|
2337 | s qsurf(:,:,islope), |
---|
2338 | s inertiesoil_tifeedback(:,:,islope)) |
---|
2339 | ENDDO |
---|
2340 | CALL soil(ngrid,nsoilmx,.false.,inertiesoil_tifeedback, |
---|
2341 | s ptimestep,tsurf,tsoil,capcal,fluxgrd) |
---|
2342 | ELSE |
---|
2343 | CALL soil(ngrid,nsoilmx,.false.,inertiesoil, |
---|
2344 | s ptimestep,tsurf,tsoil,capcal,fluxgrd) |
---|
2345 | ENDIF |
---|
2346 | ENDIF |
---|
2347 | |
---|
2348 | c To avoid negative values |
---|
2349 | IF (rdstorm) THEN |
---|
2350 | where (pq(:,:,igcm_stormdust_mass) + |
---|
2351 | & ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.) |
---|
2352 | pdq(:,:,igcm_stormdust_mass) = |
---|
2353 | & - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30 |
---|
2354 | pdq(:,:,igcm_stormdust_number) = |
---|
2355 | & - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30 |
---|
2356 | end where |
---|
2357 | where (pq(:,:,igcm_stormdust_number) + |
---|
2358 | & ptimestep*pdq(:,:,igcm_stormdust_number) < 0.) |
---|
2359 | pdq(:,:,igcm_stormdust_mass) = |
---|
2360 | & - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30 |
---|
2361 | pdq(:,:,igcm_stormdust_number) = |
---|
2362 | & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 |
---|
2363 | end where |
---|
2364 | |
---|
2365 | where (pq(:,:,igcm_dust_mass) + |
---|
2366 | & ptimestep*pdq(:,:,igcm_dust_mass) < 0.) |
---|
2367 | pdq(:,:,igcm_dust_mass) = |
---|
2368 | & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 |
---|
2369 | pdq(:,:,igcm_dust_number) = |
---|
2370 | & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 |
---|
2371 | end where |
---|
2372 | where (pq(:,:,igcm_dust_number) + |
---|
2373 | & ptimestep*pdq(:,:,igcm_dust_number) < 0.) |
---|
2374 | pdq(:,:,igcm_dust_mass) = |
---|
2375 | & - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30 |
---|
2376 | pdq(:,:,igcm_dust_number) = |
---|
2377 | & - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30 |
---|
2378 | end where |
---|
2379 | ENDIF !(rdstorm) |
---|
2380 | |
---|
2381 | c----------------------------------------------------------------------- |
---|
2382 | c J. Naar : Surface and sub-surface water ice |
---|
2383 | c----------------------------------------------------------------------- |
---|
2384 | c |
---|
2385 | c |
---|
2386 | c Increment Watercap (surface h2o reservoirs): |
---|
2387 | c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
2388 | |
---|
2389 | DO ig=1,ngrid |
---|
2390 | DO islope = 1,nslope |
---|
2391 | watercap(ig,islope)=watercap(ig,islope)+ |
---|
2392 | s ptimestep*dwatercap(ig,islope) |
---|
2393 | ENDDO |
---|
2394 | ENDDO |
---|
2395 | |
---|
2396 | IF (refill_watercap) THEN |
---|
2397 | |
---|
2398 | DO ig=1,ngrid |
---|
2399 | DO islope = 1,nslope |
---|
2400 | if (watercaptag(ig).and. (qsurf(ig,igcm_h2o_ice,islope) |
---|
2401 | & .gt.frost_metam_threshold)) then |
---|
2402 | |
---|
2403 | watercap(ig,islope)=watercap(ig,islope) |
---|
2404 | & +qsurf(ig,igcm_h2o_ice,islope) |
---|
2405 | & - frost_metam_threshold |
---|
2406 | qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold |
---|
2407 | endif ! (watercaptag(ig).and. |
---|
2408 | ENDDO |
---|
2409 | ENDDO |
---|
2410 | |
---|
2411 | ENDIF ! (refill_watercap) THEN |
---|
2412 | |
---|
2413 | |
---|
2414 | c----------------------------------------------------------------------- |
---|
2415 | c 13. Write output files |
---|
2416 | c ---------------------- |
---|
2417 | |
---|
2418 | call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf, |
---|
2419 | & albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg) |
---|
2420 | |
---|
2421 | c ------------------------------- |
---|
2422 | c Dynamical fields incrementation |
---|
2423 | c ------------------------------- |
---|
2424 | c (FOR OUTPUT ONLY : the actual model integration is performed in the dynamics) |
---|
2425 | ! temperature, zonal and meridional wind |
---|
2426 | DO l=1,nlayer |
---|
2427 | DO ig=1,ngrid |
---|
2428 | zt(ig,l)=pt(ig,l) + pdt(ig,l)*ptimestep |
---|
2429 | zu(ig,l)=pu(ig,l) + pdu(ig,l)*ptimestep |
---|
2430 | zv(ig,l)=pv(ig,l) + pdv(ig,l)*ptimestep |
---|
2431 | ENDDO |
---|
2432 | ENDDO |
---|
2433 | |
---|
2434 | ! tracers |
---|
2435 | DO iq=1, nq |
---|
2436 | DO l=1,nlayer |
---|
2437 | DO ig=1,ngrid |
---|
2438 | zq(ig,l,iq)=pq(ig,l,iq) +pdq(ig,l,iq)*ptimestep |
---|
2439 | ENDDO |
---|
2440 | ENDDO |
---|
2441 | ENDDO |
---|
2442 | |
---|
2443 | ! Density |
---|
2444 | DO l=1,nlayer |
---|
2445 | DO ig=1,ngrid |
---|
2446 | rho(ig,l) = zplay(ig,l)/(rnew(ig,l)*zt(ig,l)) |
---|
2447 | ENDDO |
---|
2448 | ENDDO |
---|
2449 | |
---|
2450 | ! Potential Temperature |
---|
2451 | |
---|
2452 | DO ig=1,ngrid |
---|
2453 | DO l=1,nlayer |
---|
2454 | zh(ig,l) = zt(ig,l)*(zplev(ig,1)/zplay(ig,l))**rcp |
---|
2455 | ENDDO |
---|
2456 | ENDDO |
---|
2457 | |
---|
2458 | c Compute surface stress : (NB: z0 is a common in surfdat.h) |
---|
2459 | c DO ig=1,ngrid |
---|
2460 | c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2 |
---|
2461 | c zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2) |
---|
2462 | c ENDDO |
---|
2463 | |
---|
2464 | c Sum of fluxes in solar spectral bands (for output only) |
---|
2465 | fluxtop_dn_sw_tot(1:ngrid)=fluxtop_dn_sw(1:ngrid,1) + |
---|
2466 | & fluxtop_dn_sw(1:ngrid,2) |
---|
2467 | fluxtop_up_sw_tot(1:ngrid)=fluxtop_up_sw(1:ngrid,1) + |
---|
2468 | & fluxtop_up_sw(1:ngrid,2) |
---|
2469 | fluxsurf_dn_sw_tot(1:ngrid,1:nslope)= |
---|
2470 | & fluxsurf_dn_sw(1:ngrid,1,1:nslope) + |
---|
2471 | & fluxsurf_dn_sw(1:ngrid,2,1:nslope) |
---|
2472 | fluxsurf_up_sw_tot(1:ngrid)=fluxsurf_up_sw(1:ngrid,1) + |
---|
2473 | & fluxsurf_up_sw(1:ngrid,2) |
---|
2474 | |
---|
2475 | c ******* TEST ****************************************************** |
---|
2476 | ztim1 = 999 |
---|
2477 | DO l=1,nlayer |
---|
2478 | DO ig=1,ngrid |
---|
2479 | if (pt(ig,l).lt.ztim1) then |
---|
2480 | ztim1 = pt(ig,l) |
---|
2481 | igmin = ig |
---|
2482 | lmin = l |
---|
2483 | end if |
---|
2484 | ENDDO |
---|
2485 | ENDDO |
---|
2486 | if(min(pt(igmin,lmin),zt(igmin,lmin)).lt.70.) then |
---|
2487 | write(*,*) 'PHYSIQ: stability WARNING :' |
---|
2488 | write(*,*) 'pt, zt Tmin = ', pt(igmin,lmin), zt(igmin,lmin), |
---|
2489 | & 'ig l =', igmin, lmin |
---|
2490 | end if |
---|
2491 | c ******************************************************************* |
---|
2492 | |
---|
2493 | c --------------------- |
---|
2494 | c Outputs to the screen |
---|
2495 | c --------------------- |
---|
2496 | |
---|
2497 | IF (lwrite) THEN |
---|
2498 | PRINT*,'Global diagnostics for the physics' |
---|
2499 | PRINT*,'Variables and their increments x and dx/dt * dt' |
---|
2500 | WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps' |
---|
2501 | WRITE(*,'(2f10.5,2f15.5)') |
---|
2502 | s tsurf(igout,:),zdtsurf(igout,:)*ptimestep, |
---|
2503 | s zplev(igout,1),pdpsrf(igout)*ptimestep |
---|
2504 | WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT' |
---|
2505 | WRITE(*,'(i4,6f10.5)') (l, |
---|
2506 | s pu(igout,l),pdu(igout,l)*ptimestep, |
---|
2507 | s pv(igout,l),pdv(igout,l)*ptimestep, |
---|
2508 | s pt(igout,l),pdt(igout,l)*ptimestep, |
---|
2509 | s l=1,nlayer) |
---|
2510 | ENDIF ! of IF (lwrite) |
---|
2511 | |
---|
2512 | c ---------------------------------------------------------- |
---|
2513 | c ---------------------------------------------------------- |
---|
2514 | c INTERPOLATIONS IN THE SURFACE-LAYER |
---|
2515 | c ---------------------------------------------------------- |
---|
2516 | c ---------------------------------------------------------- |
---|
2517 | |
---|
2518 | n_out=0 ! number of elements in the z_out array. |
---|
2519 | ! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set |
---|
2520 | ! to 5 |
---|
2521 | IF (n_out .ne. 0) THEN |
---|
2522 | |
---|
2523 | IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out)) |
---|
2524 | IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out)) |
---|
2525 | IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out)) |
---|
2526 | |
---|
2527 | z_out(:)=[3.,2.,1.,0.5,0.1] |
---|
2528 | u_out(:,:)=0. |
---|
2529 | T_out(:,:)=0. |
---|
2530 | |
---|
2531 | call pbl_parameters(ngrid,nlayer,ps,zplay,z0, |
---|
2532 | & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf(:,iflat), |
---|
2533 | & zh,z_out,n_out,T_out,u_out,ustar,tstar,L_mo,vhf,vvv) |
---|
2534 | ! pourquoi ustar recalcule ici? fait dans vdifc. |
---|
2535 | |
---|
2536 | #ifndef MESOSCALE |
---|
2537 | DO n=1,n_out |
---|
2538 | write(zstring, '(F8.6)') z_out(n) |
---|
2539 | call write_output('T_out_'//trim(zstring), |
---|
2540 | & 'potential temperature at z_out','K',T_out(:,n)) |
---|
2541 | call write_output('u_out_'//trim(zstring), |
---|
2542 | & 'horizontal velocity norm at z_out','m/s',u_out(:,n)) |
---|
2543 | ENDDO |
---|
2544 | call write_output('u_star', |
---|
2545 | & 'friction velocity','m/s',ustar) |
---|
2546 | call write_output('teta_star', |
---|
2547 | & 'friction potential temperature','K',tstar) |
---|
2548 | call write_output('vvv', |
---|
2549 | & 'Vertical velocity variance at zout','m',vvv) |
---|
2550 | call write_output('vhf', |
---|
2551 | & 'Vertical heat flux at zout','m',vhf) |
---|
2552 | #else |
---|
2553 | T_out1(:)=T_out(:,1) |
---|
2554 | u_out1(:)=u_out(:,1) |
---|
2555 | #endif |
---|
2556 | |
---|
2557 | ENDIF |
---|
2558 | |
---|
2559 | c ---------------------------------------------------------- |
---|
2560 | c ---------------------------------------------------------- |
---|
2561 | c END OF SURFACE LAYER INTERPOLATIONS |
---|
2562 | c ---------------------------------------------------------- |
---|
2563 | c ---------------------------------------------------------- |
---|
2564 | |
---|
2565 | #ifndef MESOSCALE |
---|
2566 | c ------------------------------------------------------------------- |
---|
2567 | c Writing NetCDF file "RESTARTFI" at the end of the run |
---|
2568 | c ------------------------------------------------------------------- |
---|
2569 | c Note: 'restartfi' is stored just before dynamics are stored |
---|
2570 | c in 'restart'. Between now and the writting of 'restart', |
---|
2571 | c there will have been the itau=itau+1 instruction and |
---|
2572 | c a reset of 'time' (lastacll = .true. when itau+1= itaufin) |
---|
2573 | c thus we store for time=time+dtvr |
---|
2574 | |
---|
2575 | ! default: not writing a restart file at this time step |
---|
2576 | write_restart=.false. |
---|
2577 | IF (ecritstart.GT.0) THEN |
---|
2578 | ! For when we store multiple time steps in the restart file |
---|
2579 | IF (MODULO(icount*iphysiq,ecritstart).EQ.0) THEN |
---|
2580 | write_restart=.true. |
---|
2581 | ENDIF |
---|
2582 | ENDIF |
---|
2583 | IF (lastcall) THEN |
---|
2584 | ! Always write a restart at the end of the simulation |
---|
2585 | write_restart=.true. |
---|
2586 | ENDIF |
---|
2587 | |
---|
2588 | IF (write_restart) THEN |
---|
2589 | IF (grid_type==unstructured) THEN !IF DYNAMICO |
---|
2590 | |
---|
2591 | ! When running Dynamico, no need to add a dynamics time step to ztime_fin |
---|
2592 | IF (ptime.LE. 1.E-10) THEN |
---|
2593 | ! Residual ptime occurs with Dynamico |
---|
2594 | ztime_fin = pday !+ ptime + ptimestep/(float(iphysiq)*daysec) |
---|
2595 | . - day_ini - time_phys |
---|
2596 | ELSE |
---|
2597 | ztime_fin = pday + ptime !+ ptimestep/(float(iphysiq)*daysec) |
---|
2598 | . - day_ini - time_phys |
---|
2599 | ENDIF |
---|
2600 | if (ecritstart==0) then |
---|
2601 | ztime_fin = ztime_fin-(day_end-day_ini) |
---|
2602 | endif |
---|
2603 | |
---|
2604 | ELSE ! IF LMDZ |
---|
2605 | |
---|
2606 | if (ecritstart.GT.0) then !IF MULTIPLE RESTARTS nothing change |
---|
2607 | ztime_fin = pday - day_ini + ptime |
---|
2608 | . + ptimestep/(float(iphysiq)*daysec) |
---|
2609 | else !IF ONE RESTART final time in top of day_end |
---|
2610 | ztime_fin = pday - day_ini-(day_end-day_ini) |
---|
2611 | . + ptime + ptimestep/(float(iphysiq)*daysec) |
---|
2612 | endif |
---|
2613 | |
---|
2614 | ENDIF ! of IF (grid_type==unstructured) |
---|
2615 | write(*,'(A,I7,A,F12.5)') |
---|
2616 | . 'PHYSIQ: writing in restartfi ; icount=', |
---|
2617 | . icount,' date=',ztime_fin |
---|
2618 | |
---|
2619 | call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, |
---|
2620 | . ptimestep,ztime_fin, |
---|
2621 | . tsurf,tsoil,inertiesoil,albedo, |
---|
2622 | . emis,q2,qsurf,tauscaling,totcloudfrac,wstar, |
---|
2623 | . watercap,perenial_co2ice) |
---|
2624 | ENDIF ! of IF (write_restart) |
---|
2625 | |
---|
2626 | #endif |
---|
2627 | |
---|
2628 | c IF (ngrid.NE.1) then |
---|
2629 | |
---|
2630 | c ------------------------------------------------------------------- |
---|
2631 | c Calculation of diagnostic variables written in both stats and |
---|
2632 | c diagfi files |
---|
2633 | c ------------------------------------------------------------------- |
---|
2634 | |
---|
2635 | if (water) then |
---|
2636 | do ig=1,ngrid |
---|
2637 | if(mu0(ig).le.0.01) then |
---|
2638 | fluxsurf_dir_dn_sw(ig) = 0. |
---|
2639 | else |
---|
2640 | fluxsurf_dir_dn_sw(ig) = flux_1AU/dist_sol/dist_sol*mu0(ig)* |
---|
2641 | & exp(-(tau(ig,iaer_dust_doubleq)+ |
---|
2642 | & tau(ig,iaer_h2o_ice))/mu0(ig)) |
---|
2643 | endif |
---|
2644 | enddo |
---|
2645 | endif |
---|
2646 | |
---|
2647 | ! Density-scaled opacities |
---|
2648 | do ig=1,ngrid |
---|
2649 | dsodust(ig,:) = |
---|
2650 | & dsodust(ig,:)*tauscaling(ig) |
---|
2651 | dsords(ig,:) = |
---|
2652 | & dsords(ig,:)*tauscaling(ig) |
---|
2653 | dsotop(ig,:) = |
---|
2654 | & dsotop(ig,:)*tauscaling(ig) |
---|
2655 | enddo |
---|
2656 | |
---|
2657 | if(doubleq) then |
---|
2658 | do ig=1,ngrid |
---|
2659 | IF (sedimentation) THEN |
---|
2660 | dqdustsurf(ig) = |
---|
2661 | & zdqssed(ig,igcm_dust_mass)*tauscaling(ig) |
---|
2662 | dndustsurf(ig) = |
---|
2663 | & zdqssed(ig,igcm_dust_number)*tauscaling(ig) |
---|
2664 | ENDIF |
---|
2665 | ndust(ig,:) = |
---|
2666 | & zq(ig,:,igcm_dust_number)*tauscaling(ig) |
---|
2667 | qdust(ig,:) = |
---|
2668 | & zq(ig,:,igcm_dust_mass)*tauscaling(ig) |
---|
2669 | enddo |
---|
2670 | if (scavenging) then |
---|
2671 | do ig=1,ngrid |
---|
2672 | IF (sedimentation) THEN |
---|
2673 | dqdustsurf(ig) = dqdustsurf(ig) + |
---|
2674 | & zdqssed(ig,igcm_ccn_mass)*tauscaling(ig) |
---|
2675 | dndustsurf(ig) = dndustsurf(ig) + |
---|
2676 | & zdqssed(ig,igcm_ccn_number)*tauscaling(ig) |
---|
2677 | ENDIF |
---|
2678 | nccn(ig,:) = |
---|
2679 | & zq(ig,:,igcm_ccn_number)*tauscaling(ig) |
---|
2680 | qccn(ig,:) = |
---|
2681 | & zq(ig,:,igcm_ccn_mass)*tauscaling(ig) |
---|
2682 | enddo |
---|
2683 | endif |
---|
2684 | endif ! of (doubleq) |
---|
2685 | |
---|
2686 | if (rdstorm) then ! diagnostics of stormdust tendancies for 1D and 3D |
---|
2687 | mstormdtot(:)=0 |
---|
2688 | mdusttot(:)=0 |
---|
2689 | qdusttotal(:,:)=0 |
---|
2690 | do ig=1,ngrid |
---|
2691 | rdsdqdustsurf(ig) = |
---|
2692 | & zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig) |
---|
2693 | rdsdndustsurf(ig) = |
---|
2694 | & zdqssed(ig,igcm_stormdust_number)*tauscaling(ig) |
---|
2695 | rdsndust(ig,:) = |
---|
2696 | & pq(ig,:,igcm_stormdust_number)*tauscaling(ig) |
---|
2697 | rdsqdust(ig,:) = |
---|
2698 | & pq(ig,:,igcm_stormdust_mass)*tauscaling(ig) |
---|
2699 | do l=1,nlayer |
---|
2700 | mstormdtot(ig) = mstormdtot(ig) + |
---|
2701 | & zq(ig,l,igcm_stormdust_mass) * |
---|
2702 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2703 | mdusttot(ig) = mdusttot(ig) + |
---|
2704 | & zq(ig,l,igcm_dust_mass) * |
---|
2705 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2706 | qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust |
---|
2707 | enddo |
---|
2708 | enddo |
---|
2709 | endif !(rdstorm) |
---|
2710 | |
---|
2711 | if (water) then |
---|
2712 | mtot(:)=0 |
---|
2713 | icetot(:)=0 |
---|
2714 | rave(:)=0 |
---|
2715 | tauTES(:)=0 |
---|
2716 | |
---|
2717 | IF (hdo) then |
---|
2718 | mtotD(:)=0 |
---|
2719 | icetotD(:)=0 |
---|
2720 | ENDIF !hdo |
---|
2721 | |
---|
2722 | do ig=1,ngrid |
---|
2723 | do l=1,nlayer |
---|
2724 | mtot(ig) = mtot(ig) + |
---|
2725 | & zq(ig,l,igcm_h2o_vap) * |
---|
2726 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2727 | icetot(ig) = icetot(ig) + |
---|
2728 | & zq(ig,l,igcm_h2o_ice) * |
---|
2729 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2730 | IF (hdo) then |
---|
2731 | mtotD(ig) = mtotD(ig) + |
---|
2732 | & zq(ig,l,igcm_hdo_vap) * |
---|
2733 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2734 | icetotD(ig) = icetotD(ig) + |
---|
2735 | & zq(ig,l,igcm_hdo_ice) * |
---|
2736 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2737 | ENDIF !hdo |
---|
2738 | |
---|
2739 | c Computing abs optical depth at 825 cm-1 in each |
---|
2740 | c layer to simulate NEW TES retrieval |
---|
2741 | Qabsice = min( |
---|
2742 | & max(0.4e6*rice(ig,l)*(1.+nuice_ref)-0.05 ,0.),1.2 |
---|
2743 | & ) |
---|
2744 | opTES(ig,l)= 0.75 * Qabsice * |
---|
2745 | & zq(ig,l,igcm_h2o_ice) * |
---|
2746 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2747 | & / (rho_ice * rice(ig,l) * (1.+nuice_ref)) |
---|
2748 | tauTES(ig)=tauTES(ig)+ opTES(ig,l) |
---|
2749 | enddo |
---|
2750 | c rave(ig)=rave(ig)/max(icetot(ig),1.e-30) ! mass weight |
---|
2751 | c if (icetot(ig)*1e3.lt.0.01) rave(ig)=0. |
---|
2752 | enddo |
---|
2753 | call watersat(ngrid*nlayer,zt,zplay,zqsat) |
---|
2754 | satu(:,:) = zq(:,:,igcm_h2o_vap)/zqsat(:,:) |
---|
2755 | |
---|
2756 | if (scavenging) then |
---|
2757 | Nccntot(:)= 0 |
---|
2758 | Mccntot(:)= 0 |
---|
2759 | rave(:)=0 |
---|
2760 | do ig=1,ngrid |
---|
2761 | do l=1,nlayer |
---|
2762 | Nccntot(ig) = Nccntot(ig) + |
---|
2763 | & zq(ig,l,igcm_ccn_number)*tauscaling(ig) |
---|
2764 | & *(zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2765 | Mccntot(ig) = Mccntot(ig) + |
---|
2766 | & zq(ig,l,igcm_ccn_mass)*tauscaling(ig) |
---|
2767 | & *(zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2768 | cccc Column integrated effective ice radius |
---|
2769 | cccc is weighted by total ice surface area (BETTER than total ice mass) |
---|
2770 | rave(ig) = rave(ig) + |
---|
2771 | & tauscaling(ig) * |
---|
2772 | & zq(ig,l,igcm_ccn_number) * |
---|
2773 | & (zplev(ig,l) - zplev(ig,l+1)) / g * |
---|
2774 | & rice(ig,l) * rice(ig,l)* (1.+nuice_ref) |
---|
2775 | enddo |
---|
2776 | rave(ig)=(icetot(ig)/rho_ice+Mccntot(ig)/rho_dust)*0.75 |
---|
2777 | & /max(pi*rave(ig),1.e-30) ! surface weight |
---|
2778 | if (icetot(ig)*1e3.lt.0.01) rave(ig)=0. |
---|
2779 | enddo |
---|
2780 | else ! of if (scavenging) |
---|
2781 | rave(:)=0 |
---|
2782 | do ig=1,ngrid |
---|
2783 | do l=1,nlayer |
---|
2784 | rave(ig) = rave(ig) + |
---|
2785 | & zq(ig,l,igcm_h2o_ice) * |
---|
2786 | & (zplev(ig,l) - zplev(ig,l+1)) / g * |
---|
2787 | & rice(ig,l) * (1.+nuice_ref) |
---|
2788 | enddo |
---|
2789 | rave(ig) = max(rave(ig) / |
---|
2790 | & max(icetot(ig),1.e-30),1.e-30) ! mass weight |
---|
2791 | enddo |
---|
2792 | endif ! of if (scavenging) |
---|
2793 | |
---|
2794 | !Alternative A. Pottier weighting |
---|
2795 | rave2(:) = 0. |
---|
2796 | totrave2(:) = 0. |
---|
2797 | do ig=1,ngrid |
---|
2798 | do l=1,nlayer |
---|
2799 | rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l) |
---|
2800 | totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice) |
---|
2801 | end do |
---|
2802 | rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30) |
---|
2803 | end do |
---|
2804 | |
---|
2805 | endif ! of if (water) |
---|
2806 | |
---|
2807 | if (co2clouds) then |
---|
2808 | mtotco2(1:ngrid) = 0. |
---|
2809 | icetotco2(1:ngrid) = 0. |
---|
2810 | vaptotco2(1:ngrid) = 0. |
---|
2811 | do ig=1,ngrid |
---|
2812 | do l=1,nlayer |
---|
2813 | vaptotco2(ig) = vaptotco2(ig) + |
---|
2814 | & zq(ig,l,igcm_co2) * |
---|
2815 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2816 | icetotco2(ig) = icetot(ig) + |
---|
2817 | & zq(ig,l,igcm_co2_ice) * |
---|
2818 | & (zplev(ig,l) - zplev(ig,l+1)) / g |
---|
2819 | end do |
---|
2820 | mtotco2(ig) = icetotco2(ig) + vaptotco2(ig) |
---|
2821 | end do |
---|
2822 | end if |
---|
2823 | |
---|
2824 | #ifndef MESOSCALE |
---|
2825 | c ----------------------------------------------------------------- |
---|
2826 | c WSTATS: Saving statistics |
---|
2827 | c ----------------------------------------------------------------- |
---|
2828 | c ("stats" stores and accumulates key variables in file "stats.nc" |
---|
2829 | c which can later be used to make the statistic files of the run: |
---|
2830 | c if flag "callstats" from callphys.def is .true.) |
---|
2831 | |
---|
2832 | call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) |
---|
2833 | call wstats(ngrid,"tsurf","Surface temperature","K",2 |
---|
2834 | & ,tsurf(:,iflat)) |
---|
2835 | call wstats(ngrid,"co2ice","CO2 ice cover", |
---|
2836 | & "kg.m-2",2,qsurf(:,igcm_co2,iflat)) |
---|
2837 | call wstats(ngrid,"watercap","H2O ice cover", |
---|
2838 | & "kg.m-2",2,watercap(:,iflat)) |
---|
2839 | call wstats(ngrid,"tau_pref_scenario", |
---|
2840 | & "prescribed visible dod at 610 Pa","NU", |
---|
2841 | & 2,tau_pref_scenario) |
---|
2842 | call wstats(ngrid,"tau_pref_gcm", |
---|
2843 | & "visible dod at 610 Pa in the GCM","NU", |
---|
2844 | & 2,tau_pref_gcm) |
---|
2845 | call wstats(ngrid,"fluxsurf_lw", |
---|
2846 | & "Thermal IR radiative flux to surface","W.m-2",2, |
---|
2847 | & fluxsurf_lw(:,iflat)) |
---|
2848 | call wstats(ngrid,"fluxsurf_dn_sw", |
---|
2849 | & "Incoming Solar radiative flux to surface","W.m-2",2, |
---|
2850 | & fluxsurf_dn_sw_tot(:,iflat)) |
---|
2851 | call wstats(ngrid,"fluxsurf_up_sw", |
---|
2852 | & "Reflected Solar radiative flux from surface","W.m-2",2, |
---|
2853 | & fluxsurf_up_sw_tot) |
---|
2854 | call wstats(ngrid,"fluxtop_lw", |
---|
2855 | & "Thermal IR radiative flux to space","W.m-2",2, |
---|
2856 | & fluxtop_lw) |
---|
2857 | call wstats(ngrid,"fluxtop_dn_sw", |
---|
2858 | & "Incoming Solar radiative flux from space","W.m-2",2, |
---|
2859 | & fluxtop_dn_sw_tot) |
---|
2860 | call wstats(ngrid,"fluxtop_up_sw", |
---|
2861 | & "Outgoing Solar radiative flux to space","W.m-2",2, |
---|
2862 | & fluxtop_up_sw_tot) |
---|
2863 | call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) |
---|
2864 | call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) |
---|
2865 | call wstats(ngrid,"v","Meridional (North-South) wind", |
---|
2866 | & "m.s-1",3,zv) |
---|
2867 | call wstats(ngrid,"w","Vertical (down-up) wind", |
---|
2868 | & "m.s-1",3,pw) |
---|
2869 | call wstats(ngrid,"rho","Atmospheric density","kg/m3",3,rho) |
---|
2870 | call wstats(ngrid,"pressure","Pressure","Pa",3,zplay) |
---|
2871 | call wstats(ngrid,"q2", |
---|
2872 | & "Boundary layer eddy kinetic energy", |
---|
2873 | & "m2.s-2",3,q2) |
---|
2874 | call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, |
---|
2875 | & emis(:,iflat)) |
---|
2876 | call wstats(ngrid,"fluxsurf_dir_dn_sw", |
---|
2877 | & "Direct incoming SW flux at surface", |
---|
2878 | & "W.m-2",2,fluxsurf_dir_dn_sw) |
---|
2879 | |
---|
2880 | if (calltherm) then |
---|
2881 | call wstats(ngrid,"zmax_th","Height of thermals", |
---|
2882 | & "m",2,zmax_th) |
---|
2883 | call wstats(ngrid,"hfmax_th","Max thermals heat flux", |
---|
2884 | & "K.m/s",2,hfmax_th) |
---|
2885 | call wstats(ngrid,"wstar", |
---|
2886 | & "Max vertical velocity in thermals", |
---|
2887 | & "m/s",2,wstar) |
---|
2888 | endif |
---|
2889 | |
---|
2890 | if (water) then |
---|
2891 | vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap) |
---|
2892 | & *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap) |
---|
2893 | call wstats(ngrid,"vmr_h2ovap", |
---|
2894 | & "H2O vapor volume mixing ratio","mol/mol", |
---|
2895 | & 3,vmr) |
---|
2896 | vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
2897 | & *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice) |
---|
2898 | call wstats(ngrid,"vmr_h2oice", |
---|
2899 | & "H2O ice volume mixing ratio","mol/mol", |
---|
2900 | & 3,vmr) |
---|
2901 | ! also store vmr_ice*rice for better diagnostics of rice |
---|
2902 | vmr(1:ngrid,1:nlayer)=vmr(1:ngrid,1:nlayer)* |
---|
2903 | & rice(1:ngrid,1:nlayer) |
---|
2904 | call wstats(ngrid,"vmr_h2oice_rice", |
---|
2905 | & "H2O ice mixing ratio times ice particule size", |
---|
2906 | & "(mol/mol)*m", |
---|
2907 | & 3,vmr) |
---|
2908 | vmr=zqsat(1:ngrid,1:nlayer) |
---|
2909 | & *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap) |
---|
2910 | call wstats(ngrid,"vmr_h2osat", |
---|
2911 | & "saturation volume mixing ratio","mol/mol", |
---|
2912 | & 3,vmr) |
---|
2913 | call wstats(ngrid,"h2o_ice_s", |
---|
2914 | & "surface h2o_ice","kg/m2", |
---|
2915 | & 2,qsurf(1,igcm_h2o_ice,iflat)) |
---|
2916 | call wstats(ngrid,'albedo', |
---|
2917 | & 'albedo', |
---|
2918 | & '',2,albedo(1,1,iflat)) |
---|
2919 | call wstats(ngrid,"mtot", |
---|
2920 | & "total mass of water vapor","kg/m2", |
---|
2921 | & 2,mtot) |
---|
2922 | call wstats(ngrid,"icetot", |
---|
2923 | & "total mass of water ice","kg/m2", |
---|
2924 | & 2,icetot) |
---|
2925 | call wstats(ngrid,"reffice", |
---|
2926 | & "Mean reff","m", |
---|
2927 | & 2,rave) |
---|
2928 | call wstats(ngrid,"Nccntot", |
---|
2929 | & "condensation nuclei","Nbr/m2", |
---|
2930 | & 2,Nccntot) |
---|
2931 | call wstats(ngrid,"Mccntot", |
---|
2932 | & "condensation nuclei mass","kg/m2", |
---|
2933 | & 2,Mccntot) |
---|
2934 | call wstats(ngrid,"rice", |
---|
2935 | & "Ice particle size","m", |
---|
2936 | & 3,rice) |
---|
2937 | if (.not.activice) then |
---|
2938 | call wstats(ngrid,"tauTESap", |
---|
2939 | & "tau abs 825 cm-1","", |
---|
2940 | & 2,tauTES) |
---|
2941 | else |
---|
2942 | call wstats(ngrid,'tauTES', |
---|
2943 | & 'tau abs 825 cm-1', |
---|
2944 | & '',2,taucloudtes) |
---|
2945 | endif |
---|
2946 | |
---|
2947 | endif ! of if (water) |
---|
2948 | |
---|
2949 | if (co2clouds) then |
---|
2950 | call wstats(ngrid,"mtotco2", |
---|
2951 | & "total mass atm of co2","kg/m2", |
---|
2952 | & 2,mtotco2) |
---|
2953 | call wstats(ngrid,"icetotco2", |
---|
2954 | & "total mass atm of co2 ice","kg/m2", |
---|
2955 | & 2,icetotco2) |
---|
2956 | call wstats(ngrid,"vaptotco2", |
---|
2957 | & "total mass atm of co2 vapor","kg/m2", |
---|
2958 | & 2,icetotco2) |
---|
2959 | end if |
---|
2960 | |
---|
2961 | |
---|
2962 | if (dustbin.ne.0) then |
---|
2963 | |
---|
2964 | call wstats(ngrid,'tau','taudust','SI',2,tau(1,1)) |
---|
2965 | |
---|
2966 | if (doubleq) then |
---|
2967 | call wstats(ngrid,'dqsdust', |
---|
2968 | & 'deposited surface dust mass', |
---|
2969 | & 'kg.m-2.s-1',2,dqdustsurf) |
---|
2970 | call wstats(ngrid,'dqndust', |
---|
2971 | & 'deposited surface dust number', |
---|
2972 | & 'number.m-2.s-1',2,dndustsurf) |
---|
2973 | call wstats(ngrid,'reffdust','reffdust', |
---|
2974 | & 'm',3,rdust*ref_r0) |
---|
2975 | call wstats(ngrid,'dustq','Dust mass mr', |
---|
2976 | & 'kg/kg',3,qdust) |
---|
2977 | call wstats(ngrid,'dustN','Dust number', |
---|
2978 | & 'part/kg',3,ndust) |
---|
2979 | if (rdstorm) then |
---|
2980 | call wstats(ngrid,'reffstormdust','reffdust', |
---|
2981 | & 'm',3,rstormdust*ref_r0) |
---|
2982 | call wstats(ngrid,'rdsdustq','Dust mass mr', |
---|
2983 | & 'kg/kg',3,rdsqdust) |
---|
2984 | call wstats(ngrid,'rdsdustN','Dust number', |
---|
2985 | & 'part/kg',3,rdsndust) |
---|
2986 | end if |
---|
2987 | else |
---|
2988 | do iq=1,dustbin |
---|
2989 | write(str2(1:2),'(i2.2)') iq |
---|
2990 | call wstats(ngrid,'q'//str2,'mix. ratio', |
---|
2991 | & 'kg/kg',3,zq(1,1,iq)) |
---|
2992 | call wstats(ngrid,'qsurf'//str2,'qsurf', |
---|
2993 | & 'kg.m-2',2,qsurf(1,iq,iflat)) |
---|
2994 | end do |
---|
2995 | endif ! (doubleq) |
---|
2996 | |
---|
2997 | if (scavenging) then |
---|
2998 | call wstats(ngrid,'ccnq','CCN mass mr', |
---|
2999 | & 'kg/kg',3,qccn) |
---|
3000 | call wstats(ngrid,'ccnN','CCN number', |
---|
3001 | & 'part/kg',3,nccn) |
---|
3002 | endif ! (scavenging) |
---|
3003 | |
---|
3004 | endif ! (dustbin.ne.0) |
---|
3005 | |
---|
3006 | if (photochem) then |
---|
3007 | do iq=1,nq |
---|
3008 | if (noms(iq) .ne. "dust_mass" .and. |
---|
3009 | $ noms(iq) .ne. "dust_number" .and. |
---|
3010 | $ noms(iq) .ne. "ccn_mass" .and. |
---|
3011 | $ noms(iq) .ne. "ccn_number" .and. |
---|
3012 | $ noms(iq) .ne. "ccnco2_mass" .and. |
---|
3013 | $ noms(iq) .ne. "ccnco2_number" .and. |
---|
3014 | $ noms(iq) .ne. "stormdust_mass" .and. |
---|
3015 | $ noms(iq) .ne. "stormdust_number" .and. |
---|
3016 | $ noms(iq) .ne. "topdust_mass" .and. |
---|
3017 | $ noms(iq) .ne. "topdust_number") then |
---|
3018 | ! volume mixing ratio |
---|
3019 | |
---|
3020 | vmr(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq) |
---|
3021 | & *mmean(1:ngrid,1:nlayer)/mmol(iq) |
---|
3022 | |
---|
3023 | call wstats(ngrid,"vmr_"//trim(noms(iq)), |
---|
3024 | $ "Volume mixing ratio","mol/mol",3,vmr) |
---|
3025 | if ((noms(iq).eq."o") |
---|
3026 | $ .or. (noms(iq).eq."co2") |
---|
3027 | $ .or. (noms(iq).eq."o3") |
---|
3028 | $ .or. (noms(iq).eq."ar") |
---|
3029 | $ .or. (noms(iq).eq."o2") |
---|
3030 | $ .or. (noms(iq).eq."h2o_vap") ) then |
---|
3031 | call write_output("vmr_"//trim(noms(iq)), |
---|
3032 | $ "Volume mixing ratio","mol/mol",vmr(:,:)) |
---|
3033 | end if |
---|
3034 | |
---|
3035 | ! number density (molecule.cm-3) |
---|
3036 | |
---|
3037 | rhopart(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,iq) |
---|
3038 | & *rho(1:ngrid,1:nlayer)*n_avog/ |
---|
3039 | & (1000*mmol(iq)) |
---|
3040 | |
---|
3041 | call wstats(ngrid,"num_"//trim(noms(iq)), |
---|
3042 | $ "Number density","cm-3",3,rhopart) |
---|
3043 | call write_output("num_"//trim(noms(iq)), |
---|
3044 | $ "Number density","cm-3",rhopart(:,:)) |
---|
3045 | |
---|
3046 | ! vertical column (molecule.cm-2) |
---|
3047 | |
---|
3048 | do ig = 1,ngrid |
---|
3049 | colden(ig,iq) = 0. |
---|
3050 | end do |
---|
3051 | do l=1,nlayer |
---|
3052 | do ig=1,ngrid |
---|
3053 | colden(ig,iq) = colden(ig,iq) + zq(ig,l,iq) |
---|
3054 | $ *(zplev(ig,l)-zplev(ig,l+1)) |
---|
3055 | $ *6.022e22/(mmol(iq)*g) |
---|
3056 | end do |
---|
3057 | end do |
---|
3058 | |
---|
3059 | call wstats(ngrid,"c_"//trim(noms(iq)), |
---|
3060 | $ "column","mol cm-2",2,colden(1,iq)) |
---|
3061 | call write_output("c_"//trim(noms(iq)), |
---|
3062 | $ "column","mol cm-2",colden(:,iq)) |
---|
3063 | |
---|
3064 | ! global mass (g) |
---|
3065 | |
---|
3066 | call planetwide_sumval(colden(:,iq)/6.022e23 |
---|
3067 | $ *mmol(iq)*1.e4*cell_area(:),mass(iq)) |
---|
3068 | |
---|
3069 | call write_output("mass_"//trim(noms(iq)), |
---|
3070 | $ "global mass","g",mass(iq)) |
---|
3071 | |
---|
3072 | end if ! of if (noms(iq) .ne. "dust_mass" ...) |
---|
3073 | end do ! of do iq=1,nq |
---|
3074 | end if ! of if (photochem) |
---|
3075 | |
---|
3076 | |
---|
3077 | IF(lastcall.and.callstats) THEN |
---|
3078 | write (*,*) "Writing stats..." |
---|
3079 | call mkstats(ierr) |
---|
3080 | ENDIF |
---|
3081 | |
---|
3082 | |
---|
3083 | c (Store EOF for Mars Climate database software) |
---|
3084 | IF (calleofdump) THEN |
---|
3085 | CALL eofdump(ngrid, nlayer, zu, zv, zt, rho, ps) |
---|
3086 | ENDIF |
---|
3087 | #endif |
---|
3088 | !endif of ifndef MESOSCALE |
---|
3089 | |
---|
3090 | #ifdef MESOSCALE |
---|
3091 | |
---|
3092 | !! see comm_wrf. |
---|
3093 | !! not needed when an array is already in a shared module. |
---|
3094 | !! --> example : hfmax_th, zmax_th |
---|
3095 | |
---|
3096 | CALL allocate_comm_wrf(ngrid,nlayer) |
---|
3097 | |
---|
3098 | !state real HR_SW ikj misc 1 - h "HR_SW" "HEATING RATE SW" "K/s" |
---|
3099 | comm_HR_SW(1:ngrid,1:nlayer) = zdtsw(1:ngrid,1:nlayer) |
---|
3100 | !state real HR_LW ikj misc 1 - h "HR_LW" "HEATING RATE LW" "K/s" |
---|
3101 | comm_HR_LW(1:ngrid,1:nlayer) = zdtlw(1:ngrid,1:nlayer) |
---|
3102 | !state real SWDOWNZ ij misc 1 - h "SWDOWNZ" "DOWNWARD SW FLUX AT SURFACE" "W m-2" |
---|
3103 | comm_SWDOWNZ(1:ngrid) = fluxsurf_dn_sw_tot(1:ngrid) |
---|
3104 | !state real TAU_DUST ij misc 1 - h "TAU_DUST" "REFERENCE VISIBLE DUST OPACITY" "" |
---|
3105 | comm_TAU_DUST(1:ngrid) = tau_pref_gcm(1:ngrid) |
---|
3106 | !state real RDUST ikj misc 1 - h "RDUST" "DUST RADIUS" "m" |
---|
3107 | comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer) |
---|
3108 | !state real QSURFDUST ij misc 1 - h "QSURFDUST" "DUST MASS AT SURFACE" "kg m-2" |
---|
3109 | IF (igcm_dust_mass .ne. 0) THEN |
---|
3110 | comm_QSURFDUST(1:ngrid) = qsurf(1:ngrid,igcm_dust_mass) |
---|
3111 | ELSE |
---|
3112 | comm_QSURFDUST(1:ngrid) = 0. |
---|
3113 | ENDIF |
---|
3114 | !state real MTOT ij misc 1 - h "MTOT" "TOTAL MASS WATER VAPOR in pmic" "pmic" |
---|
3115 | comm_MTOT(1:ngrid) = mtot(1:ngrid) * 1.e6 / rho_ice |
---|
3116 | !state real ICETOT ij misc 1 - h "ICETOT" "TOTAL MASS WATER ICE" "kg m-2" |
---|
3117 | comm_ICETOT(1:ngrid) = icetot(1:ngrid) * 1.e6 / rho_ice |
---|
3118 | !state real VMR_ICE ikj misc 1 - h "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm" |
---|
3119 | IF (igcm_h2o_ice .ne. 0) THEN |
---|
3120 | comm_VMR_ICE(1:ngrid,1:nlayer) = 1.e6 |
---|
3121 | . * zq(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
3122 | . * mmean(1:ngrid,1:nlayer) / mmol(igcm_h2o_ice) |
---|
3123 | ELSE |
---|
3124 | comm_VMR_ICE(1:ngrid,1:nlayer) = 0. |
---|
3125 | ENDIF |
---|
3126 | !state real TAU_ICE ij misc 1 - h "TAU_ICE" "CLOUD OD at 825 cm-1 TES" "" |
---|
3127 | if (activice) then |
---|
3128 | comm_TAU_ICE(1:ngrid) = taucloudtes(1:ngrid) |
---|
3129 | else |
---|
3130 | comm_TAU_ICE(1:ngrid) = tauTES(1:ngrid) |
---|
3131 | endif |
---|
3132 | !state real RICE ikj misc 1 - h "RICE" "ICE RADIUS" "m" |
---|
3133 | comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer) |
---|
3134 | |
---|
3135 | !! calculate sensible heat flux in W/m2 for outputs |
---|
3136 | !! -- the one computed in vdifc is not the real one |
---|
3137 | !! -- vdifc must have been called |
---|
3138 | if (.not.callrichsl) then |
---|
3139 | sensibFlux(1:ngrid) = zflubid(1:ngrid) |
---|
3140 | . - capcal(1:ngrid)*zdtsdif(1:ngrid) |
---|
3141 | else |
---|
3142 | sensibFlux(1:ngrid) = |
---|
3143 | & (pplay(1:ngrid,1)/(r*pt(1:ngrid,1)))*cpp |
---|
3144 | & *sqrt(pu(1:ngrid,1)*pu(1:ngrid,1)+pv(1:ngrid,1)*pv(1:ngrid,1) |
---|
3145 | & +(log(1.+0.7*wstar(1:ngrid) + 2.3*wstar(1:ngrid)**2))**2) |
---|
3146 | & *zcdh(1:ngrid)*(tsurf(1:ngrid)-zh(1:ngrid,1)) |
---|
3147 | endif |
---|
3148 | |
---|
3149 | #else |
---|
3150 | #ifndef MESOINI |
---|
3151 | |
---|
3152 | c ========================================================== |
---|
3153 | c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing |
---|
3154 | c any variable for diagnostic (output with period |
---|
3155 | c "ecritphy", set in "run.def") |
---|
3156 | c ========================================================== |
---|
3157 | c WRITEDIAGFI can ALSO be called from any other subroutines |
---|
3158 | c for any variables !! |
---|
3159 | call write_output("emis","Surface emissivity","", |
---|
3160 | & emis(:,iflat)) |
---|
3161 | do islope=1,nslope |
---|
3162 | write(str2(1:2),'(i2.2)') islope |
---|
3163 | call write_output("emis_slope"//str2, |
---|
3164 | & "Surface emissivity","",emis(:,islope)) |
---|
3165 | ENDDO |
---|
3166 | call write_output("zzlay","Midlayer altitude", |
---|
3167 | & "m",zzlay(:,:)) |
---|
3168 | call write_output("zzlev","Interlayer altitude", |
---|
3169 | & "m",zzlev(:,1:nlayer)) |
---|
3170 | call write_output("pphi","Geopotential","m2s-2", |
---|
3171 | & pphi(:,:)) |
---|
3172 | call write_output("phisfi","Surface geopotential", |
---|
3173 | & "m2s-2",phisfi(:)) |
---|
3174 | call write_output("tsurf","Surface temperature","K", |
---|
3175 | & tsurf(:,iflat)) |
---|
3176 | do islope=1,nslope |
---|
3177 | write(str2(1:2),'(i2.2)') islope |
---|
3178 | call write_output("tsurf_slope"//str2, |
---|
3179 | & "Surface temperature","K", |
---|
3180 | & tsurf(:,islope)) |
---|
3181 | ENDDO |
---|
3182 | call write_output("ps","surface pressure","Pa",ps(:)) |
---|
3183 | call write_output("co2ice","co2 ice thickness" |
---|
3184 | & ,"kg.m-2",qsurf(:,igcm_co2,iflat)) |
---|
3185 | do islope=1,nslope |
---|
3186 | write(str2(1:2),'(i2.2)') islope |
---|
3187 | call write_output("co2ice_slope"//str2,"co2 ice thickness" |
---|
3188 | & ,"kg.m-2",qsurf(:,igcm_co2,islope)) |
---|
3189 | ENDDO |
---|
3190 | call write_output("watercap","Perennial water ice thickness" |
---|
3191 | & ,"kg.m-2",watercap(:,iflat)) |
---|
3192 | do islope=1,nslope |
---|
3193 | write(str2(1:2),'(i2.2)') islope |
---|
3194 | call write_output("watercap_slope"//str2, |
---|
3195 | & "Perennial water ice thickness" |
---|
3196 | & ,"kg.m-2",watercap(:,islope)) |
---|
3197 | ENDDO |
---|
3198 | call write_output("temp_layer1","temperature in layer 1", |
---|
3199 | & "K",zt(:,1)) |
---|
3200 | call write_output("temp7","temperature in layer 7", |
---|
3201 | & "K",zt(:,7)) |
---|
3202 | call write_output("fluxsurf_lw","fluxsurf_lw","W.m-2", |
---|
3203 | & fluxsurf_lw(:,iflat)) |
---|
3204 | do islope=1,nslope |
---|
3205 | write(str2(1:2),'(i2.2)') islope |
---|
3206 | call write_output("fluxsurf_lw_slope"//str2, |
---|
3207 | & "fluxsurf_lw","W.m-2", |
---|
3208 | & fluxsurf_lw(:,islope)) |
---|
3209 | ENDDO |
---|
3210 | call write_output("fluxsurf_dn_sw","fluxsurf_dn_sw", |
---|
3211 | & "W.m-2",fluxsurf_dn_sw_tot(:,iflat)) |
---|
3212 | do islope=1,nslope |
---|
3213 | write(str2(1:2),'(i2.2)') islope |
---|
3214 | call write_output("fluxsurf_dn_sw_slope"//str2, |
---|
3215 | & "fluxsurf_dn_sw", |
---|
3216 | & "W.m-2",fluxsurf_dn_sw_tot(:,islope)) |
---|
3217 | ENDDO |
---|
3218 | call write_output("fluxtop_lw","fluxtop_lw","W.m-2", |
---|
3219 | & fluxtop_lw(:)) |
---|
3220 | call write_output("fluxtop_up_sw","fluxtop_up_sw","W.m-2", |
---|
3221 | & fluxtop_up_sw_tot(:)) |
---|
3222 | call write_output("temp","temperature","K",zt(:,:)) |
---|
3223 | call write_output("Sols","Time","sols",zday) |
---|
3224 | call write_output("Ls","Solar longitude","deg", |
---|
3225 | & zls*180./pi) |
---|
3226 | call write_output("u","Zonal wind","m.s-1",zu(:,:)) |
---|
3227 | call write_output("v","Meridional wind","m.s-1",zv(:,:)) |
---|
3228 | call write_output("w","Vertical wind","m.s-1",pw(:,:)) |
---|
3229 | call write_output("rho","density","kg.m-3",rho(:,:)) |
---|
3230 | call write_output("pressure","Pressure","Pa",zplay(:,:)) |
---|
3231 | call write_output("zplev","Interlayer pressure","Pa", |
---|
3232 | & zplev(:,1:nlayer)) |
---|
3233 | call write_output('sw_htrt','sw heat. rate', |
---|
3234 | & 'K/s',zdtsw(:,:)) |
---|
3235 | call write_output('lw_htrt','lw heat. rate', |
---|
3236 | & 'K/s',zdtlw(:,:)) |
---|
3237 | call write_output("local_time","Local time", |
---|
3238 | & 'sol',local_time(:)) |
---|
3239 | if (.not.activice) then |
---|
3240 | CALL write_output('tauTESap', |
---|
3241 | & 'tau abs 825 cm-1', |
---|
3242 | & '',tauTES(:)) |
---|
3243 | else |
---|
3244 | CALL write_output('tauTES', |
---|
3245 | & 'tau abs 825 cm-1', |
---|
3246 | & '',taucloudtes(:)) |
---|
3247 | endif |
---|
3248 | #else |
---|
3249 | !!! this is to ensure correct initialisation of mesoscale model |
---|
3250 | call write_output("tsurf","Surface temperature","K", |
---|
3251 | & tsurf(:,iflat)) |
---|
3252 | call write_output("ps","surface pressure","Pa",ps(:)) |
---|
3253 | call write_output("co2ice","co2 ice thickness","kg.m-2", |
---|
3254 | & qsurf(:,igcm_co2,iflat)) |
---|
3255 | call write_output("temp","temperature","K",zt(:,:)) |
---|
3256 | call write_output("u","Zonal wind","m.s-1",zu(:,:)) |
---|
3257 | call write_output("v","Meridional wind","m.s-1",zv(:,:)) |
---|
3258 | call write_output("emis","Surface emissivity","", |
---|
3259 | & emis(:,iflat)) |
---|
3260 | call write_output("tsoil","Soil temperature", |
---|
3261 | & "K",tsoil(:,:,iflat)) |
---|
3262 | call write_output("inertiedat","Soil inertia", |
---|
3263 | & "K",inertiedat(:,:)) |
---|
3264 | #endif |
---|
3265 | |
---|
3266 | c ---------------------------------------------------------- |
---|
3267 | c Outputs of the CO2 cycle |
---|
3268 | c ---------------------------------------------------------- |
---|
3269 | |
---|
3270 | if (igcm_co2.ne.0) then |
---|
3271 | call write_output("co2","co2 mass mixing ratio", |
---|
3272 | & "kg.kg-1",zq(:,:,igcm_co2)) |
---|
3273 | |
---|
3274 | if (co2clouds) then |
---|
3275 | call write_output('ccnqco2','CCNco2 mmr', |
---|
3276 | & 'kg.kg-1',zq(:,:,igcm_ccnco2_mass)) |
---|
3277 | |
---|
3278 | call write_output('ccnNco2','CCNco2 number', |
---|
3279 | & 'part.kg-1',zq(:,:,igcm_ccnco2_number)) |
---|
3280 | |
---|
3281 | call write_output('co2_ice','co2_ice mmr in atm', |
---|
3282 | & 'kg.kg-1',zq(:,:,igcm_co2_ice)) |
---|
3283 | |
---|
3284 | call write_output("mtotco2","total mass atm of co2", |
---|
3285 | & "kg.m-2",mtotco2(:)) |
---|
3286 | call write_output("icetotco2","total mass atm of co2 ice", |
---|
3287 | & "kg.m-2", icetotco2(:)) |
---|
3288 | call write_output("vaptotco2","total mass atm of co2 "// |
---|
3289 | & "vapor","kg.m-2", vaptotco2(:)) |
---|
3290 | if (co2useh2o) then |
---|
3291 | call write_output('ccnqco2_h2o_m_ice', |
---|
3292 | & 'CCNco2_h2o_mass_ice mmr', |
---|
3293 | & 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ice)) |
---|
3294 | |
---|
3295 | call write_output('ccnqco2_h2o_m_ccn', |
---|
3296 | & 'CCNco2_h2o_mass_ccn mmr', |
---|
3297 | & 'kg.kg-1',zq(:,:,igcm_ccnco2_h2o_mass_ccn)) |
---|
3298 | |
---|
3299 | call write_output('ccnNco2_h2o','CCNco2_h2o number', |
---|
3300 | & 'part.kg-1',zq(:,:,igcm_ccnco2_h2o_number)) |
---|
3301 | end if |
---|
3302 | |
---|
3303 | if (meteo_flux) then |
---|
3304 | call write_output('ccnqco2_meteor','CCNco2_meteor mmr', |
---|
3305 | & 'kg.kg-1',zq(:,:,igcm_ccnco2_meteor_mass)) |
---|
3306 | |
---|
3307 | call write_output('ccnNco2_meteor','CCNco2_meteor number', |
---|
3308 | & 'part.kg-1',zq(:,:,igcm_ccnco2_meteor_number)) |
---|
3309 | end if |
---|
3310 | |
---|
3311 | end if ! of if (co2clouds) |
---|
3312 | end if ! of if (igcm_co2.ne.0) |
---|
3313 | |
---|
3314 | ! Output He tracer, if there is one |
---|
3315 | if (igcm_he.ne.0) then |
---|
3316 | call write_output("he","helium mass mixing ratio", |
---|
3317 | & "kg/kg",zq(:,:,igcm_he)) |
---|
3318 | vmr = zq(1:ngrid,1:nlayer,igcm_he) |
---|
3319 | & * mmean(1:ngrid,1:nlayer)/mmol(igcm_he) |
---|
3320 | call write_output('vmr_he','helium vol. mixing ratio', |
---|
3321 | & 'mol/mol',vmr(:,:)) |
---|
3322 | end if |
---|
3323 | |
---|
3324 | c ---------------------------------------------------------- |
---|
3325 | c Outputs of the water cycle |
---|
3326 | c ---------------------------------------------------------- |
---|
3327 | if (water) then |
---|
3328 | #ifdef MESOINI |
---|
3329 | !!!! waterice = q01, voir readmeteo.F90 |
---|
3330 | call write_output('q01',noms(igcm_h2o_ice), |
---|
3331 | & 'kg/kg', |
---|
3332 | & zq(:,:,igcm_h2o_ice)) |
---|
3333 | !!!! watervapor = q02, voir readmeteo.F90 |
---|
3334 | call write_output('q02',noms(igcm_h2o_vap), |
---|
3335 | & 'kg/kg', |
---|
3336 | & zq(:,:,igcm_h2o_vap)) |
---|
3337 | !!!! surface waterice qsurf02 (voir readmeteo) |
---|
3338 | call write_output('qsurf02','surface tracer', |
---|
3339 | & 'kg.m-2', |
---|
3340 | & qsurf(:,igcm_h2o_ice,iflat)) |
---|
3341 | do islope=1,nslope |
---|
3342 | write(str2(1:2),'(i2.2)') islope |
---|
3343 | call write_output('qsurf02_slope'//str2, |
---|
3344 | & 'surface tracer','kg.m-2', |
---|
3345 | & qsurf(:,igcm_h2o_ice,islope)) |
---|
3346 | ENDDO |
---|
3347 | #endif |
---|
3348 | call write_output('mtot', |
---|
3349 | & 'total mass of water vapor', |
---|
3350 | & 'kg/m2',mtot(:)) |
---|
3351 | call write_output('icetot', |
---|
3352 | & 'total mass of water ice', |
---|
3353 | & 'kg/m2',icetot(:)) |
---|
3354 | vmr = zq(1:ngrid,1:nlayer,igcm_h2o_ice) |
---|
3355 | & * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice) |
---|
3356 | call write_output('vmr_h2oice','h2o ice vmr', |
---|
3357 | & 'mol/mol',vmr(:,:)) |
---|
3358 | vmr = zq(1:ngrid,1:nlayer,igcm_h2o_vap) |
---|
3359 | & * mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap) |
---|
3360 | call write_output('vmr_h2ovap','h2o vap vmr', |
---|
3361 | & 'mol/mol',vmr(:,:)) |
---|
3362 | call write_output('reffice', |
---|
3363 | & 'Mean reff', |
---|
3364 | & 'm',rave(:)) |
---|
3365 | call write_output('h2o_ice','h2o_ice','kg/kg', |
---|
3366 | & zq(:,:,igcm_h2o_ice)) |
---|
3367 | call write_output('h2o_vap','h2o_vap','kg/kg', |
---|
3368 | & zq(:,:,igcm_h2o_vap)) |
---|
3369 | |
---|
3370 | if (hdo) then |
---|
3371 | vmr=zq(1:ngrid,1:nlayer,igcm_hdo_ice) |
---|
3372 | & *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_ice) |
---|
3373 | CALL write_output('vmr_hdoice','hdo ice vmr', |
---|
3374 | & 'mol/mol',vmr(:,:)) |
---|
3375 | vmr=zq(1:ngrid,1:nlayer,igcm_hdo_vap) |
---|
3376 | & *mmean(1:ngrid,1:nlayer)/mmol(igcm_hdo_vap) |
---|
3377 | CALL write_output('vmr_hdovap','hdo vap vmr', |
---|
3378 | & 'mol/mol',vmr(:,:)) |
---|
3379 | call write_output('hdo_ice','hdo_ice','kg/kg', |
---|
3380 | & zq(:,:,igcm_hdo_ice)) |
---|
3381 | call write_output('hdo_vap','hdo_vap','kg/kg', |
---|
3382 | & zq(:,:,igcm_hdo_vap)) |
---|
3383 | |
---|
3384 | CALL write_output('mtotD', |
---|
3385 | & 'total mass of HDO vapor', |
---|
3386 | & 'kg/m2',mtotD(:)) |
---|
3387 | CALL write_output('icetotD', |
---|
3388 | & 'total mass of HDO ice', |
---|
3389 | & 'kg/m2',icetotD(:)) |
---|
3390 | |
---|
3391 | C Calculation of the D/H ratio |
---|
3392 | do l=1,nlayer |
---|
3393 | do ig=1,ngrid |
---|
3394 | if (zq(ig,l,igcm_h2o_vap).gt.qperemin) then |
---|
3395 | DoH_vap(ig,l) = ( zq(ig,l,igcm_hdo_vap)/ |
---|
3396 | & zq(ig,l,igcm_h2o_vap) )*1./(2.*155.76e-6) |
---|
3397 | else |
---|
3398 | DoH_vap(ig,l) = 0. |
---|
3399 | endif |
---|
3400 | enddo |
---|
3401 | enddo |
---|
3402 | |
---|
3403 | do l=1,nlayer |
---|
3404 | do ig=1,ngrid |
---|
3405 | if (zq(ig,l,igcm_h2o_ice).gt.qperemin) then |
---|
3406 | DoH_ice(ig,l) = ( zq(ig,l,igcm_hdo_ice)/ |
---|
3407 | & zq(ig,l,igcm_h2o_ice) )/(2.*155.76e-6) |
---|
3408 | else |
---|
3409 | DoH_ice(ig,l) = 0. |
---|
3410 | endif |
---|
3411 | enddo |
---|
3412 | enddo |
---|
3413 | |
---|
3414 | CALL write_output('DoH_vap', |
---|
3415 | & 'D/H ratio in vapor', |
---|
3416 | & ' ',DoH_vap(:,:)) |
---|
3417 | CALL write_output('DoH_ice', |
---|
3418 | & 'D/H ratio in ice', |
---|
3419 | & '',DoH_ice(:,:)) |
---|
3420 | |
---|
3421 | endif !hdo |
---|
3422 | |
---|
3423 | !A. Pottier |
---|
3424 | ! CALL write_output('rmoym', |
---|
3425 | ! & 'alternative reffice', |
---|
3426 | ! & 'm',rave2(:)) |
---|
3427 | call write_output('h2o_saturation', |
---|
3428 | & 'h2o vap saturation ratio','',satu(:,:)) |
---|
3429 | if (scavenging) then |
---|
3430 | CALL write_output("Nccntot", |
---|
3431 | & "condensation nuclei","Nbr/m2", |
---|
3432 | & Nccntot(:)) |
---|
3433 | CALL write_output("Mccntot", |
---|
3434 | & "mass condensation nuclei","kg/m2", |
---|
3435 | & Mccntot(:)) |
---|
3436 | endif |
---|
3437 | call write_output('rice','Water ice particle size', |
---|
3438 | & 'm',rice(:,:)) |
---|
3439 | call write_output('h2o_ice_s', |
---|
3440 | & 'surface h2o_ice', |
---|
3441 | & 'kg.m-2',qsurf(:,igcm_h2o_ice,iflat)) |
---|
3442 | do islope=1,nslope |
---|
3443 | write(str2(1:2),'(i2.2)') islope |
---|
3444 | call write_output('h2o_ice_s_slope'//str2, |
---|
3445 | & 'surface h2o_ice', |
---|
3446 | & 'kg.m-2',qsurf(:,igcm_h2o_ice,islope)) |
---|
3447 | ENDDO |
---|
3448 | if (hdo) then |
---|
3449 | call write_output('hdo_ice_s', |
---|
3450 | & 'surface hdo_ice', |
---|
3451 | & 'kg.m-2',qsurf(:,igcm_hdo_ice,iflat)) |
---|
3452 | do islope=1,nslope |
---|
3453 | write(str2(1:2),'(i2.2)') islope |
---|
3454 | call write_output('hdo_ice_s_slope'//str2, |
---|
3455 | & 'surface hdo_ice', |
---|
3456 | & 'kg.m-2',qsurf(:,igcm_hdo_ice,islope)) |
---|
3457 | ENDDO |
---|
3458 | |
---|
3459 | do ig=1,ngrid |
---|
3460 | if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then |
---|
3461 | DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/ |
---|
3462 | & qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6 |
---|
3463 | else |
---|
3464 | DoH_surf(ig) = 0. |
---|
3465 | endif |
---|
3466 | enddo |
---|
3467 | |
---|
3468 | call write_output('DoH_surf', |
---|
3469 | & 'surface D/H', |
---|
3470 | & '',DoH_surf(:)) |
---|
3471 | endif ! hdo |
---|
3472 | |
---|
3473 | CALL write_output('albedo', |
---|
3474 | & 'albedo', |
---|
3475 | & '',albedo(:,1,iflat)) |
---|
3476 | do islope=1,nslope |
---|
3477 | write(str2(1:2),'(i2.2)') islope |
---|
3478 | CALL write_output('albedo_slope'//str2, |
---|
3479 | & 'albedo', |
---|
3480 | & '',albedo(:,1,islope)) |
---|
3481 | ENDDO |
---|
3482 | if (tifeedback) then |
---|
3483 | call write_output("soiltemp", |
---|
3484 | & "Soil temperature","K", |
---|
3485 | & tsoil(:,:,iflat)) |
---|
3486 | call write_output('soilti', |
---|
3487 | & 'Soil Thermal Inertia', |
---|
3488 | & 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat)) |
---|
3489 | do islope=1,nslope |
---|
3490 | write(str2(1:2),'(i2.2)') islope |
---|
3491 | call write_output('soilti_slope'//str2, |
---|
3492 | & 'Soil Thermal Inertia', |
---|
3493 | & 'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,islope)) |
---|
3494 | ENDDO |
---|
3495 | endif |
---|
3496 | !A. Pottier |
---|
3497 | if (CLFvarying) then !AP14 nebulosity |
---|
3498 | call write_output('totcloudfrac', |
---|
3499 | & 'Total cloud fraction', |
---|
3500 | & ' ',totcloudfrac(:)) |
---|
3501 | end if !clf varying |
---|
3502 | end if !(water) |
---|
3503 | |
---|
3504 | c ---------------------------------------------------------- |
---|
3505 | c Outputs of the dust cycle |
---|
3506 | c ---------------------------------------------------------- |
---|
3507 | |
---|
3508 | call write_output('tau_pref_scenario', |
---|
3509 | & 'Prescribed visible dust optical depth at 610Pa', |
---|
3510 | & 'NU',tau_pref_scenario(:)) |
---|
3511 | |
---|
3512 | call write_output('tau_pref_gcm', |
---|
3513 | & 'Visible dust optical depth at 610Pa in the GCM', |
---|
3514 | & 'NU',tau_pref_gcm(:)) |
---|
3515 | |
---|
3516 | if (reff_driven_IRtoVIS_scenario) then |
---|
3517 | call write_output('IRtoVIScoef', |
---|
3518 | & 'Conversion coeff for dust tau from abs9.3um to ext0.67um', |
---|
3519 | & '/',IRtoVIScoef(:)) |
---|
3520 | endif |
---|
3521 | |
---|
3522 | if (dustbin.ne.0) then |
---|
3523 | |
---|
3524 | #ifndef MESOINI |
---|
3525 | if (doubleq) then |
---|
3526 | call write_output('dqsdust', |
---|
3527 | & 'deposited surface dust mass', |
---|
3528 | & 'kg.m-2.s-1',dqdustsurf(:)) |
---|
3529 | call write_output('dqndust', |
---|
3530 | & 'deposited surface dust number', |
---|
3531 | & 'number.m-2.s-1',dndustsurf(:)) |
---|
3532 | call write_output('reffdust','reffdust', |
---|
3533 | & 'm',rdust(:,:)*ref_r0) |
---|
3534 | call write_output('dustq','Dust mass mr', |
---|
3535 | & 'kg/kg',qdust(:,:)) |
---|
3536 | call write_output('dustN','Dust number', |
---|
3537 | & 'part/kg',ndust(:,:)) |
---|
3538 | |
---|
3539 | select case (trim(dustiropacity)) |
---|
3540 | case ("tes") |
---|
3541 | call write_output('dsodust_TES', |
---|
3542 | & 'density scaled extinction opacity of std dust at 9.3um(TES)', |
---|
3543 | & 'm2.kg-1',dsodust(:,:)) |
---|
3544 | call write_output('dso_TES', |
---|
3545 | & 'density scaled extinction opacity of all dust at 9.3um(TES)', |
---|
3546 | & 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:)) |
---|
3547 | case ("mcs") |
---|
3548 | call write_output('dsodust', |
---|
3549 | & 'density scaled extinction opacity of std dust at 21.6um(MCS)', |
---|
3550 | & 'm2.kg-1',dsodust(:,:)) |
---|
3551 | call write_output('dso', |
---|
3552 | & 'density scaled extinction opacity of all dust at 21.6um(MCS)', |
---|
3553 | & 'm2.kg-1',dsodust(:,:)+dsords(:,:)+dsotop(:,:)) |
---|
3554 | end select |
---|
3555 | else ! (doubleq=.false.) |
---|
3556 | do iq=1,dustbin |
---|
3557 | write(str2(1:2),'(i2.2)') iq |
---|
3558 | call write_output('q'//str2,'mix. ratio', |
---|
3559 | & 'kg/kg',zq(:,:,iq)) |
---|
3560 | call write_output('qsurf'//str2,'qsurf', |
---|
3561 | & 'kg.m-2',qsurf(:,iq,iflat)) |
---|
3562 | do islope=1,nslope |
---|
3563 | write(str2(1:2),'(i2.2)') islope |
---|
3564 | call write_output('qsurf_slope'//str2,'qsurf', |
---|
3565 | & 'kg.m-2',qsurf(:,iq,islope)) |
---|
3566 | ENDDO |
---|
3567 | end do |
---|
3568 | endif ! (doubleq) |
---|
3569 | |
---|
3570 | if (rdstorm) then ! writediagfi tendencies stormdust tracers |
---|
3571 | call write_output('reffstormdust','reffstormdust', |
---|
3572 | & 'm',rstormdust(:,:)*ref_r0) |
---|
3573 | call write_output('mstormdtot', |
---|
3574 | & 'total mass of stormdust only', |
---|
3575 | & 'kg.m-2',mstormdtot(:)) |
---|
3576 | call write_output('mdusttot', |
---|
3577 | & 'total mass of dust only', |
---|
3578 | & 'kg.m-2',mdusttot(:)) |
---|
3579 | call write_output('rdsdqsdust', |
---|
3580 | & 'deposited surface stormdust mass', |
---|
3581 | & 'kg.m-2.s-1',rdsdqdustsurf(:)) |
---|
3582 | call write_output('rdsdustq','storm Dust mass mr', |
---|
3583 | & 'kg/kg',rdsqdust(:,:)) |
---|
3584 | call write_output('rdsdustqmodel','storm Dust massmr', |
---|
3585 | & 'kg/kg',pq(:,:,igcm_stormdust_mass)) |
---|
3586 | call write_output('rdsdustN','storm Dust number', |
---|
3587 | & 'part/kg',rdsndust(:,:)) |
---|
3588 | call write_output("stormfract", |
---|
3589 | & "fraction of the mesh, with stormdust","none", |
---|
3590 | & totstormfract(:)) |
---|
3591 | ! call write_output('qsurf', |
---|
3592 | ! & 'stormdust injection', |
---|
3593 | ! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,iflat)) |
---|
3594 | ! do islope=1,nslope |
---|
3595 | ! write(str2(1:2),'(i2.2)') islope |
---|
3596 | ! call write_output('qsurf_slope'//str2, |
---|
3597 | ! & 'stormdust injection', |
---|
3598 | ! & 'kg.m-2',qsurf(:,igcm_stormdust_mass,islope)) |
---|
3599 | ! ENDDO |
---|
3600 | ! call write_output('pdqsurf', |
---|
3601 | ! & 'tendancy stormdust mass at surface', |
---|
3602 | ! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,iflat)) |
---|
3603 | ! do islope=1,nslope |
---|
3604 | ! write(str2(1:2),'(i2.2)') islope |
---|
3605 | ! call write_output('pdqsurf_slope'//str2, |
---|
3606 | ! & 'tendancy stormdust mass at surface', |
---|
3607 | ! & 'kg.m-2',dqsurf(:,igcm_stormdust_mass,islope)) |
---|
3608 | ! ENDDO |
---|
3609 | call write_output('wspeed_stormdust', |
---|
3610 | & 'vertical velocity of stormdust', |
---|
3611 | & 'm/s',wspeed(:,:)) |
---|
3612 | call write_output('zdqsed_dust_mass' |
---|
3613 | & ,'sedimentation tendency of background dust mmr' |
---|
3614 | & ,'kg/kg.s-1', |
---|
3615 | & zdqsed(:,:,igcm_dust_mass)) |
---|
3616 | call write_output('zdqssed_dust_mass' |
---|
3617 | & ,'sedimentation tendency of background dust on surface' |
---|
3618 | & ,'kg.m-2.s-1', |
---|
3619 | & zdqssed(:,igcm_dust_mass)) |
---|
3620 | call write_output('zdqsed_stormdust_mass' |
---|
3621 | & ,'sedimentation tendency of stormdust dust mmr' |
---|
3622 | & ,'kg/kg.s-1', |
---|
3623 | & zdqsed(:,:,igcm_stormdust_mass)) |
---|
3624 | call write_output('zdqsed_dust_number' |
---|
3625 | & ,'sedimentation tendency of background dust number' |
---|
3626 | & ,'nbr/kg.s-1', |
---|
3627 | & zdqsed(:,:,igcm_dust_number)) |
---|
3628 | call write_output('rdust','rdust', |
---|
3629 | & 'm',rdust(:,:)) |
---|
3630 | call write_output('rstormdust','rstormdust', |
---|
3631 | & 'm',rstormdust(:,:)) |
---|
3632 | |
---|
3633 | select case (trim(dustiropacity)) |
---|
3634 | case ("tes") |
---|
3635 | call write_output('dsords_TES', |
---|
3636 | & 'density scaled extinction opacity of stormdust at 9.3um(TES)', |
---|
3637 | & 'm2.kg-1',dsords(:,:)) |
---|
3638 | case ("mcs") |
---|
3639 | call write_output('dsords', |
---|
3640 | & 'density scaled extinction opacity of stormdust at 21.6um(MCS)', |
---|
3641 | & 'm2.kg-1',dsords(:,:)) |
---|
3642 | end select |
---|
3643 | endif ! (rdstorm) |
---|
3644 | |
---|
3645 | if (topflows) then |
---|
3646 | call write_output('refftopdust', |
---|
3647 | & 'Topdust dust effective radius', |
---|
3648 | & 'm',rtopdust(:,:)*ref_r0) |
---|
3649 | call write_output('topdustq','top Dust mass mr', |
---|
3650 | & 'kg/kg',pq(:,:,igcm_topdust_mass)) |
---|
3651 | call write_output('topdustN','top Dust number', |
---|
3652 | & 'part/kg',pq(:,:,igcm_topdust_number)) |
---|
3653 | select case (trim(dustiropacity)) |
---|
3654 | case ("tes") |
---|
3655 | call write_output('dsotop_TES', |
---|
3656 | & 'density scaled extinction opacity of topdust at 9.3um(TES)', |
---|
3657 | & 'm2.kg-1',dsotop(:,:)) |
---|
3658 | case ("mcs") |
---|
3659 | call write_output('dsotop', |
---|
3660 | & 'density scaled extinction opacity of topdust at 21.6um(MCS)', |
---|
3661 | & 'm2.kg-1',dsotop(:,:)) |
---|
3662 | end select |
---|
3663 | endif ! (topflows) |
---|
3664 | |
---|
3665 | if (dustscaling_mode==2) then |
---|
3666 | call write_output("dust_rad_adjust", |
---|
3667 | & "radiative adjustment coefficient for dust", |
---|
3668 | & "",dust_rad_adjust(:)) |
---|
3669 | endif |
---|
3670 | |
---|
3671 | ! if (scavenging) then ! these outputs should be in the scavenging routine |
---|
3672 | ! call write_output('ccnq','CCN mass mr', |
---|
3673 | ! & 'kg/kg',qccn(:,:)) |
---|
3674 | ! call write_output('ccnN','CCN number', |
---|
3675 | ! & 'part/kg',nccn(:,:)) |
---|
3676 | ! call write_output('surfccnq','Surf nuclei mass mr', |
---|
3677 | ! & 'kg.m-2',qsurf(:,igcm_ccn_mass,iflat)) |
---|
3678 | ! do islope=1,nslope |
---|
3679 | ! write(str2(1:2),'(i2.2)') islope |
---|
3680 | ! call write_output('surfccnq_slope'//str2, |
---|
3681 | ! & 'Surf nuclei mass mr', |
---|
3682 | ! & 'kg.m-2',qsurf(:,igcm_ccn_mass,islope)) |
---|
3683 | ! ENDDO |
---|
3684 | ! call write_output('surfccnN','Surf nuclei number', |
---|
3685 | ! & 'kg.m-2',qsurf(:,igcm_ccn_number,iflat)) |
---|
3686 | ! do islope=1,nslope |
---|
3687 | ! write(str2(1:2),'(i2.2)') islope |
---|
3688 | ! call write_output('surfccnN_slope'//str2, |
---|
3689 | ! & 'Surf nuclei number', |
---|
3690 | ! & 'kg.m-2',qsurf(:,igcm_ccn_number,islope)) |
---|
3691 | ! ENDDO |
---|
3692 | ! endif ! (scavenging) |
---|
3693 | |
---|
3694 | #else |
---|
3695 | ! !!! to initialize mesoscale we need scaled variables |
---|
3696 | ! !!! because this must correspond to starting point for tracers |
---|
3697 | ! call write_output('dustq','Dust mass mr', |
---|
3698 | ! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_mass)) |
---|
3699 | ! call write_output('dustN','Dust number', |
---|
3700 | ! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_dust_number)) |
---|
3701 | ! call write_output('ccn','Nuclei mass mr', |
---|
3702 | ! & 'kg/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_mass)) |
---|
3703 | ! call write_output('ccnN','Nuclei number', |
---|
3704 | ! & 'part/kg',3,pq(1:ngrid,1:nlayer,igcm_ccn_number)) |
---|
3705 | if (freedust) then |
---|
3706 | call write_output('dustq','Dust mass mr', |
---|
3707 | & 'kg/kg',qdust) |
---|
3708 | call write_output('dustN','Dust number', |
---|
3709 | & 'part/kg',ndust) |
---|
3710 | call write_output('ccn','CCN mass mr', |
---|
3711 | & 'kg/kg',qccn) |
---|
3712 | call write_output('ccnN','CCN number', |
---|
3713 | & 'part/kg',nccn) |
---|
3714 | else |
---|
3715 | call write_output('dustq','Dust mass mr', |
---|
3716 | & 'kg/kg',pq(:,:,igcm_dust_mass)) |
---|
3717 | call write_output('dustN','Dust number', |
---|
3718 | & 'part/kg',pq(:,:,igcm_dust_number)) |
---|
3719 | call write_output('ccn','Nuclei mass mr', |
---|
3720 | & 'kg/kg',pq(:,:,igcm_ccn_mass)) |
---|
3721 | call write_output('ccnN','Nuclei number', |
---|
3722 | & 'part/kg',pq(:,:,igcm_ccn_number)) |
---|
3723 | endif |
---|
3724 | #endif |
---|
3725 | |
---|
3726 | end if ! (dustbin.ne.0) |
---|
3727 | |
---|
3728 | c ---------------------------------------------------------- |
---|
3729 | c Thermospheric outputs |
---|
3730 | c ---------------------------------------------------------- |
---|
3731 | |
---|
3732 | if(callthermos) then |
---|
3733 | |
---|
3734 | call write_output("quv","UV heating","K/s", |
---|
3735 | $ zdteuv(:,:)) |
---|
3736 | call write_output("cond","Thermal conduction","K/s", |
---|
3737 | $ zdtconduc(:,:)) |
---|
3738 | |
---|
3739 | !H, H2 and D escape fluxes |
---|
3740 | |
---|
3741 | call write_output("PhiH","H escape flux","s-1", |
---|
3742 | $ PhiEscH) |
---|
3743 | call write_output("PhiH2","H2 escape flux","s-1", |
---|
3744 | $ PhiEscH2) |
---|
3745 | call write_output("PhiD","D escape flux","s-1", |
---|
3746 | $ PhiEscD) |
---|
3747 | |
---|
3748 | endif !(callthermos) |
---|
3749 | |
---|
3750 | call write_output("q15um","15 um cooling","K/s", |
---|
3751 | $ zdtnlte(:,:)) |
---|
3752 | call write_output("qnir","NIR heating","K/s", |
---|
3753 | $ zdtnirco2(:,:)) |
---|
3754 | |
---|
3755 | c ---------------------------------------------------------- |
---|
3756 | c ---------------------------------------------------------- |
---|
3757 | c PBL OUTPUS |
---|
3758 | c ---------------------------------------------------------- |
---|
3759 | c ---------------------------------------------------------- |
---|
3760 | |
---|
3761 | c ---------------------------------------------------------- |
---|
3762 | c Outputs of thermals |
---|
3763 | c ---------------------------------------------------------- |
---|
3764 | if (calltherm) then |
---|
3765 | call write_output('lmax_th', |
---|
3766 | & 'index of vertical extension of thermals', |
---|
3767 | & 'grid level',lmax_th_out(:)) |
---|
3768 | call write_output('zmax_th', |
---|
3769 | & 'vertical extension of thermals','m', |
---|
3770 | & zmax_th(:)) |
---|
3771 | call write_output('hfmax_th', |
---|
3772 | & 'maximum heat flux in thermals','K.m/s', |
---|
3773 | & hfmax_th(:)) |
---|
3774 | call write_output('wstar', |
---|
3775 | & 'maximum thermals vertical velocity','m/s', |
---|
3776 | & wstar(:)) |
---|
3777 | end if |
---|
3778 | |
---|
3779 | c ---------------------------------------------------------- |
---|
3780 | c ---------------------------------------------------------- |
---|
3781 | c END OF PBL OUTPUS |
---|
3782 | c ---------------------------------------------------------- |
---|
3783 | c ---------------------------------------------------------- |
---|
3784 | |
---|
3785 | |
---|
3786 | c ---------------------------------------------------------- |
---|
3787 | c Output in netcdf file "diagsoil.nc" for subterranean |
---|
3788 | c variables (output every "ecritphy", as for writediagfi) |
---|
3789 | c ---------------------------------------------------------- |
---|
3790 | |
---|
3791 | ! Write soil temperature |
---|
3792 | call write_output("soiltemp","Soil temperature","K", |
---|
3793 | & tsoil(:,:,iflat)) |
---|
3794 | do islope=1,nslope |
---|
3795 | write(str2(1:2),'(i2.2)') islope |
---|
3796 | call write_output("soiltemp_slope"//str2, |
---|
3797 | & "Soil temperature","K", |
---|
3798 | & tsoil(:,:,islope)) |
---|
3799 | ENDDO |
---|
3800 | |
---|
3801 | !PREVIOUSLY IN 1D ONLY |
---|
3802 | call write_output("dtrad","rad. heat. rate", |
---|
3803 | & "K.s-1",dtrad(:,:)) |
---|
3804 | |
---|
3805 | if (rdstorm) then |
---|
3806 | call write_output('aerosol_dust','opacity of env. dust','' |
---|
3807 | & ,aerosol(:,:,iaer_dust_doubleq)) |
---|
3808 | call write_output('aerosol_stormdust', |
---|
3809 | & 'opacity of storm dust','' |
---|
3810 | & ,aerosol(:,:,iaer_stormdust_doubleq)) |
---|
3811 | call write_output('dqsdifdustq', |
---|
3812 | &'tendency due to vertical diffusion of background dust on surface' |
---|
3813 | & ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,iflat)) |
---|
3814 | do islope=1,nslope |
---|
3815 | write(str2(1:2),'(i2.2)') islope |
---|
3816 | call write_output('dqsdifdustq_slope'//str2, |
---|
3817 | &'tendency due to vertical diffusion of background dust on surface' |
---|
3818 | & ,'kg.m-2.s-1',zdqsdif(:,igcm_dust_mass,islope)) |
---|
3819 | ENDDO |
---|
3820 | call write_output('dqsdifrdsq', |
---|
3821 | & 'tendency due to vertical diffusion of stormdust on surface', |
---|
3822 | & 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,iflat)) |
---|
3823 | do islope=1,nslope |
---|
3824 | write(str2(1:2),'(i2.2)') islope |
---|
3825 | call write_output('dqsdifrdsq_slope'//str2, |
---|
3826 | & 'tendency due to vertical diffusion of stormdust on surface', |
---|
3827 | & 'kg.m-2.s-1',zdqsdif(:,igcm_stormdust_mass,islope)) |
---|
3828 | ENDDO |
---|
3829 | endif !(rdstorm) |
---|
3830 | |
---|
3831 | if(water) then |
---|
3832 | if (.not.scavenging) then |
---|
3833 | call write_output('zdqcloud_ice','cloud ice', |
---|
3834 | & 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice)) |
---|
3835 | call write_output('zdqcloud_vap','cloud vap', |
---|
3836 | & 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_vap)) |
---|
3837 | call write_output('zdqcloud','cloud', |
---|
3838 | & 'kg.m-2.s-1',zdqcloud(:,:,igcm_h2o_ice) |
---|
3839 | & +zdqcloud(:,:,igcm_h2o_vap)) |
---|
3840 | IF (hdo) THEN |
---|
3841 | call write_output('zdqcloud_iceD','cloud ice hdo', |
---|
3842 | & 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_ice)) |
---|
3843 | call write_output('zdqcloud_vapD','cloud vap hdo', |
---|
3844 | & 'kg.m-2.s-1',zdqcloud(:,:,igcm_hdo_vap)) |
---|
3845 | ENDIF ! hdo |
---|
3846 | endif !not.scavenging |
---|
3847 | ENDIF ! of IF (water) |
---|
3848 | |
---|
3849 | ! Output needed by the PEM |
---|
3850 | DO ig = 1,ngrid |
---|
3851 | ztmp1 =(1/m_co2 - 1/m_noco2) |
---|
3852 | ztmp2=1/m_noco2 |
---|
3853 | pvap_surf(ig) = 1/(ztmp1*zq(ig,1,igcm_co2)+ztmp2) |
---|
3854 | & * zq(ig,1,igcm_h2o_vap)/(mmol(igcm_h2o_vap)*1.e-3)*ps(ig) |
---|
3855 | |
---|
3856 | DO islope = 1,nslope |
---|
3857 | ! Clapeyron law for psat (psat = exp(beta/Th2o+alpha)),following Murphy and Koop 2005 |
---|
3858 | rhowater_surf_sat(ig,islope) = |
---|
3859 | & exp(beta_clap_h2o/tsurf(ig,islope)+alpha_clap_h2o) |
---|
3860 | & / tsurf(ig,islope) |
---|
3861 | & * mmol(igcm_h2o_vap)/(mugaz*r) |
---|
3862 | |
---|
3863 | if(qsurf(ig,igcm_h2o_ice,islope).gt.(1.e-4)) then |
---|
3864 | ! we consider to be at saturation above 1.e-4 kg.m-2 |
---|
3865 | rhowater_surf(ig,islope) = rhowater_surf_sat(ig,islope) |
---|
3866 | else |
---|
3867 | ! otherwise, use vapor partial pressure |
---|
3868 | rhowater_surf(ig,islope) = pvap_surf(ig) |
---|
3869 | & / tsurf(ig,islope) |
---|
3870 | & * mmol(igcm_h2o_vap)/(mugaz*r) |
---|
3871 | endif |
---|
3872 | DO isoil = 1,nsoilmx |
---|
3873 | rhowater_soil(ig,isoil,islope) = |
---|
3874 | & exp(beta_clap_h2o/tsoil(ig,isoil,islope)+alpha_clap_h2o) |
---|
3875 | & / tsoil(ig,isoil,islope) |
---|
3876 | & * mmol(igcm_h2o_vap)/(mugaz*r) |
---|
3877 | ENDDO |
---|
3878 | ENDDO |
---|
3879 | ENDDO |
---|
3880 | |
---|
3881 | DO islope = 1,nslope |
---|
3882 | write(str2(1:2),'(i2.2)') islope |
---|
3883 | CALL write_output("Waterdensity_soil_slope"//str2, |
---|
3884 | & "rhowater_soil_slope"//str2,'kg.m-3', |
---|
3885 | & rhowater_soil(:,:,islope)) |
---|
3886 | CALL write_output("Waterdensity_surface"//str2, |
---|
3887 | & "rhowater_surface"//str2,'kg.m-3', |
---|
3888 | & rhowater_surf(:,islope)) |
---|
3889 | ENDDO |
---|
3890 | |
---|
3891 | CALL write_output("h2o_layer1","h2o mass mr in the first layer", |
---|
3892 | & 'kg/kg',zq(:,1,igcm_h2o_vap)) |
---|
3893 | CALL write_output("co2_layer1","co2 mass mr in the first layer", |
---|
3894 | & 'kg/kg',zq(:,1,igcm_co2)) |
---|
3895 | |
---|
3896 | !PREVIOUSLY IN 1D ONLY |
---|
3897 | |
---|
3898 | c ========================================================== |
---|
3899 | c END OF WRITEDIAGFI |
---|
3900 | c ========================================================== |
---|
3901 | #endif |
---|
3902 | ! of ifdef MESOSCALE |
---|
3903 | |
---|
3904 | c ELSE ! if(ngrid.eq.1) |
---|
3905 | |
---|
3906 | c#ifndef MESOSCALE |
---|
3907 | c write(*, |
---|
3908 | c & '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)') |
---|
3909 | c & zls*180./pi,odpref,tau_pref_scenario |
---|
3910 | c#endif |
---|
3911 | |
---|
3912 | c END IF ! if(ngrid.ne.1) |
---|
3913 | |
---|
3914 | |
---|
3915 | ! test for co2 conservation with co2 microphysics |
---|
3916 | if (igcm_co2_ice.ne.0) then |
---|
3917 | co2totB = 0. ! added by C.M. |
---|
3918 | do ig=1,ngrid |
---|
3919 | do l=1,nlayer |
---|
3920 | co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g* |
---|
3921 | & (pq(ig,l,igcm_co2)+pq(ig,l,igcm_co2_ice) |
---|
3922 | & +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep) |
---|
3923 | enddo |
---|
3924 | co2totB = co2totB + qsurf(ig,igcm_co2,iflat) |
---|
3925 | enddo |
---|
3926 | else |
---|
3927 | co2totB = 0. ! added by C.M. |
---|
3928 | do ig=1,ngrid |
---|
3929 | do l=1,nlayer |
---|
3930 | co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g* |
---|
3931 | & (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep) |
---|
3932 | enddo |
---|
3933 | co2totB = co2totB + qsurf(ig,igcm_co2,iflat) |
---|
3934 | enddo |
---|
3935 | endif ! of if (igcm_co2_ice.ne.0) |
---|
3936 | co2conservation = (co2totA-co2totB)/co2totA |
---|
3937 | call write_output( 'co2conservation', |
---|
3938 | & 'Total CO2 mass conservation in physic', |
---|
3939 | & 'kg', co2conservation) |
---|
3940 | ! XIOS outputs |
---|
3941 | #ifdef CPP_XIOS |
---|
3942 | ! Send fields to XIOS: (NB these fields must also be defined as |
---|
3943 | ! <field id="..." /> in context_lmdz_physics.xml to be correctly used) |
---|
3944 | |
---|
3945 | CALL send_xios_field("controle",tab_cntrl_mod,1) |
---|
3946 | |
---|
3947 | CALL send_xios_field("ap",ap,1) |
---|
3948 | CALL send_xios_field("bp",bp,1) |
---|
3949 | CALL send_xios_field("aps",aps,1) |
---|
3950 | CALL send_xios_field("bps",bps,1) |
---|
3951 | |
---|
3952 | if (lastcall.and.is_omp_master) then |
---|
3953 | write(*,*) "physiq lastcall: call xios_context_finalize" |
---|
3954 | call xios_context_finalize |
---|
3955 | endif |
---|
3956 | #endif |
---|
3957 | |
---|
3958 | if (check_physics_outputs) then |
---|
3959 | ! Check the validity of updated fields at the end of the physics step |
---|
3960 | call check_physics_fields("end of physiq:",zt,zu,zv,zplev,zq) |
---|
3961 | endif |
---|
3962 | |
---|
3963 | icount=icount+1 |
---|
3964 | |
---|
3965 | END SUBROUTINE physiq |
---|
3966 | |
---|
3967 | END MODULE physiq_mod |
---|