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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 154 - (show annotations)
Tue Jul 7 17:49:23 2015 UTC (8 years, 10 months ago) by guez
File size: 61058 byte(s)
Removed argument dtphys of physiq. Use it directly from comconst in
physiq instead.

Donwgraded variables eignfnu, eignfnv of module inifgn_m to dummy
arguments of SUBROUTINE inifgn. They were not used elsewhere than in
the calling procedure inifilr. Renamed argument dv of inifgn to eignval_v.

Made alboc and alboc_cd independent of the size of arguments. Now we
can call them only at indices knindex in interfsurf_hq, where we need
them. Fixed a bug in alboc_cd: rmu0 was modified, and the
corresponding actual argument in interfsurf_hq is an intent(in)
argument of interfsurf_hq.

Variables of size knon instead of klon in interfsur_lim and interfsurf_hq.

Removed argument alb_new of interfsurf_hq because it was the same than
alblw. Simplified test on cycle_diurne, following LMDZ.

Moved tests on nbapp_rad from physiq to read_clesphys2. No need for
separate counter itaprad, we can use itap. Define lmt_pas and radpas
from integer input parameters instead of real-type computed values.

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

  ViewVC Help
Powered by ViewVC 1.1.21