/[lmdze]/trunk/libf/phylmd/physiq.f90
ViewVC logotype

Contents of /trunk/libf/phylmd/physiq.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 52 - (show annotations)
Fri Sep 23 12:28:01 2011 UTC (12 years, 7 months ago) by guez
File size: 71804 byte(s)
Split "conflx.f" into single-procedure files in directory "Conflx".

Split "cv_routines.f" into single-procedure files in directory
"CV_routines". Made module "cvparam" from included file
"cvparam.h". No included file other than "netcdf.inc" left in LMDZE.

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

  ViewVC Help
Powered by ViewVC 1.1.21