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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 227 - (show annotations)
Thu Nov 2 15:47:03 2017 UTC (6 years, 6 months ago) by guez
File size: 38905 byte(s)
Rename phisinit to phis in restart.nc: clearer, same name as Fortran variable.

In aaam_bud, use rlat and rlon from phyetat0_m instead of having these
module variables associated to actual arguments in physiq.

In clmain, too many wind variables make the procedure hard to
understand. Use yu(:knon, 1) and yv(:knon, 1) instead of u1lay(:knon)
and v1lay(:knon). Note that when yu(:knon, 1) and yv(:knon, 1) are
used as actual arguments, they are probably copied to new arrays since
the elements are not contiguous. Rename yu10m to wind10m because this
is the norm of wind vector, not its zonal component. Rename yustar to
ustar. Rename uzon and vmer to u1 and v1 since these are wind
components at first layer and u1 and v1 are the names of corresponding
dummy arguments in stdlevvar.

In clmain, rename yzlev to zlev.

In clmain, screenc, stdlevvar and coefcdrag, remove the code
corresponding to zxli true (not used in LMDZ either).

Subroutine ustarhb becomes a function. Simplifications using the fact
that zx_alf2 = 0 and zx_alf1 = 1 (discarding the possibility to change
this).

In procedure vdif_kcay, remove unused dummy argument plev. Remove
useless computations of sss and sssq.

In clouds_gno, exp(100.) would overflow in single precision. Set
maximum to exp(80.) instead.

In physiq, use u(:, 1) and v(:, 1) as arguments to phytrac instead of
creating ad hoc variables yu1 and yv1.

In stdlevvar, rename dummy argument u_10m to wind10m, following the
corresponding modification in clmain. Simplifications using the fact
that ok_pred = 0 and ok_corr = 1 (discarding the possibility to change
this).

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

  ViewVC Help
Powered by ViewVC 1.1.21