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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 76 - (show annotations)
Fri Nov 15 18:45:49 2013 UTC (10 years, 6 months ago) by guez
Original Path: trunk/phylmd/physiq.f90
File size: 68698 byte(s)
Moved everything out of libf.
1 module physiq_m
2
3 IMPLICIT none
4
5 contains
6
7 SUBROUTINE physiq(lafin, rdayvrai, time, dtphys, paprs, play, pphi, pphis, &
8 u, v, t, qx, omega, d_u, d_v, d_t, d_qx, d_ps, dudyn, PVteta)
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 calendar, ONLY: ymds2ju
22 use calltherm_m, only: calltherm
23 USE clesphys, ONLY: cdhmax, cdmmax, co2_ppm, ecrit_hf, ecrit_ins, &
24 ecrit_mth, ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin
25 USE clesphys2, ONLY: cycle_diurne, iflag_con, nbapp_rad, new_oliq, &
26 ok_orodr, ok_orolf, soil_model
27 USE clmain_m, ONLY: clmain
28 use clouds_gno_m, only: clouds_gno
29 USE comgeomphy, ONLY: airephy, cuphy, cvphy
30 USE concvl_m, ONLY: concvl
31 USE conf_gcm_m, ONLY: offline, raz_date
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: iim, jjm, llm, nqmx
39 USE dimphy, ONLY: klon, nbtr
40 USE dimsoil, ONLY: nsoilmx
41 use drag_noro_m, only: drag_noro
42 USE fcttre, ONLY: foeew, qsatl, qsats, thermcep
43 use fisrtilp_m, only: fisrtilp
44 USE hgardfou_m, ONLY: hgardfou
45 USE histsync_m, ONLY: histsync
46 USE histwrite_m, ONLY: histwrite
47 USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &
48 nbsrf
49 USE ini_histhf_m, ONLY: ini_histhf
50 USE ini_histday_m, ONLY: ini_histday
51 USE ini_histins_m, ONLY: ini_histins
52 use newmicro_m, only: newmicro
53 USE oasis_m, ONLY: ok_oasis
54 USE orbite_m, ONLY: orbite, zenang
55 USE ozonecm_m, ONLY: ozonecm
56 USE phyetat0_m, ONLY: phyetat0, rlat, rlon
57 USE phyredem_m, ONLY: phyredem
58 USE phystokenc_m, ONLY: phystokenc
59 USE phytrac_m, ONLY: phytrac
60 USE qcheck_m, ONLY: qcheck
61 use radlwsw_m, only: radlwsw
62 use readsulfate_m, only: readsulfate
63 use sugwd_m, only: sugwd
64 USE suphec_m, ONLY: ra, rcpd, retv, rg, rlvtt, romega, rsigma, rtt
65 USE temps, ONLY: annee_ref, day_ref, itau_phy
66 use unit_nml_m, only: unit_nml
67 USE yoethf_m, ONLY: r2es, rvtmp2
68
69 ! Arguments:
70
71 REAL, intent(in):: rdayvrai
72 ! (elapsed time since January 1st 0h of the starting year, in days)
73
74 REAL, intent(in):: time ! heure de la journée en fraction de jour
75 REAL, intent(in):: dtphys ! pas d'integration pour la physique (seconde)
76 logical, intent(in):: lafin ! dernier passage
77
78 REAL, intent(in):: paprs(klon, llm + 1)
79 ! (pression pour chaque inter-couche, en Pa)
80
81 REAL, intent(in):: play(klon, llm)
82 ! (input pression pour le mileu de chaque couche (en Pa))
83
84 REAL, intent(in):: pphi(klon, llm)
85 ! (input geopotentiel de chaque couche (g z) (reference sol))
86
87 REAL, intent(in):: pphis(klon) ! input geopotentiel du sol
88
89 REAL, intent(in):: u(klon, llm)
90 ! vitesse dans la direction X (de O a E) en m/s
91
92 REAL, intent(in):: v(klon, llm) ! vitesse Y (de S a N) en m/s
93 REAL, intent(in):: t(klon, llm) ! input temperature (K)
94
95 REAL, intent(in):: qx(klon, llm, nqmx)
96 ! (humidité spécifique et fractions massiques des autres traceurs)
97
98 REAL omega(klon, llm) ! input vitesse verticale en Pa/s
99 REAL, intent(out):: d_u(klon, llm) ! tendance physique de "u" (m/s/s)
100 REAL, intent(out):: d_v(klon, llm) ! tendance physique de "v" (m/s/s)
101 REAL, intent(out):: d_t(klon, llm) ! tendance physique de "t" (K/s)
102 REAL d_qx(klon, llm, nqmx) ! output tendance physique de "qx" (kg/kg/s)
103 REAL d_ps(klon) ! output tendance physique de la pression au sol
104
105 LOGICAL:: firstcal = .true.
106
107 INTEGER nbteta
108 PARAMETER(nbteta = 3)
109
110 REAL PVteta(klon, nbteta)
111 ! (output vorticite potentielle a des thetas constantes)
112
113 LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
114 PARAMETER (ok_gust = .FALSE.)
115
116 LOGICAL check ! Verifier la conservation du modele en eau
117 PARAMETER (check = .FALSE.)
118
119 LOGICAL, PARAMETER:: ok_stratus = .FALSE.
120 ! Ajouter artificiellement les stratus
121
122 ! Parametres lies au coupleur OASIS:
123 INTEGER, SAVE:: npas, nexca
124 logical rnpb
125 parameter(rnpb = .true.)
126
127 character(len = 6):: ocean = 'force '
128 ! (type de modèle océan à utiliser: "force" ou "slab" mais pas "couple")
129
130 ! "slab" ocean
131 REAL, save:: tslab(klon) ! temperature of ocean slab
132 REAL, save:: seaice(klon) ! glace de mer (kg/m2)
133 REAL fluxo(klon) ! flux turbulents ocean-glace de mer
134 REAL fluxg(klon) ! flux turbulents ocean-atmosphere
135
136 ! Modele thermique du sol, a activer pour le cycle diurne:
137 logical:: ok_veget = .false. ! type de modele de vegetation utilise
138
139 logical:: ok_journe = .false., ok_mensuel = .true., ok_instan = .false.
140 ! sorties journalieres, mensuelles et instantanees dans les
141 ! fichiers histday, histmth et histins
142
143 LOGICAL ok_region ! sortir le fichier regional
144 PARAMETER (ok_region = .FALSE.)
145
146 ! pour phsystoke avec thermiques
147 REAL fm_therm(klon, llm + 1)
148 REAL entr_therm(klon, llm)
149 real, save:: q2(klon, llm + 1, nbsrf)
150
151 INTEGER ivap ! indice de traceurs pour vapeur d'eau
152 PARAMETER (ivap = 1)
153 INTEGER iliq ! indice de traceurs pour eau liquide
154 PARAMETER (iliq = 2)
155
156 REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
157 LOGICAL, save:: ancien_ok
158
159 REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)
160 REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg/kg/s)
161
162 real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
163
164 !IM Amip2 PV a theta constante
165
166 CHARACTER(LEN = 3) ctetaSTD(nbteta)
167 DATA ctetaSTD/'350', '380', '405'/
168 REAL rtetaSTD(nbteta)
169 DATA rtetaSTD/350., 380., 405./
170
171 !MI Amip2 PV a theta constante
172
173 REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)
174 REAL swup0(klon, llm + 1), swup(klon, llm + 1)
175 SAVE swdn0, swdn, swup0, swup
176
177 REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
178 REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)
179 SAVE lwdn0, lwdn, lwup0, lwup
180
181 !IM Amip2
182 ! variables a une pression donnee
183
184 integer nlevSTD
185 PARAMETER(nlevSTD = 17)
186 real rlevSTD(nlevSTD)
187 DATA rlevSTD/100000., 92500., 85000., 70000., &
188 60000., 50000., 40000., 30000., 25000., 20000., &
189 15000., 10000., 7000., 5000., 3000., 2000., 1000./
190 CHARACTER(LEN = 4) clevSTD(nlevSTD)
191 DATA clevSTD/'1000', '925 ', '850 ', '700 ', '600 ', &
192 '500 ', '400 ', '300 ', '250 ', '200 ', '150 ', '100 ', &
193 '70 ', '50 ', '30 ', '20 ', '10 '/
194
195 ! prw: precipitable water
196 real prw(klon)
197
198 ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
199 ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
200 REAL flwp(klon), fiwp(klon)
201 REAL flwc(klon, llm), fiwc(klon, llm)
202
203 INTEGER kmax, lmax
204 PARAMETER(kmax = 8, lmax = 8)
205 INTEGER kmaxm1, lmaxm1
206 PARAMETER(kmaxm1 = kmax-1, lmaxm1 = lmax-1)
207
208 REAL zx_tau(kmaxm1), zx_pc(lmaxm1)
209 DATA zx_tau/0., 0.3, 1.3, 3.6, 9.4, 23., 60./
210 DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
211
212 ! cldtopres pression au sommet des nuages
213 REAL cldtopres(lmaxm1)
214 DATA cldtopres/50., 180., 310., 440., 560., 680., 800./
215
216 ! taulev: numero du niveau de tau dans les sorties ISCCP
217 CHARACTER(LEN = 4) taulev(kmaxm1)
218
219 DATA taulev/'tau0', 'tau1', 'tau2', 'tau3', 'tau4', 'tau5', 'tau6'/
220 CHARACTER(LEN = 3) pclev(lmaxm1)
221 DATA pclev/'pc1', 'pc2', 'pc3', 'pc4', 'pc5', 'pc6', 'pc7'/
222
223 CHARACTER(LEN = 28) cnameisccp(lmaxm1, kmaxm1)
224 DATA cnameisccp/'pc< 50hPa, tau< 0.3', 'pc= 50-180hPa, tau< 0.3', &
225 'pc= 180-310hPa, tau< 0.3', 'pc= 310-440hPa, tau< 0.3', &
226 'pc= 440-560hPa, tau< 0.3', 'pc= 560-680hPa, tau< 0.3', &
227 'pc= 680-800hPa, tau< 0.3', 'pc< 50hPa, tau= 0.3-1.3', &
228 'pc= 50-180hPa, tau= 0.3-1.3', 'pc= 180-310hPa, tau= 0.3-1.3', &
229 'pc= 310-440hPa, tau= 0.3-1.3', 'pc= 440-560hPa, tau= 0.3-1.3', &
230 'pc= 560-680hPa, tau= 0.3-1.3', 'pc= 680-800hPa, tau= 0.3-1.3', &
231 'pc< 50hPa, tau= 1.3-3.6', 'pc= 50-180hPa, tau= 1.3-3.6', &
232 'pc= 180-310hPa, tau= 1.3-3.6', 'pc= 310-440hPa, tau= 1.3-3.6', &
233 'pc= 440-560hPa, tau= 1.3-3.6', 'pc= 560-680hPa, tau= 1.3-3.6', &
234 'pc= 680-800hPa, tau= 1.3-3.6', 'pc< 50hPa, tau= 3.6-9.4', &
235 'pc= 50-180hPa, tau= 3.6-9.4', 'pc= 180-310hPa, tau= 3.6-9.4', &
236 'pc= 310-440hPa, tau= 3.6-9.4', 'pc= 440-560hPa, tau= 3.6-9.4', &
237 'pc= 560-680hPa, tau= 3.6-9.4', 'pc= 680-800hPa, tau= 3.6-9.4', &
238 'pc< 50hPa, tau= 9.4-23', 'pc= 50-180hPa, tau= 9.4-23', &
239 'pc= 180-310hPa, tau= 9.4-23', 'pc= 310-440hPa, tau= 9.4-23', &
240 'pc= 440-560hPa, tau= 9.4-23', 'pc= 560-680hPa, tau= 9.4-23', &
241 'pc= 680-800hPa, tau= 9.4-23', 'pc< 50hPa, tau= 23-60', &
242 'pc= 50-180hPa, tau= 23-60', 'pc= 180-310hPa, tau= 23-60', &
243 'pc= 310-440hPa, tau= 23-60', 'pc= 440-560hPa, tau= 23-60', &
244 'pc= 560-680hPa, tau= 23-60', 'pc= 680-800hPa, tau= 23-60', &
245 'pc< 50hPa, tau> 60.', 'pc= 50-180hPa, tau> 60.', &
246 'pc= 180-310hPa, tau> 60.', 'pc= 310-440hPa, tau> 60.', &
247 'pc= 440-560hPa, tau> 60.', 'pc= 560-680hPa, tau> 60.', &
248 'pc= 680-800hPa, tau> 60.'/
249
250 !IM ISCCP simulator v3.4
251
252 integer nid_hf, nid_hf3d
253 save nid_hf, nid_hf3d
254
255 ! Variables propres a la physique
256
257 INTEGER, save:: radpas
258 ! (Radiative transfer computations are made every "radpas" call to
259 ! "physiq".)
260
261 REAL radsol(klon)
262 SAVE radsol ! bilan radiatif au sol calcule par code radiatif
263
264 INTEGER, SAVE:: itap ! number of calls to "physiq"
265
266 REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
267
268 REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
269 ! soil temperature of surface fraction
270
271 REAL, save:: fevap(klon, nbsrf) ! evaporation
272 REAL fluxlat(klon, nbsrf)
273 SAVE fluxlat
274
275 REAL fqsurf(klon, nbsrf)
276 SAVE fqsurf ! humidite de l'air au contact de la surface
277
278 REAL, save:: qsol(klon) ! hauteur d'eau dans le sol
279
280 REAL fsnow(klon, nbsrf)
281 SAVE fsnow ! epaisseur neigeuse
282
283 REAL falbe(klon, nbsrf)
284 SAVE falbe ! albedo par type de surface
285 REAL falblw(klon, nbsrf)
286 SAVE falblw ! albedo par type de surface
287
288 ! Paramètres de l'orographie à l'échelle sous-maille (OESM) :
289 REAL, save:: zmea(klon) ! orographie moyenne
290 REAL, save:: zstd(klon) ! deviation standard de l'OESM
291 REAL, save:: zsig(klon) ! pente de l'OESM
292 REAL, save:: zgam(klon) ! anisotropie de l'OESM
293 REAL, save:: zthe(klon) ! orientation de l'OESM
294 REAL, save:: zpic(klon) ! Maximum de l'OESM
295 REAL, save:: zval(klon) ! Minimum de l'OESM
296 REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
297
298 REAL zulow(klon), zvlow(klon)
299
300 INTEGER igwd, idx(klon), itest(klon)
301
302 REAL agesno(klon, nbsrf)
303 SAVE agesno ! age de la neige
304
305 REAL run_off_lic_0(klon)
306 SAVE run_off_lic_0
307 !KE43
308 ! Variables liees a la convection de K. Emanuel (sb):
309
310 REAL bas, top ! cloud base and top levels
311 SAVE bas
312 SAVE top
313
314 REAL Ma(klon, llm) ! undilute upward mass flux
315 SAVE Ma
316 REAL qcondc(klon, llm) ! in-cld water content from convect
317 SAVE qcondc
318 REAL, save:: sig1(klon, llm), w01(klon, llm)
319 REAL, save:: wd(klon)
320
321 ! Variables locales pour la couche limite (al1):
322
323 ! Variables locales:
324
325 REAL cdragh(klon) ! drag coefficient pour T and Q
326 REAL cdragm(klon) ! drag coefficient pour vent
327
328 ! Pour phytrac :
329 REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
330 REAL yu1(klon) ! vents dans la premiere couche U
331 REAL yv1(klon) ! vents dans la premiere couche V
332 REAL ffonte(klon, nbsrf) !Flux thermique utilise pour fondre la neige
333 REAL fqcalving(klon, nbsrf) !Flux d'eau "perdue" par la surface
334 ! !et necessaire pour limiter la
335 ! !hauteur de neige, en kg/m2/s
336 REAL zxffonte(klon), zxfqcalving(klon)
337
338 REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
339 save pfrac_impa
340 REAL pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation
341 save pfrac_nucl
342 REAL pfrac_1nucl(klon, llm)! Produits des coefs lessi nucl (alpha = 1)
343 save pfrac_1nucl
344 REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)
345 REAL frac_nucl(klon, llm) ! idem (nucleation)
346
347 REAL, save:: rain_fall(klon) ! pluie
348 REAL, save:: snow_fall(klon) ! neige
349
350 REAL rain_tiedtke(klon), snow_tiedtke(klon)
351
352 REAL evap(klon), devap(klon) ! evaporation and its derivative
353 REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
354 REAL dlw(klon) ! derivee infra rouge
355 SAVE dlw
356 REAL bils(klon) ! bilan de chaleur au sol
357 REAL fder(klon) ! Derive de flux (sensible et latente)
358 save fder
359 REAL ve(klon) ! integr. verticale du transport meri. de l'energie
360 REAL vq(klon) ! integr. verticale du transport meri. de l'eau
361 REAL ue(klon) ! integr. verticale du transport zonal de l'energie
362 REAL uq(klon) ! integr. verticale du transport zonal de l'eau
363
364 REAL frugs(klon, nbsrf) ! longueur de rugosite
365 save frugs
366 REAL zxrugs(klon) ! longueur de rugosite
367
368 ! Conditions aux limites
369
370 INTEGER julien
371
372 INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
373 REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
374 REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
375
376 REAL albsol(klon)
377 SAVE albsol ! albedo du sol total
378 REAL albsollw(klon)
379 SAVE albsollw ! albedo du sol total
380
381 REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
382
383 ! Declaration des procedures appelees
384
385 EXTERNAL alboc ! calculer l'albedo sur ocean
386 !KE43
387 EXTERNAL conema3 ! convect4.3
388 EXTERNAL nuage ! calculer les proprietes radiatives
389 EXTERNAL transp ! transport total de l'eau et de l'energie
390
391 ! Variables locales
392
393 real, save:: clwcon(klon, llm), rnebcon(klon, llm)
394 real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
395
396 REAL rhcl(klon, llm) ! humiditi relative ciel clair
397 REAL dialiq(klon, llm) ! eau liquide nuageuse
398 REAL diafra(klon, llm) ! fraction nuageuse
399 REAL cldliq(klon, llm) ! eau liquide nuageuse
400 REAL cldfra(klon, llm) ! fraction nuageuse
401 REAL cldtau(klon, llm) ! epaisseur optique
402 REAL cldemi(klon, llm) ! emissivite infrarouge
403
404 REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite
405 REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur
406 REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u
407 REAL fluxv(klon, llm, nbsrf) ! flux turbulent de vitesse v
408
409 REAL zxfluxt(klon, llm)
410 REAL zxfluxq(klon, llm)
411 REAL zxfluxu(klon, llm)
412 REAL zxfluxv(klon, llm)
413
414 ! Le rayonnement n'est pas calculé tous les pas, il faut donc que
415 ! les variables soient rémanentes.
416 REAL, save:: heat(klon, llm) ! chauffage solaire
417 REAL heat0(klon, llm) ! chauffage solaire ciel clair
418 REAL, save:: cool(klon, llm) ! refroidissement infrarouge
419 REAL cool0(klon, llm) ! refroidissement infrarouge ciel clair
420 REAL, save:: topsw(klon), toplw(klon), solsw(klon)
421 REAL, save:: sollw(klon) ! rayonnement infrarouge montant à la surface
422 real, save:: sollwdown(klon) ! downward LW flux at surface
423 REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
424 REAL albpla(klon)
425 REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
426 REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
427 SAVE albpla
428 SAVE heat0, cool0
429
430 INTEGER itaprad
431 SAVE itaprad
432
433 REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
434 REAL conv_t(klon, llm) ! convergence of temperature (K/s)
435
436 REAL cldl(klon), cldm(klon), cldh(klon) !nuages bas, moyen et haut
437 REAL cldt(klon), cldq(klon) !nuage total, eau liquide integree
438
439 REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
440
441 REAL dist, rmu0(klon), fract(klon)
442 REAL zdtime ! pas de temps du rayonnement (s)
443 real zlongi
444 REAL z_avant(klon), z_apres(klon), z_factor(klon)
445 REAL za, zb
446 REAL zx_t, zx_qs, zdelta, zcor
447 real zqsat(klon, llm)
448 INTEGER i, k, iq, nsrf
449 REAL, PARAMETER:: t_coup = 234.
450 REAL zphi(klon, llm)
451
452 !IM cf. AM Variables locales pour la CLA (hbtm2)
453
454 REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
455 REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
456 REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
457 REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
458 REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
459 REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite
460 REAL, SAVE:: therm(klon, nbsrf)
461 REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
462 REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
463 REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
464 ! Grdeurs de sorties
465 REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
466 REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
467 REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
468 REAL s_trmb3(klon)
469
470 ! Variables locales pour la convection de K. Emanuel :
471
472 REAL upwd(klon, llm) ! saturated updraft mass flux
473 REAL dnwd(klon, llm) ! saturated downdraft mass flux
474 REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
475 REAL tvp(klon, llm) ! virtual temp of lifted parcel
476 REAL cape(klon) ! CAPE
477 SAVE cape
478
479 REAL pbase(klon) ! cloud base pressure
480 SAVE pbase
481 REAL bbase(klon) ! cloud base buoyancy
482 SAVE bbase
483 REAL rflag(klon) ! flag fonctionnement de convect
484 INTEGER iflagctrl(klon) ! flag fonctionnement de convect
485 ! -- convect43:
486 REAL dtvpdt1(klon, llm), dtvpdq1(klon, llm)
487 REAL dplcldt(klon), dplcldr(klon)
488
489 ! Variables du changement
490
491 ! con: convection
492 ! lsc: large scale condensation
493 ! ajs: ajustement sec
494 ! eva: évaporation de l'eau liquide nuageuse
495 ! vdf: vertical diffusion in boundary layer
496 REAL d_t_con(klon, llm), d_q_con(klon, llm)
497 REAL d_u_con(klon, llm), d_v_con(klon, llm)
498 REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
499 REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
500 REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
501 REAL rneb(klon, llm)
502
503 REAL mfu(klon, llm), mfd(klon, llm)
504 REAL pen_u(klon, llm), pen_d(klon, llm)
505 REAL pde_u(klon, llm), pde_d(klon, llm)
506 INTEGER kcbot(klon), kctop(klon), kdtop(klon)
507 REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
508 REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
509
510 INTEGER, save:: ibas_con(klon), itop_con(klon)
511
512 REAL rain_con(klon), rain_lsc(klon)
513 REAL snow_con(klon), snow_lsc(klon)
514 REAL d_ts(klon, nbsrf)
515
516 REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
517 REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
518
519 REAL d_u_oro(klon, llm), d_v_oro(klon, llm)
520 REAL d_t_oro(klon, llm)
521 REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
522 REAL d_t_lif(klon, llm)
523
524 REAL, save:: ratqs(klon, llm)
525 real ratqss(klon, llm), ratqsc(klon, llm)
526 real:: ratqsbas = 0.01, ratqshaut = 0.3
527
528 ! Parametres lies au nouveau schema de nuages (SB, PDF)
529 real:: fact_cldcon = 0.375
530 real:: facttemps = 1.e-4
531 logical:: ok_newmicro = .true.
532 real facteur
533
534 integer:: iflag_cldcon = 1
535 logical ptconv(klon, llm)
536
537 ! Variables locales pour effectuer les appels en série :
538
539 REAL t_seri(klon, llm), q_seri(klon, llm)
540 REAL ql_seri(klon, llm), qs_seri(klon, llm)
541 REAL u_seri(klon, llm), v_seri(klon, llm)
542
543 REAL tr_seri(klon, llm, nbtr)
544 REAL d_tr(klon, llm, nbtr)
545
546 REAL zx_rh(klon, llm)
547
548 REAL zustrdr(klon), zvstrdr(klon)
549 REAL zustrli(klon), zvstrli(klon)
550 REAL zustrph(klon), zvstrph(klon)
551 REAL aam, torsfc
552
553 REAL dudyn(iim + 1, jjm + 1, llm)
554
555 REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
556 REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
557
558 INTEGER, SAVE:: nid_day, nid_ins
559
560 REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
561 REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
562 REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
563 REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
564
565 REAL zsto
566
567 logical ok_sync
568 real date0
569
570 ! Variables liées au bilan d'énergie et d'enthalpie :
571 REAL ztsol(klon)
572 REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
573 REAL, SAVE:: d_h_vcol_phy
574 REAL fs_bound, fq_bound
575 REAL zero_v(klon)
576 CHARACTER(LEN = 15) tit
577 INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
578 INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
579
580 REAL d_t_ec(klon, llm) ! tendance due à la conversion Ec -> E thermique
581 REAL ZRCPD
582
583 REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m
584 REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
585 REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille
586 REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes s/1 maille
587
588 ! Aerosol effects:
589
590 REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
591
592 REAL, save:: sulfate_pi(klon, llm)
593 ! SO4 aerosol concentration, in micro g/m3, pre-industrial value
594
595 REAL cldtaupi(klon, llm)
596 ! cloud optical thickness for pre-industrial (pi) aerosols
597
598 REAL re(klon, llm) ! Cloud droplet effective radius
599 REAL fl(klon, llm) ! denominator of re
600
601 ! Aerosol optical properties
602 REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
603 REAL, save:: cg_ae(klon, llm, 2)
604
605 REAL topswad(klon), solswad(klon) ! aerosol direct effect
606 REAL topswai(klon), solswai(klon) ! aerosol indirect effect
607
608 REAL aerindex(klon) ! POLDER aerosol index
609
610 LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
611 LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
612
613 REAL:: bl95_b0 = 2., bl95_b1 = 0.2
614 ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
615 ! B). They link cloud droplet number concentration to aerosol mass
616 ! concentration.
617
618 SAVE u10m
619 SAVE v10m
620 SAVE t2m
621 SAVE q2m
622 SAVE ffonte
623 SAVE fqcalving
624 SAVE rain_con
625 SAVE snow_con
626 SAVE topswai
627 SAVE topswad
628 SAVE solswai
629 SAVE solswad
630 SAVE d_u_con
631 SAVE d_v_con
632
633 real zmasse(klon, llm)
634 ! (column-density of mass of air in a cell, in kg m-2)
635
636 real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
637
638 namelist /physiq_nml/ ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, &
639 fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, ratqsbas, &
640 ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, &
641 nsplit_thermals
642
643 !----------------------------------------------------------------
644
645 IF (if_ebil >= 1) zero_v = 0.
646 ok_sync = .TRUE.
647 IF (nqmx < 2) CALL abort_gcm('physiq', &
648 'eaux vapeur et liquide sont indispensables', 1)
649
650 test_firstcal: IF (firstcal) THEN
651 ! initialiser
652 u10m = 0.
653 v10m = 0.
654 t2m = 0.
655 q2m = 0.
656 ffonte = 0.
657 fqcalving = 0.
658 piz_ae = 0.
659 tau_ae = 0.
660 cg_ae = 0.
661 rain_con(:) = 0.
662 snow_con(:) = 0.
663 topswai(:) = 0.
664 topswad(:) = 0.
665 solswai(:) = 0.
666 solswad(:) = 0.
667
668 d_u_con = 0.
669 d_v_con = 0.
670 rnebcon0 = 0.
671 clwcon0 = 0.
672 rnebcon = 0.
673 clwcon = 0.
674
675 pblh =0. ! Hauteur de couche limite
676 plcl =0. ! Niveau de condensation de la CLA
677 capCL =0. ! CAPE de couche limite
678 oliqCL =0. ! eau_liqu integree de couche limite
679 cteiCL =0. ! cloud top instab. crit. couche limite
680 pblt =0. ! T a la Hauteur de couche limite
681 therm =0.
682 trmb1 =0. ! deep_cape
683 trmb2 =0. ! inhibition
684 trmb3 =0. ! Point Omega
685
686 IF (if_ebil >= 1) d_h_vcol_phy = 0.
687
688 iflag_thermals = 0
689 nsplit_thermals = 1
690 print *, "Enter namelist 'physiq_nml'."
691 read(unit=*, nml=physiq_nml)
692 write(unit_nml, nml=physiq_nml)
693
694 call conf_phys
695
696 ! Initialiser les compteurs:
697
698 frugs = 0.
699 itap = 0
700 itaprad = 0
701 CALL phyetat0("startphy.nc", pctsrf, ftsol, ftsoil, ocean, tslab, &
702 seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, &
703 snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, zmea, &
704 zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &
705 ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
706
707 ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
708 q2 = 1e-8
709
710 radpas = NINT(86400. / dtphys / nbapp_rad)
711
712 ! on remet le calendrier a zero
713 IF (raz_date) itau_phy = 0
714
715 PRINT *, 'cycle_diurne = ', cycle_diurne
716 CALL printflag(radpas, ocean /= 'force', ok_oasis, ok_journe, &
717 ok_instan, ok_region)
718
719 IF (dtphys * REAL(radpas) > 21600. .AND. cycle_diurne) THEN
720 print *, "Au minimum 4 appels par jour si cycle diurne"
721 call abort_gcm('physiq', &
722 "Nombre d'appels au rayonnement insuffisant", 1)
723 ENDIF
724
725 ! Initialisation pour le schéma de convection d'Emanuel :
726 IF (iflag_con >= 3) THEN
727 ibas_con = 1
728 itop_con = 1
729 ENDIF
730
731 IF (ok_orodr) THEN
732 rugoro = MAX(1e-5, zstd * zsig / 2)
733 CALL SUGWD(paprs, play)
734 else
735 rugoro = 0.
736 ENDIF
737
738 lmt_pas = NINT(86400. / dtphys) ! tous les jours
739 print *, 'Number of time steps of "physics" per day: ', lmt_pas
740
741 ecrit_ins = NINT(ecrit_ins/dtphys)
742 ecrit_hf = NINT(ecrit_hf/dtphys)
743 ecrit_mth = NINT(ecrit_mth/dtphys)
744 ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
745 ecrit_reg = NINT(ecrit_reg/dtphys)
746
747 ! Initialiser le couplage si necessaire
748
749 npas = 0
750 nexca = 0
751
752 ! Initialisation des sorties
753
754 call ini_histhf(dtphys, nid_hf, nid_hf3d)
755 call ini_histday(dtphys, ok_journe, nid_day, nqmx)
756 call ini_histins(dtphys, ok_instan, nid_ins)
757 CALL ymds2ju(annee_ref, 1, int(day_ref), 0., date0)
758 ! Positionner date0 pour initialisation de ORCHIDEE
759 print *, 'physiq date0: ', date0
760 ENDIF test_firstcal
761
762 ! Mettre a zero des variables de sortie (pour securite)
763
764 DO i = 1, klon
765 d_ps(i) = 0.
766 ENDDO
767 DO iq = 1, nqmx
768 DO k = 1, llm
769 DO i = 1, klon
770 d_qx(i, k, iq) = 0.
771 ENDDO
772 ENDDO
773 ENDDO
774 da = 0.
775 mp = 0.
776 phi = 0.
777
778 ! Ne pas affecter les valeurs entrées de u, v, h, et q :
779
780 DO k = 1, llm
781 DO i = 1, klon
782 t_seri(i, k) = t(i, k)
783 u_seri(i, k) = u(i, k)
784 v_seri(i, k) = v(i, k)
785 q_seri(i, k) = qx(i, k, ivap)
786 ql_seri(i, k) = qx(i, k, iliq)
787 qs_seri(i, k) = 0.
788 ENDDO
789 ENDDO
790 IF (nqmx >= 3) THEN
791 tr_seri(:, :, :nqmx-2) = qx(:, :, 3:nqmx)
792 ELSE
793 tr_seri(:, :, 1) = 0.
794 ENDIF
795
796 DO i = 1, klon
797 ztsol(i) = 0.
798 ENDDO
799 DO nsrf = 1, nbsrf
800 DO i = 1, klon
801 ztsol(i) = ztsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
802 ENDDO
803 ENDDO
804
805 IF (if_ebil >= 1) THEN
806 tit = 'after dynamics'
807 CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
808 ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
809 d_ql, d_qs, d_ec)
810 ! Comme les tendances de la physique sont ajoutés dans la
811 ! dynamique, la variation d'enthalpie par la dynamique devrait
812 ! être égale à la variation de la physique au pas de temps
813 ! précédent. Donc la somme de ces 2 variations devrait être
814 ! nulle.
815 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
816 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
817 d_qt, 0., fs_bound, fq_bound)
818 END IF
819
820 ! Diagnostic de la tendance dynamique :
821 IF (ancien_ok) THEN
822 DO k = 1, llm
823 DO i = 1, klon
824 d_t_dyn(i, k) = (t_seri(i, k) - t_ancien(i, k)) / dtphys
825 d_q_dyn(i, k) = (q_seri(i, k) - q_ancien(i, k)) / dtphys
826 ENDDO
827 ENDDO
828 ELSE
829 DO k = 1, llm
830 DO i = 1, klon
831 d_t_dyn(i, k) = 0.
832 d_q_dyn(i, k) = 0.
833 ENDDO
834 ENDDO
835 ancien_ok = .TRUE.
836 ENDIF
837
838 ! Ajouter le geopotentiel du sol:
839 DO k = 1, llm
840 DO i = 1, klon
841 zphi(i, k) = pphi(i, k) + pphis(i)
842 ENDDO
843 ENDDO
844
845 ! Check temperatures:
846 CALL hgardfou(t_seri, ftsol)
847
848 ! Incrementer le compteur de la physique
849 itap = itap + 1
850 julien = MOD(NINT(rdayvrai), 360)
851 if (julien == 0) julien = 360
852
853 forall (k = 1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k + 1)) / rg
854
855 ! Mettre en action les conditions aux limites (albedo, sst etc.).
856
857 ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
858 wo = ozonecm(REAL(julien), paprs)
859
860 ! Évaporation de l'eau liquide nuageuse :
861 DO k = 1, llm
862 DO i = 1, klon
863 zb = MAX(0., ql_seri(i, k))
864 t_seri(i, k) = t_seri(i, k) &
865 - zb * RLVTT / RCPD / (1. + RVTMP2 * q_seri(i, k))
866 q_seri(i, k) = q_seri(i, k) + zb
867 ENDDO
868 ENDDO
869 ql_seri = 0.
870
871 IF (if_ebil >= 2) THEN
872 tit = 'after reevap'
873 CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
874 ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
875 d_ql, d_qs, d_ec)
876 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
877 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
878 fs_bound, fq_bound)
879
880 END IF
881
882 ! Appeler la diffusion verticale (programme de couche limite)
883
884 DO i = 1, klon
885 zxrugs(i) = 0.
886 ENDDO
887 DO nsrf = 1, nbsrf
888 DO i = 1, klon
889 frugs(i, nsrf) = MAX(frugs(i, nsrf), 0.000015)
890 ENDDO
891 ENDDO
892 DO nsrf = 1, nbsrf
893 DO i = 1, klon
894 zxrugs(i) = zxrugs(i) + frugs(i, nsrf)*pctsrf(i, nsrf)
895 ENDDO
896 ENDDO
897
898 ! calculs necessaires au calcul de l'albedo dans l'interface
899
900 CALL orbite(REAL(julien), zlongi, dist)
901 IF (cycle_diurne) THEN
902 zdtime = dtphys * REAL(radpas)
903 CALL zenang(zlongi, time, zdtime, rmu0, fract)
904 ELSE
905 rmu0 = -999.999
906 ENDIF
907
908 ! Calcul de l'abedo moyen par maille
909 albsol(:) = 0.
910 albsollw(:) = 0.
911 DO nsrf = 1, nbsrf
912 DO i = 1, klon
913 albsol(i) = albsol(i) + falbe(i, nsrf) * pctsrf(i, nsrf)
914 albsollw(i) = albsollw(i) + falblw(i, nsrf) * pctsrf(i, nsrf)
915 ENDDO
916 ENDDO
917
918 ! Répartition sous maille des flux longwave et shortwave
919 ! Répartition du longwave par sous-surface linéarisée
920
921 DO nsrf = 1, nbsrf
922 DO i = 1, klon
923 fsollw(i, nsrf) = sollw(i) &
924 + 4. * RSIGMA * ztsol(i)**3 * (ztsol(i) - ftsol(i, nsrf))
925 fsolsw(i, nsrf) = solsw(i) * (1. - falbe(i, nsrf)) / (1. - albsol(i))
926 ENDDO
927 ENDDO
928
929 fder = dlw
930
931 ! Couche limite:
932
933 CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, &
934 u_seri, v_seri, julien, rmu0, co2_ppm, ok_veget, ocean, &
935 ftsol, soil_model, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, &
936 qsol, paprs, play, fsnow, fqsurf, fevap, falbe, falblw, fluxlat, &
937 rain_fall, snow_fall, fsolsw, fsollw, fder, rlon, rlat, &
938 frugs, firstcal, agesno, rugoro, d_t_vdf, &
939 d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, &
940 cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, &
941 pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &
942 fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab, seaice)
943
944 ! Incrémentation des flux
945
946 zxfluxt = 0.
947 zxfluxq = 0.
948 zxfluxu = 0.
949 zxfluxv = 0.
950 DO nsrf = 1, nbsrf
951 DO k = 1, llm
952 DO i = 1, klon
953 zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
954 zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
955 zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
956 zxfluxv(i, k) = zxfluxv(i, k) + fluxv(i, k, nsrf) * pctsrf(i, nsrf)
957 END DO
958 END DO
959 END DO
960 DO i = 1, klon
961 sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
962 evap(i) = - zxfluxq(i, 1) ! flux d'évaporation au sol
963 fder(i) = dlw(i) + dsens(i) + devap(i)
964 ENDDO
965
966 DO k = 1, llm
967 DO i = 1, klon
968 t_seri(i, k) = t_seri(i, k) + d_t_vdf(i, k)
969 q_seri(i, k) = q_seri(i, k) + d_q_vdf(i, k)
970 u_seri(i, k) = u_seri(i, k) + d_u_vdf(i, k)
971 v_seri(i, k) = v_seri(i, k) + d_v_vdf(i, k)
972 ENDDO
973 ENDDO
974
975 IF (if_ebil >= 2) THEN
976 tit = 'after clmain'
977 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
978 ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
979 d_ql, d_qs, d_ec)
980 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
981 sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
982 fs_bound, fq_bound)
983 END IF
984
985 ! Update surface temperature:
986
987 DO i = 1, klon
988 zxtsol(i) = 0.
989 zxfluxlat(i) = 0.
990
991 zt2m(i) = 0.
992 zq2m(i) = 0.
993 zu10m(i) = 0.
994 zv10m(i) = 0.
995 zxffonte(i) = 0.
996 zxfqcalving(i) = 0.
997
998 s_pblh(i) = 0.
999 s_lcl(i) = 0.
1000 s_capCL(i) = 0.
1001 s_oliqCL(i) = 0.
1002 s_cteiCL(i) = 0.
1003 s_pblT(i) = 0.
1004 s_therm(i) = 0.
1005 s_trmb1(i) = 0.
1006 s_trmb2(i) = 0.
1007 s_trmb3(i) = 0.
1008
1009 IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &
1010 + pctsrf(i, is_sic) - 1.) > EPSFRA) print *, &
1011 'physiq : problème sous surface au point ', i, pctsrf(i, 1 : nbsrf)
1012 ENDDO
1013 DO nsrf = 1, nbsrf
1014 DO i = 1, klon
1015 ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)
1016 zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
1017 zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)
1018
1019 zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)
1020 zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)
1021 zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)
1022 zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)
1023 zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)
1024 zxfqcalving(i) = zxfqcalving(i) + &
1025 fqcalving(i, nsrf)*pctsrf(i, nsrf)
1026 s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)
1027 s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)
1028 s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)
1029 s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)
1030 s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)
1031 s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)
1032 s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)
1033 s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)
1034 s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)
1035 s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)
1036 ENDDO
1037 ENDDO
1038
1039 ! Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1040
1041 DO nsrf = 1, nbsrf
1042 DO i = 1, klon
1043 IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
1044
1045 IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)
1046 IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)
1047 IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)
1048 IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)
1049 IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)
1050 IF (pctsrf(i, nsrf) < epsfra) &
1051 fqcalving(i, nsrf) = zxfqcalving(i)
1052 IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf) = s_pblh(i)
1053 IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf) = s_lcl(i)
1054 IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf) = s_capCL(i)
1055 IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf) = s_oliqCL(i)
1056 IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf) = s_cteiCL(i)
1057 IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf) = s_pblT(i)
1058 IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf) = s_therm(i)
1059 IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf) = s_trmb1(i)
1060 IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf) = s_trmb2(i)
1061 IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf) = s_trmb3(i)
1062 ENDDO
1063 ENDDO
1064
1065 ! Calculer la derive du flux infrarouge
1066
1067 DO i = 1, klon
1068 dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
1069 ENDDO
1070
1071 ! Appeler la convection (au choix)
1072
1073 DO k = 1, llm
1074 DO i = 1, klon
1075 conv_q(i, k) = d_q_dyn(i, k) + d_q_vdf(i, k)/dtphys
1076 conv_t(i, k) = d_t_dyn(i, k) + d_t_vdf(i, k)/dtphys
1077 ENDDO
1078 ENDDO
1079
1080 IF (check) THEN
1081 za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1082 print *, "avantcon = ", za
1083 ENDIF
1084
1085 if (iflag_con == 2) then
1086 z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
1087 CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:-1), &
1088 q_seri(:, llm:1:-1), conv_t, conv_q, zxfluxq(:, 1), omega, &
1089 d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:-1), &
1090 mfd(:, llm:1:-1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
1091 kdtop, pmflxr, pmflxs)
1092 WHERE (rain_con < 0.) rain_con = 0.
1093 WHERE (snow_con < 0.) snow_con = 0.
1094 ibas_con = llm + 1 - kcbot
1095 itop_con = llm + 1 - kctop
1096 else
1097 ! iflag_con >= 3
1098
1099 CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, &
1100 v_seri, tr_seri, sig1, w01, d_t_con, d_q_con, &
1101 d_u_con, d_v_con, d_tr, rain_con, snow_con, ibas_con, &
1102 itop_con, upwd, dnwd, dnwd0, Ma, cape, tvp, iflagctrl, &
1103 pbase, bbase, dtvpdt1, dtvpdq1, dplcldt, dplcldr, qcondc, &
1104 wd, pmflxr, pmflxs, da, phi, mp, ntra=1)
1105 ! (number of tracers for the convection scheme of Kerry Emanuel:
1106 ! la partie traceurs est faite dans phytrac
1107 ! on met ntra = 1 pour limiter les appels mais on peut
1108 ! supprimer les calculs / ftra.)
1109
1110 clwcon0 = qcondc
1111 mfu = upwd + dnwd
1112 IF (.NOT. ok_gust) wd = 0.
1113
1114 ! Calcul des propriétés des nuages convectifs
1115
1116 DO k = 1, llm
1117 DO i = 1, klon
1118 zx_t = t_seri(i, k)
1119 IF (thermcep) THEN
1120 zdelta = MAX(0., SIGN(1., rtt-zx_t))
1121 zx_qs = r2es * FOEEW(zx_t, zdelta) / play(i, k)
1122 zx_qs = MIN(0.5, zx_qs)
1123 zcor = 1./(1.-retv*zx_qs)
1124 zx_qs = zx_qs*zcor
1125 ELSE
1126 IF (zx_t < t_coup) THEN
1127 zx_qs = qsats(zx_t)/play(i, k)
1128 ELSE
1129 zx_qs = qsatl(zx_t)/play(i, k)
1130 ENDIF
1131 ENDIF
1132 zqsat(i, k) = zx_qs
1133 ENDDO
1134 ENDDO
1135
1136 ! calcul des proprietes des nuages convectifs
1137 clwcon0 = fact_cldcon * clwcon0
1138 call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
1139 rnebcon0)
1140
1141 mfd = 0.
1142 pen_u = 0.
1143 pen_d = 0.
1144 pde_d = 0.
1145 pde_u = 0.
1146 END if
1147
1148 DO k = 1, llm
1149 DO i = 1, klon
1150 t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
1151 q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
1152 u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
1153 v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
1154 ENDDO
1155 ENDDO
1156
1157 IF (if_ebil >= 2) THEN
1158 tit = 'after convect'
1159 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1160 ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1161 d_ql, d_qs, d_ec)
1162 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1163 zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec, &
1164 fs_bound, fq_bound)
1165 END IF
1166
1167 IF (check) THEN
1168 za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1169 print *, "aprescon = ", za
1170 zx_t = 0.
1171 za = 0.
1172 DO i = 1, klon
1173 za = za + airephy(i)/REAL(klon)
1174 zx_t = zx_t + (rain_con(i)+ &
1175 snow_con(i))*airephy(i)/REAL(klon)
1176 ENDDO
1177 zx_t = zx_t/za*dtphys
1178 print *, "Precip = ", zx_t
1179 ENDIF
1180
1181 IF (iflag_con == 2) THEN
1182 z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
1183 z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
1184 DO k = 1, llm
1185 DO i = 1, klon
1186 IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
1187 q_seri(i, k) = q_seri(i, k) * z_factor(i)
1188 ENDIF
1189 ENDDO
1190 ENDDO
1191 ENDIF
1192
1193 ! Convection sèche (thermiques ou ajustement)
1194
1195 d_t_ajs = 0.
1196 d_u_ajs = 0.
1197 d_v_ajs = 0.
1198 d_q_ajs = 0.
1199 fm_therm = 0.
1200 entr_therm = 0.
1201
1202 if (iflag_thermals == 0) then
1203 ! Ajustement sec
1204 CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)
1205 t_seri = t_seri + d_t_ajs
1206 q_seri = q_seri + d_q_ajs
1207 else
1208 ! Thermiques
1209 call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
1210 q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
1211 endif
1212
1213 IF (if_ebil >= 2) THEN
1214 tit = 'after dry_adjust'
1215 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1216 ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1217 d_ql, d_qs, d_ec)
1218 END IF
1219
1220 ! Caclul des ratqs
1221
1222 ! ratqs convectifs à l'ancienne en fonction de (q(z = 0) - q) / q
1223 ! on écrase le tableau ratqsc calculé par clouds_gno
1224 if (iflag_cldcon == 1) then
1225 do k = 1, llm
1226 do i = 1, klon
1227 if(ptconv(i, k)) then
1228 ratqsc(i, k) = ratqsbas + fact_cldcon &
1229 * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
1230 else
1231 ratqsc(i, k) = 0.
1232 endif
1233 enddo
1234 enddo
1235 endif
1236
1237 ! ratqs stables
1238 do k = 1, llm
1239 do i = 1, klon
1240 ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
1241 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
1242 enddo
1243 enddo
1244
1245 ! ratqs final
1246 if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
1247 ! les ratqs sont une conbinaison de ratqss et ratqsc
1248 ! ratqs final
1249 ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1250 ! relaxation des ratqs
1251 ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
1252 ratqs = max(ratqs, ratqsc)
1253 else
1254 ! on ne prend que le ratqs stable pour fisrtilp
1255 ratqs = ratqss
1256 endif
1257
1258 CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
1259 d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
1260 pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
1261 psfl, rhcl)
1262
1263 WHERE (rain_lsc < 0) rain_lsc = 0.
1264 WHERE (snow_lsc < 0) snow_lsc = 0.
1265 DO k = 1, llm
1266 DO i = 1, klon
1267 t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
1268 q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
1269 ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
1270 cldfra(i, k) = rneb(i, k)
1271 IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
1272 ENDDO
1273 ENDDO
1274 IF (check) THEN
1275 za = qcheck(klon, llm, paprs, q_seri, ql_seri, airephy)
1276 print *, "apresilp = ", za
1277 zx_t = 0.
1278 za = 0.
1279 DO i = 1, klon
1280 za = za + airephy(i)/REAL(klon)
1281 zx_t = zx_t + (rain_lsc(i) &
1282 + snow_lsc(i))*airephy(i)/REAL(klon)
1283 ENDDO
1284 zx_t = zx_t/za*dtphys
1285 print *, "Precip = ", zx_t
1286 ENDIF
1287
1288 IF (if_ebil >= 2) THEN
1289 tit = 'after fisrt'
1290 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1291 ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1292 d_ql, d_qs, d_ec)
1293 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1294 zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec, &
1295 fs_bound, fq_bound)
1296 END IF
1297
1298 ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1299
1300 ! 1. NUAGES CONVECTIFS
1301
1302 IF (iflag_cldcon <= -1) THEN
1303 ! seulement pour Tiedtke
1304 snow_tiedtke = 0.
1305 if (iflag_cldcon == -1) then
1306 rain_tiedtke = rain_con
1307 else
1308 rain_tiedtke = 0.
1309 do k = 1, llm
1310 do i = 1, klon
1311 if (d_q_con(i, k) < 0.) then
1312 rain_tiedtke(i) = rain_tiedtke(i)-d_q_con(i, k)/dtphys &
1313 *zmasse(i, k)
1314 endif
1315 enddo
1316 enddo
1317 endif
1318
1319 ! Nuages diagnostiques pour Tiedtke
1320 CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1321 itop_con, diafra, dialiq)
1322 DO k = 1, llm
1323 DO i = 1, klon
1324 IF (diafra(i, k) > cldfra(i, k)) THEN
1325 cldliq(i, k) = dialiq(i, k)
1326 cldfra(i, k) = diafra(i, k)
1327 ENDIF
1328 ENDDO
1329 ENDDO
1330 ELSE IF (iflag_cldcon == 3) THEN
1331 ! On prend pour les nuages convectifs le maximum du calcul de
1332 ! la convection et du calcul du pas de temps précédent diminué
1333 ! d'un facteur facttemps.
1334 facteur = dtphys * facttemps
1335 do k = 1, llm
1336 do i = 1, klon
1337 rnebcon(i, k) = rnebcon(i, k) * facteur
1338 if (rnebcon0(i, k) * clwcon0(i, k) &
1339 > rnebcon(i, k) * clwcon(i, k)) then
1340 rnebcon(i, k) = rnebcon0(i, k)
1341 clwcon(i, k) = clwcon0(i, k)
1342 endif
1343 enddo
1344 enddo
1345
1346 ! On prend la somme des fractions nuageuses et des contenus en eau
1347 cldfra = min(max(cldfra, rnebcon), 1.)
1348 cldliq = cldliq + rnebcon*clwcon
1349 ENDIF
1350
1351 ! 2. Nuages stratiformes
1352
1353 IF (ok_stratus) THEN
1354 CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
1355 DO k = 1, llm
1356 DO i = 1, klon
1357 IF (diafra(i, k) > cldfra(i, k)) THEN
1358 cldliq(i, k) = dialiq(i, k)
1359 cldfra(i, k) = diafra(i, k)
1360 ENDIF
1361 ENDDO
1362 ENDDO
1363 ENDIF
1364
1365 ! Precipitation totale
1366 DO i = 1, klon
1367 rain_fall(i) = rain_con(i) + rain_lsc(i)
1368 snow_fall(i) = snow_con(i) + snow_lsc(i)
1369 ENDDO
1370
1371 IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1372 dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &
1373 d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1374
1375 ! Humidité relative pour diagnostic :
1376 DO k = 1, llm
1377 DO i = 1, klon
1378 zx_t = t_seri(i, k)
1379 IF (thermcep) THEN
1380 zdelta = MAX(0., SIGN(1., rtt-zx_t))
1381 zx_qs = r2es * FOEEW(zx_t, zdelta)/play(i, k)
1382 zx_qs = MIN(0.5, zx_qs)
1383 zcor = 1./(1.-retv*zx_qs)
1384 zx_qs = zx_qs*zcor
1385 ELSE
1386 IF (zx_t < t_coup) THEN
1387 zx_qs = qsats(zx_t)/play(i, k)
1388 ELSE
1389 zx_qs = qsatl(zx_t)/play(i, k)
1390 ENDIF
1391 ENDIF
1392 zx_rh(i, k) = q_seri(i, k)/zx_qs
1393 zqsat(i, k) = zx_qs
1394 ENDDO
1395 ENDDO
1396
1397 ! Introduce the aerosol direct and first indirect radiative forcings:
1398 IF (ok_ade .OR. ok_aie) THEN
1399 ! Get sulfate aerosol distribution :
1400 CALL readsulfate(rdayvrai, firstcal, sulfate)
1401 CALL readsulfate_preind(rdayvrai, firstcal, sulfate_pi)
1402
1403 CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &
1404 aerindex)
1405 ELSE
1406 tau_ae = 0.
1407 piz_ae = 0.
1408 cg_ae = 0.
1409 ENDIF
1410
1411 ! Paramètres optiques des nuages et quelques paramètres pour diagnostics :
1412 if (ok_newmicro) then
1413 CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1414 cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1415 sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
1416 else
1417 CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1418 cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
1419 bl95_b1, cldtaupi, re, fl)
1420 endif
1421
1422 ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1423 IF (MOD(itaprad, radpas) == 0) THEN
1424 DO i = 1, klon
1425 albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &
1426 + falbe(i, is_lic) * pctsrf(i, is_lic) &
1427 + falbe(i, is_ter) * pctsrf(i, is_ter) &
1428 + falbe(i, is_sic) * pctsrf(i, is_sic)
1429 albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &
1430 + falblw(i, is_lic) * pctsrf(i, is_lic) &
1431 + falblw(i, is_ter) * pctsrf(i, is_ter) &
1432 + falblw(i, is_sic) * pctsrf(i, is_sic)
1433 ENDDO
1434 ! Rayonnement (compatible Arpege-IFS) :
1435 CALL radlwsw(dist, rmu0, fract, paprs, play, zxtsol, albsol, &
1436 albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &
1437 heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &
1438 sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &
1439 lwup, swdn0, swdn, swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, &
1440 cg_ae, topswad, solswad, cldtaupi, topswai, solswai)
1441 itaprad = 0
1442 ENDIF
1443 itaprad = itaprad + 1
1444
1445 ! Ajouter la tendance des rayonnements (tous les pas)
1446
1447 DO k = 1, llm
1448 DO i = 1, klon
1449 t_seri(i, k) = t_seri(i, k) + (heat(i, k)-cool(i, k)) * dtphys/86400.
1450 ENDDO
1451 ENDDO
1452
1453 IF (if_ebil >= 2) THEN
1454 tit = 'after rad'
1455 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1456 ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1457 d_ql, d_qs, d_ec)
1458 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1459 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec, &
1460 fs_bound, fq_bound)
1461 END IF
1462
1463 ! Calculer l'hydrologie de la surface
1464 DO i = 1, klon
1465 zxqsurf(i) = 0.
1466 zxsnow(i) = 0.
1467 ENDDO
1468 DO nsrf = 1, nbsrf
1469 DO i = 1, klon
1470 zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)
1471 zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)
1472 ENDDO
1473 ENDDO
1474
1475 ! Calculer le bilan du sol et la dérive de température (couplage)
1476
1477 DO i = 1, klon
1478 bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1479 ENDDO
1480
1481 ! Paramétrisation de l'orographie à l'échelle sous-maille :
1482
1483 IF (ok_orodr) THEN
1484 ! selection des points pour lesquels le shema est actif:
1485 igwd = 0
1486 DO i = 1, klon
1487 itest(i) = 0
1488 IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.)) THEN
1489 itest(i) = 1
1490 igwd = igwd + 1
1491 idx(igwd) = i
1492 ENDIF
1493 ENDDO
1494
1495 CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1496 zthe, zpic, zval, igwd, idx, itest, t_seri, u_seri, v_seri, &
1497 zulow, zvlow, zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1498
1499 ! ajout des tendances
1500 DO k = 1, llm
1501 DO i = 1, klon
1502 t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
1503 u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
1504 v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
1505 ENDDO
1506 ENDDO
1507 ENDIF
1508
1509 IF (ok_orolf) THEN
1510 ! Sélection des points pour lesquels le schéma est actif :
1511 igwd = 0
1512 DO i = 1, klon
1513 itest(i) = 0
1514 IF ((zpic(i) - zmea(i)) > 100.) THEN
1515 itest(i) = 1
1516 igwd = igwd + 1
1517 idx(igwd) = i
1518 ENDIF
1519 ENDDO
1520
1521 CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
1522 itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &
1523 d_t_lif, d_u_lif, d_v_lif)
1524
1525 ! Ajout des tendances :
1526 DO k = 1, llm
1527 DO i = 1, klon
1528 t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
1529 u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
1530 v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
1531 ENDDO
1532 ENDDO
1533 ENDIF
1534
1535 ! Stress nécessaires : toute la physique
1536
1537 DO i = 1, klon
1538 zustrph(i) = 0.
1539 zvstrph(i) = 0.
1540 ENDDO
1541 DO k = 1, llm
1542 DO i = 1, klon
1543 zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1544 * zmasse(i, k)
1545 zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1546 * zmasse(i, k)
1547 ENDDO
1548 ENDDO
1549
1550 CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &
1551 zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1552
1553 IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1554 2, dtphys, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, paprs, &
1555 d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1556
1557 ! Calcul des tendances traceurs
1558 call phytrac(rnpb, itap, lmt_pas, julien, time, firstcal, lafin, nqmx-2, &
1559 dtphys, u, t, paprs, play, mfu, mfd, pen_u, pde_u, pen_d, pde_d, &
1560 ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, pctsrf, frac_impa, &
1561 frac_nucl, pphis, albsol, rhcl, cldfra, rneb, diafra, cldliq, &
1562 pmflxr, pmflxs, prfl, psfl, da, phi, mp, upwd, dnwd, tr_seri, zmasse)
1563
1564 IF (offline) THEN
1565 call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, pde_u, &
1566 pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
1567 pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
1568 ENDIF
1569
1570 ! Calculer le transport de l'eau et de l'energie (diagnostique)
1571 CALL transp(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, &
1572 ue, uq)
1573
1574 ! diag. bilKP
1575
1576 CALL transp_lay(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
1577 ve_lay, vq_lay, ue_lay, uq_lay)
1578
1579 ! Accumuler les variables a stocker dans les fichiers histoire:
1580
1581 ! conversion Ec -> E thermique
1582 DO k = 1, llm
1583 DO i = 1, klon
1584 ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))
1585 d_t_ec(i, k) = 0.5 / ZRCPD &
1586 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
1587 t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
1588 d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1589 END DO
1590 END DO
1591
1592 IF (if_ebil >= 1) THEN
1593 tit = 'after physic'
1594 CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1595 ql_seri, qs_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_qw, &
1596 d_ql, d_qs, d_ec)
1597 ! Comme les tendances de la physique sont ajoute dans la dynamique,
1598 ! on devrait avoir que la variation d'entalpie par la dynamique
1599 ! est egale a la variation de la physique au pas de temps precedent.
1600 ! Donc la somme de ces 2 variations devrait etre nulle.
1601 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1602 evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec, &
1603 fs_bound, fq_bound)
1604
1605 d_h_vcol_phy = d_h_vcol
1606
1607 END IF
1608
1609 ! SORTIES
1610
1611 ! prw = eau precipitable
1612 DO i = 1, klon
1613 prw(i) = 0.
1614 DO k = 1, llm
1615 prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)
1616 ENDDO
1617 ENDDO
1618
1619 ! Convertir les incrementations en tendances
1620
1621 DO k = 1, llm
1622 DO i = 1, klon
1623 d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1624 d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1625 d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1626 d_qx(i, k, ivap) = (q_seri(i, k) - qx(i, k, ivap)) / dtphys
1627 d_qx(i, k, iliq) = (ql_seri(i, k) - qx(i, k, iliq)) / dtphys
1628 ENDDO
1629 ENDDO
1630
1631 IF (nqmx >= 3) THEN
1632 DO iq = 3, nqmx
1633 DO k = 1, llm
1634 DO i = 1, klon
1635 d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / dtphys
1636 ENDDO
1637 ENDDO
1638 ENDDO
1639 ENDIF
1640
1641 ! Sauvegarder les valeurs de t et q a la fin de la physique:
1642 DO k = 1, llm
1643 DO i = 1, klon
1644 t_ancien(i, k) = t_seri(i, k)
1645 q_ancien(i, k) = q_seri(i, k)
1646 ENDDO
1647 ENDDO
1648
1649 ! Ecriture des sorties
1650 call write_histhf
1651 call write_histday
1652 call write_histins
1653
1654 ! Si c'est la fin, il faut conserver l'etat de redemarrage
1655 IF (lafin) THEN
1656 itau_phy = itau_phy + itap
1657 CALL phyredem("restartphy.nc", rlat, rlon, pctsrf, ftsol, ftsoil, &
1658 tslab, seaice, fqsurf, qsol, fsnow, falbe, falblw, fevap, &
1659 rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
1660 agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
1661 q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
1662 ENDIF
1663
1664 firstcal = .FALSE.
1665
1666 contains
1667
1668 subroutine write_histday
1669
1670 use gr_phy_write_3d_m, only: gr_phy_write_3d
1671 integer itau_w ! pas de temps ecriture
1672
1673 !------------------------------------------------
1674
1675 if (ok_journe) THEN
1676 itau_w = itau_phy + itap
1677 if (nqmx <= 4) then
1678 call histwrite(nid_day, "Sigma_O3_Royer", itau_w, &
1679 gr_phy_write_3d(wo) * 1e3)
1680 ! (convert "wo" from kDU to DU)
1681 end if
1682 if (ok_sync) then
1683 call histsync(nid_day)
1684 endif
1685 ENDIF
1686
1687 End subroutine write_histday
1688
1689 !****************************
1690
1691 subroutine write_histhf
1692
1693 ! From phylmd/write_histhf.h, version 1.5 2005/05/25 13:10:09
1694
1695 !------------------------------------------------
1696
1697 call write_histhf3d
1698
1699 IF (ok_sync) THEN
1700 call histsync(nid_hf)
1701 ENDIF
1702
1703 end subroutine write_histhf
1704
1705 !***************************************************************
1706
1707 subroutine write_histins
1708
1709 ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09
1710
1711 real zout
1712 integer itau_w ! pas de temps ecriture
1713
1714 !--------------------------------------------------
1715
1716 IF (ok_instan) THEN
1717 ! Champs 2D:
1718
1719 zsto = dtphys * ecrit_ins
1720 zout = dtphys * ecrit_ins
1721 itau_w = itau_phy + itap
1722
1723 i = NINT(zout/zsto)
1724 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)
1725 CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)
1726
1727 i = NINT(zout/zsto)
1728 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)
1729 CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)
1730
1731 DO i = 1, klon
1732 zx_tmp_fi2d(i) = paprs(i, 1)
1733 ENDDO
1734 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1735 CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d)
1736
1737 DO i = 1, klon
1738 zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
1739 ENDDO
1740 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1741 CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)
1742
1743 DO i = 1, klon
1744 zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
1745 ENDDO
1746 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1747 CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)
1748
1749 DO i = 1, klon
1750 zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
1751 ENDDO
1752 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1753 CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)
1754
1755 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxtsol, zx_tmp_2d)
1756 CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)
1757 !ccIM
1758 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zt2m, zx_tmp_2d)
1759 CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)
1760
1761 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zq2m, zx_tmp_2d)
1762 CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)
1763
1764 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zu10m, zx_tmp_2d)
1765 CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)
1766
1767 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zv10m, zx_tmp_2d)
1768 CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)
1769
1770 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, snow_fall, zx_tmp_2d)
1771 CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)
1772
1773 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragm, zx_tmp_2d)
1774 CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)
1775
1776 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragh, zx_tmp_2d)
1777 CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)
1778
1779 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, toplw, zx_tmp_2d)
1780 CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)
1781
1782 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, evap, zx_tmp_2d)
1783 CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)
1784
1785 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, solsw, zx_tmp_2d)
1786 CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)
1787
1788 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollw, zx_tmp_2d)
1789 CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)
1790
1791 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollwdown, zx_tmp_2d)
1792 CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)
1793
1794 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)
1795 CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)
1796
1797 zx_tmp_fi2d(1:klon) = -1*sens(1:klon)
1798 ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)
1799 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1800 CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)
1801
1802 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, fder, zx_tmp_2d)
1803 CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)
1804
1805 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_oce), zx_tmp_2d)
1806 CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)
1807
1808 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_ter), zx_tmp_2d)
1809 CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)
1810
1811 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_lic), zx_tmp_2d)
1812 CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)
1813
1814 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_sic), zx_tmp_2d)
1815 CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)
1816
1817 DO nsrf = 1, nbsrf
1818 !XXX
1819 zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)*100.
1820 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1821 CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &
1822 zx_tmp_2d)
1823
1824 zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)
1825 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1826 CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &
1827 zx_tmp_2d)
1828
1829 zx_tmp_fi2d(1 : klon) = fluxt(1 : klon, 1, nsrf)
1830 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1831 CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &
1832 zx_tmp_2d)
1833
1834 zx_tmp_fi2d(1 : klon) = fluxlat(1 : klon, nsrf)
1835 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1836 CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &
1837 zx_tmp_2d)
1838
1839 zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, nsrf)
1840 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1841 CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &
1842 zx_tmp_2d)
1843
1844 zx_tmp_fi2d(1 : klon) = fluxu(1 : klon, 1, nsrf)
1845 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1846 CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &
1847 zx_tmp_2d)
1848
1849 zx_tmp_fi2d(1 : klon) = fluxv(1 : klon, 1, nsrf)
1850 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1851 CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &
1852 zx_tmp_2d)
1853
1854 zx_tmp_fi2d(1 : klon) = frugs(1 : klon, nsrf)
1855 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1856 CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &
1857 zx_tmp_2d)
1858
1859 zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)
1860 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1861 CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &
1862 zx_tmp_2d)
1863
1864 END DO
1865 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)
1866 CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)
1867 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsollw, zx_tmp_2d)
1868 CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)
1869
1870 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)
1871 CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)
1872
1873 !HBTM2
1874
1875 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblh, zx_tmp_2d)
1876 CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d)
1877
1878 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblt, zx_tmp_2d)
1879 CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d)
1880
1881 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_lcl, zx_tmp_2d)
1882 CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d)
1883
1884 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_capCL, zx_tmp_2d)
1885 CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d)
1886
1887 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_oliqCL, zx_tmp_2d)
1888 CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d)
1889
1890 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_cteiCL, zx_tmp_2d)
1891 CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d)
1892
1893 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_therm, zx_tmp_2d)
1894 CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d)
1895
1896 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb1, zx_tmp_2d)
1897 CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d)
1898
1899 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb2, zx_tmp_2d)
1900 CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d)
1901
1902 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb3, zx_tmp_2d)
1903 CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)
1904
1905 ! Champs 3D:
1906
1907 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)
1908 CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d)
1909
1910 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)
1911 CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d)
1912
1913 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)
1914 CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d)
1915
1916 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, zphi, zx_tmp_3d)
1917 CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d)
1918
1919 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, play, zx_tmp_3d)
1920 CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d)
1921
1922 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_t_vdf, zx_tmp_3d)
1923 CALL histwrite(nid_ins, "dtvdf", itau_w, zx_tmp_3d)
1924
1925 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)
1926 CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
1927
1928 if (ok_sync) then
1929 call histsync(nid_ins)
1930 endif
1931 ENDIF
1932
1933 end subroutine write_histins
1934
1935 !****************************************************
1936
1937 subroutine write_histhf3d
1938
1939 ! From phylmd/write_histhf3d.h, version 1.2 2005/05/25 13:10:09
1940
1941 integer itau_w ! pas de temps ecriture
1942
1943 !-------------------------------------------------------
1944
1945 itau_w = itau_phy + itap
1946
1947 ! Champs 3D:
1948
1949 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)
1950 CALL histwrite(nid_hf3d, "temp", itau_w, zx_tmp_3d)
1951
1952 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, qx(1, 1, ivap), zx_tmp_3d)
1953 CALL histwrite(nid_hf3d, "ovap", itau_w, zx_tmp_3d)
1954
1955 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)
1956 CALL histwrite(nid_hf3d, "vitu", itau_w, zx_tmp_3d)
1957
1958 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)
1959 CALL histwrite(nid_hf3d, "vitv", itau_w, zx_tmp_3d)
1960
1961 if (nbtr >= 3) then
1962 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, tr_seri(1, 1, 3), &
1963 zx_tmp_3d)
1964 CALL histwrite(nid_hf3d, "O3", itau_w, zx_tmp_3d)
1965 end if
1966
1967 if (ok_sync) then
1968 call histsync(nid_hf3d)
1969 endif
1970
1971 end subroutine write_histhf3d
1972
1973 END SUBROUTINE physiq
1974
1975 end module physiq_m

  ViewVC Help
Powered by ViewVC 1.1.21