/[lmdze]/trunk/phylmd/physiq.f
ViewVC logotype

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 191 - (show annotations)
Mon May 9 19:56:28 2016 UTC (8 years ago) by guez
Original Path: trunk/Sources/phylmd/physiq.f
File size: 49430 byte(s)
Extracted the call to read_comdissnew out of conf_gcm.

Made ok_instan a variable of module clesphys, itau_phy a variable of
module phyetat0_m, nid_ins a variable of module ini_histins_m, itap a
variable of new module time_phylmdz, so that histwrite_phy can be
called from any procedure without the need to cascade those variables
into that procedure. Made itau_w a variable of module time_phylmdz so
that it is computed only once per time step of physics.

Extracted variables of module clesphys which were in namelist
conf_phys_nml into their own namelist, clesphys_nml, and created
procedure read_clesphys reading clesphys_nml, to avoid side effect.

No need for double precision in procedure getso4fromfile. Assume there
is a single variable for the whole year in the NetCDF file instead of
one variable per month.

Created generic procedure histwrite_phy and removed procedure
write_histins, following LMDZ. histwrite_phy has only two arguments,
can be called from anywhere, and should manage the logic of writing or
not writing into various history files with various operations. So the
test on ok_instan goes inside histwrite_phy.

Test for raz_date in phyetat0 instead of physiq to avoid side effect.

Created procedure increment_itap to avoid side effect.

Removed unnecessary differences between procedures readsulfate and
readsulfate_pi.

