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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 221 - (hide annotations)
Thu Apr 20 14:44:47 2017 UTC (7 years, 1 month ago) by guez
File size: 39184 byte(s) 
clcdrag is no longer used in LMDZ. Replaced by cdrag in LMDZ. In cdrag
in LMDZ, zxli is a symbolic constant, false. So removed case zxli true
in LMDZE.

read_sst is called zero (if no ocean point on the whole planet) time or
once per call of physiq. If mod(itap - 1, lmt_pas) == 0 then we have
advanced in time of lmt_pas and deja_lu is necessarily false.

qsat[sl] and dqsat[sl] were never called.

Added output of qsurf in histins, following LMDZ.

Last dummy argument dtime of phystokenc is always the same as first
dummy argument pdtphys, removed dtime.

Removed make rules for nag_xref95, since it does not exist any longer.
1 guez 3 module physiq_m
2    
3     IMPLICIT none
4    
5     contains
6    
7 guez 154 SUBROUTINE physiq(lafin, dayvrai, time, paprs, play, pphi, pphis, u, v, t, &
8     qx, omega, d_u, d_v, d_t, d_qx)
9 guez 3
10 guez 72 ! From phylmd/physiq.F, version 1.22 2006/02/20 09:38:28
11     ! (subversion revision 678)
12    
13 guez 154 ! Author: Z. X. Li (LMD/CNRS) 1993
14 guez 3
15 guez 49 ! This is the main procedure for the "physics" part of the program.
16 guez 3
17 guez 56 use aaam_bud_m, only: aaam_bud
18 guez 51 USE abort_gcm_m, ONLY: abort_gcm
19 guez 53 use ajsec_m, only: ajsec
20 guez 52 use calltherm_m, only: calltherm
21 guez 202 USE clesphys, ONLY: cdhmax, cdmmax, ecrit_ins, ksta, ksta_ter, ok_kzmin, &
22     ok_instan
23 guez 209 USE clesphys2, ONLY: conv_emanuel, nbapp_rad, new_oliq, ok_orodr, ok_orolf
24 guez 51 USE clmain_m, ONLY: clmain
25 guez 72 use clouds_gno_m, only: clouds_gno
26 guez 154 use comconst, only: dtphys
27 guez 137 USE comgeomphy, ONLY: airephy
28 guez 51 USE concvl_m, ONLY: concvl
29 guez 207 USE conf_gcm_m, ONLY: offline, lmt_pas
30 guez 51 USE conf_phys_m, ONLY: conf_phys
31 guez 62 use conflx_m, only: conflx
32 guez 51 USE ctherm, ONLY: iflag_thermals, nsplit_thermals
33 guez 52 use diagcld2_m, only: diagcld2
34 guez 90 USE dimens_m, ONLY: llm, nqmx
35 guez 98 USE dimphy, ONLY: klon
36 guez 51 USE dimsoil, ONLY: nsoilmx
37 guez 52 use drag_noro_m, only: drag_noro
38 guez 129 use dynetat0_m, only: day_ref, annee_ref
39 guez 221 USE fcttre, ONLY: foeew
40 guez 68 use fisrtilp_m, only: fisrtilp
41 guez 51 USE hgardfou_m, ONLY: hgardfou
42 guez 191 USE histsync_m, ONLY: histsync
43     USE histwrite_phy_m, ONLY: histwrite_phy
44 guez 51 USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &
45     nbsrf
46 guez 191 USE ini_histins_m, ONLY: ini_histins, nid_ins
47 guez 157 use netcdf95, only: NF95_CLOSE
48 guez 68 use newmicro_m, only: newmicro
49 guez 191 use nr_util, only: assert
50 guez 175 use nuage_m, only: nuage
51 guez 98 USE orbite_m, ONLY: orbite
52 guez 51 USE ozonecm_m, ONLY: ozonecm
53     USE phyetat0_m, ONLY: phyetat0, rlat, rlon
54     USE phyredem_m, ONLY: phyredem
55 guez 157 USE phyredem0_m, ONLY: phyredem0
56 guez 51 USE phystokenc_m, ONLY: phystokenc
57     USE phytrac_m, ONLY: phytrac
58 guez 53 use radlwsw_m, only: radlwsw
59 guez 158 use yoegwd, only: sugwd
60 guez 212 USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt, rmo3, md
61 guez 191 use time_phylmdz, only: itap, increment_itap
62 guez 169 use transp_m, only: transp
63 guez 178 use transp_lay_m, only: transp_lay
64 guez 68 use unit_nml_m, only: unit_nml
65 guez 92 USE ymds2ju_m, ONLY: ymds2ju
66 guez 51 USE yoethf_m, ONLY: r2es, rvtmp2
67 guez 98 use zenang_m, only: zenang
68 guez 3
69 guez 91 logical, intent(in):: lafin ! dernier passage
70 guez 3
71 guez 130 integer, intent(in):: dayvrai
72     ! current day number, based at value 1 on January 1st of annee_ref
73 guez 15
74 guez 90 REAL, intent(in):: time ! heure de la journ\'ee en fraction de jour
75 guez 3
76 guez 98 REAL, intent(in):: paprs(:, :) ! (klon, llm + 1)
77     ! pression pour chaque inter-couche, en Pa
78 guez 15
79 guez 98 REAL, intent(in):: play(:, :) ! (klon, llm)
80     ! pression pour le mileu de chaque couche (en Pa)
81 guez 3
82 guez 190 REAL, intent(in):: pphi(:, :) ! (klon, llm)
83 guez 97 ! géopotentiel de chaque couche (référence sol)
84 guez 3
85 guez 98 REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
86 guez 3
87 guez 98 REAL, intent(in):: u(:, :) ! (klon, llm)
88 guez 202 ! vitesse dans la direction X (de O a E) en m / s
89 guez 51
90 guez 202 REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m / s
91 guez 98 REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
92 guez 3
93 guez 98 REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
94 guez 90 ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
95 guez 3
96 guez 202 REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa / s
97 guez 98 REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)
98     REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
99 guez 202 REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K / s)
100 guez 3
101 guez 98 REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
102     ! tendance physique de "qx" (s-1)
103    
104 guez 91 ! Local:
105    
106 guez 35 LOGICAL:: firstcal = .true.
107    
108 guez 51 LOGICAL, PARAMETER:: ok_stratus = .FALSE.
109 guez 49 ! Ajouter artificiellement les stratus
110    
111 guez 200 ! pour phystoke avec thermiques
112 guez 51 REAL fm_therm(klon, llm + 1)
113 guez 3 REAL entr_therm(klon, llm)
114 guez 51 real, save:: q2(klon, llm + 1, nbsrf)
115 guez 3
116 guez 98 INTEGER, PARAMETER:: ivap = 1 ! indice de traceur pour vapeur d'eau
117     INTEGER, PARAMETER:: iliq = 2 ! indice de traceur pour eau liquide
118 guez 3
119 guez 49 REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
120     LOGICAL, save:: ancien_ok
121 guez 3
122 guez 202 REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K / s)
123     REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg / kg / s)
124 guez 3
125     real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
126    
127 guez 205 REAL, save:: swdn0(klon, llm + 1), swdn(klon, llm + 1)
128     REAL, save:: swup0(klon, llm + 1), swup(klon, llm + 1)
129 guez 3
130 guez 205 REAL, save:: lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
131     REAL, save:: lwup0(klon, llm + 1), lwup(klon, llm + 1)
132 guez 3
133     ! prw: precipitable water
134     real prw(klon)
135    
136 guez 202 ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg / m2)
137     ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg / kg)
138 guez 3 REAL flwp(klon), fiwp(klon)
139     REAL flwc(klon, llm), fiwc(klon, llm)
140    
141     ! Variables propres a la physique
142    
143     INTEGER, save:: radpas
144 guez 125 ! Radiative transfer computations are made every "radpas" call to
145     ! "physiq".
146 guez 3
147 guez 205 REAL, save:: radsol(klon) ! bilan radiatif au sol calcule par code radiatif
148 guez 49 REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
149 guez 3
150 guez 49 REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
151     ! soil temperature of surface fraction
152 guez 3
153 guez 71 REAL, save:: fevap(klon, nbsrf) ! evaporation
154 guez 214 REAL fluxlat(klon, nbsrf)
155 guez 3
156 guez 98 REAL, save:: fqsurf(klon, nbsrf)
157     ! humidite de l'air au contact de la surface
158 guez 3
159 guez 215 REAL, save:: qsol(klon) ! column-density of water in soil, in kg m-2
160     REAL, save:: fsnow(klon, nbsrf) ! \'epaisseur neigeuse
161 guez 155 REAL, save:: falbe(klon, nbsrf) ! albedo visible par type de surface
162 guez 3
163 guez 90 ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
164 guez 13 REAL, save:: zmea(klon) ! orographie moyenne
165     REAL, save:: zstd(klon) ! deviation standard de l'OESM
166     REAL, save:: zsig(klon) ! pente de l'OESM
167     REAL, save:: zgam(klon) ! anisotropie de l'OESM
168     REAL, save:: zthe(klon) ! orientation de l'OESM
169     REAL, save:: zpic(klon) ! Maximum de l'OESM
170     REAL, save:: zval(klon) ! Minimum de l'OESM
171     REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
172 guez 3 REAL zulow(klon), zvlow(klon)
173 guez 178 INTEGER igwd, itest(klon)
174 guez 3
175 guez 189 REAL, save:: agesno(klon, nbsrf) ! age de la neige
176     REAL, save:: run_off_lic_0(klon)
177 guez 3
178 guez 189 ! Variables li\'ees \`a la convection d'Emanuel :
179     REAL, save:: Ma(klon, llm) ! undilute upward mass flux
180     REAL, save:: qcondc(klon, llm) ! in-cld water content from convect
181 guez 72 REAL, save:: sig1(klon, llm), w01(klon, llm)
182 guez 3
183 guez 189 ! Variables pour la couche limite (Alain Lahellec) :
184 guez 3 REAL cdragh(klon) ! drag coefficient pour T and Q
185     REAL cdragm(klon) ! drag coefficient pour vent
186    
187 guez 69 ! Pour phytrac :
188 guez 47 REAL ycoefh(klon, llm) ! coef d'echange pour phytrac
189     REAL yu1(klon) ! vents dans la premiere couche U
190     REAL yv1(klon) ! vents dans la premiere couche V
191 guez 191
192 guez 205 REAL, save:: ffonte(klon, nbsrf)
193     ! flux thermique utilise pour fondre la neige
194    
195     REAL, save:: fqcalving(klon, nbsrf)
196 guez 191 ! flux d'eau "perdue" par la surface et necessaire pour limiter la
197 guez 202 ! hauteur de neige, en kg / m2 / s
198 guez 191
199 guez 3 REAL zxffonte(klon), zxfqcalving(klon)
200    
201 guez 205 REAL, save:: pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
202     REAL, save:: pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation
203    
204     REAL, save:: pfrac_1nucl(klon, llm)
205     ! Produits des coefs lessi nucl (alpha = 1)
206    
207 guez 213 REAL frac_impa(klon, llm) ! fraction d'a\'erosols lessiv\'es (impaction)
208 guez 3 REAL frac_nucl(klon, llm) ! idem (nucleation)
209    
210 guez 101 REAL, save:: rain_fall(klon)
211 guez 202 ! liquid water mass flux (kg / m2 / s), positive down
212 guez 62
213 guez 101 REAL, save:: snow_fall(klon)
214 guez 202 ! solid water mass flux (kg / m2 / s), positive down
215 guez 101
216 guez 3 REAL rain_tiedtke(klon), snow_tiedtke(klon)
217    
218 guez 206 REAL evap(klon) ! flux d'\'evaporation au sol
219     real devap(klon) ! derivative of the evaporation flux at the surface
220     REAL sens(klon) ! flux de chaleur sensible au sol
221     real dsens(klon) ! derivee du flux de chaleur sensible au sol
222 guez 205 REAL, save:: dlw(klon) ! derivee infra rouge
223 guez 3 REAL bils(klon) ! bilan de chaleur au sol
224 guez 190 REAL, save:: fder(klon) ! Derive de flux (sensible et latente)
225 guez 3 REAL ve(klon) ! integr. verticale du transport meri. de l'energie
226     REAL vq(klon) ! integr. verticale du transport meri. de l'eau
227     REAL ue(klon) ! integr. verticale du transport zonal de l'energie
228     REAL uq(klon) ! integr. verticale du transport zonal de l'eau
229    
230 guez 98 REAL, save:: frugs(klon, nbsrf) ! longueur de rugosite
231 guez 3 REAL zxrugs(klon) ! longueur de rugosite
232    
233     ! Conditions aux limites
234    
235     INTEGER julien
236 guez 70 REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
237 guez 217 REAL, save:: albsol(klon) ! albedo du sol total, visible, moyen par maille
238 guez 17 REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
239 guez 212 real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
240 guez 3
241 guez 72 real, save:: clwcon(klon, llm), rnebcon(klon, llm)
242     real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
243 guez 3
244 guez 47 REAL rhcl(klon, llm) ! humiditi relative ciel clair
245     REAL dialiq(klon, llm) ! eau liquide nuageuse
246     REAL diafra(klon, llm) ! fraction nuageuse
247     REAL cldliq(klon, llm) ! eau liquide nuageuse
248     REAL cldfra(klon, llm) ! fraction nuageuse
249     REAL cldtau(klon, llm) ! epaisseur optique
250     REAL cldemi(klon, llm) ! emissivite infrarouge
251 guez 3
252 guez 206 REAL flux_q(klon, nbsrf) ! flux turbulent d'humidite à la surface
253     REAL flux_t(klon, nbsrf) ! flux turbulent de chaleur à la surface
254     REAL flux_u(klon, nbsrf) ! flux turbulent de vitesse u à la surface
255     REAL flux_v(klon, nbsrf) ! flux turbulent de vitesse v à la surface
256 guez 3
257 guez 90 ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
258     ! les variables soient r\'emanentes.
259 guez 53 REAL, save:: heat(klon, llm) ! chauffage solaire
260 guez 154 REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
261 guez 62 REAL, save:: cool(klon, llm) ! refroidissement infrarouge
262 guez 154 REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
263 guez 72 REAL, save:: topsw(klon), toplw(klon), solsw(klon)
264 guez 90 REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
265 guez 72 real, save:: sollwdown(klon) ! downward LW flux at surface
266 guez 62 REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
267 guez 154 REAL, save:: albpla(klon)
268 guez 191 REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
269     REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
270 guez 3
271 guez 202 REAL conv_q(klon, llm) ! convergence de l'humidite (kg / kg / s)
272     REAL conv_t(klon, llm) ! convergence of temperature (K / s)
273 guez 3
274 guez 191 REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut
275     REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree
276 guez 3
277 guez 221 REAL zxfluxlat(klon)
278 guez 118 REAL dist, mu0(klon), fract(klon)
279     real longi
280 guez 3 REAL z_avant(klon), z_apres(klon), z_factor(klon)
281 guez 205 REAL zb
282 guez 103 REAL zx_t, zx_qs, zcor
283 guez 3 real zqsat(klon, llm)
284     INTEGER i, k, iq, nsrf
285     REAL zphi(klon, llm)
286    
287 guez 200 ! cf. Anne Mathieu, variables pour la couche limite atmosphérique (hbtm)
288 guez 3
289 guez 49 REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
290     REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
291     REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
292     REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
293     REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
294 guez 207 REAL, SAVE:: pblt(klon, nbsrf) ! T \`a la hauteur de couche limite
295 guez 49 REAL, SAVE:: therm(klon, nbsrf)
296     REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
297 guez 190 REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
298 guez 49 REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
299 guez 186 ! Grandeurs de sorties
300 guez 3 REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
301     REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
302     REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
303     REAL s_trmb3(klon)
304    
305 guez 175 ! Variables pour la convection de K. Emanuel :
306 guez 3
307 guez 47 REAL upwd(klon, llm) ! saturated updraft mass flux
308     REAL dnwd(klon, llm) ! saturated downdraft mass flux
309 guez 205 REAL, save:: cape(klon)
310 guez 3
311 guez 47 INTEGER iflagctrl(klon) ! flag fonctionnement de convect
312 guez 3
313     ! Variables du changement
314    
315     ! con: convection
316 guez 51 ! lsc: large scale condensation
317 guez 3 ! ajs: ajustement sec
318 guez 90 ! eva: \'evaporation de l'eau liquide nuageuse
319 guez 51 ! vdf: vertical diffusion in boundary layer
320 guez 3 REAL d_t_con(klon, llm), d_q_con(klon, llm)
321 guez 205 REAL, save:: d_u_con(klon, llm), d_v_con(klon, llm)
322 guez 3 REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
323     REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
324     REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
325     REAL rneb(klon, llm)
326    
327 guez 71 REAL mfu(klon, llm), mfd(klon, llm)
328 guez 3 REAL pen_u(klon, llm), pen_d(klon, llm)
329     REAL pde_u(klon, llm), pde_d(klon, llm)
330     INTEGER kcbot(klon), kctop(klon), kdtop(klon)
331 guez 51 REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
332     REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
333 guez 3
334 guez 62 INTEGER, save:: ibas_con(klon), itop_con(klon)
335 guez 183 real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
336 guez 3
337 guez 205 REAL, save:: rain_con(klon)
338     real rain_lsc(klon)
339 guez 183 REAL, save:: snow_con(klon) ! neige (mm / s)
340 guez 180 real snow_lsc(klon)
341 guez 221 REAL d_ts(klon, nbsrf) ! variation of ftsol
342 guez 3
343     REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
344     REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
345    
346     REAL d_u_oro(klon, llm), d_v_oro(klon, llm)
347     REAL d_t_oro(klon, llm)
348     REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
349     REAL d_t_lif(klon, llm)
350    
351 guez 68 REAL, save:: ratqs(klon, llm)
352     real ratqss(klon, llm), ratqsc(klon, llm)
353     real:: ratqsbas = 0.01, ratqshaut = 0.3
354 guez 3
355     ! Parametres lies au nouveau schema de nuages (SB, PDF)
356 guez 68 real:: fact_cldcon = 0.375
357     real:: facttemps = 1.e-4
358     logical:: ok_newmicro = .true.
359 guez 3 real facteur
360    
361 guez 68 integer:: iflag_cldcon = 1
362 guez 3 logical ptconv(klon, llm)
363    
364 guez 175 ! Variables pour effectuer les appels en s\'erie :
365 guez 3
366     REAL t_seri(klon, llm), q_seri(klon, llm)
367 guez 98 REAL ql_seri(klon, llm)
368 guez 3 REAL u_seri(klon, llm), v_seri(klon, llm)
369 guez 98 REAL tr_seri(klon, llm, nqmx - 2)
370 guez 3
371     REAL zx_rh(klon, llm)
372    
373     REAL zustrdr(klon), zvstrdr(klon)
374     REAL zustrli(klon), zvstrli(klon)
375     REAL zustrph(klon), zvstrph(klon)
376     REAL aam, torsfc
377    
378     REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
379     REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
380     REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
381     REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
382    
383     real date0
384 guez 221 REAL tsol(klon)
385 guez 51
386 guez 202 REAL d_t_ec(klon, llm)
387 guez 213 ! tendance due \`a la conversion d'\'energie cin\'etique en
388     ! énergie thermique
389 guez 200
390 guez 205 REAL, save:: t2m(klon, nbsrf), q2m(klon, nbsrf)
391     ! temperature and humidity at 2 m
392    
393     REAL, save:: u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
394 guez 207 REAL zt2m(klon), zq2m(klon) ! température, humidité 2 m moyenne sur 1 maille
395     REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes sur 1 maille
396 guez 3
397 guez 69 ! Aerosol effects:
398    
399 guez 205 REAL, save:: topswad(klon), solswad(klon) ! aerosol direct effect
400 guez 68 LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
401 guez 3
402 guez 68 REAL:: bl95_b0 = 2., bl95_b1 = 0.2
403 guez 69 ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
404     ! B). They link cloud droplet number concentration to aerosol mass
405     ! concentration.
406 guez 68
407 guez 190 real zmasse(klon, llm)
408 guez 17 ! (column-density of mass of air in a cell, in kg m-2)
409    
410 guez 191 integer, save:: ncid_startphy
411 guez 17
412 guez 204 namelist /physiq_nml/ fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, &
413 guez 217 ratqsbas, ratqshaut, ok_ade, bl95_b0, bl95_b1, iflag_thermals, &
414     nsplit_thermals
415 guez 68
416 guez 3 !----------------------------------------------------------------
417    
418 guez 69 IF (nqmx < 2) CALL abort_gcm('physiq', &
419 guez 171 'eaux vapeur et liquide sont indispensables')
420 guez 3
421 guez 7 test_firstcal: IF (firstcal) THEN
422 guez 47 ! initialiser
423 guez 51 u10m = 0.
424     v10m = 0.
425     t2m = 0.
426     q2m = 0.
427     ffonte = 0.
428     fqcalving = 0.
429 guez 98 rain_con = 0.
430     snow_con = 0.
431 guez 72 d_u_con = 0.
432     d_v_con = 0.
433     rnebcon0 = 0.
434     clwcon0 = 0.
435     rnebcon = 0.
436     clwcon = 0.
437 guez 47 pblh =0. ! Hauteur de couche limite
438     plcl =0. ! Niveau de condensation de la CLA
439     capCL =0. ! CAPE de couche limite
440     oliqCL =0. ! eau_liqu integree de couche limite
441     cteiCL =0. ! cloud top instab. crit. couche limite
442 guez 207 pblt =0.
443 guez 47 therm =0.
444     trmb1 =0. ! deep_cape
445 guez 190 trmb2 =0. ! inhibition
446 guez 47 trmb3 =0. ! Point Omega
447 guez 3
448 guez 68 iflag_thermals = 0
449     nsplit_thermals = 1
450     print *, "Enter namelist 'physiq_nml'."
451     read(unit=*, nml=physiq_nml)
452     write(unit_nml, nml=physiq_nml)
453    
454     call conf_phys
455 guez 3
456     ! Initialiser les compteurs:
457    
458     frugs = 0.
459 guez 191 CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
460     fevap, rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
461     agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
462     q_ancien, ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
463     w01, ncid_startphy)
464 guez 3
465 guez 47 ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
466 guez 69 q2 = 1e-8
467 guez 3
468 guez 154 radpas = lmt_pas / nbapp_rad
469 guez 191 print *, "radpas = ", radpas
470 guez 154
471 guez 90 ! Initialisation pour le sch\'ema de convection d'Emanuel :
472 guez 182 IF (conv_emanuel) THEN
473 guez 69 ibas_con = 1
474     itop_con = 1
475 guez 3 ENDIF
476    
477     IF (ok_orodr) THEN
478 guez 13 rugoro = MAX(1e-5, zstd * zsig / 2)
479 guez 54 CALL SUGWD(paprs, play)
480 guez 13 else
481     rugoro = 0.
482 guez 3 ENDIF
483    
484 guez 202 ecrit_ins = NINT(ecrit_ins / dtphys)
485 guez 3
486 guez 47 ! Initialisation des sorties
487 guez 3
488 guez 217 call ini_histins(dtphys, ok_newmicro)
489 guez 129 CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
490 guez 69 ! Positionner date0 pour initialisation de ORCHIDEE
491     print *, 'physiq date0: ', date0
492 guez 202 CALL phyredem0
493 guez 7 ENDIF test_firstcal
494 guez 3
495 guez 91 ! We will modify variables *_seri and we will not touch variables
496 guez 98 ! u, v, t, qx:
497     t_seri = t
498     u_seri = u
499     v_seri = v
500     q_seri = qx(:, :, ivap)
501     ql_seri = qx(:, :, iliq)
502 guez 157 tr_seri = qx(:, :, 3:nqmx)
503 guez 3
504 guez 221 tsol = sum(ftsol * pctsrf, dim = 2)
505 guez 3
506 guez 51 ! Diagnostic de la tendance dynamique :
507 guez 3 IF (ancien_ok) THEN
508     DO k = 1, llm
509     DO i = 1, klon
510 guez 49 d_t_dyn(i, k) = (t_seri(i, k) - t_ancien(i, k)) / dtphys
511     d_q_dyn(i, k) = (q_seri(i, k) - q_ancien(i, k)) / dtphys
512 guez 3 ENDDO
513     ENDDO
514     ELSE
515     DO k = 1, llm
516     DO i = 1, klon
517 guez 72 d_t_dyn(i, k) = 0.
518     d_q_dyn(i, k) = 0.
519 guez 3 ENDDO
520     ENDDO
521     ancien_ok = .TRUE.
522     ENDIF
523    
524     ! Ajouter le geopotentiel du sol:
525     DO k = 1, llm
526     DO i = 1, klon
527     zphi(i, k) = pphi(i, k) + pphis(i)
528     ENDDO
529     ENDDO
530    
531 guez 49 ! Check temperatures:
532 guez 3 CALL hgardfou(t_seri, ftsol)
533    
534 guez 191 call increment_itap
535 guez 130 julien = MOD(dayvrai, 360)
536 guez 3 if (julien == 0) julien = 360
537    
538 guez 103 forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
539 guez 17
540 guez 90 ! \'Evaporation de l'eau liquide nuageuse :
541 guez 51 DO k = 1, llm
542 guez 3 DO i = 1, klon
543 guez 51 zb = MAX(0., ql_seri(i, k))
544     t_seri(i, k) = t_seri(i, k) &
545     - zb * RLVTT / RCPD / (1. + RVTMP2 * q_seri(i, k))
546 guez 3 q_seri(i, k) = q_seri(i, k) + zb
547     ENDDO
548     ENDDO
549 guez 51 ql_seri = 0.
550 guez 3
551 guez 98 frugs = MAX(frugs, 0.000015)
552     zxrugs = sum(frugs * pctsrf, dim = 2)
553 guez 3
554 guez 191 ! Calculs n\'ecessaires au calcul de l'albedo dans l'interface avec
555 guez 118 ! la surface.
556 guez 3
557 guez 118 CALL orbite(REAL(julien), longi, dist)
558 guez 209 CALL zenang(longi, time, dtphys * radpas, mu0, fract)
559 guez 98 albsol = sum(falbe * pctsrf, dim = 2)
560 guez 3
561 guez 90 ! R\'epartition sous maille des flux longwave et shortwave
562     ! R\'epartition du longwave par sous-surface lin\'earis\'ee
563 guez 3
564 guez 98 forall (nsrf = 1: nbsrf)
565 guez 221 fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
566     * (tsol - ftsol(:, nsrf))
567 guez 98 fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
568     END forall
569 guez 3
570     fder = dlw
571    
572 guez 202 CALL clmain(dtphys, pctsrf, t_seri, q_seri, u_seri, v_seri, julien, mu0, &
573     ftsol, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &
574     paprs, play, fsnow, fqsurf, fevap, falbe, fluxlat, rain_fall, &
575 guez 209 snow_fall, fsolsw, fsollw, fder, frugs, agesno, rugoro, d_t_vdf, &
576     d_q_vdf, d_u_vdf, d_v_vdf, d_ts, flux_t, flux_q, flux_u, flux_v, &
577     cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, &
578     v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, &
579     plcl, fqcalving, ffonte, run_off_lic_0)
580 guez 3
581 guez 90 ! Incr\'ementation des flux
582 guez 40
583 guez 206 sens = - sum(flux_t * pctsrf, dim = 2)
584     evap = - sum(flux_q * pctsrf, dim = 2)
585     fder = dlw + dsens + devap
586 guez 3
587     DO k = 1, llm
588     DO i = 1, klon
589     t_seri(i, k) = t_seri(i, k) + d_t_vdf(i, k)
590     q_seri(i, k) = q_seri(i, k) + d_q_vdf(i, k)
591     u_seri(i, k) = u_seri(i, k) + d_u_vdf(i, k)
592     v_seri(i, k) = v_seri(i, k) + d_v_vdf(i, k)
593     ENDDO
594     ENDDO
595    
596 guez 49 ! Update surface temperature:
597 guez 3
598 guez 191 call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
599 guez 202 ftsol = ftsol + d_ts
600 guez 221 tsol = sum(ftsol * pctsrf, dim = 2)
601 guez 208 zxfluxlat = sum(fluxlat * pctsrf, dim = 2)
602     zt2m = sum(t2m * pctsrf, dim = 2)
603     zq2m = sum(q2m * pctsrf, dim = 2)
604     zu10m = sum(u10m * pctsrf, dim = 2)
605     zv10m = sum(v10m * pctsrf, dim = 2)
606     zxffonte = sum(ffonte * pctsrf, dim = 2)
607     zxfqcalving = sum(fqcalving * pctsrf, dim = 2)
608     s_pblh = sum(pblh * pctsrf, dim = 2)
609     s_lcl = sum(plcl * pctsrf, dim = 2)
610     s_capCL = sum(capCL * pctsrf, dim = 2)
611     s_oliqCL = sum(oliqCL * pctsrf, dim = 2)
612     s_cteiCL = sum(cteiCL * pctsrf, dim = 2)
613     s_pblT = sum(pblT * pctsrf, dim = 2)
614     s_therm = sum(therm * pctsrf, dim = 2)
615     s_trmb1 = sum(trmb1 * pctsrf, dim = 2)
616     s_trmb2 = sum(trmb2 * pctsrf, dim = 2)
617     s_trmb3 = sum(trmb3 * pctsrf, dim = 2)
618 guez 3
619 guez 205 ! Si une sous-fraction n'existe pas, elle prend la valeur moyenne :
620 guez 3 DO nsrf = 1, nbsrf
621     DO i = 1, klon
622 guez 205 IF (pctsrf(i, nsrf) < epsfra) then
623 guez 221 ftsol(i, nsrf) = tsol(i)
624 guez 205 t2m(i, nsrf) = zt2m(i)
625     q2m(i, nsrf) = zq2m(i)
626     u10m(i, nsrf) = zu10m(i)
627     v10m(i, nsrf) = zv10m(i)
628     ffonte(i, nsrf) = zxffonte(i)
629     fqcalving(i, nsrf) = zxfqcalving(i)
630     pblh(i, nsrf) = s_pblh(i)
631     plcl(i, nsrf) = s_lcl(i)
632     capCL(i, nsrf) = s_capCL(i)
633     oliqCL(i, nsrf) = s_oliqCL(i)
634     cteiCL(i, nsrf) = s_cteiCL(i)
635     pblT(i, nsrf) = s_pblT(i)
636     therm(i, nsrf) = s_therm(i)
637     trmb1(i, nsrf) = s_trmb1(i)
638     trmb2(i, nsrf) = s_trmb2(i)
639     trmb3(i, nsrf) = s_trmb3(i)
640     end IF
641 guez 3 ENDDO
642     ENDDO
643    
644 guez 98 ! Calculer la dérive du flux infrarouge
645 guez 3
646     DO i = 1, klon
647 guez 221 dlw(i) = - 4. * RSIGMA * tsol(i)**3
648 guez 3 ENDDO
649    
650 guez 190 ! Appeler la convection
651 guez 3
652 guez 182 if (conv_emanuel) then
653 guez 195 CALL concvl(paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, w01, &
654     d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, itop_con, &
655 guez 206 upwd, dnwd, Ma, cape, iflagctrl, qcondc, pmflxr, da, phi, mp)
656 guez 183 snow_con = 0.
657 guez 62 clwcon0 = qcondc
658 guez 71 mfu = upwd + dnwd
659 guez 3
660 guez 207 zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
661     zqsat = zqsat / (1. - retv * zqsat)
662 guez 3
663 guez 103 ! Properties of convective clouds
664 guez 71 clwcon0 = fact_cldcon * clwcon0
665 guez 62 call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
666     rnebcon0)
667 guez 72
668 guez 190 forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
669 guez 72 mfd = 0.
670     pen_u = 0.
671     pen_d = 0.
672     pde_d = 0.
673     pde_u = 0.
674 guez 182 else
675     conv_q = d_q_dyn + d_q_vdf / dtphys
676     conv_t = d_t_dyn + d_t_vdf / dtphys
677     z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
678     CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
679 guez 206 q_seri(:, llm:1:- 1), conv_t, conv_q, - evap, omega, &
680 guez 182 d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
681     mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
682     kdtop, pmflxr, pmflxs)
683     WHERE (rain_con < 0.) rain_con = 0.
684     WHERE (snow_con < 0.) snow_con = 0.
685     ibas_con = llm + 1 - kcbot
686     itop_con = llm + 1 - kctop
687 guez 69 END if
688 guez 3
689     DO k = 1, llm
690     DO i = 1, klon
691     t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
692     q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
693     u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
694     v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
695     ENDDO
696     ENDDO
697    
698 guez 182 IF (.not. conv_emanuel) THEN
699 guez 69 z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
700     z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
701 guez 3 DO k = 1, llm
702     DO i = 1, klon
703 guez 52 IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
704 guez 3 q_seri(i, k) = q_seri(i, k) * z_factor(i)
705     ENDIF
706     ENDDO
707     ENDDO
708     ENDIF
709    
710 guez 90 ! Convection s\`eche (thermiques ou ajustement)
711 guez 3
712 guez 51 d_t_ajs = 0.
713     d_u_ajs = 0.
714     d_v_ajs = 0.
715     d_q_ajs = 0.
716     fm_therm = 0.
717     entr_therm = 0.
718 guez 3
719 guez 47 if (iflag_thermals == 0) then
720     ! Ajustement sec
721     CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)
722 guez 13 t_seri = t_seri + d_t_ajs
723     q_seri = q_seri + d_q_ajs
724 guez 3 else
725 guez 47 call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
726     q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
727 guez 3 endif
728    
729 guez 47 ! Caclul des ratqs
730 guez 3
731 guez 90 ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
732     ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
733 guez 3 if (iflag_cldcon == 1) then
734 guez 51 do k = 1, llm
735     do i = 1, klon
736 guez 3 if(ptconv(i, k)) then
737 guez 70 ratqsc(i, k) = ratqsbas + fact_cldcon &
738     * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
739 guez 3 else
740 guez 51 ratqsc(i, k) = 0.
741 guez 3 endif
742     enddo
743     enddo
744     endif
745    
746 guez 47 ! ratqs stables
747 guez 51 do k = 1, llm
748     do i = 1, klon
749 guez 70 ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
750 guez 190 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
751 guez 3 enddo
752     enddo
753    
754 guez 47 ! ratqs final
755 guez 69 if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
756 guez 47 ! les ratqs sont une conbinaison de ratqss et ratqsc
757     ! ratqs final
758     ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
759     ! relaxation des ratqs
760 guez 70 ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
761 guez 51 ratqs = max(ratqs, ratqsc)
762 guez 3 else
763 guez 47 ! on ne prend que le ratqs stable pour fisrtilp
764 guez 51 ratqs = ratqss
765 guez 3 endif
766    
767 guez 51 CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
768     d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
769     pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
770     psfl, rhcl)
771 guez 3
772     WHERE (rain_lsc < 0) rain_lsc = 0.
773     WHERE (snow_lsc < 0) snow_lsc = 0.
774     DO k = 1, llm
775     DO i = 1, klon
776     t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
777     q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
778     ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
779     cldfra(i, k) = rneb(i, k)
780     IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
781     ENDDO
782     ENDDO
783    
784 guez 47 ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
785 guez 3
786     ! 1. NUAGES CONVECTIFS
787    
788 guez 174 IF (iflag_cldcon <= - 1) THEN
789 guez 62 ! seulement pour Tiedtke
790 guez 51 snow_tiedtke = 0.
791 guez 174 if (iflag_cldcon == - 1) then
792 guez 51 rain_tiedtke = rain_con
793 guez 3 else
794 guez 51 rain_tiedtke = 0.
795     do k = 1, llm
796     do i = 1, klon
797 guez 7 if (d_q_con(i, k) < 0.) then
798 guez 202 rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k) / dtphys &
799     * zmasse(i, k)
800 guez 3 endif
801     enddo
802     enddo
803     endif
804    
805     ! Nuages diagnostiques pour Tiedtke
806 guez 69 CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
807     itop_con, diafra, dialiq)
808 guez 3 DO k = 1, llm
809     DO i = 1, klon
810 guez 51 IF (diafra(i, k) > cldfra(i, k)) THEN
811 guez 3 cldliq(i, k) = dialiq(i, k)
812     cldfra(i, k) = diafra(i, k)
813     ENDIF
814     ENDDO
815     ENDDO
816     ELSE IF (iflag_cldcon == 3) THEN
817 guez 72 ! On prend pour les nuages convectifs le maximum du calcul de
818 guez 90 ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
819 guez 72 ! d'un facteur facttemps.
820     facteur = dtphys * facttemps
821 guez 51 do k = 1, llm
822     do i = 1, klon
823 guez 70 rnebcon(i, k) = rnebcon(i, k) * facteur
824 guez 72 if (rnebcon0(i, k) * clwcon0(i, k) &
825     > rnebcon(i, k) * clwcon(i, k)) then
826 guez 51 rnebcon(i, k) = rnebcon0(i, k)
827     clwcon(i, k) = clwcon0(i, k)
828 guez 3 endif
829     enddo
830     enddo
831    
832 guez 47 ! On prend la somme des fractions nuageuses et des contenus en eau
833 guez 51 cldfra = min(max(cldfra, rnebcon), 1.)
834 guez 202 cldliq = cldliq + rnebcon * clwcon
835 guez 3 ENDIF
836    
837 guez 51 ! 2. Nuages stratiformes
838 guez 3
839     IF (ok_stratus) THEN
840 guez 47 CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
841 guez 3 DO k = 1, llm
842     DO i = 1, klon
843 guez 51 IF (diafra(i, k) > cldfra(i, k)) THEN
844 guez 3 cldliq(i, k) = dialiq(i, k)
845     cldfra(i, k) = diafra(i, k)
846     ENDIF
847     ENDDO
848     ENDDO
849     ENDIF
850    
851     ! Precipitation totale
852     DO i = 1, klon
853     rain_fall(i) = rain_con(i) + rain_lsc(i)
854     snow_fall(i) = snow_con(i) + snow_lsc(i)
855     ENDDO
856    
857 guez 90 ! Humidit\'e relative pour diagnostic :
858 guez 3 DO k = 1, llm
859     DO i = 1, klon
860     zx_t = t_seri(i, k)
861 guez 207 zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t) / play(i, k)
862     zx_qs = MIN(0.5, zx_qs)
863     zcor = 1. / (1. - retv * zx_qs)
864     zx_qs = zx_qs * zcor
865 guez 202 zx_rh(i, k) = q_seri(i, k) / zx_qs
866 guez 51 zqsat(i, k) = zx_qs
867 guez 3 ENDDO
868     ENDDO
869 guez 52
870 guez 97 ! Param\`etres optiques des nuages et quelques param\`etres pour
871     ! diagnostics :
872 guez 3 if (ok_newmicro) then
873 guez 69 CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
874 guez 217 cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc)
875 guez 3 else
876 guez 52 CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
877 guez 217 cldl, cldm, cldt, cldq)
878 guez 3 endif
879    
880 guez 154 IF (MOD(itap - 1, radpas) == 0) THEN
881 guez 212 wo = ozonecm(REAL(julien), paprs)
882 guez 155 albsol = sum(falbe * pctsrf, dim = 2)
883 guez 221 CALL radlwsw(dist, mu0, fract, paprs, play, tsol, albsol, t_seri, &
884 guez 155 q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
885     radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
886     toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
887 guez 217 swup0, swup, ok_ade, topswad, solswad)
888 guez 3 ENDIF
889 guez 118
890 guez 3 ! Ajouter la tendance des rayonnements (tous les pas)
891     DO k = 1, llm
892     DO i = 1, klon
893 guez 202 t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys &
894     / 86400.
895 guez 3 ENDDO
896     ENDDO
897    
898 guez 90 ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
899 guez 3 DO i = 1, klon
900     bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
901     ENDDO
902    
903 guez 90 ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
904 guez 3
905     IF (ok_orodr) THEN
906 guez 174 ! S\'election des points pour lesquels le sch\'ema est actif :
907 guez 51 igwd = 0
908     DO i = 1, klon
909     itest(i) = 0
910 guez 174 IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
911 guez 51 itest(i) = 1
912     igwd = igwd + 1
913 guez 3 ENDIF
914     ENDDO
915    
916 guez 51 CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
917 guez 150 zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
918     zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
919 guez 3
920 guez 47 ! ajout des tendances
921 guez 3 DO k = 1, llm
922     DO i = 1, klon
923     t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
924     u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
925     v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
926     ENDDO
927     ENDDO
928 guez 13 ENDIF
929 guez 3
930     IF (ok_orolf) THEN
931 guez 90 ! S\'election des points pour lesquels le sch\'ema est actif :
932 guez 51 igwd = 0
933     DO i = 1, klon
934     itest(i) = 0
935 guez 174 IF (zpic(i) - zmea(i) > 100.) THEN
936 guez 51 itest(i) = 1
937     igwd = igwd + 1
938 guez 3 ENDIF
939     ENDDO
940    
941 guez 47 CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
942     itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &
943 guez 3 d_t_lif, d_u_lif, d_v_lif)
944    
945 guez 51 ! Ajout des tendances :
946 guez 3 DO k = 1, llm
947     DO i = 1, klon
948     t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
949     u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
950     v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
951     ENDDO
952     ENDDO
953 guez 49 ENDIF
954 guez 3
955 guez 90 ! Stress n\'ecessaires : toute la physique
956 guez 3
957     DO i = 1, klon
958 guez 51 zustrph(i) = 0.
959     zvstrph(i) = 0.
960 guez 3 ENDDO
961     DO k = 1, llm
962     DO i = 1, klon
963 guez 62 zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
964     * zmasse(i, k)
965     zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
966     * zmasse(i, k)
967 guez 3 ENDDO
968     ENDDO
969    
970 guez 171 CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
971     zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
972 guez 3
973 guez 47 ! Calcul des tendances traceurs
974 guez 202 call phytrac(julien, time, firstcal, lafin, dtphys, t, paprs, play, mfu, &
975     mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, &
976     pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, tr_seri, &
977     zmasse, ncid_startphy)
978 guez 3
979 guez 190 IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &
980     pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &
981 guez 221 frac_impa, frac_nucl, pphis, airephy)
982 guez 3
983     ! Calculer le transport de l'eau et de l'energie (diagnostique)
984 guez 171 CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
985 guez 3
986 guez 31 ! diag. bilKP
987 guez 3
988 guez 178 CALL transp_lay(paprs, t_seri, q_seri, u_seri, v_seri, zphi, &
989 guez 3 ve_lay, vq_lay, ue_lay, uq_lay)
990    
991     ! Accumuler les variables a stocker dans les fichiers histoire:
992    
993 guez 200 ! conversion Ec en énergie thermique
994 guez 3 DO k = 1, llm
995     DO i = 1, klon
996 guez 213 d_t_ec(i, k) = 0.5 / (RCPD * (1. + RVTMP2 * q_seri(i, k))) &
997 guez 51 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
998     t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
999     d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1000 guez 3 END DO
1001     END DO
1002 guez 51
1003 guez 47 ! SORTIES
1004 guez 3
1005 guez 69 ! prw = eau precipitable
1006 guez 3 DO i = 1, klon
1007     prw(i) = 0.
1008     DO k = 1, llm
1009 guez 202 prw(i) = prw(i) + q_seri(i, k) * zmasse(i, k)
1010 guez 3 ENDDO
1011     ENDDO
1012    
1013     ! Convertir les incrementations en tendances
1014    
1015     DO k = 1, llm
1016     DO i = 1, klon
1017 guez 49 d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1018     d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1019     d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1020     d_qx(i, k, ivap) = (q_seri(i, k) - qx(i, k, ivap)) / dtphys
1021     d_qx(i, k, iliq) = (ql_seri(i, k) - qx(i, k, iliq)) / dtphys
1022 guez 3 ENDDO
1023     ENDDO
1024    
1025 guez 98 DO iq = 3, nqmx
1026     DO k = 1, llm
1027     DO i = 1, klon
1028 guez 174 d_qx(i, k, iq) = (tr_seri(i, k, iq - 2) - qx(i, k, iq)) / dtphys
1029 guez 3 ENDDO
1030     ENDDO
1031 guez 98 ENDDO
1032 guez 3
1033     ! Sauvegarder les valeurs de t et q a la fin de la physique:
1034     DO k = 1, llm
1035     DO i = 1, klon
1036     t_ancien(i, k) = t_seri(i, k)
1037     q_ancien(i, k) = q_seri(i, k)
1038     ENDDO
1039     ENDDO
1040    
1041 guez 191 CALL histwrite_phy("phis", pphis)
1042     CALL histwrite_phy("aire", airephy)
1043     CALL histwrite_phy("psol", paprs(:, 1))
1044     CALL histwrite_phy("precip", rain_fall + snow_fall)
1045     CALL histwrite_phy("plul", rain_lsc + snow_lsc)
1046     CALL histwrite_phy("pluc", rain_con + snow_con)
1047 guez 221 CALL histwrite_phy("tsol", tsol)
1048 guez 191 CALL histwrite_phy("t2m", zt2m)
1049     CALL histwrite_phy("q2m", zq2m)
1050     CALL histwrite_phy("u10m", zu10m)
1051     CALL histwrite_phy("v10m", zv10m)
1052     CALL histwrite_phy("snow", snow_fall)
1053     CALL histwrite_phy("cdrm", cdragm)
1054     CALL histwrite_phy("cdrh", cdragh)
1055     CALL histwrite_phy("topl", toplw)
1056     CALL histwrite_phy("evap", evap)
1057     CALL histwrite_phy("sols", solsw)
1058     CALL histwrite_phy("soll", sollw)
1059     CALL histwrite_phy("solldown", sollwdown)
1060     CALL histwrite_phy("bils", bils)
1061     CALL histwrite_phy("sens", - sens)
1062     CALL histwrite_phy("fder", fder)
1063     CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
1064     CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1065     CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1066     CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
1067 guez 3
1068 guez 191 DO nsrf = 1, nbsrf
1069 guez 202 CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf) * 100.)
1070 guez 191 CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1071 guez 206 CALL histwrite_phy("sens_"//clnsurf(nsrf), flux_t(:, nsrf))
1072 guez 191 CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1073     CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1074 guez 206 CALL histwrite_phy("taux_"//clnsurf(nsrf), flux_u(:, nsrf))
1075     CALL histwrite_phy("tauy_"//clnsurf(nsrf), flux_v(:, nsrf))
1076 guez 191 CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1077     CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1078     END DO
1079    
1080     CALL histwrite_phy("albs", albsol)
1081 guez 212 CALL histwrite_phy("tro3", wo * dobson_u * 1e3 / zmasse / rmo3 * md)
1082 guez 191 CALL histwrite_phy("rugs", zxrugs)
1083     CALL histwrite_phy("s_pblh", s_pblh)
1084     CALL histwrite_phy("s_pblt", s_pblt)
1085     CALL histwrite_phy("s_lcl", s_lcl)
1086     CALL histwrite_phy("s_capCL", s_capCL)
1087     CALL histwrite_phy("s_oliqCL", s_oliqCL)
1088     CALL histwrite_phy("s_cteiCL", s_cteiCL)
1089     CALL histwrite_phy("s_therm", s_therm)
1090     CALL histwrite_phy("s_trmb1", s_trmb1)
1091     CALL histwrite_phy("s_trmb2", s_trmb2)
1092     CALL histwrite_phy("s_trmb3", s_trmb3)
1093 guez 206
1094     if (conv_emanuel) then
1095     CALL histwrite_phy("ptop", ema_pct)
1096     CALL histwrite_phy("dnwd0", - mp)
1097     end if
1098    
1099 guez 191 CALL histwrite_phy("temp", t_seri)
1100     CALL histwrite_phy("vitu", u_seri)
1101     CALL histwrite_phy("vitv", v_seri)
1102     CALL histwrite_phy("geop", zphi)
1103     CALL histwrite_phy("pres", play)
1104     CALL histwrite_phy("dtvdf", d_t_vdf)
1105     CALL histwrite_phy("dqvdf", d_q_vdf)
1106     CALL histwrite_phy("rhum", zx_rh)
1107 guez 213 CALL histwrite_phy("d_t_ec", d_t_ec)
1108     CALL histwrite_phy("dtsw0", heat0 / 86400.)
1109     CALL histwrite_phy("dtlw0", - cool0 / 86400.)
1110 guez 215 CALL histwrite_phy("msnow", sum(fsnow * pctsrf, dim = 2))
1111 guez 221 call histwrite_phy("qsurf", sum(fqsurf * pctsrf, dim = 2))
1112 guez 191
1113     if (ok_instan) call histsync(nid_ins)
1114    
1115 guez 157 IF (lafin) then
1116     call NF95_CLOSE(ncid_startphy)
1117 guez 175 CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1118 guez 157 fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1119     radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1120     t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1121     w01)
1122     end IF
1123 guez 3
1124 guez 35 firstcal = .FALSE.
1125    
1126 guez 3 END SUBROUTINE physiq
1127    
1128     end module physiq_m
  ViewVC Help
Powered by ViewVC 1.1.21