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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 79 - (show annotations)
Fri Feb 28 17:52:47 2014 UTC (10 years, 2 months ago) by guez
Original Path: trunk/phylmd/physiq.f90
File size: 68555 byte(s)
Moved procedure iniconst inside module comconst. Removed useless
variables of module comconst: im, jm, lllm, imp1, jmp1, lllmm1,
lllmp1, lcl, cotot, unsim. Move definition of dtvr that was in
dynetat0 and etat0 to iniconst. Moved comparison of dtvr from day_step
and start.nc that was in gcm to dynetat0. Moved call to disvert out of
iniconst. Moved call to iniconst in gcm before call to dynetat0.

Removed unused argument pvteta of physiq (not used either in LMDZ).

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

  ViewVC Help
Powered by ViewVC 1.1.21