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

Contents of /trunk/Sources/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 205 - (show annotations)
Tue Jun 21 15:16:03 2016 UTC (7 years, 10 months ago) by guez
File size: 42581 byte(s)
dnwd0 is just - mp. Compute it simply in concvl.

da, phi and mp were set to 0 in physiq before the call to
concvl. Clearer to set da1, phi1 and mp1 to 0 in cv_driver so they are
intent out.

qcheck was debugging, printed to standard output and was called
several times per time step of physics.

zxtsol was a duplicate of ztsol.

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

  ViewVC Help
Powered by ViewVC 1.1.21