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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 7 - (show annotations)
Mon Mar 31 12:24:17 2008 UTC (16 years, 1 month ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 93342 byte(s)
This revision is not in working order. Pending some moving of files.

Important changes. In the program "etat0_lim": ozone coefficients from
Mobidic are regridded in time instead of pressure ; consequences in
"etat0". In the program "gcm", ozone coefficients from Mobidic are
read once per day only for the current day and regridded in pressure ;
consequences in "o3_chem_m", "regr_pr_coefoz", "phytrac" and
"regr_pr_comb_coefoz_m".

NetCDF95 is a library and does not export NetCDF.

New variables "nag_gl_options", "nag_fcalls_options" and
"nag_cross_options" in "nag_tools.mk".

"check_coefoz.jnl" rewritten entirely for new version of
"coefoz_LMDZ.nc".

Target "obj_etat0_lim" moved from "GNUmakefile" to "nag_rules.mk".

Added some "intent" attributes in "calfis", "clmain", "clqh",
"cltrac", "cltracrn", "cvltr", "ini_undefSTD", "moy_undefSTD",
"nflxtr", "phystokenc", "phytrac", "readsulfate", "readsulfate_preind"
and "undefSTD".

In "dynetat0", "dynredem0" and "gcm", "phis" has rank 2 instead of
1. "phis" has assumed shape in "dynredem0".

Added module containing "dynredem0". Changed some calls with NetCDF
Fortran 77 interface to calls with NetCDF95 interface.

Replaced calls to "ssum" by calls to "sum" in "inigeom".

In "make.sh", new option "-c" to change compiler.

In "aaam_bud", argument "rjour" deleted.

In "physiq": renamed some variables; deleted variable "xjour".

In "phytrac": renamed some variables; new argument "lmt_pas".

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

  ViewVC Help
Powered by ViewVC 1.1.21