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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 54 - (show annotations)
Tue Dec 6 15:07:04 2011 UTC (12 years, 5 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 71789 byte(s)
Removed Numerical Recipes procedure "ran1". Replaced calls to "ran1"
in "inidissip" by calls to intrinsic procedures.

Split file "interface_surf.f90" into a file with a module containing
only variables, "interface_surf", and single-procedure files. Gathered
files into directory "Interface_surf".

Added argument "cdivu" to "gradiv" and "gradiv2", "cdivh" to
"divgrad2" and "divgrad", and "crot" to "nxgraro2" and
"nxgrarot". "dissip" now uses variables "cdivu", "cdivh" and "crot"
from module "inidissip_m", so it can pass them to "gradiv2",
etc. Thanks to this modification, we avoid a circular dependency
betwwen "inidissip.f90" and "gradiv2.f90", etc. The value -1. used by
"gradiv2", for instance, during computation of eigenvalues is not the
value "cdivu" computed by "inidissip".

Extracted procedure "start_inter_3d" from module "startdyn", to its
own module.

In "inidissip", unrolled loop on "ii". I find it clearer now.

Moved variables "matriceun", "matriceus", "matricevn", "matricevs",
"matrinvn" and "matrinvs" from module "parafilt" to module
"inifilr_m". Moved variables "jfiltnu", "jfiltnv", "jfiltsu",
"jfiltsv" from module "coefils" to module "inifilr_m".

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

  ViewVC Help
Powered by ViewVC 1.1.21