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

Contents of /trunk/Sources/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 23 - (show annotations)
Mon Dec 14 15:25:16 2009 UTC (14 years, 4 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 75335 byte(s)
Split "orografi.f": one file for each procedure. Put the created files
in new directory "Orography".

Removed argument "vcov" of procedure "sortvarc". Removed arguments
"itau" and "time" of procedure "caldyn0". Removed arguments "itau",
"time" and "vcov" of procedure "sortvarc0".

Removed argument "time" of procedure "dynredem1". Removed NetCDF
variable "temps" in files "start.nc" and "restart.nc", because its
value is always 0.

Removed argument "nq" of procedures "iniadvtrac" and "leapfrog". The
number of "tracers read in "traceur.def" must now be equal to "nqmx",
or "nqmx" must equal 4 if there is no file "traceur.def". Replaced
variable "nq" by constant "nqmx" in "leapfrog".

NetCDF variable for ozone field in "coefoz.nc" must now be called
"tro3" instead of "r".

Fixed bug in "zenang".

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

  ViewVC Help
Powered by ViewVC 1.1.21