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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 15 - (show annotations)
Fri Aug 1 15:24:12 2008 UTC (15 years, 9 months ago) by guez
File size: 78386 byte(s)
-- Minor modification of input/output:

Added variable "Sigma_O3_Royer" to "histday.nc". "ecrit_day" is not
modified in "physiq". Removed variables "pyu1", "pyv1", "ftsol1",
"ftsol2", "ftsol3", "ftsol4", "psrf1", "psrf2", "psrf3", "psrf4"
"mfu", "mfd", "en_u", "en_d", "de_d", "de_u", "coefh" from
"histrac.nc".

Variable "raz_date" of module "conf_gcm_m" has logical type instead of
integer type.

-- Should not change any result at run time:

Modified calls to "IOIPSL_Lionel" procedures because the interfaces of
these procedures have been simplified.

Changed name of variable in module "start_init_orog_m": "masque" to
"mask".

Created a module containing procedure "phyredem".

Removed arguments "punjours", "pdayref" and "ptimestep" of procedure
"iniphysiq".

Renamed procedure "gr_phy_write" to "gr_phy_write_2d". Created
procedure "gr_phy_write_3d".

Removed procedures "ini_undefstd", "moy_undefSTD", "calcul_STDlev",
"calcul_divers".

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

  ViewVC Help
Powered by ViewVC 1.1.21