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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 99 - (show annotations)
Wed Jul 2 18:39:15 2014 UTC (9 years, 10 months ago) by guez
Original Path: trunk/phylmd/physiq.f
File size: 62266 byte(s) 
Created procedure test_disvert (following LMDZ). Added procedures
hybrid and funcd in module disvert_m. Upgraded compute_ab from
internal procedure of disvert to module procedure. Added variables y,
ya in module disvert_m. Upgraded s from local variable of procedure
disvert to module variable.

Renamed allowed value of variable vert_sampling in procedure disvert
from "read" to "read_hybrid". Added possibility to read pressure
values, value "read_pressure". Replaced vertical distribution for
value "param" by the distribution "strato_correct" from LMDZ (but kept
the value "param"). In case "tropo", replaced 1 by dsigmin (following
LMDZ). In case "strato", replaced 0.3 by dsigmin (following LMDZ).

Changed computation of bp in procedure compute_ab.

Removed debugindex case in clmain. Removed useless argument rlon of
procedure clmain. Removed useless variables ytaux, ytauy of procedure
clmain.

Removed intermediary variables tsol, qsol, tsolsrf, tslab in procedure
etat0.

Removed variable ok_veget:. coupling with the model Orchid is not
possible. Removed variable ocean: modeling an ocean slab is not
possible.

Removed useless variables tmp_rriv and tmp_rcoa from module
interface_surf.

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