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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Jul 8 17:03:45 2015 UTC (8 years, 10 months ago) by guez
Original Path: trunk/Sources/phylmd/physiq.f
File size: 60331 byte(s)
Do not write any longer to startphy.nc nor read from restartphy.nc the
NetCDF variable ALBLW: it was the same than ALBE. ALBE was for the
visible and ALBLW for the near infrared. In physiq, use only variables
falbe and albsol, removed falblw and albsollw. See revision 888 of
LMDZ.

Removed unused arguments pdp of SUBROUTINE lwbv, ptave of SUBROUTINE
lwv, kuaer of SUBROUTINE lwvd, nq of SUBROUTINE initphysto.

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

  ViewVC Help
Powered by ViewVC 1.1.21