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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 72 - (show annotations)
Tue Jul 23 13:00:07 2013 UTC (10 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 68785 byte(s)
NaN to signalling NaN in gfortran_debug.mk.

Removed unused procedures in getincom and getincom2. In procedure
conf_interface, replaced call to getincom by new namelist. Moved
procedure conf_interface into module interface_surf.

Added variables sig1 and w01 to startphy.nc and restartphy.nc, for
procedure cv_driver. Renamed (ema_)?work1 and (ema_)?work2 to sig1 and
w01 in concvl and physiq.

Deleted unused arguments of clmain, clqh and intersurf_hq, among which
(y)?sollwdown. Following LMDZ, in physiq, read sollw instead of
sollwdown from startphy.nc, write sollw instead of sollwdown to
restartphy.nc.

In procedure sw, initialized zfs[ud][pn]a[di], for runs where ok_ade
and ok_aie are false. (Following LMDZ.)

Added dimension klev to startphy.nc and restartphy.nc, and deleted
dimension horizon_vertical. Made t_ancien and q_ancien two-dimensional
NetCDF variables. Bug fix: in phyetat0, define ratqs, clwcon and
rnebcon for vertical levels >=2.

Bug fix: set mfg, p[de]n_[ud] to 0. when iflag_con >= 3. (Following LMDZ.)

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

  ViewVC Help
Powered by ViewVC 1.1.21