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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 209 - (hide annotations)
Wed Dec 7 17:37:21 2016 UTC (7 years, 5 months ago) by guez
File size: 40294 byte(s) 
The program did not work with cycle_diurne set to false. mu0 in
physiq, which is supposed to be a cosine, was set to -999.999. So prmu
in swu had a value of the order of 1e3. So zrmum1 in sw2s had a value
of the order of 1e3. So zrayl in sw2s had a value of the order of
1e15. So ztray and ptauaz in swclr also had a large value. So zcorae
at line 138 in swclr had a large negative value, which resulted in
overflow at line 138.

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