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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 13 - (show annotations)
Fri Jul 25 19:59:34 2008 UTC (15 years, 10 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 91675 byte(s)
-- Minor change of behaviour:

"etat0" does not compute "rugsrel" nor "radpas". Deleted arguments
"radpas" and "rugsrel" of "phyredem". Deleted argument "rugsrel" of
"phyetat0". "startphy.nc" does not contain the variable "RUGSREL". In
"physiq", "rugoro" is set to 0 if not "ok_orodr". The whole program
"etat0_lim" does not use "clesphys2".

-- Minor modification of input/output:

Created subroutine "read_clesphys2". Variables of "clesphys2" are read
in "read_clesphys2" instead of "conf_gcm". "printflag" does not print
variables of "clesphys2".

-- Should not change any result at run time:

References to module "numer_rec" instead of individual modules of
"Numer_rec_Lionel".

Deleted argument "clesphy0" of "calfis", "physiq", "conf_gcm",
"leapfrog", "phyetat0". Deleted variable "clesphy0" in
"gcm". "phyetat0" does not modify variables of "clesphys2".

The program unit "gcm" does not modify "itau_phy".

Added some "intent" attributes.

"regr11_lint" does not call "polint".

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

  ViewVC Help
Powered by ViewVC 1.1.21