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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 51 - (show annotations)
Tue Sep 20 09:14:34 2011 UTC (12 years, 7 months ago) by guez
File size: 72007 byte(s)
Split "getincom.f90" into "getincom.f90" and "getincom2.f90". Split
"nuage.f" into "nuage.f90", "diagcld1.f90" and "diagcld2.f90". Created
module "chem" from included file "chem.h". Moved "YOEGWD.f90" to
directory "Orography".

In "physiq", for evaporation of water, "zlsdcp" was equal to
"zlvdc". Removed useless variables.

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

  ViewVC Help
Powered by ViewVC 1.1.21