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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 68 - (show annotations)
Wed Nov 14 16:59:30 2012 UTC (11 years, 6 months ago) by guez
Original Path: trunk/libf/phylmd/physiq.f90
File size: 70590 byte(s)
Split "flincom.f90" into "flinclo.f90", "flinfindcood.f90",
"flininfo.f90" and "flinopen_nozoom.f90", in directory
"IOIPSL/Flincom".

Renamed "etat0_lim" to "ce0l", as in LMDZ.

Split "readsulfate.f" into "readsulfate.f90", "readsulfate_preind.f90"
and "getso4fromfile.f90".

In etat0, renamed variable q3d to q, as in "dynredem1". Replaced calls
to Flicom procedures by calls to NetCDF95.

In leapfrog, added call to writehist.

Extracted ASCII art from "grid_noro" into a file
"grid_noro.txt". Transformed explicit-shape local arrays into
automatic arrays, so that test on values of iim and jjm is no longer
needed. Test on weight:
          IF (weight(ii, jj) /= 0.) THEN
is useless. There is already a test before:
    if (any(weight == 0.)) stop "zero weight in grid_noro"

In "aeropt", replaced duplicated lines with different values of inu by
a loop on inu.

Removed arguments of "conf_phys". Corresponding variables are now
defined in "physiq", in a namelist. In "conf_phys", read a namelist
instead of using getin.

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

  ViewVC Help
Powered by ViewVC 1.1.21