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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 137 - (show annotations)
Wed May 6 15:51:03 2015 UTC (9 years ago) by guez
Original Path: trunk/Sources/phylmd/physiq.f
File size: 61669 byte(s)
Removed unused argument missval in ma_fucoll_r[1-3]1, ma_fufill_r[1-3]1.

Split filtreg into two procedures: filtreg_scal and filtreg_v. I did
not like the test on the extent of the argument and there was no
common code between the two cases: jjm and jjm + 1. Also, it is
simpler now to just remove the argument "direct" from filtreg_v instead
of allowing it and then stopping the program if it is false.

Removed the computation of pkf in reanalyse2nat, was not used.

As a consequence of the split of filtreg, had to extract the
computation of pkf out of exner_hyb. This is clearer anyway because we
want to be able to call exner_hyb with any size in the first two
dimensions (as in test_disvert). But at the same time exner_hyb
required particular sizes for the computation of pkf. It was
awkward. The only computation of pkf is now in leapfrog.

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

  ViewVC Help
Powered by ViewVC 1.1.21