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

Annotation of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 212 - (hide annotations)
Thu Jan 12 12:31:31 2017 UTC (7 years, 4 months ago) by guez
Original Path: trunk/Sources/phylmd/physiq.f
File size: 40452 byte(s) 
Moved variables from module com_io_dyn to module inithist_m, where
they are defined.

Split grid_atob.f into grille_m.f and dist_sphe.f. Extracted ASCCI art
to documentation. In grille_m, use automatic arrays instead of maximum
size. In grille_m, instead of printing data for every problematic
point, print a single diagnostic message.

Removed variables top_height, overlap, lev_histhf, lev_histday,
lev_histmth, type_run, ok_isccp, ok_regdyn, lonmin_ins, lonmax_ins,
latmin_ins, latmax_ins of module clesphys, not used.

Removed variable itap of module histwrite_phy_m, not used. There is a
variable itap in module time_phylmdz.

Added output of tro3.

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