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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 62 - (show annotations)
Thu Jul 26 14:37:37 2012 UTC (11 years, 9 months ago) by guez
File size: 70376 byte(s)
Changed handling of compiler in compilation system.

Removed the prefix letters "y", "p", "t" or "z" in some names of variables.

Replaced calls to NetCDF by calls to NetCDF95.

Extracted "ioget_calendar" procedures from "calendar.f90" into a
separate file.

Extracted to a separate file, "mathop2.f90", procedures that were not
part of the generic interface "mathop" in "mathop.f90".

Removed computation of "dq" in "bilan_dyn", which was not used.

In "iniadvtrac", removed schemes 20 Slopes and 30 Prather. Was not
compatible with declarations of array sizes.

In "clcdrag", "ustarhb", "vdif_kcay", "yamada4" and "coefkz", changed
the size of some arrays from "klon" to "knon".

Removed possible call to "conema3" in "physiq".

Removed unused argument "cd" in "yamada".

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

  ViewVC Help
Powered by ViewVC 1.1.21