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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 151 - (show annotations)
Tue Jun 23 15:14:20 2015 UTC (8 years, 10 months ago) by guez
File size: 61491 byte(s)
In procedure inifilr, only a part of the arrays modfrstu and modfrstv
were defined. Split these into 4 arrays that are fully defined and
used: modfrst[ns][uv].

Clarified the logic for the computation of jfilt[ns][uv]. Changed the
initial value of the search so that the initial values for northern
hemisphere and southern hemisphere cannot be the same.

Clarified the logic for the computation of modfrst[ns][uv]: removed
the cycle and exit instructions.

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 date0
507
508 ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
509 REAL ztsol(klon)
510 REAL d_h_vcol, d_qt, d_ec
511 REAL, SAVE:: d_h_vcol_phy
512 REAL zero_v(klon)
513 CHARACTER(LEN = 20) tit
514 INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
515 INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
516
517 REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique
518 REAL ZRCPD
519
520 REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m
521 REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
522 REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille
523 REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes s/1 maille
524
525 ! Aerosol effects:
526
527 REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
528
529 REAL, save:: sulfate_pi(klon, llm)
530 ! SO4 aerosol concentration, in micro g/m3, pre-industrial value
531
532 REAL cldtaupi(klon, llm)
533 ! cloud optical thickness for pre-industrial (pi) aerosols
534
535 REAL re(klon, llm) ! Cloud droplet effective radius
536 REAL fl(klon, llm) ! denominator of re
537
538 ! Aerosol optical properties
539 REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
540 REAL, save:: cg_ae(klon, llm, 2)
541
542 REAL topswad(klon), solswad(klon) ! aerosol direct effect
543 REAL topswai(klon), solswai(klon) ! aerosol indirect effect
544
545 REAL aerindex(klon) ! POLDER aerosol index
546
547 LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
548 LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
549
550 REAL:: bl95_b0 = 2., bl95_b1 = 0.2
551 ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
552 ! B). They link cloud droplet number concentration to aerosol mass
553 ! concentration.
554
555 SAVE u10m
556 SAVE v10m
557 SAVE t2m
558 SAVE q2m
559 SAVE ffonte
560 SAVE fqcalving
561 SAVE rain_con
562 SAVE snow_con
563 SAVE topswai
564 SAVE topswad
565 SAVE solswai
566 SAVE solswad
567 SAVE d_u_con
568 SAVE d_v_con
569
570 real zmasse(klon, llm)
571 ! (column-density of mass of air in a cell, in kg m-2)
572
573 real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
574
575 namelist /physiq_nml/ ok_journe, ok_mensuel, ok_instan, fact_cldcon, &
576 facttemps, ok_newmicro, iflag_cldcon, ratqsbas, ratqshaut, if_ebil, &
577 ok_ade, ok_aie, bl95_b0, bl95_b1, iflag_thermals, nsplit_thermals
578
579 !----------------------------------------------------------------
580
581 IF (if_ebil >= 1) zero_v = 0.
582 IF (nqmx < 2) CALL abort_gcm('physiq', &
583 'eaux vapeur et liquide sont indispensables', 1)
584
585 test_firstcal: IF (firstcal) THEN
586 ! initialiser
587 u10m = 0.
588 v10m = 0.
589 t2m = 0.
590 q2m = 0.
591 ffonte = 0.
592 fqcalving = 0.
593 piz_ae = 0.
594 tau_ae = 0.
595 cg_ae = 0.
596 rain_con = 0.
597 snow_con = 0.
598 topswai = 0.
599 topswad = 0.
600 solswai = 0.
601 solswad = 0.
602
603 d_u_con = 0.
604 d_v_con = 0.
605 rnebcon0 = 0.
606 clwcon0 = 0.
607 rnebcon = 0.
608 clwcon = 0.
609
610 pblh =0. ! Hauteur de couche limite
611 plcl =0. ! Niveau de condensation de la CLA
612 capCL =0. ! CAPE de couche limite
613 oliqCL =0. ! eau_liqu integree de couche limite
614 cteiCL =0. ! cloud top instab. crit. couche limite
615 pblt =0. ! T a la Hauteur de couche limite
616 therm =0.
617 trmb1 =0. ! deep_cape
618 trmb2 =0. ! inhibition
619 trmb3 =0. ! Point Omega
620
621 IF (if_ebil >= 1) d_h_vcol_phy = 0.
622
623 iflag_thermals = 0
624 nsplit_thermals = 1
625 print *, "Enter namelist 'physiq_nml'."
626 read(unit=*, nml=physiq_nml)
627 write(unit_nml, nml=physiq_nml)
628
629 call conf_phys
630
631 ! Initialiser les compteurs:
632
633 frugs = 0.
634 itap = 0
635 itaprad = 0
636 CALL phyetat0(pctsrf, ftsol, ftsoil, tslab, seaice, fqsurf, qsol, &
637 fsnow, falbe, falblw, fevap, rain_fall, snow_fall, solsw, sollw, &
638 dlw, radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, &
639 zval, t_ancien, q_ancien, ancien_ok, rnebcon, ratqs, clwcon, &
640 run_off_lic_0, sig1, w01)
641
642 ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
643 q2 = 1e-8
644
645 radpas = NINT(86400. / dtphys / nbapp_rad)
646
647 ! on remet le calendrier a zero
648 IF (raz_date) itau_phy = 0
649
650 CALL printflag(radpas, ok_journe, ok_instan, ok_region)
651
652 IF (dtphys * radpas > 21600. .AND. cycle_diurne) THEN
653 print *, "Au minimum 4 appels par jour si cycle diurne"
654 call abort_gcm('physiq', &
655 "Nombre d'appels au rayonnement insuffisant", 1)
656 ENDIF
657
658 ! Initialisation pour le sch\'ema de convection d'Emanuel :
659 IF (iflag_con >= 3) THEN
660 ibas_con = 1
661 itop_con = 1
662 ENDIF
663
664 IF (ok_orodr) THEN
665 rugoro = MAX(1e-5, zstd * zsig / 2)
666 CALL SUGWD(paprs, play)
667 else
668 rugoro = 0.
669 ENDIF
670
671 lmt_pas = NINT(86400. / dtphys) ! tous les jours
672 print *, 'Number of time steps of "physics" per day: ', lmt_pas
673
674 ecrit_ins = NINT(ecrit_ins/dtphys)
675 ecrit_hf = NINT(ecrit_hf/dtphys)
676 ecrit_mth = NINT(ecrit_mth/dtphys)
677 ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
678 ecrit_reg = NINT(ecrit_reg/dtphys)
679
680 ! Initialisation des sorties
681
682 call ini_histins(dtphys, ok_instan, nid_ins)
683 CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
684 ! Positionner date0 pour initialisation de ORCHIDEE
685 print *, 'physiq date0: ', date0
686 ENDIF test_firstcal
687
688 ! We will modify variables *_seri and we will not touch variables
689 ! u, v, t, qx:
690 t_seri = t
691 u_seri = u
692 v_seri = v
693 q_seri = qx(:, :, ivap)
694 ql_seri = qx(:, :, iliq)
695 tr_seri = qx(:, :, 3: nqmx)
696
697 ztsol = sum(ftsol * pctsrf, dim = 2)
698
699 IF (if_ebil >= 1) THEN
700 tit = 'after dynamics'
701 CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
702 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
703 ! Comme les tendances de la physique sont ajout\'es dans la
704 ! dynamique, la variation d'enthalpie par la dynamique devrait
705 ! \^etre \'egale \`a la variation de la physique au pas de temps
706 ! pr\'ec\'edent. Donc la somme de ces 2 variations devrait \^etre
707 ! nulle.
708 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
709 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
710 d_qt, 0.)
711 END IF
712
713 ! Diagnostic de la tendance dynamique :
714 IF (ancien_ok) THEN
715 DO k = 1, llm
716 DO i = 1, klon
717 d_t_dyn(i, k) = (t_seri(i, k) - t_ancien(i, k)) / dtphys
718 d_q_dyn(i, k) = (q_seri(i, k) - q_ancien(i, k)) / dtphys
719 ENDDO
720 ENDDO
721 ELSE
722 DO k = 1, llm
723 DO i = 1, klon
724 d_t_dyn(i, k) = 0.
725 d_q_dyn(i, k) = 0.
726 ENDDO
727 ENDDO
728 ancien_ok = .TRUE.
729 ENDIF
730
731 ! Ajouter le geopotentiel du sol:
732 DO k = 1, llm
733 DO i = 1, klon
734 zphi(i, k) = pphi(i, k) + pphis(i)
735 ENDDO
736 ENDDO
737
738 ! Check temperatures:
739 CALL hgardfou(t_seri, ftsol)
740
741 ! Incrémenter le compteur de la physique
742 itap = itap + 1
743 julien = MOD(dayvrai, 360)
744 if (julien == 0) julien = 360
745
746 forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
747
748 ! Prescrire l'ozone :
749 wo = ozonecm(REAL(julien), paprs)
750
751 ! \'Evaporation de l'eau liquide nuageuse :
752 DO k = 1, llm
753 DO i = 1, klon
754 zb = MAX(0., ql_seri(i, k))
755 t_seri(i, k) = t_seri(i, k) &
756 - zb * RLVTT / RCPD / (1. + RVTMP2 * q_seri(i, k))
757 q_seri(i, k) = q_seri(i, k) + zb
758 ENDDO
759 ENDDO
760 ql_seri = 0.
761
762 IF (if_ebil >= 2) THEN
763 tit = 'after reevap'
764 CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
765 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
766 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
767 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
768 END IF
769
770 frugs = MAX(frugs, 0.000015)
771 zxrugs = sum(frugs * pctsrf, dim = 2)
772
773 ! Calculs nécessaires au calcul de l'albedo dans l'interface avec
774 ! la surface.
775
776 CALL orbite(REAL(julien), longi, dist)
777 IF (cycle_diurne) THEN
778 CALL zenang(longi, time, dtphys * radpas, mu0, fract)
779 ELSE
780 mu0 = -999.999
781 ENDIF
782
783 ! Calcul de l'abedo moyen par maille
784 albsol = sum(falbe * pctsrf, dim = 2)
785 albsollw = sum(falblw * pctsrf, dim = 2)
786
787 ! R\'epartition sous maille des flux longwave et shortwave
788 ! R\'epartition du longwave par sous-surface lin\'earis\'ee
789
790 forall (nsrf = 1: nbsrf)
791 fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
792 * (ztsol - ftsol(:, nsrf))
793 fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
794 END forall
795
796 fder = dlw
797
798 ! Couche limite:
799
800 CALL clmain(dtphys, itap, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, &
801 v_seri, julien, mu0, co2_ppm, ftsol, cdmmax, cdhmax, &
802 ksta, ksta_ter, ok_kzmin, ftsoil, qsol, paprs, play, fsnow, fqsurf, &
803 fevap, falbe, falblw, fluxlat, rain_fall, snow_fall, fsolsw, fsollw, &
804 fder, rlat, frugs, firstcal, agesno, rugoro, d_t_vdf, d_q_vdf, &
805 d_u_vdf, d_v_vdf, d_ts, fluxt, fluxq, fluxu, fluxv, cdragh, cdragm, &
806 q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, v10m, pblh, &
807 capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, plcl, &
808 fqcalving, ffonte, run_off_lic_0, fluxo, fluxg, tslab)
809
810 ! Incr\'ementation des flux
811
812 zxfluxt = 0.
813 zxfluxq = 0.
814 zxfluxu = 0.
815 zxfluxv = 0.
816 DO nsrf = 1, nbsrf
817 DO k = 1, llm
818 DO i = 1, klon
819 zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
820 zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
821 zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
822 zxfluxv(i, k) = zxfluxv(i, k) + fluxv(i, k, nsrf) * pctsrf(i, nsrf)
823 END DO
824 END DO
825 END DO
826 DO i = 1, klon
827 sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
828 evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
829 fder(i) = dlw(i) + dsens(i) + devap(i)
830 ENDDO
831
832 DO k = 1, llm
833 DO i = 1, klon
834 t_seri(i, k) = t_seri(i, k) + d_t_vdf(i, k)
835 q_seri(i, k) = q_seri(i, k) + d_q_vdf(i, k)
836 u_seri(i, k) = u_seri(i, k) + d_u_vdf(i, k)
837 v_seri(i, k) = v_seri(i, k) + d_v_vdf(i, k)
838 ENDDO
839 ENDDO
840
841 IF (if_ebil >= 2) THEN
842 tit = 'after clmain'
843 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
844 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
845 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
846 sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
847 END IF
848
849 ! Update surface temperature:
850
851 DO i = 1, klon
852 zxtsol(i) = 0.
853 zxfluxlat(i) = 0.
854
855 zt2m(i) = 0.
856 zq2m(i) = 0.
857 zu10m(i) = 0.
858 zv10m(i) = 0.
859 zxffonte(i) = 0.
860 zxfqcalving(i) = 0.
861
862 s_pblh(i) = 0.
863 s_lcl(i) = 0.
864 s_capCL(i) = 0.
865 s_oliqCL(i) = 0.
866 s_cteiCL(i) = 0.
867 s_pblT(i) = 0.
868 s_therm(i) = 0.
869 s_trmb1(i) = 0.
870 s_trmb2(i) = 0.
871 s_trmb3(i) = 0.
872
873 IF (abs(pctsrf(i, is_ter) + pctsrf(i, is_lic) + pctsrf(i, is_oce) &
874 + pctsrf(i, is_sic) - 1.) > EPSFRA) print *, &
875 'physiq : probl\`eme sous surface au point ', i, &
876 pctsrf(i, 1 : nbsrf)
877 ENDDO
878 DO nsrf = 1, nbsrf
879 DO i = 1, klon
880 ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)
881 zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
882 zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)
883
884 zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)
885 zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)
886 zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)
887 zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)
888 zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)
889 zxfqcalving(i) = zxfqcalving(i) + &
890 fqcalving(i, nsrf)*pctsrf(i, nsrf)
891 s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)
892 s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)
893 s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)
894 s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)
895 s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)
896 s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)
897 s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)
898 s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)
899 s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)
900 s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)
901 ENDDO
902 ENDDO
903
904 ! Si une sous-fraction n'existe pas, elle prend la température moyenne :
905 DO nsrf = 1, nbsrf
906 DO i = 1, klon
907 IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
908
909 IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)
910 IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)
911 IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)
912 IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)
913 IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)
914 IF (pctsrf(i, nsrf) < epsfra) &
915 fqcalving(i, nsrf) = zxfqcalving(i)
916 IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf) = s_pblh(i)
917 IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf) = s_lcl(i)
918 IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf) = s_capCL(i)
919 IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf) = s_oliqCL(i)
920 IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf) = s_cteiCL(i)
921 IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf) = s_pblT(i)
922 IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf) = s_therm(i)
923 IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf) = s_trmb1(i)
924 IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf) = s_trmb2(i)
925 IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf) = s_trmb3(i)
926 ENDDO
927 ENDDO
928
929 ! Calculer la dérive du flux infrarouge
930
931 DO i = 1, klon
932 dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
933 ENDDO
934
935 IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
936
937 ! Appeler la convection (au choix)
938
939 if (iflag_con == 2) then
940 conv_q = d_q_dyn + d_q_vdf / dtphys
941 conv_t = d_t_dyn + d_t_vdf / dtphys
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 IF (thermcep) THEN
967 zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
968 zqsat = zqsat / (1. - retv * zqsat)
969 ELSE
970 zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
971 ENDIF
972
973 ! Properties of convective clouds
974 clwcon0 = fact_cldcon * clwcon0
975 call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
976 rnebcon0)
977
978 mfd = 0.
979 pen_u = 0.
980 pen_d = 0.
981 pde_d = 0.
982 pde_u = 0.
983 END if
984
985 DO k = 1, llm
986 DO i = 1, klon
987 t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
988 q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
989 u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
990 v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
991 ENDDO
992 ENDDO
993
994 IF (if_ebil >= 2) THEN
995 tit = 'after convect'
996 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
997 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
998 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
999 zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec)
1000 END IF
1001
1002 IF (check) THEN
1003 za = qcheck(paprs, q_seri, ql_seri)
1004 print *, "aprescon = ", za
1005 zx_t = 0.
1006 za = 0.
1007 DO i = 1, klon
1008 za = za + airephy(i)/REAL(klon)
1009 zx_t = zx_t + (rain_con(i)+ &
1010 snow_con(i))*airephy(i)/REAL(klon)
1011 ENDDO
1012 zx_t = zx_t/za*dtphys
1013 print *, "Precip = ", zx_t
1014 ENDIF
1015
1016 IF (iflag_con == 2) THEN
1017 z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
1018 z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
1019 DO k = 1, llm
1020 DO i = 1, klon
1021 IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
1022 q_seri(i, k) = q_seri(i, k) * z_factor(i)
1023 ENDIF
1024 ENDDO
1025 ENDDO
1026 ENDIF
1027
1028 ! Convection s\`eche (thermiques ou ajustement)
1029
1030 d_t_ajs = 0.
1031 d_u_ajs = 0.
1032 d_v_ajs = 0.
1033 d_q_ajs = 0.
1034 fm_therm = 0.
1035 entr_therm = 0.
1036
1037 if (iflag_thermals == 0) then
1038 ! Ajustement sec
1039 CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)
1040 t_seri = t_seri + d_t_ajs
1041 q_seri = q_seri + d_q_ajs
1042 else
1043 ! Thermiques
1044 call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
1045 q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
1046 endif
1047
1048 IF (if_ebil >= 2) THEN
1049 tit = 'after dry_adjust'
1050 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1051 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1052 END IF
1053
1054 ! Caclul des ratqs
1055
1056 ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
1057 ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
1058 if (iflag_cldcon == 1) then
1059 do k = 1, llm
1060 do i = 1, klon
1061 if(ptconv(i, k)) then
1062 ratqsc(i, k) = ratqsbas + fact_cldcon &
1063 * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
1064 else
1065 ratqsc(i, k) = 0.
1066 endif
1067 enddo
1068 enddo
1069 endif
1070
1071 ! ratqs stables
1072 do k = 1, llm
1073 do i = 1, klon
1074 ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
1075 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
1076 enddo
1077 enddo
1078
1079 ! ratqs final
1080 if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
1081 ! les ratqs sont une conbinaison de ratqss et ratqsc
1082 ! ratqs final
1083 ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1084 ! relaxation des ratqs
1085 ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
1086 ratqs = max(ratqs, ratqsc)
1087 else
1088 ! on ne prend que le ratqs stable pour fisrtilp
1089 ratqs = ratqss
1090 endif
1091
1092 CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
1093 d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
1094 pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
1095 psfl, rhcl)
1096
1097 WHERE (rain_lsc < 0) rain_lsc = 0.
1098 WHERE (snow_lsc < 0) snow_lsc = 0.
1099 DO k = 1, llm
1100 DO i = 1, klon
1101 t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
1102 q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
1103 ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
1104 cldfra(i, k) = rneb(i, k)
1105 IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
1106 ENDDO
1107 ENDDO
1108 IF (check) THEN
1109 za = qcheck(paprs, q_seri, ql_seri)
1110 print *, "apresilp = ", za
1111 zx_t = 0.
1112 za = 0.
1113 DO i = 1, klon
1114 za = za + airephy(i)/REAL(klon)
1115 zx_t = zx_t + (rain_lsc(i) &
1116 + snow_lsc(i))*airephy(i)/REAL(klon)
1117 ENDDO
1118 zx_t = zx_t/za*dtphys
1119 print *, "Precip = ", zx_t
1120 ENDIF
1121
1122 IF (if_ebil >= 2) THEN
1123 tit = 'after fisrt'
1124 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1125 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1126 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1127 zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec)
1128 END IF
1129
1130 ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1131
1132 ! 1. NUAGES CONVECTIFS
1133
1134 IF (iflag_cldcon <= -1) THEN
1135 ! seulement pour Tiedtke
1136 snow_tiedtke = 0.
1137 if (iflag_cldcon == -1) then
1138 rain_tiedtke = rain_con
1139 else
1140 rain_tiedtke = 0.
1141 do k = 1, llm
1142 do i = 1, klon
1143 if (d_q_con(i, k) < 0.) then
1144 rain_tiedtke(i) = rain_tiedtke(i)-d_q_con(i, k)/dtphys &
1145 *zmasse(i, k)
1146 endif
1147 enddo
1148 enddo
1149 endif
1150
1151 ! Nuages diagnostiques pour Tiedtke
1152 CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1153 itop_con, diafra, dialiq)
1154 DO k = 1, llm
1155 DO i = 1, klon
1156 IF (diafra(i, k) > cldfra(i, k)) THEN
1157 cldliq(i, k) = dialiq(i, k)
1158 cldfra(i, k) = diafra(i, k)
1159 ENDIF
1160 ENDDO
1161 ENDDO
1162 ELSE IF (iflag_cldcon == 3) THEN
1163 ! On prend pour les nuages convectifs le maximum du calcul de
1164 ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
1165 ! d'un facteur facttemps.
1166 facteur = dtphys * facttemps
1167 do k = 1, llm
1168 do i = 1, klon
1169 rnebcon(i, k) = rnebcon(i, k) * facteur
1170 if (rnebcon0(i, k) * clwcon0(i, k) &
1171 > rnebcon(i, k) * clwcon(i, k)) then
1172 rnebcon(i, k) = rnebcon0(i, k)
1173 clwcon(i, k) = clwcon0(i, k)
1174 endif
1175 enddo
1176 enddo
1177
1178 ! On prend la somme des fractions nuageuses et des contenus en eau
1179 cldfra = min(max(cldfra, rnebcon), 1.)
1180 cldliq = cldliq + rnebcon*clwcon
1181 ENDIF
1182
1183 ! 2. Nuages stratiformes
1184
1185 IF (ok_stratus) THEN
1186 CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
1187 DO k = 1, llm
1188 DO i = 1, klon
1189 IF (diafra(i, k) > cldfra(i, k)) THEN
1190 cldliq(i, k) = dialiq(i, k)
1191 cldfra(i, k) = diafra(i, k)
1192 ENDIF
1193 ENDDO
1194 ENDDO
1195 ENDIF
1196
1197 ! Precipitation totale
1198 DO i = 1, klon
1199 rain_fall(i) = rain_con(i) + rain_lsc(i)
1200 snow_fall(i) = snow_con(i) + snow_lsc(i)
1201 ENDDO
1202
1203 IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1204 dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1205 d_qt, d_ec)
1206
1207 ! Humidit\'e relative pour diagnostic :
1208 DO k = 1, llm
1209 DO i = 1, klon
1210 zx_t = t_seri(i, k)
1211 IF (thermcep) THEN
1212 zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)
1213 zx_qs = MIN(0.5, zx_qs)
1214 zcor = 1./(1.-retv*zx_qs)
1215 zx_qs = zx_qs*zcor
1216 ELSE
1217 IF (zx_t < t_coup) THEN
1218 zx_qs = qsats(zx_t)/play(i, k)
1219 ELSE
1220 zx_qs = qsatl(zx_t)/play(i, k)
1221 ENDIF
1222 ENDIF
1223 zx_rh(i, k) = q_seri(i, k)/zx_qs
1224 zqsat(i, k) = zx_qs
1225 ENDDO
1226 ENDDO
1227
1228 ! Introduce the aerosol direct and first indirect radiative forcings:
1229 IF (ok_ade .OR. ok_aie) THEN
1230 ! Get sulfate aerosol distribution :
1231 CALL readsulfate(dayvrai, time, firstcal, sulfate)
1232 CALL readsulfate_preind(dayvrai, time, firstcal, sulfate_pi)
1233
1234 CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &
1235 aerindex)
1236 ELSE
1237 tau_ae = 0.
1238 piz_ae = 0.
1239 cg_ae = 0.
1240 ENDIF
1241
1242 ! Param\`etres optiques des nuages et quelques param\`etres pour
1243 ! diagnostics :
1244 if (ok_newmicro) then
1245 CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1246 cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1247 sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
1248 else
1249 CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1250 cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
1251 bl95_b1, cldtaupi, re, fl)
1252 endif
1253
1254 IF (MOD(itaprad, radpas) == 0) THEN
1255 ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1256 DO i = 1, klon
1257 albsol(i) = falbe(i, is_oce) * pctsrf(i, is_oce) &
1258 + falbe(i, is_lic) * pctsrf(i, is_lic) &
1259 + falbe(i, is_ter) * pctsrf(i, is_ter) &
1260 + falbe(i, is_sic) * pctsrf(i, is_sic)
1261 albsollw(i) = falblw(i, is_oce) * pctsrf(i, is_oce) &
1262 + falblw(i, is_lic) * pctsrf(i, is_lic) &
1263 + falblw(i, is_ter) * pctsrf(i, is_ter) &
1264 + falblw(i, is_sic) * pctsrf(i, is_sic)
1265 ENDDO
1266 ! Rayonnement (compatible Arpege-IFS) :
1267 CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, &
1268 albsollw, t_seri, q_seri, wo, cldfra, cldemi, cldtau, heat, &
1269 heat0, cool, cool0, radsol, albpla, topsw, toplw, solsw, sollw, &
1270 sollwdown, topsw0, toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, &
1271 lwup, swdn0, swdn, swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, &
1272 cg_ae, topswad, solswad, cldtaupi, topswai, solswai)
1273 itaprad = 0
1274 ENDIF
1275
1276 itaprad = itaprad + 1
1277
1278 ! Ajouter la tendance des rayonnements (tous les pas)
1279
1280 DO k = 1, llm
1281 DO i = 1, klon
1282 t_seri(i, k) = t_seri(i, k) + (heat(i, k)-cool(i, k)) * dtphys/86400.
1283 ENDDO
1284 ENDDO
1285
1286 IF (if_ebil >= 2) THEN
1287 tit = 'after rad'
1288 CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1289 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1290 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1291 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
1292 END IF
1293
1294 ! Calculer l'hydrologie de la surface
1295 DO i = 1, klon
1296 zxqsurf(i) = 0.
1297 zxsnow(i) = 0.
1298 ENDDO
1299 DO nsrf = 1, nbsrf
1300 DO i = 1, klon
1301 zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)
1302 zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)
1303 ENDDO
1304 ENDDO
1305
1306 ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1307
1308 DO i = 1, klon
1309 bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1310 ENDDO
1311
1312 ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1313
1314 IF (ok_orodr) THEN
1315 ! selection des points pour lesquels le shema est actif:
1316 igwd = 0
1317 DO i = 1, klon
1318 itest(i) = 0
1319 IF (((zpic(i)-zmea(i)) > 100.).AND.(zstd(i) > 10.)) THEN
1320 itest(i) = 1
1321 igwd = igwd + 1
1322 idx(igwd) = i
1323 ENDIF
1324 ENDDO
1325
1326 CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1327 zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1328 zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1329
1330 ! ajout des tendances
1331 DO k = 1, llm
1332 DO i = 1, klon
1333 t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
1334 u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
1335 v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
1336 ENDDO
1337 ENDDO
1338 ENDIF
1339
1340 IF (ok_orolf) THEN
1341 ! S\'election des points pour lesquels le sch\'ema est actif :
1342 igwd = 0
1343 DO i = 1, klon
1344 itest(i) = 0
1345 IF ((zpic(i) - zmea(i)) > 100.) THEN
1346 itest(i) = 1
1347 igwd = igwd + 1
1348 idx(igwd) = i
1349 ENDIF
1350 ENDDO
1351
1352 CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
1353 itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &
1354 d_t_lif, d_u_lif, d_v_lif)
1355
1356 ! Ajout des tendances :
1357 DO k = 1, llm
1358 DO i = 1, klon
1359 t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
1360 u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
1361 v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
1362 ENDDO
1363 ENDDO
1364 ENDIF
1365
1366 ! Stress n\'ecessaires : toute la physique
1367
1368 DO i = 1, klon
1369 zustrph(i) = 0.
1370 zvstrph(i) = 0.
1371 ENDDO
1372 DO k = 1, llm
1373 DO i = 1, klon
1374 zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1375 * zmasse(i, k)
1376 zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1377 * zmasse(i, k)
1378 ENDDO
1379 ENDDO
1380
1381 CALL aaam_bud(ra, rg, romega, rlat, rlon, pphis, zustrdr, zustrli, &
1382 zustrph, zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1383
1384 IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1385 2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1386 d_qt, d_ec)
1387
1388 ! Calcul des tendances traceurs
1389 call phytrac(itap, lmt_pas, julien, time, firstcal, lafin, dtphys, t, &
1390 paprs, play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, &
1391 yu1, yv1, ftsol, pctsrf, frac_impa, frac_nucl, pphis, da, phi, mp, &
1392 upwd, dnwd, tr_seri, zmasse)
1393
1394 IF (offline) call phystokenc(dtphys, rlon, rlat, t, mfu, mfd, pen_u, &
1395 pde_u, pen_d, pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, &
1396 pctsrf, frac_impa, frac_nucl, pphis, airephy, dtphys, itap)
1397
1398 ! Calculer le transport de l'eau et de l'energie (diagnostique)
1399 CALL transp(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, &
1400 ue, uq)
1401
1402 ! diag. bilKP
1403
1404 CALL transp_lay(paprs, zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
1405 ve_lay, vq_lay, ue_lay, uq_lay)
1406
1407 ! Accumuler les variables a stocker dans les fichiers histoire:
1408
1409 ! conversion Ec -> E thermique
1410 DO k = 1, llm
1411 DO i = 1, klon
1412 ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))
1413 d_t_ec(i, k) = 0.5 / ZRCPD &
1414 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
1415 t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
1416 d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1417 END DO
1418 END DO
1419
1420 IF (if_ebil >= 1) THEN
1421 tit = 'after physic'
1422 CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1423 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1424 ! Comme les tendances de la physique sont ajoute dans la dynamique,
1425 ! on devrait avoir que la variation d'entalpie par la dynamique
1426 ! est egale a la variation de la physique au pas de temps precedent.
1427 ! Donc la somme de ces 2 variations devrait etre nulle.
1428 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1429 evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec)
1430 d_h_vcol_phy = d_h_vcol
1431 END IF
1432
1433 ! SORTIES
1434
1435 ! prw = eau precipitable
1436 DO i = 1, klon
1437 prw(i) = 0.
1438 DO k = 1, llm
1439 prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)
1440 ENDDO
1441 ENDDO
1442
1443 ! Convertir les incrementations en tendances
1444
1445 DO k = 1, llm
1446 DO i = 1, klon
1447 d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1448 d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1449 d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1450 d_qx(i, k, ivap) = (q_seri(i, k) - qx(i, k, ivap)) / dtphys
1451 d_qx(i, k, iliq) = (ql_seri(i, k) - qx(i, k, iliq)) / dtphys
1452 ENDDO
1453 ENDDO
1454
1455 DO iq = 3, nqmx
1456 DO k = 1, llm
1457 DO i = 1, klon
1458 d_qx(i, k, iq) = (tr_seri(i, k, iq-2) - qx(i, k, iq)) / dtphys
1459 ENDDO
1460 ENDDO
1461 ENDDO
1462
1463 ! Sauvegarder les valeurs de t et q a la fin de la physique:
1464 DO k = 1, llm
1465 DO i = 1, klon
1466 t_ancien(i, k) = t_seri(i, k)
1467 q_ancien(i, k) = q_seri(i, k)
1468 ENDDO
1469 ENDDO
1470
1471 ! Ecriture des sorties
1472 call write_histins
1473
1474 ! Si c'est la fin, il faut conserver l'etat de redemarrage
1475 IF (lafin) THEN
1476 itau_phy = itau_phy + itap
1477 CALL phyredem("restartphy.nc", pctsrf, ftsol, ftsoil, tslab, seaice, &
1478 fqsurf, qsol, fsnow, falbe, falblw, fevap, rain_fall, snow_fall, &
1479 solsw, sollw, dlw, radsol, frugs, agesno, zmea, zstd, zsig, zgam, &
1480 zthe, zpic, zval, t_ancien, q_ancien, rnebcon, ratqs, clwcon, &
1481 run_off_lic_0, sig1, w01)
1482 ENDIF
1483
1484 firstcal = .FALSE.
1485
1486 contains
1487
1488 subroutine write_histins
1489
1490 ! From phylmd/write_histins.h, version 1.2 2005/05/25 13:10:09
1491
1492 use dimens_m, only: iim, jjm
1493 USE histsync_m, ONLY: histsync
1494 USE histwrite_m, ONLY: histwrite
1495
1496 integer i, itau_w ! pas de temps ecriture
1497 REAL zx_tmp_2d(iim, jjm + 1), zx_tmp_3d(iim, jjm + 1, llm)
1498
1499 !--------------------------------------------------
1500
1501 IF (ok_instan) THEN
1502 ! Champs 2D:
1503
1504 itau_w = itau_phy + itap
1505
1506 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, pphis, zx_tmp_2d)
1507 CALL histwrite(nid_ins, "phis", itau_w, zx_tmp_2d)
1508
1509 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, airephy, zx_tmp_2d)
1510 CALL histwrite(nid_ins, "aire", itau_w, zx_tmp_2d)
1511
1512 DO i = 1, klon
1513 zx_tmp_fi2d(i) = paprs(i, 1)
1514 ENDDO
1515 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1516 CALL histwrite(nid_ins, "psol", itau_w, zx_tmp_2d)
1517
1518 DO i = 1, klon
1519 zx_tmp_fi2d(i) = rain_fall(i) + snow_fall(i)
1520 ENDDO
1521 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1522 CALL histwrite(nid_ins, "precip", itau_w, zx_tmp_2d)
1523
1524 DO i = 1, klon
1525 zx_tmp_fi2d(i) = rain_lsc(i) + snow_lsc(i)
1526 ENDDO
1527 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1528 CALL histwrite(nid_ins, "plul", itau_w, zx_tmp_2d)
1529
1530 DO i = 1, klon
1531 zx_tmp_fi2d(i) = rain_con(i) + snow_con(i)
1532 ENDDO
1533 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1534 CALL histwrite(nid_ins, "pluc", itau_w, zx_tmp_2d)
1535
1536 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxtsol, zx_tmp_2d)
1537 CALL histwrite(nid_ins, "tsol", itau_w, zx_tmp_2d)
1538 !ccIM
1539 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zt2m, zx_tmp_2d)
1540 CALL histwrite(nid_ins, "t2m", itau_w, zx_tmp_2d)
1541
1542 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zq2m, zx_tmp_2d)
1543 CALL histwrite(nid_ins, "q2m", itau_w, zx_tmp_2d)
1544
1545 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zu10m, zx_tmp_2d)
1546 CALL histwrite(nid_ins, "u10m", itau_w, zx_tmp_2d)
1547
1548 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zv10m, zx_tmp_2d)
1549 CALL histwrite(nid_ins, "v10m", itau_w, zx_tmp_2d)
1550
1551 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, snow_fall, zx_tmp_2d)
1552 CALL histwrite(nid_ins, "snow", itau_w, zx_tmp_2d)
1553
1554 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragm, zx_tmp_2d)
1555 CALL histwrite(nid_ins, "cdrm", itau_w, zx_tmp_2d)
1556
1557 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, cdragh, zx_tmp_2d)
1558 CALL histwrite(nid_ins, "cdrh", itau_w, zx_tmp_2d)
1559
1560 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, toplw, zx_tmp_2d)
1561 CALL histwrite(nid_ins, "topl", itau_w, zx_tmp_2d)
1562
1563 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, evap, zx_tmp_2d)
1564 CALL histwrite(nid_ins, "evap", itau_w, zx_tmp_2d)
1565
1566 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, solsw, zx_tmp_2d)
1567 CALL histwrite(nid_ins, "sols", itau_w, zx_tmp_2d)
1568
1569 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollw, zx_tmp_2d)
1570 CALL histwrite(nid_ins, "soll", itau_w, zx_tmp_2d)
1571
1572 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sollwdown, zx_tmp_2d)
1573 CALL histwrite(nid_ins, "solldown", itau_w, zx_tmp_2d)
1574
1575 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, bils, zx_tmp_2d)
1576 CALL histwrite(nid_ins, "bils", itau_w, zx_tmp_2d)
1577
1578 zx_tmp_fi2d(1:klon) = -1*sens(1:klon)
1579 ! CALL gr_fi_ecrit(1, klon, iim, jjm + 1, sens, zx_tmp_2d)
1580 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1581 CALL histwrite(nid_ins, "sens", itau_w, zx_tmp_2d)
1582
1583 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, fder, zx_tmp_2d)
1584 CALL histwrite(nid_ins, "fder", itau_w, zx_tmp_2d)
1585
1586 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_oce), zx_tmp_2d)
1587 CALL histwrite(nid_ins, "dtsvdfo", itau_w, zx_tmp_2d)
1588
1589 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_ter), zx_tmp_2d)
1590 CALL histwrite(nid_ins, "dtsvdft", itau_w, zx_tmp_2d)
1591
1592 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_lic), zx_tmp_2d)
1593 CALL histwrite(nid_ins, "dtsvdfg", itau_w, zx_tmp_2d)
1594
1595 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, d_ts(1, is_sic), zx_tmp_2d)
1596 CALL histwrite(nid_ins, "dtsvdfi", itau_w, zx_tmp_2d)
1597
1598 DO nsrf = 1, nbsrf
1599 !XXX
1600 zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)*100.
1601 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1602 CALL histwrite(nid_ins, "pourc_"//clnsurf(nsrf), itau_w, &
1603 zx_tmp_2d)
1604
1605 zx_tmp_fi2d(1 : klon) = pctsrf(1 : klon, nsrf)
1606 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1607 CALL histwrite(nid_ins, "fract_"//clnsurf(nsrf), itau_w, &
1608 zx_tmp_2d)
1609
1610 zx_tmp_fi2d(1 : klon) = fluxt(1 : klon, 1, nsrf)
1611 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1612 CALL histwrite(nid_ins, "sens_"//clnsurf(nsrf), itau_w, &
1613 zx_tmp_2d)
1614
1615 zx_tmp_fi2d(1 : klon) = fluxlat(1 : klon, nsrf)
1616 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1617 CALL histwrite(nid_ins, "lat_"//clnsurf(nsrf), itau_w, &
1618 zx_tmp_2d)
1619
1620 zx_tmp_fi2d(1 : klon) = ftsol(1 : klon, nsrf)
1621 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1622 CALL histwrite(nid_ins, "tsol_"//clnsurf(nsrf), itau_w, &
1623 zx_tmp_2d)
1624
1625 zx_tmp_fi2d(1 : klon) = fluxu(1 : klon, 1, nsrf)
1626 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1627 CALL histwrite(nid_ins, "taux_"//clnsurf(nsrf), itau_w, &
1628 zx_tmp_2d)
1629
1630 zx_tmp_fi2d(1 : klon) = fluxv(1 : klon, 1, nsrf)
1631 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1632 CALL histwrite(nid_ins, "tauy_"//clnsurf(nsrf), itau_w, &
1633 zx_tmp_2d)
1634
1635 zx_tmp_fi2d(1 : klon) = frugs(1 : klon, nsrf)
1636 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1637 CALL histwrite(nid_ins, "rugs_"//clnsurf(nsrf), itau_w, &
1638 zx_tmp_2d)
1639
1640 zx_tmp_fi2d(1 : klon) = falbe(1 : klon, nsrf)
1641 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zx_tmp_fi2d, zx_tmp_2d)
1642 CALL histwrite(nid_ins, "albe_"//clnsurf(nsrf), itau_w, &
1643 zx_tmp_2d)
1644
1645 END DO
1646 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsol, zx_tmp_2d)
1647 CALL histwrite(nid_ins, "albs", itau_w, zx_tmp_2d)
1648 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, albsollw, zx_tmp_2d)
1649 CALL histwrite(nid_ins, "albslw", itau_w, zx_tmp_2d)
1650
1651 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, zxrugs, zx_tmp_2d)
1652 CALL histwrite(nid_ins, "rugs", itau_w, zx_tmp_2d)
1653
1654 !HBTM2
1655
1656 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblh, zx_tmp_2d)
1657 CALL histwrite(nid_ins, "s_pblh", itau_w, zx_tmp_2d)
1658
1659 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_pblt, zx_tmp_2d)
1660 CALL histwrite(nid_ins, "s_pblt", itau_w, zx_tmp_2d)
1661
1662 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_lcl, zx_tmp_2d)
1663 CALL histwrite(nid_ins, "s_lcl", itau_w, zx_tmp_2d)
1664
1665 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_capCL, zx_tmp_2d)
1666 CALL histwrite(nid_ins, "s_capCL", itau_w, zx_tmp_2d)
1667
1668 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_oliqCL, zx_tmp_2d)
1669 CALL histwrite(nid_ins, "s_oliqCL", itau_w, zx_tmp_2d)
1670
1671 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_cteiCL, zx_tmp_2d)
1672 CALL histwrite(nid_ins, "s_cteiCL", itau_w, zx_tmp_2d)
1673
1674 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_therm, zx_tmp_2d)
1675 CALL histwrite(nid_ins, "s_therm", itau_w, zx_tmp_2d)
1676
1677 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb1, zx_tmp_2d)
1678 CALL histwrite(nid_ins, "s_trmb1", itau_w, zx_tmp_2d)
1679
1680 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb2, zx_tmp_2d)
1681 CALL histwrite(nid_ins, "s_trmb2", itau_w, zx_tmp_2d)
1682
1683 CALL gr_fi_ecrit(1, klon, iim, jjm + 1, s_trmb3, zx_tmp_2d)
1684 CALL histwrite(nid_ins, "s_trmb3", itau_w, zx_tmp_2d)
1685
1686 ! Champs 3D:
1687
1688 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, t_seri, zx_tmp_3d)
1689 CALL histwrite(nid_ins, "temp", itau_w, zx_tmp_3d)
1690
1691 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, u_seri, zx_tmp_3d)
1692 CALL histwrite(nid_ins, "vitu", itau_w, zx_tmp_3d)
1693
1694 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, v_seri, zx_tmp_3d)
1695 CALL histwrite(nid_ins, "vitv", itau_w, zx_tmp_3d)
1696
1697 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, zphi, zx_tmp_3d)
1698 CALL histwrite(nid_ins, "geop", itau_w, zx_tmp_3d)
1699
1700 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, play, zx_tmp_3d)
1701 CALL histwrite(nid_ins, "pres", itau_w, zx_tmp_3d)
1702
1703 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_t_vdf, zx_tmp_3d)
1704 CALL histwrite(nid_ins, "dtvdf", itau_w, zx_tmp_3d)
1705
1706 CALL gr_fi_ecrit(llm, klon, iim, jjm + 1, d_q_vdf, zx_tmp_3d)
1707 CALL histwrite(nid_ins, "dqvdf", itau_w, zx_tmp_3d)
1708
1709 call histsync(nid_ins)
1710 ENDIF
1711
1712 end subroutine write_histins
1713
1714 END SUBROUTINE physiq
1715
1716 end module physiq_m

  ViewVC Help
Powered by ViewVC 1.1.21