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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 38 - (show annotations)
Thu Jan 6 17:52:19 2011 UTC (13 years, 4 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 75225 byte(s)
Extracted ASCII art from "inigeom" into a separate text file in the
documentation.

"test_disvert" now creates a separate file for layer thicknesses.

Moved variables from module "yomcst" to module "suphec_m" because this
is where those variables are defined. Kept in "yomcst" only parameters
of Earth orbit. Gave the attribute "parameter" to some variables of
module "suphec_m".

Variables of module "yoethf" were defined in procedure "suphec". Moved
these definitions to a new procedure "yoethf" in module "yoethf_m".

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

  ViewVC Help
Powered by ViewVC 1.1.21