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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 226 - (hide annotations)
Mon Oct 16 13:04:05 2017 UTC (6 years, 5 months ago) by guez
File size: 39032 byte(s)
In clmain, u1lay and v1lay do not depend on the subsurface. So output
values of zu1 and zv1 are simply u(:, 1) and v(:, 1). Just remove the
corresponding computations from clmain and define yu1 and yv1 in physiq.

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

  ViewVC Help
Powered by ViewVC 1.1.21