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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 18 - (show annotations)
Thu Aug 7 12:29:13 2008 UTC (15 years, 9 months ago) by guez
File size: 76785 byte(s)
In module "regr_pr", rewrote scanning of horizontal positions as a
single set of loops, using a mask.

Added some "intent" attributes.

In "dynredem0", replaced calls to Fortran 77 interface of NetCDF by
calls to NetCDF95. Removed calls to "nf_redef", regrouped all writing
operations. In "dynredem1", replaced some calls to Fortran 77
interface of NetCDF by calls to Fortran 90 interface.

Renamed variable "nqmax" to "nq_phys".

In "physiq", if "nq >= 5" then "wo" is computed from the
parameterization of "Cariolle".

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

  ViewVC Help
Powered by ViewVC 1.1.21