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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 98 - (show annotations)
Tue May 13 17:23:16 2014 UTC (10 years ago) by guez
File size: 62705 byte(s)
Split inter_barxy.f : one procedure per module, one module per
file. Grouped the files into a directory.

Split orbite.f.

Value of raz_date read from the namelist is taken into account
(resetting the step counter) even if annee_ref == anneeref and day_ref
== dayref. raz_date is no longer modified by gcm main unit. (Following
LMDZ.)

Removed argument klon of interfsur_lim. Renamed arguments lmt_alb,
lmt_rug to alb_new, z0_new (same name as corresponding actual
arguments in interfsurf_hq).

Removed argument klon of interfsurf_hq.

Removed arguments qs and d_qs of diagetpq. Were always
zero. Downgraded arguments d_qw, d_ql of diagetpq to local variables,
they were not used in physiq. Removed all computations for solid water
in diagetpq, was just zero.


Downgraded arguments fs_bound, fq_bound of diagphy to local variables,
they were not used in physiq. Encapsulated in a test on iprt all
computations in diagphy.

Removed parameter nbtr of module dimphy. Replaced it everywhere in the
program by nqmx - 2.

Removed parameter rnpb of procedure physiq. Kept the true case in
physiq and phytrac. Could not work with false case anyway.

Removed arguments klon, llm, airephy of qcheck. Removed argument ftsol
of initrrnpb, was not used.

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

  ViewVC Help
Powered by ViewVC 1.1.21