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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 56 - (show annotations)
Tue Jan 10 19:02:02 2012 UTC (12 years, 4 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 71712 byte(s)
Imported "writehist.f" from LMDZ.

Moved module variable "histaveid" from "com_io_dyn" to "initdynav_m".

In "inithist", access directly module variables from "com_io_dyn"
instead of going through the arguments. Copying from LMDZ, write "u"
and scalar variables to separate files. Create a new variable for the
new file in "com_io_dyn". Copying from LMDZ, change the vertical axes
of the three files.

Removed some useless initializations in "dissip".

In "bilan_dyn", removed useless variable "time". Avoiding the
approximate test on "dt_cum" being a multiple of "dt_app", just
compute "ncum" from known usage of "bilan_dyn" and compute "dt_cum"
from "ncum". Change "periodav" from real to integer in
"conf_gcm_m". Since "day_step" is required to be a multiple of
"iperiod", so is "ncum".

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

  ViewVC Help
Powered by ViewVC 1.1.21