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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 47 - (show annotations)
Fri Jul 1 15:00:48 2011 UTC (12 years, 10 months ago) by guez
File size: 72649 byte(s)
Split "thermcell.f" and "cv3_routines.f".
Removed copies of files that are now in "L_util".
Moved "mva9" and "diagetpq" to their own files.
Unified variable names across procedures.

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

  ViewVC Help
Powered by ViewVC 1.1.21