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

Annotation of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 202 - (hide annotations)
Wed Jun 8 12:23:41 2016 UTC (7 years, 11 months ago) by guez
Original Path: trunk/Sources/phylmd/physiq.f
File size: 48454 byte(s) 
Promoted lmt_pas from local variable of physiq to variable of module
conf_gcm_m.

Removed variable run_off of module interface_surf. Was not
used. Called run_off_ter in LMDZ, but not used nor printed there
either.

Simplified logic in interfoce_lim. The way it was convoluted with
interfsurf_hq and clmain was quite a mess. Extracted reading of SST
into a separate procedure: read_sst. We do not need SST and pctsrf_new
at the same time: SST is not needed for sea-ice surface. I did not
like this programming: going through the procedure repeatedly for
different purposes and testing inside whether there was something to
do or it was already done. Reading is now only controlled by itap and
lmt_pas, instead of debut, jour, jour_lu and deja_lu. Now we do not
copy from pct_tmp to pctsrf_new every time step.

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