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

Annotation of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 215 - (hide annotations)
Tue Mar 28 12:46:28 2017 UTC (7 years, 1 month ago) by guez
Original Path: trunk/Sources/phylmd/physiq.f
File size: 40577 byte(s) 
size(snow) is now knon in interfsurf_hq.

Renamed snow to fsnow in clmain, same name as corresponding actual
argument. We can then rename ysnow to simply snow in clmain, same name
as corresponding dummy argument of clqh. No need to initialize local
snow to 0 since it is only used with indices 1:knon and already
initialized from fsnow for each type of surface. If there is no point
for a given type of surface, fsnow should be reset to 0 for this
type. We need to give a valid value to fsnow in this case even if it
will be multiplied by pctsrf = 0 in physiq.

In physiq, no need for intermediate zxsnow for output.

Removed unused arguments tsurf, p1lay, beta, coef1lay, ps, t1lay,
q1lay, u1lay, v1lay, petAcoef, peqAcoef, petBcoef, peqBcoef of
fonte_neige, with unused computations of zx_qs and zcor. (Same was
done in LMDZ.)
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 207 USE fcttre, ONLY: foeew, qsatl, qsats
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 155 REAL, save:: albsol(klon) ! albedo du sol total visible
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 215 REAL zxqsurf(klon), zxfluxlat(klon)
278 guez 3
279 guez 118 REAL dist, mu0(klon), fract(klon)
280     real longi
281 guez 3 REAL z_avant(klon), z_apres(klon), z_factor(klon)
282 guez 205 REAL zb
283 guez 103 REAL zx_t, zx_qs, zcor
284 guez 3 real zqsat(klon, llm)
285     INTEGER i, k, iq, nsrf
286     REAL zphi(klon, llm)
287    
288 guez 200 ! cf. Anne Mathieu, variables pour la couche limite atmosphérique (hbtm)
289 guez 3
290 guez 49 REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
291     REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
292     REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
293     REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
294     REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
295 guez 207 REAL, SAVE:: pblt(klon, nbsrf) ! T \`a la hauteur de couche limite
296 guez 49 REAL, SAVE:: therm(klon, nbsrf)
297     REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
298 guez 190 REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
299 guez 49 REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
300 guez 186 ! Grandeurs de sorties
301 guez 3 REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
302     REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
303     REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
304     REAL s_trmb3(klon)
305    
306 guez 175 ! Variables pour la convection de K. Emanuel :
307 guez 3
308 guez 47 REAL upwd(klon, llm) ! saturated updraft mass flux
309     REAL dnwd(klon, llm) ! saturated downdraft mass flux
310 guez 205 REAL, save:: cape(klon)
311 guez 3
312 guez 47 INTEGER iflagctrl(klon) ! flag fonctionnement de convect
313 guez 3
314     ! Variables du changement
315    
316     ! con: convection
317 guez 51 ! lsc: large scale condensation
318 guez 3 ! ajs: ajustement sec
319 guez 90 ! eva: \'evaporation de l'eau liquide nuageuse
320 guez 51 ! vdf: vertical diffusion in boundary layer
321 guez 3 REAL d_t_con(klon, llm), d_q_con(klon, llm)
322 guez 205 REAL, save:: d_u_con(klon, llm), d_v_con(klon, llm)
323 guez 3 REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
324     REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
325     REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
326     REAL rneb(klon, llm)
327    
328 guez 71 REAL mfu(klon, llm), mfd(klon, llm)
329 guez 3 REAL pen_u(klon, llm), pen_d(klon, llm)
330     REAL pde_u(klon, llm), pde_d(klon, llm)
331     INTEGER kcbot(klon), kctop(klon), kdtop(klon)
332 guez 51 REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
333     REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
334 guez 3
335 guez 62 INTEGER, save:: ibas_con(klon), itop_con(klon)
336 guez 183 real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
337 guez 3
338 guez 205 REAL, save:: rain_con(klon)
339     real rain_lsc(klon)
340 guez 183 REAL, save:: snow_con(klon) ! neige (mm / s)
341 guez 180 real snow_lsc(klon)
342 guez 3 REAL d_ts(klon, nbsrf)
343    
344     REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
345     REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
346    
347     REAL d_u_oro(klon, llm), d_v_oro(klon, llm)
348     REAL d_t_oro(klon, llm)
349     REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
350     REAL d_t_lif(klon, llm)
351    
352 guez 68 REAL, save:: ratqs(klon, llm)
353     real ratqss(klon, llm), ratqsc(klon, llm)
354     real:: ratqsbas = 0.01, ratqshaut = 0.3
355 guez 3
356     ! Parametres lies au nouveau schema de nuages (SB, PDF)
357 guez 68 real:: fact_cldcon = 0.375
358     real:: facttemps = 1.e-4
359     logical:: ok_newmicro = .true.
360 guez 3 real facteur
361    
362 guez 68 integer:: iflag_cldcon = 1
363 guez 3 logical ptconv(klon, llm)
364    
365 guez 175 ! Variables pour effectuer les appels en s\'erie :
366 guez 3
367     REAL t_seri(klon, llm), q_seri(klon, llm)
368 guez 98 REAL ql_seri(klon, llm)
369 guez 3 REAL u_seri(klon, llm), v_seri(klon, llm)
370 guez 98 REAL tr_seri(klon, llm, nqmx - 2)
371 guez 3
372     REAL zx_rh(klon, llm)
373    
374     REAL zustrdr(klon), zvstrdr(klon)
375     REAL zustrli(klon), zvstrli(klon)
376     REAL zustrph(klon), zvstrph(klon)
377     REAL aam, torsfc
378    
379     REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
380     REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
381     REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
382     REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
383    
384     real date0
385     REAL ztsol(klon)
386 guez 51
387 guez 202 REAL d_t_ec(klon, llm)
388 guez 213 ! tendance due \`a la conversion d'\'energie cin\'etique en
389     ! énergie thermique
390 guez 200
391 guez 205 REAL, save:: t2m(klon, nbsrf), q2m(klon, nbsrf)
392     ! temperature and humidity at 2 m
393    
394     REAL, save:: u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m
395 guez 207 REAL zt2m(klon), zq2m(klon) ! température, humidité 2 m moyenne sur 1 maille
396     REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes sur 1 maille
397 guez 3
398 guez 69 ! Aerosol effects:
399    
400 guez 202 REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g / m3)
401 guez 69
402 guez 49 REAL, save:: sulfate_pi(klon, llm)
403 guez 202 ! SO4 aerosol concentration, in \mu g / m3, pre-industrial value
404 guez 3
405     REAL cldtaupi(klon, llm)
406 guez 191 ! cloud optical thickness for pre-industrial aerosols
407 guez 3
408 guez 47 REAL re(klon, llm) ! Cloud droplet effective radius
409     REAL fl(klon, llm) ! denominator of re
410 guez 3
411     ! Aerosol optical properties
412 guez 68 REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)
413     REAL, save:: cg_ae(klon, llm, 2)
414 guez 3
415 guez 205 REAL, save:: topswad(klon), solswad(klon) ! aerosol direct effect
416     REAL, save:: topswai(klon), solswai(klon) ! aerosol indirect effect
417 guez 3
418 guez 68 LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
419     LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect
420 guez 3
421 guez 68 REAL:: bl95_b0 = 2., bl95_b1 = 0.2
422 guez 69 ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
423     ! B). They link cloud droplet number concentration to aerosol mass
424     ! concentration.
425 guez 68
426 guez 190 real zmasse(klon, llm)
427 guez 17 ! (column-density of mass of air in a cell, in kg m-2)
428    
429 guez 191 integer, save:: ncid_startphy
430 guez 17
431 guez 204 namelist /physiq_nml/ fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, &
432     ratqsbas, ratqshaut, ok_ade, ok_aie, bl95_b0, bl95_b1, &
433     iflag_thermals, nsplit_thermals
434 guez 68
435 guez 3 !----------------------------------------------------------------
436    
437 guez 69 IF (nqmx < 2) CALL abort_gcm('physiq', &
438 guez 171 'eaux vapeur et liquide sont indispensables')
439 guez 3
440 guez 7 test_firstcal: IF (firstcal) THEN
441 guez 47 ! initialiser
442 guez 51 u10m = 0.
443     v10m = 0.
444     t2m = 0.
445     q2m = 0.
446     ffonte = 0.
447     fqcalving = 0.
448     piz_ae = 0.
449     tau_ae = 0.
450     cg_ae = 0.
451 guez 98 rain_con = 0.
452     snow_con = 0.
453     topswai = 0.
454     topswad = 0.
455     solswai = 0.
456     solswad = 0.
457 guez 3
458 guez 72 d_u_con = 0.
459     d_v_con = 0.
460     rnebcon0 = 0.
461     clwcon0 = 0.
462     rnebcon = 0.
463     clwcon = 0.
464 guez 3
465 guez 47 pblh =0. ! Hauteur de couche limite
466     plcl =0. ! Niveau de condensation de la CLA
467     capCL =0. ! CAPE de couche limite
468     oliqCL =0. ! eau_liqu integree de couche limite
469     cteiCL =0. ! cloud top instab. crit. couche limite
470 guez 207 pblt =0.
471 guez 47 therm =0.
472     trmb1 =0. ! deep_cape
473 guez 190 trmb2 =0. ! inhibition
474 guez 47 trmb3 =0. ! Point Omega
475 guez 3
476 guez 68 iflag_thermals = 0
477     nsplit_thermals = 1
478     print *, "Enter namelist 'physiq_nml'."
479     read(unit=*, nml=physiq_nml)
480     write(unit_nml, nml=physiq_nml)
481    
482     call conf_phys
483 guez 3
484     ! Initialiser les compteurs:
485    
486     frugs = 0.
487 guez 191 CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
488     fevap, rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &
489     agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
490     q_ancien, ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
491     w01, ncid_startphy)
492 guez 3
493 guez 47 ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
494 guez 69 q2 = 1e-8
495 guez 3
496 guez 154 radpas = lmt_pas / nbapp_rad
497 guez 191 print *, "radpas = ", radpas
498 guez 154
499 guez 90 ! Initialisation pour le sch\'ema de convection d'Emanuel :
500 guez 182 IF (conv_emanuel) THEN
501 guez 69 ibas_con = 1
502     itop_con = 1
503 guez 3 ENDIF
504    
505     IF (ok_orodr) THEN
506 guez 13 rugoro = MAX(1e-5, zstd * zsig / 2)
507 guez 54 CALL SUGWD(paprs, play)
508 guez 13 else
509     rugoro = 0.
510 guez 3 ENDIF
511    
512 guez 202 ecrit_ins = NINT(ecrit_ins / dtphys)
513 guez 3
514 guez 47 ! Initialisation des sorties
515 guez 3
516 guez 191 call ini_histins(dtphys)
517 guez 129 CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
518 guez 69 ! Positionner date0 pour initialisation de ORCHIDEE
519     print *, 'physiq date0: ', date0
520 guez 202 CALL phyredem0
521 guez 7 ENDIF test_firstcal
522 guez 3
523 guez 91 ! We will modify variables *_seri and we will not touch variables
524 guez 98 ! u, v, t, qx:
525     t_seri = t
526     u_seri = u
527     v_seri = v
528     q_seri = qx(:, :, ivap)
529     ql_seri = qx(:, :, iliq)
530 guez 157 tr_seri = qx(:, :, 3:nqmx)
531 guez 3
532 guez 98 ztsol = sum(ftsol * pctsrf, dim = 2)
533 guez 3
534 guez 51 ! Diagnostic de la tendance dynamique :
535 guez 3 IF (ancien_ok) THEN
536     DO k = 1, llm
537     DO i = 1, klon
538 guez 49 d_t_dyn(i, k) = (t_seri(i, k) - t_ancien(i, k)) / dtphys
539     d_q_dyn(i, k) = (q_seri(i, k) - q_ancien(i, k)) / dtphys
540 guez 3 ENDDO
541     ENDDO
542     ELSE
543     DO k = 1, llm
544     DO i = 1, klon
545 guez 72 d_t_dyn(i, k) = 0.
546     d_q_dyn(i, k) = 0.
547 guez 3 ENDDO
548     ENDDO
549     ancien_ok = .TRUE.
550     ENDIF
551    
552     ! Ajouter le geopotentiel du sol:
553     DO k = 1, llm
554     DO i = 1, klon
555     zphi(i, k) = pphi(i, k) + pphis(i)
556     ENDDO
557     ENDDO
558    
559 guez 49 ! Check temperatures:
560 guez 3 CALL hgardfou(t_seri, ftsol)
561    
562 guez 191 call increment_itap
563 guez 130 julien = MOD(dayvrai, 360)
564 guez 3 if (julien == 0) julien = 360
565    
566 guez 103 forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
567 guez 17
568 guez 90 ! \'Evaporation de l'eau liquide nuageuse :
569 guez 51 DO k = 1, llm
570 guez 3 DO i = 1, klon
571 guez 51 zb = MAX(0., ql_seri(i, k))
572     t_seri(i, k) = t_seri(i, k) &
573     - zb * RLVTT / RCPD / (1. + RVTMP2 * q_seri(i, k))
574 guez 3 q_seri(i, k) = q_seri(i, k) + zb
575     ENDDO
576     ENDDO
577 guez 51 ql_seri = 0.
578 guez 3
579 guez 98 frugs = MAX(frugs, 0.000015)
580     zxrugs = sum(frugs * pctsrf, dim = 2)
581 guez 3
582 guez 191 ! Calculs n\'ecessaires au calcul de l'albedo dans l'interface avec
583 guez 118 ! la surface.
584 guez 3
585 guez 118 CALL orbite(REAL(julien), longi, dist)
586 guez 209 CALL zenang(longi, time, dtphys * radpas, mu0, fract)
587 guez 3
588 guez 47 ! Calcul de l'abedo moyen par maille
589 guez 98 albsol = sum(falbe * pctsrf, dim = 2)
590 guez 3
591 guez 90 ! R\'epartition sous maille des flux longwave et shortwave
592     ! R\'epartition du longwave par sous-surface lin\'earis\'ee
593 guez 3
594 guez 98 forall (nsrf = 1: nbsrf)
595     fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
596     * (ztsol - ftsol(:, nsrf))
597     fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
598     END forall
599 guez 3
600     fder = dlw
601    
602 guez 202 CALL clmain(dtphys, pctsrf, t_seri, q_seri, u_seri, v_seri, julien, mu0, &
603     ftsol, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, ftsoil, qsol, &
604     paprs, play, fsnow, fqsurf, fevap, falbe, fluxlat, rain_fall, &
605 guez 209 snow_fall, fsolsw, fsollw, fder, frugs, agesno, rugoro, d_t_vdf, &
606     d_q_vdf, d_u_vdf, d_v_vdf, d_ts, flux_t, flux_q, flux_u, flux_v, &
607     cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, yv1, t2m, q2m, u10m, &
608     v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, trmb1, trmb2, trmb3, &
609     plcl, fqcalving, ffonte, run_off_lic_0)
610 guez 3
611 guez 90 ! Incr\'ementation des flux
612 guez 40
613 guez 206 sens = - sum(flux_t * pctsrf, dim = 2)
614     evap = - sum(flux_q * pctsrf, dim = 2)
615     fder = dlw + dsens + devap
616 guez 3
617     DO k = 1, llm
618     DO i = 1, klon
619     t_seri(i, k) = t_seri(i, k) + d_t_vdf(i, k)
620     q_seri(i, k) = q_seri(i, k) + d_q_vdf(i, k)
621     u_seri(i, k) = u_seri(i, k) + d_u_vdf(i, k)
622     v_seri(i, k) = v_seri(i, k) + d_v_vdf(i, k)
623     ENDDO
624     ENDDO
625    
626 guez 49 ! Update surface temperature:
627 guez 3
628 guez 191 call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
629 guez 202 ftsol = ftsol + d_ts
630 guez 205 ztsol = sum(ftsol * pctsrf, dim = 2)
631 guez 208 zxfluxlat = sum(fluxlat * pctsrf, dim = 2)
632     zt2m = sum(t2m * pctsrf, dim = 2)
633     zq2m = sum(q2m * pctsrf, dim = 2)
634     zu10m = sum(u10m * pctsrf, dim = 2)
635     zv10m = sum(v10m * pctsrf, dim = 2)
636     zxffonte = sum(ffonte * pctsrf, dim = 2)
637     zxfqcalving = sum(fqcalving * pctsrf, dim = 2)
638     s_pblh = sum(pblh * pctsrf, dim = 2)
639     s_lcl = sum(plcl * pctsrf, dim = 2)
640     s_capCL = sum(capCL * pctsrf, dim = 2)
641     s_oliqCL = sum(oliqCL * pctsrf, dim = 2)
642     s_cteiCL = sum(cteiCL * pctsrf, dim = 2)
643     s_pblT = sum(pblT * pctsrf, dim = 2)
644     s_therm = sum(therm * pctsrf, dim = 2)
645     s_trmb1 = sum(trmb1 * pctsrf, dim = 2)
646     s_trmb2 = sum(trmb2 * pctsrf, dim = 2)
647     s_trmb3 = sum(trmb3 * pctsrf, dim = 2)
648 guez 3
649 guez 205 ! Si une sous-fraction n'existe pas, elle prend la valeur moyenne :
650 guez 3 DO nsrf = 1, nbsrf
651     DO i = 1, klon
652 guez 205 IF (pctsrf(i, nsrf) < epsfra) then
653     ftsol(i, nsrf) = ztsol(i)
654     t2m(i, nsrf) = zt2m(i)
655     q2m(i, nsrf) = zq2m(i)
656     u10m(i, nsrf) = zu10m(i)
657     v10m(i, nsrf) = zv10m(i)
658     ffonte(i, nsrf) = zxffonte(i)
659     fqcalving(i, nsrf) = zxfqcalving(i)
660     pblh(i, nsrf) = s_pblh(i)
661     plcl(i, nsrf) = s_lcl(i)
662     capCL(i, nsrf) = s_capCL(i)
663     oliqCL(i, nsrf) = s_oliqCL(i)
664     cteiCL(i, nsrf) = s_cteiCL(i)
665     pblT(i, nsrf) = s_pblT(i)
666     therm(i, nsrf) = s_therm(i)
667     trmb1(i, nsrf) = s_trmb1(i)
668     trmb2(i, nsrf) = s_trmb2(i)
669     trmb3(i, nsrf) = s_trmb3(i)
670     end IF
671 guez 3 ENDDO
672     ENDDO
673    
674 guez 98 ! Calculer la dérive du flux infrarouge
675 guez 3
676     DO i = 1, klon
677 guez 205 dlw(i) = - 4. * RSIGMA * ztsol(i)**3
678 guez 3 ENDDO
679    
680 guez 190 ! Appeler la convection
681 guez 3
682 guez 182 if (conv_emanuel) then
683 guez 195 CALL concvl(paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, w01, &
684     d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, itop_con, &
685 guez 206 upwd, dnwd, Ma, cape, iflagctrl, qcondc, pmflxr, da, phi, mp)
686 guez 183 snow_con = 0.
687 guez 62 clwcon0 = qcondc
688 guez 71 mfu = upwd + dnwd
689 guez 3
690 guez 207 zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
691     zqsat = zqsat / (1. - retv * zqsat)
692 guez 3
693 guez 103 ! Properties of convective clouds
694 guez 71 clwcon0 = fact_cldcon * clwcon0
695 guez 62 call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
696     rnebcon0)
697 guez 72
698 guez 190 forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
699 guez 72 mfd = 0.
700     pen_u = 0.
701     pen_d = 0.
702     pde_d = 0.
703     pde_u = 0.
704 guez 182 else
705     conv_q = d_q_dyn + d_q_vdf / dtphys
706     conv_t = d_t_dyn + d_t_vdf / dtphys
707     z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
708     CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
709 guez 206 q_seri(:, llm:1:- 1), conv_t, conv_q, - evap, omega, &
710 guez 182 d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
711     mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
712     kdtop, pmflxr, pmflxs)
713     WHERE (rain_con < 0.) rain_con = 0.
714     WHERE (snow_con < 0.) snow_con = 0.
715     ibas_con = llm + 1 - kcbot
716     itop_con = llm + 1 - kctop
717 guez 69 END if
718 guez 3
719     DO k = 1, llm
720     DO i = 1, klon
721     t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
722     q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
723     u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
724     v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
725     ENDDO
726     ENDDO
727    
728 guez 182 IF (.not. conv_emanuel) THEN
729 guez 69 z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
730     z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
731 guez 3 DO k = 1, llm
732     DO i = 1, klon
733 guez 52 IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
734 guez 3 q_seri(i, k) = q_seri(i, k) * z_factor(i)
735     ENDIF
736     ENDDO
737     ENDDO
738     ENDIF
739    
740 guez 90 ! Convection s\`eche (thermiques ou ajustement)
741 guez 3
742 guez 51 d_t_ajs = 0.
743     d_u_ajs = 0.
744     d_v_ajs = 0.
745     d_q_ajs = 0.
746     fm_therm = 0.
747     entr_therm = 0.
748 guez 3
749 guez 47 if (iflag_thermals == 0) then
750     ! Ajustement sec
751     CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)
752 guez 13 t_seri = t_seri + d_t_ajs
753     q_seri = q_seri + d_q_ajs
754 guez 3 else
755 guez 47 call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
756     q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
757 guez 3 endif
758    
759 guez 47 ! Caclul des ratqs
760 guez 3
761 guez 90 ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
762     ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
763 guez 3 if (iflag_cldcon == 1) then
764 guez 51 do k = 1, llm
765     do i = 1, klon
766 guez 3 if(ptconv(i, k)) then
767 guez 70 ratqsc(i, k) = ratqsbas + fact_cldcon &
768     * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
769 guez 3 else
770 guez 51 ratqsc(i, k) = 0.
771 guez 3 endif
772     enddo
773     enddo
774     endif
775    
776 guez 47 ! ratqs stables
777 guez 51 do k = 1, llm
778     do i = 1, klon
779 guez 70 ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
780 guez 190 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
781 guez 3 enddo
782     enddo
783    
784 guez 47 ! ratqs final
785 guez 69 if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
786 guez 47 ! les ratqs sont une conbinaison de ratqss et ratqsc
787     ! ratqs final
788     ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
789     ! relaxation des ratqs
790 guez 70 ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
791 guez 51 ratqs = max(ratqs, ratqsc)
792 guez 3 else
793 guez 47 ! on ne prend que le ratqs stable pour fisrtilp
794 guez 51 ratqs = ratqss
795 guez 3 endif
796    
797 guez 51 CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
798     d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
799     pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
800     psfl, rhcl)
801 guez 3
802     WHERE (rain_lsc < 0) rain_lsc = 0.
803     WHERE (snow_lsc < 0) snow_lsc = 0.
804     DO k = 1, llm
805     DO i = 1, klon
806     t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
807     q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
808     ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
809     cldfra(i, k) = rneb(i, k)
810     IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
811     ENDDO
812     ENDDO
813    
814 guez 47 ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
815 guez 3
816     ! 1. NUAGES CONVECTIFS
817    
818 guez 174 IF (iflag_cldcon <= - 1) THEN
819 guez 62 ! seulement pour Tiedtke
820 guez 51 snow_tiedtke = 0.
821 guez 174 if (iflag_cldcon == - 1) then
822 guez 51 rain_tiedtke = rain_con
823 guez 3 else
824 guez 51 rain_tiedtke = 0.
825     do k = 1, llm
826     do i = 1, klon
827 guez 7 if (d_q_con(i, k) < 0.) then
828 guez 202 rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k) / dtphys &
829     * zmasse(i, k)
830 guez 3 endif
831     enddo
832     enddo
833     endif
834    
835     ! Nuages diagnostiques pour Tiedtke
836 guez 69 CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
837     itop_con, diafra, dialiq)
838 guez 3 DO k = 1, llm
839     DO i = 1, klon
840 guez 51 IF (diafra(i, k) > cldfra(i, k)) THEN
841 guez 3 cldliq(i, k) = dialiq(i, k)
842     cldfra(i, k) = diafra(i, k)
843     ENDIF
844     ENDDO
845     ENDDO
846     ELSE IF (iflag_cldcon == 3) THEN
847 guez 72 ! On prend pour les nuages convectifs le maximum du calcul de
848 guez 90 ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
849 guez 72 ! d'un facteur facttemps.
850     facteur = dtphys * facttemps
851 guez 51 do k = 1, llm
852     do i = 1, klon
853 guez 70 rnebcon(i, k) = rnebcon(i, k) * facteur
854 guez 72 if (rnebcon0(i, k) * clwcon0(i, k) &
855     > rnebcon(i, k) * clwcon(i, k)) then
856 guez 51 rnebcon(i, k) = rnebcon0(i, k)
857     clwcon(i, k) = clwcon0(i, k)
858 guez 3 endif
859     enddo
860     enddo
861    
862 guez 47 ! On prend la somme des fractions nuageuses et des contenus en eau
863 guez 51 cldfra = min(max(cldfra, rnebcon), 1.)
864 guez 202 cldliq = cldliq + rnebcon * clwcon
865 guez 3 ENDIF
866    
867 guez 51 ! 2. Nuages stratiformes
868 guez 3
869     IF (ok_stratus) THEN
870 guez 47 CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
871 guez 3 DO k = 1, llm
872     DO i = 1, klon
873 guez 51 IF (diafra(i, k) > cldfra(i, k)) THEN
874 guez 3 cldliq(i, k) = dialiq(i, k)
875     cldfra(i, k) = diafra(i, k)
876     ENDIF
877     ENDDO
878     ENDDO
879     ENDIF
880    
881     ! Precipitation totale
882     DO i = 1, klon
883     rain_fall(i) = rain_con(i) + rain_lsc(i)
884     snow_fall(i) = snow_con(i) + snow_lsc(i)
885     ENDDO
886    
887 guez 90 ! Humidit\'e relative pour diagnostic :
888 guez 3 DO k = 1, llm
889     DO i = 1, klon
890     zx_t = t_seri(i, k)
891 guez 207 zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t) / play(i, k)
892     zx_qs = MIN(0.5, zx_qs)
893     zcor = 1. / (1. - retv * zx_qs)
894     zx_qs = zx_qs * zcor
895 guez 202 zx_rh(i, k) = q_seri(i, k) / zx_qs
896 guez 51 zqsat(i, k) = zx_qs
897 guez 3 ENDDO
898     ENDDO
899 guez 52
900     ! Introduce the aerosol direct and first indirect radiative forcings:
901 guez 192 tau_ae = 0.
902     piz_ae = 0.
903     cg_ae = 0.
904 guez 3
905 guez 97 ! Param\`etres optiques des nuages et quelques param\`etres pour
906     ! diagnostics :
907 guez 3 if (ok_newmicro) then
908 guez 69 CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
909     cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
910     sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
911 guez 3 else
912 guez 52 CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
913     cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
914     bl95_b1, cldtaupi, re, fl)
915 guez 3 endif
916    
917 guez 154 IF (MOD(itap - 1, radpas) == 0) THEN
918 guez 212 ! Prescrire l'ozone :
919     wo = ozonecm(REAL(julien), paprs)
920    
921 guez 118 ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
922 guez 155 ! Calcul de l'abedo moyen par maille
923     albsol = sum(falbe * pctsrf, dim = 2)
924    
925 guez 62 ! Rayonnement (compatible Arpege-IFS) :
926 guez 205 CALL radlwsw(dist, mu0, fract, paprs, play, ztsol, albsol, t_seri, &
927 guez 155 q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
928     radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
929     toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
930     swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &
931     solswad, cldtaupi, topswai, solswai)
932 guez 3 ENDIF
933 guez 118
934 guez 3 ! Ajouter la tendance des rayonnements (tous les pas)
935     DO k = 1, llm
936     DO i = 1, klon
937 guez 202 t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys &
938     / 86400.
939 guez 3 ENDDO
940     ENDDO
941    
942     ! Calculer l'hydrologie de la surface
943 guez 208 zxqsurf = sum(fqsurf * pctsrf, dim = 2)
944 guez 3
945 guez 90 ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
946 guez 3 DO i = 1, klon
947     bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
948     ENDDO
949    
950 guez 90 ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
951 guez 3
952     IF (ok_orodr) THEN
953 guez 174 ! S\'election des points pour lesquels le sch\'ema est actif :
954 guez 51 igwd = 0
955     DO i = 1, klon
956     itest(i) = 0
957 guez 174 IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
958 guez 51 itest(i) = 1
959     igwd = igwd + 1
960 guez 3 ENDIF
961     ENDDO
962    
963 guez 51 CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
964 guez 150 zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
965     zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
966 guez 3
967 guez 47 ! ajout des tendances
968 guez 3 DO k = 1, llm
969     DO i = 1, klon
970     t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
971     u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
972     v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
973     ENDDO
974     ENDDO
975 guez 13 ENDIF
976 guez 3
977     IF (ok_orolf) THEN
978 guez 90 ! S\'election des points pour lesquels le sch\'ema est actif :
979 guez 51 igwd = 0
980     DO i = 1, klon
981     itest(i) = 0
982 guez 174 IF (zpic(i) - zmea(i) > 100.) THEN
983 guez 51 itest(i) = 1
984     igwd = igwd + 1
985 guez 3 ENDIF
986     ENDDO
987    
988 guez 47 CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
989     itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &
990 guez 3 d_t_lif, d_u_lif, d_v_lif)
991    
992 guez 51 ! Ajout des tendances :
993 guez 3 DO k = 1, llm
994     DO i = 1, klon
995     t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
996     u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
997     v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
998     ENDDO
999     ENDDO
1000 guez 49 ENDIF
1001 guez 3
1002 guez 90 ! Stress n\'ecessaires : toute la physique
1003 guez 3
1004     DO i = 1, klon
1005 guez 51 zustrph(i) = 0.
1006     zvstrph(i) = 0.
1007 guez 3 ENDDO
1008     DO k = 1, llm
1009     DO i = 1, klon
1010 guez 62 zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1011     * zmasse(i, k)
1012     zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1013     * zmasse(i, k)
1014 guez 3 ENDDO
1015     ENDDO
1016    
1017 guez 171 CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1018     zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1019 guez 3
1020 guez 47 ! Calcul des tendances traceurs
1021 guez 202 call phytrac(julien, time, firstcal, lafin, dtphys, t, paprs, play, mfu, &
1022     mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, yv1, ftsol, &
1023     pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, tr_seri, &
1024     zmasse, ncid_startphy)
1025 guez 3
1026 guez 190 IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &
1027     pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &
1028 guez 191 frac_impa, frac_nucl, pphis, airephy, dtphys)
1029 guez 3
1030     ! Calculer le transport de l'eau et de l'energie (diagnostique)
1031 guez 171 CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
1032 guez 3
1033 guez 31 ! diag. bilKP
1034 guez 3
1035 guez 178 CALL transp_lay(paprs, t_seri, q_seri, u_seri, v_seri, zphi, &
1036 guez 3 ve_lay, vq_lay, ue_lay, uq_lay)
1037    
1038     ! Accumuler les variables a stocker dans les fichiers histoire:
1039    
1040 guez 200 ! conversion Ec en énergie thermique
1041 guez 3 DO k = 1, llm
1042     DO i = 1, klon
1043 guez 213 d_t_ec(i, k) = 0.5 / (RCPD * (1. + RVTMP2 * q_seri(i, k))) &
1044 guez 51 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
1045     t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
1046     d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1047 guez 3 END DO
1048     END DO
1049 guez 51
1050 guez 47 ! SORTIES
1051 guez 3
1052 guez 69 ! prw = eau precipitable
1053 guez 3 DO i = 1, klon
1054     prw(i) = 0.
1055     DO k = 1, llm
1056 guez 202 prw(i) = prw(i) + q_seri(i, k) * zmasse(i, k)
1057 guez 3 ENDDO
1058     ENDDO
1059    
1060     ! Convertir les incrementations en tendances
1061    
1062     DO k = 1, llm
1063     DO i = 1, klon
1064 guez 49 d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1065     d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1066     d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1067     d_qx(i, k, ivap) = (q_seri(i, k) - qx(i, k, ivap)) / dtphys
1068     d_qx(i, k, iliq) = (ql_seri(i, k) - qx(i, k, iliq)) / dtphys
1069 guez 3 ENDDO
1070     ENDDO
1071    
1072 guez 98 DO iq = 3, nqmx
1073     DO k = 1, llm
1074     DO i = 1, klon
1075 guez 174 d_qx(i, k, iq) = (tr_seri(i, k, iq - 2) - qx(i, k, iq)) / dtphys
1076 guez 3 ENDDO
1077     ENDDO
1078 guez 98 ENDDO
1079 guez 3
1080     ! Sauvegarder les valeurs de t et q a la fin de la physique:
1081     DO k = 1, llm
1082     DO i = 1, klon
1083     t_ancien(i, k) = t_seri(i, k)
1084     q_ancien(i, k) = q_seri(i, k)
1085     ENDDO
1086     ENDDO
1087    
1088 guez 191 CALL histwrite_phy("phis", pphis)
1089     CALL histwrite_phy("aire", airephy)
1090     CALL histwrite_phy("psol", paprs(:, 1))
1091     CALL histwrite_phy("precip", rain_fall + snow_fall)
1092     CALL histwrite_phy("plul", rain_lsc + snow_lsc)
1093     CALL histwrite_phy("pluc", rain_con + snow_con)
1094 guez 205 CALL histwrite_phy("tsol", ztsol)
1095 guez 191 CALL histwrite_phy("t2m", zt2m)
1096     CALL histwrite_phy("q2m", zq2m)
1097     CALL histwrite_phy("u10m", zu10m)
1098     CALL histwrite_phy("v10m", zv10m)
1099     CALL histwrite_phy("snow", snow_fall)
1100     CALL histwrite_phy("cdrm", cdragm)
1101     CALL histwrite_phy("cdrh", cdragh)
1102     CALL histwrite_phy("topl", toplw)
1103     CALL histwrite_phy("evap", evap)
1104     CALL histwrite_phy("sols", solsw)
1105     CALL histwrite_phy("soll", sollw)
1106     CALL histwrite_phy("solldown", sollwdown)
1107     CALL histwrite_phy("bils", bils)
1108     CALL histwrite_phy("sens", - sens)
1109     CALL histwrite_phy("fder", fder)
1110     CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
1111     CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1112     CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1113     CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
1114 guez 3
1115 guez 191 DO nsrf = 1, nbsrf
1116 guez 202 CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf) * 100.)
1117 guez 191 CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1118 guez 206 CALL histwrite_phy("sens_"//clnsurf(nsrf), flux_t(:, nsrf))
1119 guez 191 CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1120     CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1121 guez 206 CALL histwrite_phy("taux_"//clnsurf(nsrf), flux_u(:, nsrf))
1122     CALL histwrite_phy("tauy_"//clnsurf(nsrf), flux_v(:, nsrf))
1123 guez 191 CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1124     CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1125     END DO
1126    
1127     CALL histwrite_phy("albs", albsol)
1128 guez 212 CALL histwrite_phy("tro3", wo * dobson_u * 1e3 / zmasse / rmo3 * md)
1129 guez 191 CALL histwrite_phy("rugs", zxrugs)
1130     CALL histwrite_phy("s_pblh", s_pblh)
1131     CALL histwrite_phy("s_pblt", s_pblt)
1132     CALL histwrite_phy("s_lcl", s_lcl)
1133     CALL histwrite_phy("s_capCL", s_capCL)
1134     CALL histwrite_phy("s_oliqCL", s_oliqCL)
1135     CALL histwrite_phy("s_cteiCL", s_cteiCL)
1136     CALL histwrite_phy("s_therm", s_therm)
1137     CALL histwrite_phy("s_trmb1", s_trmb1)
1138     CALL histwrite_phy("s_trmb2", s_trmb2)
1139     CALL histwrite_phy("s_trmb3", s_trmb3)
1140 guez 206
1141     if (conv_emanuel) then
1142     CALL histwrite_phy("ptop", ema_pct)
1143     CALL histwrite_phy("dnwd0", - mp)
1144     end if
1145    
1146 guez 191 CALL histwrite_phy("temp", t_seri)
1147     CALL histwrite_phy("vitu", u_seri)
1148     CALL histwrite_phy("vitv", v_seri)
1149     CALL histwrite_phy("geop", zphi)
1150     CALL histwrite_phy("pres", play)
1151     CALL histwrite_phy("dtvdf", d_t_vdf)
1152     CALL histwrite_phy("dqvdf", d_q_vdf)
1153     CALL histwrite_phy("rhum", zx_rh)
1154 guez 213 CALL histwrite_phy("d_t_ec", d_t_ec)
1155     CALL histwrite_phy("dtsw0", heat0 / 86400.)
1156     CALL histwrite_phy("dtlw0", - cool0 / 86400.)
1157 guez 215 CALL histwrite_phy("msnow", sum(fsnow * pctsrf, dim = 2))
1158 guez 191
1159     if (ok_instan) call histsync(nid_ins)
1160    
1161 guez 157 IF (lafin) then
1162     call NF95_CLOSE(ncid_startphy)
1163 guez 175 CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1164 guez 157 fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1165     radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1166     t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1167     w01)
1168     end IF
1169 guez 3
1170 guez 35 firstcal = .FALSE.
1171    
1172 guez 3 END SUBROUTINE physiq
1173    
1174     end module physiq_m
  ViewVC Help
Powered by ViewVC 1.1.21