1 module physiq_m
2
3 IMPLICIT none
4
5 contains
6
7 SUBROUTINE physiq(lafin, dayvrai, time, paprs, play, pphi, pphis, u, v, t, &
8 qx, omega, d_u, d_v, d_t, d_qx)
9
10 ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
11 ! (subversion revision 678)
12
13 ! Author: Z. X. Li (LMD/CNRS) 1993
14
15 ! This is the main procedure for the "physics" part of the program.
16
17 use aaam_bud_m, only: aaam_bud
18 USE abort_gcm_m, ONLY: abort_gcm
19 use aeropt_m, only: aeropt
20 use ajsec_m, only: ajsec
21 use calltherm_m, only: calltherm
22 USE clesphys, ONLY: cdhmax, cdmmax, ecrit_hf, ecrit_ins, ecrit_mth, &
23 ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin, ok_instan
24 USE clesphys2, ONLY: cycle_diurne, conv_emanuel, nbapp_rad, new_oliq, &
25 ok_orodr, ok_orolf
26 USE clmain_m, ONLY: clmain
27 use clouds_gno_m, only: clouds_gno
28 use comconst, only: dtphys
29 USE comgeomphy, ONLY: airephy
30 USE concvl_m, ONLY: concvl
31 USE conf_gcm_m, ONLY: offline, day_step, iphysiq
32 USE conf_phys_m, ONLY: conf_phys
33 use conflx_m, only: conflx
34 USE ctherm, ONLY: iflag_thermals, nsplit_thermals
35 use diagcld2_m, only: diagcld2
36 use diagetpq_m, only: diagetpq
37 use diagphy_m, only: diagphy
38 USE dimens_m, ONLY: llm, nqmx
39 USE dimphy, ONLY: klon
40 USE dimsoil, ONLY: nsoilmx
41 use drag_noro_m, only: drag_noro
42 use dynetat0_m, only: day_ref, annee_ref
43 USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
44 use fisrtilp_m, only: fisrtilp
45 USE hgardfou_m, ONLY: hgardfou
46 USE histsync_m, ONLY: histsync
47 USE histwrite_phy_m, ONLY: histwrite_phy
48 USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &
49 nbsrf
50 USE ini_histins_m, ONLY: ini_histins, nid_ins
51 use netcdf95, only: NF95_CLOSE
52 use newmicro_m, only: newmicro
53 use nr_util, only: assert
54 use nuage_m, only: nuage
55 USE orbite_m, ONLY: orbite
56 USE ozonecm_m, ONLY: ozonecm
57 USE phyetat0_m, ONLY: phyetat0, rlat, rlon
58 USE phyredem_m, ONLY: phyredem
59 USE phyredem0_m, ONLY: phyredem0
60 USE phystokenc_m, ONLY: phystokenc
61 USE phytrac_m, ONLY: phytrac
62 USE qcheck_m, ONLY: qcheck
63 use radlwsw_m, only: radlwsw
64 use readsulfate_m, only: readsulfate
65 use readsulfate_preind_m, only: readsulfate_preind
66 use yoegwd, only: sugwd
67 USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt
68 use time_phylmdz, only: itap, increment_itap
69 use transp_m, only: transp
70 use transp_lay_m, only: transp_lay
71 use unit_nml_m, only: unit_nml
72 USE ymds2ju_m, ONLY: ymds2ju
73 USE yoethf_m, ONLY: r2es, rvtmp2
74 use zenang_m, only: zenang
75
76 logical, intent(in):: lafin ! dernier passage
77
78 integer, intent(in):: dayvrai
79 ! current day number, based at value 1 on January 1st of annee_ref
80
81 REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
82
83 REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
84 ! pression pour chaque inter-couche, en Pa
85
86 REAL, intent(in):: play(:, :) ! (klon, llm)
87 ! pression pour le mileu de chaque couche (en Pa)
88
89 REAL, intent(in):: pphi(:, :) ! (klon, llm)
90 ! géopotentiel de chaque couche (référence sol)
91
92 REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
93
94 REAL, intent(in):: u(:, :) ! (klon, llm)
95 ! vitesse dans la direction X (de O a E) en m/s
96
97 REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m/s
98 REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
99
100 REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
101 ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
102
103 REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa/s
104 REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)
105 REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
106 REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K/s)
107
108 REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
109 ! tendance physique de "qx" (s-1)
110
111 ! Local:
112
113 LOGICAL:: firstcal = .true.
114
115 LOGICAL, PARAMETER:: check = .FALSE.
116 ! Verifier la conservation du modele en eau
117
118 LOGICAL, PARAMETER:: ok_stratus = .FALSE.
119 ! Ajouter artificiellement les stratus
120
121 ! pour phsystoke avec thermiques
122 REAL fm_therm(klon, llm + 1)
123 REAL entr_therm(klon, llm)
124 real, save:: q2(klon, llm + 1, nbsrf)
125
126 INTEGER, PARAMETER:: ivap = 1 ! indice de traceur pour vapeur d'eau
127 INTEGER, PARAMETER:: iliq = 2 ! indice de traceur pour eau liquide
128
129 REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
130 LOGICAL, save:: ancien_ok
131
132 REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)
133 REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg/kg/s)
134
135 real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
136
137 REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
138 REAL swup0(klon, llm + 1), swup(klon, llm + 1)
139 SAVE swdn0, swdn, swup0, swup
140
141 REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
142 REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
143 SAVE lwdn0, lwdn, lwup0, lwup
144
145 ! prw: precipitable water
146 real prw(klon)
147
148 ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
149 ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
150 REAL flwp(klon), fiwp(klon)
151 REAL flwc(klon, llm), fiwc(klon, llm)
152
153 ! Variables propres a la physique
154
155 INTEGER, save:: radpas
156 ! Radiative transfer computations are made every "radpas" call to
157 ! "physiq".
158
159 REAL radsol(klon)
160 SAVE radsol ! bilan radiatif au sol calcule par code radiatif
161
162 REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
163
164 REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
165 ! soil temperature of surface fraction
166
167 REAL, save:: fevap(klon, nbsrf) ! evaporation
168 REAL fluxlat(klon, nbsrf)
169 SAVE fluxlat
170
171 REAL, save:: fqsurf(klon, nbsrf)
172 ! humidite de l'air au contact de la surface
173
174 REAL, save:: qsol(klon)
175 ! column-density of water in soil, in kg m-2
176
177 REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse
178 REAL, save:: falbe(klon, nbsrf) ! albedo visible par type de surface
179
180 ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
181 REAL, save:: zmea(klon) ! orographie moyenne
182 REAL, save:: zstd(klon) ! deviation standard de l'OESM
183 REAL, save:: zsig(klon) ! pente de l'OESM
184 REAL, save:: zgam(klon) ! anisotropie de l'OESM
185 REAL, save:: zthe(klon) ! orientation de l'OESM
186 REAL, save:: zpic(klon) ! Maximum de l'OESM
187 REAL, save:: zval(klon) ! Minimum de l'OESM
188 REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
189 REAL zulow(klon), zvlow(klon)
190 INTEGER igwd, itest(klon)
191
192 REAL, save:: agesno(klon, nbsrf) ! age de la neige
193 REAL, save:: run_off_lic_0(klon)
194
195 ! Variables li\'ees \`a la convection d'Emanuel :
196 REAL, save:: Ma(klon, llm) ! undilute upward mass flux
197 REAL, save:: qcondc(klon, llm) ! in-cld water content from convect
198 REAL, save:: sig1(klon, llm), w01(klon, llm)
199
200 ! Variables pour la couche limite (Alain Lahellec) :
201 REAL cdragh(klon) ! drag coefficient pour T and Q
202 REAL cdragm(klon) ! drag coefficient pour vent
203
204 ! Pour phytrac :
205 REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
206 REAL yu1(klon) ! vents dans la premiere couche U
207 REAL yv1(klon) ! vents dans la premiere couche V
208 REAL ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige
209
210 REAL fqcalving(klon, nbsrf)
211 ! flux d'eau "perdue" par la surface et necessaire pour limiter la
212 ! hauteur de neige, en kg/m2/s
213
214 REAL zxffonte(klon), zxfqcalving(klon)
215
216 REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
217 save pfrac_impa
218 REAL pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation
219 save pfrac_nucl
220 REAL pfrac_1nucl(klon, llm)! Produits des coefs lessi nucl (alpha = 1)
221 save pfrac_1nucl
222 REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
223 REAL frac_nucl(klon, llm) ! idem (nucleation)
224
225 REAL, save:: rain_fall(klon)
226 ! liquid water mass flux (kg/m2/s), positive down
227
228 REAL, save:: snow_fall(klon)
229 ! solid water mass flux (kg/m2/s), positive down
230
231 REAL rain_tiedtke(klon), snow_tiedtke(klon)
232
233 REAL evap(klon), devap(klon) ! evaporation and its derivative
234 REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
235 REAL dlw(klon) ! derivee infra rouge
236 SAVE dlw
237 REAL bils(klon) ! bilan de chaleur au sol
238 REAL, save:: fder(klon) ! Derive de flux (sensible et latente)
239 REAL ve(klon) ! integr. verticale du transport meri. de l'energie
240 REAL vq(klon) ! integr. verticale du transport meri. de l'eau
241 REAL ue(klon) ! integr. verticale du transport zonal de l'energie
242 REAL uq(klon) ! integr. verticale du transport zonal de l'eau
243
244 REAL, save:: frugs(klon, nbsrf) ! longueur de rugosite
245 REAL zxrugs(klon) ! longueur de rugosite
246
247 ! Conditions aux limites
248
249 INTEGER julien
250 INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
251 REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
252 REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
253 REAL, save:: albsol(klon) ! albedo du sol total visible
254 REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
255
256 real, save:: clwcon(klon, llm), rnebcon(klon, llm)
257 real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
258
259 REAL rhcl(klon, llm) ! humiditi relative ciel clair
260 REAL dialiq(klon, llm) ! eau liquide nuageuse
261 REAL diafra(klon, llm) ! fraction nuageuse
262 REAL cldliq(klon, llm) ! eau liquide nuageuse
263 REAL cldfra(klon, llm) ! fraction nuageuse
264 REAL cldtau(klon, llm) ! epaisseur optique
265 REAL cldemi(klon, llm) ! emissivite infrarouge
266
267 REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite
268 REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur
269 REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u
270 REAL fluxv(klon, llm, nbsrf) ! flux turbulent de vitesse v
271
272 REAL zxfluxt(klon, llm)
273 REAL zxfluxq(klon, llm)
274 REAL zxfluxu(klon, llm)
275 REAL zxfluxv(klon, llm)
276
277 ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
278 ! les variables soient r\'emanentes.
279 REAL, save:: heat(klon, llm) ! chauffage solaire
280 REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
281 REAL, save:: cool(klon, llm) ! refroidissement infrarouge
282 REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
283 REAL, save:: topsw(klon), toplw(klon), solsw(klon)
284 REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
285 real, save:: sollwdown(klon) ! downward LW flux at surface
286 REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
287 REAL, save:: albpla(klon)
288 REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
289 REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
290
291 REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
292 REAL conv_t(klon, llm) ! convergence of temperature (K/s)
293
294 REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut
295 REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree
296
297 REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
298
299 REAL dist, mu0(klon), fract(klon)
300 real longi
301 REAL z_avant(klon), z_apres(klon), z_factor(klon)
302 REAL za, zb
303 REAL zx_t, zx_qs, zcor
304 real zqsat(klon, llm)
305 INTEGER i, k, iq, nsrf
306 REAL, PARAMETER:: t_coup = 234.
307 REAL zphi(klon, llm)
308
309 ! cf. Anne Mathieu variables pour la couche limite atmosphérique (hbtm)
310
311 REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
312 REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
313 REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
314 REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
315 REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
316 REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite
317 REAL, SAVE:: therm(klon, nbsrf)
318 REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
319 REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
320 REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
321 ! Grandeurs de sorties
322 REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
323 REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
324 REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
325 REAL s_trmb3(klon)
326
327 ! Variables pour la convection de K. Emanuel :
328
329 REAL upwd(klon, llm) ! saturated updraft mass flux
330 REAL dnwd(klon, llm) ! saturated downdraft mass flux
331 REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
332 REAL cape(klon) ! CAPE
333 SAVE cape
334
335 INTEGER iflagctrl(klon) ! flag fonctionnement de convect
336
337 ! Variables du changement
338
339 ! con: convection
340 ! lsc: large scale condensation
341 ! ajs: ajustement sec
342 ! eva: \'evaporation de l'eau liquide nuageuse
343 ! vdf: vertical diffusion in boundary layer
344 REAL d_t_con(klon, llm), d_q_con(klon, llm)
345 REAL d_u_con(klon, llm), d_v_con(klon, llm)
346 REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
347 REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
348 REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
349 REAL rneb(klon, llm)
350
351 REAL mfu(klon, llm), mfd(klon, llm)
352 REAL pen_u(klon, llm), pen_d(klon, llm)
353 REAL pde_u(klon, llm), pde_d(klon, llm)
354 INTEGER kcbot(klon), kctop(klon), kdtop(klon)
355 REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
356 REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
357
358 INTEGER, save:: ibas_con(klon), itop_con(klon)
359 real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
360
361 REAL rain_con(klon), rain_lsc(klon)
362 REAL, save:: snow_con(klon) ! neige (mm / s)
363 real snow_lsc(klon)
364 REAL d_ts(klon, nbsrf)
365
366 REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
367 REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
368
369 REAL d_u_oro(klon, llm), d_v_oro(klon, llm)
370 REAL d_t_oro(klon, llm)
371 REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
372 REAL d_t_lif(klon, llm)
373
374 REAL, save:: ratqs(klon, llm)
375 real ratqss(klon, llm), ratqsc(klon, llm)
376 real:: ratqsbas = 0.01, ratqshaut = 0.3
377
378 ! Parametres lies au nouveau schema de nuages (SB, PDF)
379 real:: fact_cldcon = 0.375
380 real:: facttemps = 1.e-4
381 logical:: ok_newmicro = .true.
382 real facteur
383
384 integer:: iflag_cldcon = 1
385 logical ptconv(klon, llm)
386
387 ! Variables pour effectuer les appels en s\'erie :
388
389 REAL t_seri(klon, llm), q_seri(klon, llm)
390 REAL ql_seri(klon, llm)
391 REAL u_seri(klon, llm), v_seri(klon, llm)
392 REAL tr_seri(klon, llm, nqmx - 2)
393
394 REAL zx_rh(klon, llm)
395
396 REAL zustrdr(klon), zvstrdr(klon)
397 REAL zustrli(klon), zvstrli(klon)
398 REAL zustrph(klon), zvstrph(klon)
399 REAL aam, torsfc
400
401 REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
402 REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
403 REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
404 REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
405
406 real date0
407
408 ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
409 REAL ztsol(klon)
410 REAL d_h_vcol, d_qt, d_ec
411 REAL, SAVE:: d_h_vcol_phy
412 REAL zero_v(klon)
413 CHARACTER(LEN = 20) tit
414 INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
415 INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
416
417 REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique
418 REAL ZRCPD
419
420 REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m
421 REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
422 REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille
423 REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes s/1 maille
424
425 ! Aerosol effects:
426
427 REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
428
429 REAL, save:: sulfate_pi(klon, llm)
430 ! SO4 aerosol concentration, in \mu g/m3, pre-industrial value
431
432 REAL cldtaupi(klon, llm)
433 ! cloud optical thickness for pre-industrial aerosols
434
435 REAL re(klon, llm) ! Cloud droplet effective radius
436 REAL fl(klon, llm) ! denominator of re
437
438 ! Aerosol optical properties
439 REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
440 REAL, save:: cg_ae(klon, llm, 2)
441
442 REAL topswad(klon), solswad(klon) ! aerosol direct effect
443 REAL topswai(klon), solswai(klon) ! aerosol indirect effect
444
445 REAL aerindex(klon) ! POLDER aerosol index
446
447 LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
448 LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
449
450 REAL:: bl95_b0 = 2., bl95_b1 = 0.2
451 ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
452 ! B). They link cloud droplet number concentration to aerosol mass
453 ! concentration.
454
455 SAVE u10m
456 SAVE v10m
457 SAVE t2m
458 SAVE q2m
459 SAVE ffonte
460 SAVE fqcalving
461 SAVE rain_con
462 SAVE topswai
463 SAVE topswad
464 SAVE solswai
465 SAVE solswad
466 SAVE d_u_con
467 SAVE d_v_con
468
469 real zmasse(klon, llm)
470 ! (column-density of mass of air in a cell, in kg m-2)
471
472 integer, save:: ncid_startphy
473
474 namelist /physiq_nml/ fact_cldcon, facttemps, ok_newmicro, &
475 iflag_cldcon, ratqsbas, ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, &
476 bl95_b1, iflag_thermals, nsplit_thermals
477
478 !----------------------------------------------------------------
479
480 IF (if_ebil >= 1) zero_v = 0.
481 IF (nqmx < 2) CALL abort_gcm('physiq', &
482 'eaux vapeur et liquide sont indispensables')
483
484 test_firstcal: IF (firstcal) THEN
485 ! initialiser
486 u10m = 0.
487 v10m = 0.
488 t2m = 0.
489 q2m = 0.
490 ffonte = 0.
491 fqcalving = 0.
492 piz_ae = 0.
493 tau_ae = 0.
494 cg_ae = 0.
495 rain_con = 0.
496 snow_con = 0.
497 topswai = 0.
498 topswad = 0.
499 solswai = 0.
500 solswad = 0.
501
502 d_u_con = 0.
503 d_v_con = 0.
504 rnebcon0 = 0.
505 clwcon0 = 0.
506 rnebcon = 0.
507 clwcon = 0.
508
509 pblh =0. ! Hauteur de couche limite
510 plcl =0. ! Niveau de condensation de la CLA
511 capCL =0. ! CAPE de couche limite
512 oliqCL =0. ! eau_liqu integree de couche limite
513 cteiCL =0. ! cloud top instab. crit. couche limite
514 pblt =0. ! T a la Hauteur de couche limite
515 therm =0.
516 trmb1 =0. ! deep_cape
517 trmb2 =0. ! inhibition
518 trmb3 =0. ! Point Omega
519
520 IF (if_ebil >= 1) d_h_vcol_phy = 0.
521
522 iflag_thermals = 0
523 nsplit_thermals = 1
524 print *, "Enter namelist 'physiq_nml'."
525 read(unit=*, nml=physiq_nml)
526 write(unit_nml, nml=physiq_nml)
527
528 call conf_phys
529
530 ! Initialiser les compteurs:
531
532 frugs = 0.
533 CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
534 fevap, rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
535 agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
536 q_ancien, ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
537 w01, ncid_startphy)
538
539 ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
540 q2 = 1e-8
541
542 lmt_pas = day_step / iphysiq
543 print *, 'Number of time steps of "physics" per day: ', lmt_pas
544
545 radpas = lmt_pas / nbapp_rad
546 print *, "radpas = ", radpas
547
548 ! Initialisation pour le sch\'ema de convection d'Emanuel :
549 IF (conv_emanuel) THEN
550 ibas_con = 1
551 itop_con = 1
552 ENDIF
553
554 IF (ok_orodr) THEN
555 rugoro = MAX(1e-5, zstd * zsig / 2)
556 CALL SUGWD(paprs, play)
557 else
558 rugoro = 0.
559 ENDIF
560
561 ecrit_ins = NINT(ecrit_ins/dtphys)
562 ecrit_hf = NINT(ecrit_hf/dtphys)
563 ecrit_mth = NINT(ecrit_mth/dtphys)
564 ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
565 ecrit_reg = NINT(ecrit_reg/dtphys)
566
567 ! Initialisation des sorties
568
569 call ini_histins(dtphys)
570 CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
571 ! Positionner date0 pour initialisation de ORCHIDEE
572 print *, 'physiq date0: ', date0
573 CALL phyredem0(lmt_pas)
574 ENDIF test_firstcal
575
576 ! We will modify variables *_seri and we will not touch variables
577 ! u, v, t, qx:
578 t_seri = t
579 u_seri = u
580 v_seri = v
581 q_seri = qx(:, :, ivap)
582 ql_seri = qx(:, :, iliq)
583 tr_seri = qx(:, :, 3:nqmx)
584
585 ztsol = sum(ftsol * pctsrf, dim = 2)
586
587 IF (if_ebil >= 1) THEN
588 tit = 'after dynamics'
589 CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
590 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
591 ! Comme les tendances de la physique sont ajout\'es dans la
592 ! dynamique, la variation d'enthalpie par la dynamique devrait
593 ! \^etre \'egale \`a la variation de la physique au pas de temps
594 ! pr\'ec\'edent. Donc la somme de ces 2 variations devrait \^etre
595 ! nulle.
596 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
597 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
598 d_qt, 0.)
599 END IF
600
601 ! Diagnostic de la tendance dynamique :
602 IF (ancien_ok) THEN
603 DO k = 1, llm
604 DO i = 1, klon
605 d_t_dyn(i, k) = (t_seri(i, k) - t_ancien(i, k)) / dtphys
606 d_q_dyn(i, k) = (q_seri(i, k) - q_ancien(i, k)) / dtphys
607 ENDDO
608 ENDDO
609 ELSE
610 DO k = 1, llm
611 DO i = 1, klon
612 d_t_dyn(i, k) = 0.
613 d_q_dyn(i, k) = 0.
614 ENDDO
615 ENDDO
616 ancien_ok = .TRUE.
617 ENDIF
618
619 ! Ajouter le geopotentiel du sol:
620 DO k = 1, llm
621 DO i = 1, klon
622 zphi(i, k) = pphi(i, k) + pphis(i)
623 ENDDO
624 ENDDO
625
626 ! Check temperatures:
627 CALL hgardfou(t_seri, ftsol)
628
629 call increment_itap
630 julien = MOD(dayvrai, 360)
631 if (julien == 0) julien = 360
632
633 forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
634
635 ! Prescrire l'ozone :
636 wo = ozonecm(REAL(julien), paprs)
637
638 ! \'Evaporation de l'eau liquide nuageuse :
639 DO k = 1, llm
640 DO i = 1, klon
641 zb = MAX(0., ql_seri(i, k))
642 t_seri(i, k) = t_seri(i, k) &
643 - zb * RLVTT / RCPD / (1. + RVTMP2 * q_seri(i, k))
644 q_seri(i, k) = q_seri(i, k) + zb
645 ENDDO
646 ENDDO
647 ql_seri = 0.
648
649 IF (if_ebil >= 2) THEN
650 tit = 'after reevap'
651 CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
652 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
653 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
654 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
655 END IF
656
657 frugs = MAX(frugs, 0.000015)
658 zxrugs = sum(frugs * pctsrf, dim = 2)
659
660 ! Calculs n\'ecessaires au calcul de l'albedo dans l'interface avec
661 ! la surface.
662
663 CALL orbite(REAL(julien), longi, dist)
664 IF (cycle_diurne) THEN
665 CALL zenang(longi, time, dtphys * radpas, mu0, fract)
666 ELSE
667 mu0 = - 999.999
668 ENDIF
669
670 ! Calcul de l'abedo moyen par maille
671 albsol = sum(falbe * pctsrf, dim = 2)
672
673 ! R\'epartition sous maille des flux longwave et shortwave
674 ! R\'epartition du longwave par sous-surface lin\'earis\'ee
675
676 forall (nsrf = 1: nbsrf)
677 fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
678 * (ztsol - ftsol(:, nsrf))
679 fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
680 END forall
681
682 fder = dlw
683
684 ! Couche limite:
685
686 CALL clmain(dtphys, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, v_seri, &
687 julien, mu0, ftsol, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, &
688 ftsoil, qsol, paprs, play, fsnow, fqsurf, fevap, falbe, fluxlat, &
689 rain_fall, snow_fall, fsolsw, fsollw, fder, rlat, frugs, firstcal, &
690 agesno, rugoro, d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, &
691 fluxq, fluxu, fluxv, cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, &
692 yv1, t2m, q2m, u10m, v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, &
693 trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)
694
695 ! Incr\'ementation des flux
696
697 zxfluxt = 0.
698 zxfluxq = 0.
699 zxfluxu = 0.
700 zxfluxv = 0.
701 DO nsrf = 1, nbsrf
702 DO k = 1, llm
703 DO i = 1, klon
704 zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
705 zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
706 zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
707 zxfluxv(i, k) = zxfluxv(i, k) + fluxv(i, k, nsrf) * pctsrf(i, nsrf)
708 END DO
709 END DO
710 END DO
711 DO i = 1, klon
712 sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
713 evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
714 fder(i) = dlw(i) + dsens(i) + devap(i)
715 ENDDO
716
717 DO k = 1, llm
718 DO i = 1, klon
719 t_seri(i, k) = t_seri(i, k) + d_t_vdf(i, k)
720 q_seri(i, k) = q_seri(i, k) + d_q_vdf(i, k)
721 u_seri(i, k) = u_seri(i, k) + d_u_vdf(i, k)
722 v_seri(i, k) = v_seri(i, k) + d_v_vdf(i, k)
723 ENDDO
724 ENDDO
725
726 IF (if_ebil >= 2) THEN
727 tit = 'after clmain'
728 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
729 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
730 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
731 sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
732 END IF
733
734 ! Update surface temperature:
735
736 DO i = 1, klon
737 zxtsol(i) = 0.
738 zxfluxlat(i) = 0.
739
740 zt2m(i) = 0.
741 zq2m(i) = 0.
742 zu10m(i) = 0.
743 zv10m(i) = 0.
744 zxffonte(i) = 0.
745 zxfqcalving(i) = 0.
746
747 s_pblh(i) = 0.
748 s_lcl(i) = 0.
749 s_capCL(i) = 0.
750 s_oliqCL(i) = 0.
751 s_cteiCL(i) = 0.
752 s_pblT(i) = 0.
753 s_therm(i) = 0.
754 s_trmb1(i) = 0.
755 s_trmb2(i) = 0.
756 s_trmb3(i) = 0.
757 ENDDO
758
759 call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
760
761 DO nsrf = 1, nbsrf
762 DO i = 1, klon
763 ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)
764 zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
765 zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)
766
767 zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)
768 zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)
769 zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)
770 zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)
771 zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)
772 zxfqcalving(i) = zxfqcalving(i) + &
773 fqcalving(i, nsrf)*pctsrf(i, nsrf)
774 s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)
775 s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)
776 s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)
777 s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)
778 s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)
779 s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)
780 s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)
781 s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)
782 s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)
783 s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)
784 ENDDO
785 ENDDO
786
787 ! Si une sous-fraction n'existe pas, elle prend la température moyenne :
788 DO nsrf = 1, nbsrf
789 DO i = 1, klon
790 IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
791
792 IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)
793 IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)
794 IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)
795 IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)
796 IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)
797 IF (pctsrf(i, nsrf) < epsfra) &
798 fqcalving(i, nsrf) = zxfqcalving(i)
799 IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf) = s_pblh(i)
800 IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf) = s_lcl(i)
801 IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf) = s_capCL(i)
802 IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf) = s_oliqCL(i)
803 IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf) = s_cteiCL(i)
804 IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf) = s_pblT(i)
805 IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf) = s_therm(i)
806 IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf) = s_trmb1(i)
807 IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf) = s_trmb2(i)
808 IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf) = s_trmb3(i)
809 ENDDO
810 ENDDO
811
812 ! Calculer la dérive du flux infrarouge
813
814 DO i = 1, klon
815 dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
816 ENDDO
817
818 IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
819
820 ! Appeler la convection
821
822 if (conv_emanuel) then
823 da = 0.
824 mp = 0.
825 phi = 0.
826 CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
827 w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, &
828 itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, qcondc, pmflxr, &
829 da, phi, mp)
830 snow_con = 0.
831 clwcon0 = qcondc
832 mfu = upwd + dnwd
833
834 IF (thermcep) THEN
835 zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
836 zqsat = zqsat / (1. - retv * zqsat)
837 ELSE
838 zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
839 ENDIF
840
841 ! Properties of convective clouds
842 clwcon0 = fact_cldcon * clwcon0
843 call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
844 rnebcon0)
845
846 forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
847 mfd = 0.
848 pen_u = 0.
849 pen_d = 0.
850 pde_d = 0.
851 pde_u = 0.
852 else
853 conv_q = d_q_dyn + d_q_vdf / dtphys
854 conv_t = d_t_dyn + d_t_vdf / dtphys
855 z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
856 CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
857 q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
858 d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
859 mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
860 kdtop, pmflxr, pmflxs)
861 WHERE (rain_con < 0.) rain_con = 0.
862 WHERE (snow_con < 0.) snow_con = 0.
863 ibas_con = llm + 1 - kcbot
864 itop_con = llm + 1 - kctop
865 END if
866
867 DO k = 1, llm
868 DO i = 1, klon
869 t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
870 q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
871 u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
872 v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
873 ENDDO
874 ENDDO
875
876 IF (if_ebil >= 2) THEN
877 tit = 'after convect'
878 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
879 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
880 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
881 zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec)
882 END IF
883
884 IF (check) THEN
885 za = qcheck(paprs, q_seri, ql_seri)
886 print *, "aprescon = ", za
887 zx_t = 0.
888 za = 0.
889 DO i = 1, klon
890 za = za + airephy(i)/REAL(klon)
891 zx_t = zx_t + (rain_con(i)+ &
892 snow_con(i))*airephy(i)/REAL(klon)
893 ENDDO
894 zx_t = zx_t/za*dtphys
895 print *, "Precip = ", zx_t
896 ENDIF
897
898 IF (.not. conv_emanuel) THEN
899 z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
900 z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
901 DO k = 1, llm
902 DO i = 1, klon
903 IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
904 q_seri(i, k) = q_seri(i, k) * z_factor(i)
905 ENDIF
906 ENDDO
907 ENDDO
908 ENDIF
909
910 ! Convection s\`eche (thermiques ou ajustement)
911
912 d_t_ajs = 0.
913 d_u_ajs = 0.
914 d_v_ajs = 0.
915 d_q_ajs = 0.
916 fm_therm = 0.
917 entr_therm = 0.
918
919 if (iflag_thermals == 0) then
920 ! Ajustement sec
921 CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)
922 t_seri = t_seri + d_t_ajs
923 q_seri = q_seri + d_q_ajs
924 else
925 ! Thermiques
926 call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
927 q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
928 endif
929
930 IF (if_ebil >= 2) THEN
931 tit = 'after dry_adjust'
932 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
933 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
934 END IF
935
936 ! Caclul des ratqs
937
938 ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
939 ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
940 if (iflag_cldcon == 1) then
941 do k = 1, llm
942 do i = 1, klon
943 if(ptconv(i, k)) then
944 ratqsc(i, k) = ratqsbas + fact_cldcon &
945 * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
946 else
947 ratqsc(i, k) = 0.
948 endif
949 enddo
950 enddo
951 endif
952
953 ! ratqs stables
954 do k = 1, llm
955 do i = 1, klon
956 ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
957 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
958 enddo
959 enddo
960
961 ! ratqs final
962 if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
963 ! les ratqs sont une conbinaison de ratqss et ratqsc
964 ! ratqs final
965 ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
966 ! relaxation des ratqs
967 ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
968 ratqs = max(ratqs, ratqsc)
969 else
970 ! on ne prend que le ratqs stable pour fisrtilp
971 ratqs = ratqss
972 endif
973
974 CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
975 d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
976 pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
977 psfl, rhcl)
978
979 WHERE (rain_lsc < 0) rain_lsc = 0.
980 WHERE (snow_lsc < 0) snow_lsc = 0.
981 DO k = 1, llm
982 DO i = 1, klon
983 t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
984 q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
985 ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
986 cldfra(i, k) = rneb(i, k)
987 IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
988 ENDDO
989 ENDDO
990 IF (check) THEN
991 za = qcheck(paprs, q_seri, ql_seri)
992 print *, "apresilp = ", za
993 zx_t = 0.
994 za = 0.
995 DO i = 1, klon
996 za = za + airephy(i)/REAL(klon)
997 zx_t = zx_t + (rain_lsc(i) &
998 + snow_lsc(i))*airephy(i)/REAL(klon)
999 ENDDO
1000 zx_t = zx_t/za*dtphys
1001 print *, "Precip = ", zx_t
1002 ENDIF
1003
1004 IF (if_ebil >= 2) THEN
1005 tit = 'after fisrt'
1006 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1007 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1008 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1009 zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec)
1010 END IF
1011
1012 ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1013
1014 ! 1. NUAGES CONVECTIFS
1015
1016 IF (iflag_cldcon <= - 1) THEN
1017 ! seulement pour Tiedtke
1018 snow_tiedtke = 0.
1019 if (iflag_cldcon == - 1) then
1020 rain_tiedtke = rain_con
1021 else
1022 rain_tiedtke = 0.
1023 do k = 1, llm
1024 do i = 1, klon
1025 if (d_q_con(i, k) < 0.) then
1026 rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k)/dtphys &
1027 *zmasse(i, k)
1028 endif
1029 enddo
1030 enddo
1031 endif
1032
1033 ! Nuages diagnostiques pour Tiedtke
1034 CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1035 itop_con, diafra, dialiq)
1036 DO k = 1, llm
1037 DO i = 1, klon
1038 IF (diafra(i, k) > cldfra(i, k)) THEN
1039 cldliq(i, k) = dialiq(i, k)
1040 cldfra(i, k) = diafra(i, k)
1041 ENDIF
1042 ENDDO
1043 ENDDO
1044 ELSE IF (iflag_cldcon == 3) THEN
1045 ! On prend pour les nuages convectifs le maximum du calcul de
1046 ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
1047 ! d'un facteur facttemps.
1048 facteur = dtphys * facttemps
1049 do k = 1, llm
1050 do i = 1, klon
1051 rnebcon(i, k) = rnebcon(i, k) * facteur
1052 if (rnebcon0(i, k) * clwcon0(i, k) &
1053 > rnebcon(i, k) * clwcon(i, k)) then
1054 rnebcon(i, k) = rnebcon0(i, k)
1055 clwcon(i, k) = clwcon0(i, k)
1056 endif
1057 enddo
1058 enddo
1059
1060 ! On prend la somme des fractions nuageuses et des contenus en eau
1061 cldfra = min(max(cldfra, rnebcon), 1.)
1062 cldliq = cldliq + rnebcon*clwcon
1063 ENDIF
1064
1065 ! 2. Nuages stratiformes
1066
1067 IF (ok_stratus) THEN
1068 CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
1069 DO k = 1, llm
1070 DO i = 1, klon
1071 IF (diafra(i, k) > cldfra(i, k)) THEN
1072 cldliq(i, k) = dialiq(i, k)
1073 cldfra(i, k) = diafra(i, k)
1074 ENDIF
1075 ENDDO
1076 ENDDO
1077 ENDIF
1078
1079 ! Precipitation totale
1080 DO i = 1, klon
1081 rain_fall(i) = rain_con(i) + rain_lsc(i)
1082 snow_fall(i) = snow_con(i) + snow_lsc(i)
1083 ENDDO
1084
1085 IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1086 dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1087 d_qt, d_ec)
1088
1089 ! Humidit\'e relative pour diagnostic :
1090 DO k = 1, llm
1091 DO i = 1, klon
1092 zx_t = t_seri(i, k)
1093 IF (thermcep) THEN
1094 zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)
1095 zx_qs = MIN(0.5, zx_qs)
1096 zcor = 1./(1. - retv*zx_qs)
1097 zx_qs = zx_qs*zcor
1098 ELSE
1099 IF (zx_t < t_coup) THEN
1100 zx_qs = qsats(zx_t)/play(i, k)
1101 ELSE
1102 zx_qs = qsatl(zx_t)/play(i, k)
1103 ENDIF
1104 ENDIF
1105 zx_rh(i, k) = q_seri(i, k)/zx_qs
1106 zqsat(i, k) = zx_qs
1107 ENDDO
1108 ENDDO
1109
1110 ! Introduce the aerosol direct and first indirect radiative forcings:
1111 IF (ok_ade .OR. ok_aie) THEN
1112 ! Get sulfate aerosol distribution :
1113 CALL readsulfate(dayvrai, time, firstcal, sulfate)
1114 CALL readsulfate_preind(dayvrai, time, firstcal, sulfate_pi)
1115
1116 CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &
1117 aerindex)
1118 ELSE
1119 tau_ae = 0.
1120 piz_ae = 0.
1121 cg_ae = 0.
1122 ENDIF
1123
1124 ! Param\`etres optiques des nuages et quelques param\`etres pour
1125 ! diagnostics :
1126 if (ok_newmicro) then
1127 CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1128 cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1129 sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
1130 else
1131 CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1132 cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
1133 bl95_b1, cldtaupi, re, fl)
1134 endif
1135
1136 IF (MOD(itap - 1, radpas) == 0) THEN
1137 ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1138 ! Calcul de l'abedo moyen par maille
1139 albsol = sum(falbe * pctsrf, dim = 2)
1140
1141 ! Rayonnement (compatible Arpege-IFS) :
1142 CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &
1143 q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
1144 radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
1145 toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
1146 swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &
1147 solswad, cldtaupi, topswai, solswai)
1148 ENDIF
1149
1150 ! Ajouter la tendance des rayonnements (tous les pas)
1151
1152 DO k = 1, llm
1153 DO i = 1, klon
1154 t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys/86400.
1155 ENDDO
1156 ENDDO
1157
1158 IF (if_ebil >= 2) THEN
1159 tit = 'after rad'
1160 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1161 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1162 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1163 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
1164 END IF
1165
1166 ! Calculer l'hydrologie de la surface
1167 DO i = 1, klon
1168 zxqsurf(i) = 0.
1169 zxsnow(i) = 0.
1170 ENDDO
1171 DO nsrf = 1, nbsrf
1172 DO i = 1, klon
1173 zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)
1174 zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)
1175 ENDDO
1176 ENDDO
1177
1178 ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1179
1180 DO i = 1, klon
1181 bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1182 ENDDO
1183
1184 ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1185
1186 IF (ok_orodr) THEN
1187 ! S\'election des points pour lesquels le sch\'ema est actif :
1188 igwd = 0
1189 DO i = 1, klon
1190 itest(i) = 0
1191 IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
1192 itest(i) = 1
1193 igwd = igwd + 1
1194 ENDIF
1195 ENDDO
1196
1197 CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1198 zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1199 zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1200
1201 ! ajout des tendances
1202 DO k = 1, llm
1203 DO i = 1, klon
1204 t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
1205 u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
1206 v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
1207 ENDDO
1208 ENDDO
1209 ENDIF
1210
1211 IF (ok_orolf) THEN
1212 ! S\'election des points pour lesquels le sch\'ema est actif :
1213 igwd = 0
1214 DO i = 1, klon
1215 itest(i) = 0
1216 IF (zpic(i) - zmea(i) > 100.) THEN
1217 itest(i) = 1
1218 igwd = igwd + 1
1219 ENDIF
1220 ENDDO
1221
1222 CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
1223 itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &
1224 d_t_lif, d_u_lif, d_v_lif)
1225
1226 ! Ajout des tendances :
1227 DO k = 1, llm
1228 DO i = 1, klon
1229 t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
1230 u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
1231 v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
1232 ENDDO
1233 ENDDO
1234 ENDIF
1235
1236 ! Stress n\'ecessaires : toute la physique
1237
1238 DO i = 1, klon
1239 zustrph(i) = 0.
1240 zvstrph(i) = 0.
1241 ENDDO
1242 DO k = 1, llm
1243 DO i = 1, klon
1244 zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1245 * zmasse(i, k)
1246 zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1247 * zmasse(i, k)
1248 ENDDO
1249 ENDDO
1250
1251 CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1252 zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1253
1254 IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1255 2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1256 d_qt, d_ec)
1257
1258 ! Calcul des tendances traceurs
1259 call phytrac(lmt_pas, julien, time, firstcal, lafin, dtphys, t, paprs, &
1260 play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, &
1261 yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
1262 tr_seri, zmasse, ncid_startphy)
1263
1264 IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &
1265 pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &
1266 frac_impa, frac_nucl, pphis, airephy, dtphys)
1267
1268 ! Calculer le transport de l'eau et de l'energie (diagnostique)
1269 CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
1270
1271 ! diag. bilKP
1272
1273 CALL transp_lay(paprs, t_seri, q_seri, u_seri, v_seri, zphi, &
1274 ve_lay, vq_lay, ue_lay, uq_lay)
1275
1276 ! Accumuler les variables a stocker dans les fichiers histoire:
1277
1278 ! conversion Ec -> E thermique
1279 DO k = 1, llm
1280 DO i = 1, klon
1281 ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))
1282 d_t_ec(i, k) = 0.5 / ZRCPD &
1283 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
1284 t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
1285 d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1286 END DO
1287 END DO
1288
1289 IF (if_ebil >= 1) THEN
1290 tit = 'after physic'
1291 CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1292 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1293 ! Comme les tendances de la physique sont ajoute dans la dynamique,
1294 ! on devrait avoir que la variation d'entalpie par la dynamique
1295 ! est egale a la variation de la physique au pas de temps precedent.
1296 ! Donc la somme de ces 2 variations devrait etre nulle.
1297 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1298 evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec)
1299 d_h_vcol_phy = d_h_vcol
1300 END IF
1301
1302 ! SORTIES
1303
1304 ! prw = eau precipitable
1305 DO i = 1, klon
1306 prw(i) = 0.
1307 DO k = 1, llm
1308 prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)
1309 ENDDO
1310 ENDDO
1311
1312 ! Convertir les incrementations en tendances
1313
1314 DO k = 1, llm
1315 DO i = 1, klon
1316 d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1317 d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1318 d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1319 d_qx(i, k, ivap) = (q_seri(i, k) - qx(i, k, ivap)) / dtphys
1320 d_qx(i, k, iliq) = (ql_seri(i, k) - qx(i, k, iliq)) / dtphys
1321 ENDDO
1322 ENDDO
1323
1324 DO iq = 3, nqmx
1325 DO k = 1, llm
1326 DO i = 1, klon
1327 d_qx(i, k, iq) = (tr_seri(i, k, iq - 2) - qx(i, k, iq)) / dtphys
1328 ENDDO
1329 ENDDO
1330 ENDDO
1331
1332 ! Sauvegarder les valeurs de t et q a la fin de la physique:
1333 DO k = 1, llm
1334 DO i = 1, klon
1335 t_ancien(i, k) = t_seri(i, k)
1336 q_ancien(i, k) = q_seri(i, k)
1337 ENDDO
1338 ENDDO
1339
1340 CALL histwrite_phy("phis", pphis)
1341 CALL histwrite_phy("aire", airephy)
1342 CALL histwrite_phy("psol", paprs(:, 1))
1343 CALL histwrite_phy("precip", rain_fall + snow_fall)
1344 CALL histwrite_phy("plul", rain_lsc + snow_lsc)
1345 CALL histwrite_phy("pluc", rain_con + snow_con)
1346 CALL histwrite_phy("tsol", zxtsol)
1347 CALL histwrite_phy("t2m", zt2m)
1348 CALL histwrite_phy("q2m", zq2m)
1349 CALL histwrite_phy("u10m", zu10m)
1350 CALL histwrite_phy("v10m", zv10m)
1351 CALL histwrite_phy("snow", snow_fall)
1352 CALL histwrite_phy("cdrm", cdragm)
1353 CALL histwrite_phy("cdrh", cdragh)
1354 CALL histwrite_phy("topl", toplw)
1355 CALL histwrite_phy("evap", evap)
1356 CALL histwrite_phy("sols", solsw)
1357 CALL histwrite_phy("soll", sollw)
1358 CALL histwrite_phy("solldown", sollwdown)
1359 CALL histwrite_phy("bils", bils)
1360 CALL histwrite_phy("sens", - sens)
1361 CALL histwrite_phy("fder", fder)
1362 CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
1363 CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1364 CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1365 CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
1366
1367 DO nsrf = 1, nbsrf
1368 CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf)*100.)
1369 CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1370 CALL histwrite_phy("sens_"//clnsurf(nsrf), fluxt(:, 1, nsrf))
1371 CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1372 CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1373 CALL histwrite_phy("taux_"//clnsurf(nsrf), fluxu(:, 1, nsrf))
1374 CALL histwrite_phy("tauy_"//clnsurf(nsrf), fluxv(:, 1, nsrf))
1375 CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1376 CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1377 END DO
1378
1379 CALL histwrite_phy("albs", albsol)
1380 CALL histwrite_phy("rugs", zxrugs)
1381 CALL histwrite_phy("s_pblh", s_pblh)
1382 CALL histwrite_phy("s_pblt", s_pblt)
1383 CALL histwrite_phy("s_lcl", s_lcl)
1384 CALL histwrite_phy("s_capCL", s_capCL)
1385 CALL histwrite_phy("s_oliqCL", s_oliqCL)
1386 CALL histwrite_phy("s_cteiCL", s_cteiCL)
1387 CALL histwrite_phy("s_therm", s_therm)
1388 CALL histwrite_phy("s_trmb1", s_trmb1)
1389 CALL histwrite_phy("s_trmb2", s_trmb2)
1390 CALL histwrite_phy("s_trmb3", s_trmb3)
1391 if (conv_emanuel) CALL histwrite_phy("ptop", ema_pct)
1392 CALL histwrite_phy("temp", t_seri)
1393 CALL histwrite_phy("vitu", u_seri)
1394 CALL histwrite_phy("vitv", v_seri)
1395 CALL histwrite_phy("geop", zphi)
1396 CALL histwrite_phy("pres", play)
1397 CALL histwrite_phy("dtvdf", d_t_vdf)
1398 CALL histwrite_phy("dqvdf", d_q_vdf)
1399 CALL histwrite_phy("rhum", zx_rh)
1400
1401 if (ok_instan) call histsync(nid_ins)
1402
1403 IF (lafin) then
1404 call NF95_CLOSE(ncid_startphy)
1405 CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1406 fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1407 radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1408 t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1409 w01)
1410 end IF
1411
1412 firstcal = .FALSE.
1413
1414 END SUBROUTINE physiq
1415
1416 end module physiq_m

  ViewVC Help
Powered by ViewVC 1.1.21