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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 16 - (show annotations)
Fri Aug 1 15:37:00 2008 UTC (15 years, 9 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 78217 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.21