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

Contents of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 101 - (show annotations)
Mon Jul 7 17:45:21 2014 UTC (9 years, 10 months ago) by guez
File size: 62337 byte(s)
Removed unused files "interfoce_slab.f" and "gath2cpl.f". Removed
unused variables coastalflow and riverflow of module
interface_surf. Removed unused arguments cal, radsol, dif_grnd,
fluxlat, fluxsens, dflux_s, dflux_l of procedure fonte_neige. Removed
unused arguments tslab, seaice of procedure interfsurf_hq and
clqh. Removed unused arguments seaice of procedure clmain.

In interfsurf_hq, used variable soil_model of module clesphys2 instead
of cascading it as an argument from physiq.

In phyetat0, stop if masque not found.

Variable TS instead of "TS[0-9][0-9]" in "(re)startphy.nc", with
additional dimension nbsrf.

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

  ViewVC Help
Powered by ViewVC 1.1.21