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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 20 - (show annotations)
Wed Oct 15 16:19:57 2008 UTC (15 years, 6 months ago) by guez
File size: 76646 byte(s)
Deleted argument "presnivs" of "physiq", "ini_histhf", "ini_histhf3d",
"ini_histday", "ini_histins", "ini_histrac", "phytrac". Access it from
"comvert" instead.

Replaced calls to NetCDF Fortran 77 interface by calls to Fortran 90
interface or to NetCDF95.

Procedure "gr_phy_write_3d" now works with a variable of arbitrary
size in the second dimension.

Annotated use statements with "only" clause.

Replaced calls to NetCDF interface version 2 by calls to Fortran 90
interface in "guide.f90" and "read_reanalyse.f".

In "write_histrac", replaced calls to "gr_fi_ecrit" by calls to
"gr_phy_write_2d" and "gr_phy_write_3d".

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

  ViewVC Help
Powered by ViewVC 1.1.21