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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 195 - (hide annotations)
Wed May 18 17:56:44 2016 UTC (8 years ago) by guez
File size: 48899 byte(s) 
In cv30_feed, iflag1 is 0 on entry so we can simplify the test for
iflag1 = 7.

In cv30_feed, for the computation of icb, replaced sequential search
(with a useless end of loop on k) by a call to locate.

In CV30 routines, replaced len, nloc, nd, na by klon or
klev. Philosophy: no more generality than actually necessary.

Converted as many variables as possible to named constants in
cv30_param_m and downgraded pbcrit, ptcrit, dtovsh, dpbase, dttrig,
tau, delta to local objects in procedures. spfac, betad and omtrain
are useless and removed.

Instead of filling the array sigp with the constant spfac in
cv30_undilute2, just made sigp a constant in cv30_unsat.

In cv_driver, define as allocatable variables that are only
used on the range (ncum, nl).
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 178 USE clesphys, ONLY: cdhmax, cdmmax, ecrit_hf, ecrit_ins, ecrit_mth, &
22 guez 191 ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin, 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 191 USE conf_gcm_m, ONLY: offline, day_step, iphysiq
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 47 ! vitesse dans la direction X (de O a E) en m/s
93 guez 51
94 guez 98 REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m/s
95     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 98 REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa/s
101     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     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 47 ! pour phsystoke 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     REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)
130 guez 47 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     ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
146     ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
147     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     ! hauteur de neige, en kg/m2/s
210    
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     ! liquid water mass flux (kg/m2/s), positive down
224 guez 62
225 guez 101 REAL, save:: snow_fall(klon)
226     ! solid water mass flux (kg/m2/s), positive down
227    
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 7 INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day
248 guez 70 REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
249     REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE
250 guez 155 REAL, save:: albsol(klon) ! albedo du sol total visible
251 guez 17 REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
252 guez 3
253 guez 72 real, save:: clwcon(klon, llm), rnebcon(klon, llm)
254     real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
255 guez 3
256 guez 47 REAL rhcl(klon, llm) ! humiditi relative ciel clair
257     REAL dialiq(klon, llm) ! eau liquide nuageuse
258     REAL diafra(klon, llm) ! fraction nuageuse
259     REAL cldliq(klon, llm) ! eau liquide nuageuse
260     REAL cldfra(klon, llm) ! fraction nuageuse
261     REAL cldtau(klon, llm) ! epaisseur optique
262     REAL cldemi(klon, llm) ! emissivite infrarouge
263 guez 3
264 guez 47 REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite
265     REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur
266     REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u
267     REAL fluxv(klon, llm, nbsrf) ! flux turbulent de vitesse v
268 guez 3
269     REAL zxfluxt(klon, llm)
270     REAL zxfluxq(klon, llm)
271     REAL zxfluxu(klon, llm)
272     REAL zxfluxv(klon, llm)
273    
274 guez 90 ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
275     ! les variables soient r\'emanentes.
276 guez 53 REAL, save:: heat(klon, llm) ! chauffage solaire
277 guez 154 REAL, save:: heat0(klon, llm) ! chauffage solaire ciel clair
278 guez 62 REAL, save:: cool(klon, llm) ! refroidissement infrarouge
279 guez 154 REAL, save:: cool0(klon, llm) ! refroidissement infrarouge ciel clair
280 guez 72 REAL, save:: topsw(klon), toplw(klon), solsw(klon)
281 guez 90 REAL, save:: sollw(klon) ! rayonnement infrarouge montant \`a la surface
282 guez 72 real, save:: sollwdown(klon) ! downward LW flux at surface
283 guez 62 REAL, save:: topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
284 guez 154 REAL, save:: albpla(klon)
285 guez 191 REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
286     REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
287 guez 3
288     REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)
289 guez 49 REAL conv_t(klon, llm) ! convergence of temperature (K/s)
290 guez 3
291 guez 191 REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut
292     REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree
293 guez 3
294     REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
295    
296 guez 118 REAL dist, mu0(klon), fract(klon)
297     real longi
298 guez 3 REAL z_avant(klon), z_apres(klon), z_factor(klon)
299     REAL za, zb
300 guez 103 REAL zx_t, zx_qs, zcor
301 guez 3 real zqsat(klon, llm)
302     INTEGER i, k, iq, nsrf
303 guez 69 REAL, PARAMETER:: t_coup = 234.
304 guez 3 REAL zphi(klon, llm)
305    
306 guez 186 ! cf. Anne Mathieu variables pour la couche limite atmosphérique (hbtm)
307 guez 3
308 guez 49 REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
309     REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
310     REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
311     REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
312     REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
313     REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite
314     REAL, SAVE:: therm(klon, nbsrf)
315     REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape
316 guez 190 REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition
317 guez 49 REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega
318 guez 186 ! Grandeurs de sorties
319 guez 3 REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
320     REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
321     REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
322     REAL s_trmb3(klon)
323    
324 guez 175 ! Variables pour la convection de K. Emanuel :
325 guez 3
326 guez 47 REAL upwd(klon, llm) ! saturated updraft mass flux
327     REAL dnwd(klon, llm) ! saturated downdraft mass flux
328     REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux
329     REAL cape(klon) ! CAPE
330 guez 3 SAVE cape
331    
332 guez 47 INTEGER iflagctrl(klon) ! flag fonctionnement de convect
333 guez 3
334     ! Variables du changement
335    
336     ! con: convection
337 guez 51 ! lsc: large scale condensation
338 guez 3 ! ajs: ajustement sec
339 guez 90 ! eva: \'evaporation de l'eau liquide nuageuse
340 guez 51 ! vdf: vertical diffusion in boundary layer
341 guez 3 REAL d_t_con(klon, llm), d_q_con(klon, llm)
342     REAL d_u_con(klon, llm), d_v_con(klon, llm)
343     REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
344     REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
345     REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
346     REAL rneb(klon, llm)
347    
348 guez 71 REAL mfu(klon, llm), mfd(klon, llm)
349 guez 3 REAL pen_u(klon, llm), pen_d(klon, llm)
350     REAL pde_u(klon, llm), pde_d(klon, llm)
351     INTEGER kcbot(klon), kctop(klon), kdtop(klon)
352 guez 51 REAL pmflxr(klon, llm + 1), pmflxs(klon, llm + 1)
353     REAL prfl(klon, llm + 1), psfl(klon, llm + 1)
354 guez 3
355 guez 62 INTEGER, save:: ibas_con(klon), itop_con(klon)
356 guez 183 real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
357 guez 3
358     REAL rain_con(klon), rain_lsc(klon)
359 guez 183 REAL, save:: snow_con(klon) ! neige (mm / s)
360 guez 180 real snow_lsc(klon)
361 guez 3 REAL d_ts(klon, nbsrf)
362    
363     REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
364     REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
365    
366     REAL d_u_oro(klon, llm), d_v_oro(klon, llm)
367     REAL d_t_oro(klon, llm)
368     REAL d_u_lif(klon, llm), d_v_lif(klon, llm)
369     REAL d_t_lif(klon, llm)
370    
371 guez 68 REAL, save:: ratqs(klon, llm)
372     real ratqss(klon, llm), ratqsc(klon, llm)
373     real:: ratqsbas = 0.01, ratqshaut = 0.3
374 guez 3
375     ! Parametres lies au nouveau schema de nuages (SB, PDF)
376 guez 68 real:: fact_cldcon = 0.375
377     real:: facttemps = 1.e-4
378     logical:: ok_newmicro = .true.
379 guez 3 real facteur
380    
381 guez 68 integer:: iflag_cldcon = 1
382 guez 3 logical ptconv(klon, llm)
383    
384 guez 175 ! Variables pour effectuer les appels en s\'erie :
385 guez 3
386     REAL t_seri(klon, llm), q_seri(klon, llm)
387 guez 98 REAL ql_seri(klon, llm)
388 guez 3 REAL u_seri(klon, llm), v_seri(klon, llm)
389 guez 98 REAL tr_seri(klon, llm, nqmx - 2)
390 guez 3
391     REAL zx_rh(klon, llm)
392    
393     REAL zustrdr(klon), zvstrdr(klon)
394     REAL zustrli(klon), zvstrli(klon)
395     REAL zustrph(klon), zvstrph(klon)
396     REAL aam, torsfc
397    
398     REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
399     REAL vq_lay(klon, llm) ! transport meri. de l'eau a chaque niveau vert.
400     REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
401     REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
402    
403     real date0
404    
405 guez 90 ! Variables li\'ees au bilan d'\'energie et d'enthalpie :
406 guez 3 REAL ztsol(klon)
407 guez 98 REAL d_h_vcol, d_qt, d_ec
408 guez 51 REAL, SAVE:: d_h_vcol_phy
409 guez 47 REAL zero_v(klon)
410 guez 97 CHARACTER(LEN = 20) tit
411 guez 51 INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics
412 guez 190 INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation
413 guez 51
414 guez 90 REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique
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     REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)
425    
426 guez 49 REAL, save:: sulfate_pi(klon, llm)
427 guez 175 ! 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 lmt_pas = day_step / iphysiq
538     print *, 'Number of time steps of "physics" per day: ', lmt_pas
539 guez 3
540 guez 154 radpas = lmt_pas / nbapp_rad
541 guez 191 print *, "radpas = ", radpas
542 guez 154
543 guez 90 ! Initialisation pour le sch\'ema de convection d'Emanuel :
544 guez 182 IF (conv_emanuel) THEN
545 guez 69 ibas_con = 1
546     itop_con = 1
547 guez 3 ENDIF
548    
549     IF (ok_orodr) THEN
550 guez 13 rugoro = MAX(1e-5, zstd * zsig / 2)
551 guez 54 CALL SUGWD(paprs, play)
552 guez 13 else
553     rugoro = 0.
554 guez 3 ENDIF
555    
556 guez 47 ecrit_ins = NINT(ecrit_ins/dtphys)
557     ecrit_hf = NINT(ecrit_hf/dtphys)
558     ecrit_mth = NINT(ecrit_mth/dtphys)
559     ecrit_tra = NINT(86400.*ecrit_tra/dtphys)
560     ecrit_reg = NINT(ecrit_reg/dtphys)
561 guez 3
562 guez 47 ! Initialisation des sorties
563 guez 3
564 guez 191 call ini_histins(dtphys)
565 guez 129 CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)
566 guez 69 ! Positionner date0 pour initialisation de ORCHIDEE
567     print *, 'physiq date0: ', date0
568 guez 191 CALL phyredem0(lmt_pas)
569 guez 7 ENDIF test_firstcal
570 guez 3
571 guez 91 ! We will modify variables *_seri and we will not touch variables
572 guez 98 ! u, v, t, qx:
573     t_seri = t
574     u_seri = u
575     v_seri = v
576     q_seri = qx(:, :, ivap)
577     ql_seri = qx(:, :, iliq)
578 guez 157 tr_seri = qx(:, :, 3:nqmx)
579 guez 3
580 guez 98 ztsol = sum(ftsol * pctsrf, dim = 2)
581 guez 3
582 guez 190 IF (if_ebil >= 1) THEN
583 guez 62 tit = 'after dynamics'
584     CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
585 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
586 guez 90 ! Comme les tendances de la physique sont ajout\'es dans la
587 guez 190 ! dynamique, la variation d'enthalpie par la dynamique devrait
588     ! \^etre \'egale \`a la variation de la physique au pas de temps
589     ! pr\'ec\'edent. Donc la somme de ces 2 variations devrait \^etre
590     ! nulle.
591 guez 62 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
592 guez 51 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &
593 guez 98 d_qt, 0.)
594 guez 3 END IF
595    
596 guez 51 ! Diagnostic de la tendance dynamique :
597 guez 3 IF (ancien_ok) THEN
598     DO k = 1, llm
599     DO i = 1, klon
600 guez 49 d_t_dyn(i, k) = (t_seri(i, k) - t_ancien(i, k)) / dtphys
601     d_q_dyn(i, k) = (q_seri(i, k) - q_ancien(i, k)) / dtphys
602 guez 3 ENDDO
603     ENDDO
604     ELSE
605     DO k = 1, llm
606     DO i = 1, klon
607 guez 72 d_t_dyn(i, k) = 0.
608     d_q_dyn(i, k) = 0.
609 guez 3 ENDDO
610     ENDDO
611     ancien_ok = .TRUE.
612     ENDIF
613    
614     ! Ajouter le geopotentiel du sol:
615     DO k = 1, llm
616     DO i = 1, klon
617     zphi(i, k) = pphi(i, k) + pphis(i)
618     ENDDO
619     ENDDO
620    
621 guez 49 ! Check temperatures:
622 guez 3 CALL hgardfou(t_seri, ftsol)
623    
624 guez 191 call increment_itap
625 guez 130 julien = MOD(dayvrai, 360)
626 guez 3 if (julien == 0) julien = 360
627    
628 guez 103 forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
629 guez 17
630 guez 98 ! Prescrire l'ozone :
631 guez 57 wo = ozonecm(REAL(julien), paprs)
632 guez 3
633 guez 90 ! \'Evaporation de l'eau liquide nuageuse :
634 guez 51 DO k = 1, llm
635 guez 3 DO i = 1, klon
636 guez 51 zb = MAX(0., ql_seri(i, k))
637     t_seri(i, k) = t_seri(i, k) &
638     - zb * RLVTT / RCPD / (1. + RVTMP2 * q_seri(i, k))
639 guez 3 q_seri(i, k) = q_seri(i, k) + zb
640     ENDDO
641     ENDDO
642 guez 51 ql_seri = 0.
643 guez 3
644 guez 190 IF (if_ebil >= 2) THEN
645 guez 62 tit = 'after reevap'
646     CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &
647 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
648 guez 62 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
649 guez 98 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
650 guez 3 END IF
651    
652 guez 98 frugs = MAX(frugs, 0.000015)
653     zxrugs = sum(frugs * pctsrf, dim = 2)
654 guez 3
655 guez 191 ! Calculs n\'ecessaires au calcul de l'albedo dans l'interface avec
656 guez 118 ! la surface.
657 guez 3
658 guez 118 CALL orbite(REAL(julien), longi, dist)
659 guez 3 IF (cycle_diurne) THEN
660 guez 125 CALL zenang(longi, time, dtphys * radpas, mu0, fract)
661 guez 3 ELSE
662 guez 174 mu0 = - 999.999
663 guez 3 ENDIF
664    
665 guez 47 ! Calcul de l'abedo moyen par maille
666 guez 98 albsol = sum(falbe * pctsrf, dim = 2)
667 guez 3
668 guez 90 ! R\'epartition sous maille des flux longwave et shortwave
669     ! R\'epartition du longwave par sous-surface lin\'earis\'ee
670 guez 3
671 guez 98 forall (nsrf = 1: nbsrf)
672     fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &
673     * (ztsol - ftsol(:, nsrf))
674     fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
675     END forall
676 guez 3
677     fder = dlw
678    
679 guez 38 ! Couche limite:
680 guez 3
681 guez 191 CALL clmain(dtphys, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, v_seri, &
682     julien, mu0, ftsol, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, &
683     ftsoil, qsol, paprs, play, fsnow, fqsurf, fevap, falbe, fluxlat, &
684     rain_fall, snow_fall, fsolsw, fsollw, fder, rlat, frugs, firstcal, &
685     agesno, rugoro, d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, &
686     fluxq, fluxu, fluxv, cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, &
687     yv1, t2m, q2m, u10m, v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, &
688     trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)
689 guez 3
690 guez 90 ! Incr\'ementation des flux
691 guez 40
692 guez 51 zxfluxt = 0.
693     zxfluxq = 0.
694     zxfluxu = 0.
695     zxfluxv = 0.
696 guez 3 DO nsrf = 1, nbsrf
697     DO k = 1, llm
698     DO i = 1, klon
699 guez 70 zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)
700     zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)
701     zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)
702     zxfluxv(i, k) = zxfluxv(i, k) + fluxv(i, k, nsrf) * pctsrf(i, nsrf)
703 guez 3 END DO
704     END DO
705     END DO
706     DO i = 1, klon
707     sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol
708 guez 90 evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol
709 guez 3 fder(i) = dlw(i) + dsens(i) + devap(i)
710     ENDDO
711    
712     DO k = 1, llm
713     DO i = 1, klon
714     t_seri(i, k) = t_seri(i, k) + d_t_vdf(i, k)
715     q_seri(i, k) = q_seri(i, k) + d_q_vdf(i, k)
716     u_seri(i, k) = u_seri(i, k) + d_u_vdf(i, k)
717     v_seri(i, k) = v_seri(i, k) + d_v_vdf(i, k)
718     ENDDO
719     ENDDO
720    
721 guez 190 IF (if_ebil >= 2) THEN
722 guez 62 tit = 'after clmain'
723     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
724 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
725 guez 62 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
726 guez 98 sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
727 guez 3 END IF
728    
729 guez 49 ! Update surface temperature:
730 guez 3
731     DO i = 1, klon
732 guez 72 zxtsol(i) = 0.
733     zxfluxlat(i) = 0.
734 guez 3
735 guez 72 zt2m(i) = 0.
736     zq2m(i) = 0.
737     zu10m(i) = 0.
738     zv10m(i) = 0.
739     zxffonte(i) = 0.
740     zxfqcalving(i) = 0.
741 guez 3
742 guez 190 s_pblh(i) = 0.
743     s_lcl(i) = 0.
744 guez 72 s_capCL(i) = 0.
745     s_oliqCL(i) = 0.
746     s_cteiCL(i) = 0.
747     s_pblT(i) = 0.
748     s_therm(i) = 0.
749     s_trmb1(i) = 0.
750     s_trmb2(i) = 0.
751     s_trmb3(i) = 0.
752 guez 191 ENDDO
753 guez 3
754 guez 191 call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
755    
756 guez 3 DO nsrf = 1, nbsrf
757     DO i = 1, klon
758     ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)
759     zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)
760     zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)
761    
762     zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)
763     zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)
764     zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)
765     zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)
766     zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)
767 guez 47 zxfqcalving(i) = zxfqcalving(i) + &
768 guez 3 fqcalving(i, nsrf)*pctsrf(i, nsrf)
769     s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)
770     s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)
771     s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)
772     s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)
773     s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)
774     s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)
775     s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)
776     s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)
777     s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)
778     s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)
779     ENDDO
780     ENDDO
781    
782 guez 97 ! Si une sous-fraction n'existe pas, elle prend la température moyenne :
783 guez 3 DO nsrf = 1, nbsrf
784     DO i = 1, klon
785 guez 47 IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)
786 guez 3
787 guez 47 IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)
788     IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)
789     IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)
790     IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)
791     IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)
792     IF (pctsrf(i, nsrf) < epsfra) &
793 guez 3 fqcalving(i, nsrf) = zxfqcalving(i)
794 guez 51 IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf) = s_pblh(i)
795     IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf) = s_lcl(i)
796     IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf) = s_capCL(i)
797     IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf) = s_oliqCL(i)
798     IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf) = s_cteiCL(i)
799     IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf) = s_pblT(i)
800     IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf) = s_therm(i)
801     IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf) = s_trmb1(i)
802     IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf) = s_trmb2(i)
803     IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf) = s_trmb3(i)
804 guez 3 ENDDO
805     ENDDO
806    
807 guez 98 ! Calculer la dérive du flux infrarouge
808 guez 3
809     DO i = 1, klon
810 guez 190 dlw(i) = - 4. * RSIGMA * zxtsol(i)**3
811 guez 3 ENDDO
812    
813 guez 103 IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)
814    
815 guez 190 ! Appeler la convection
816 guez 3
817 guez 182 if (conv_emanuel) then
818 guez 99 da = 0.
819     mp = 0.
820     phi = 0.
821 guez 195 CALL concvl(paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, w01, &
822     d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, itop_con, &
823     upwd, dnwd, dnwd0, Ma, cape, iflagctrl, qcondc, pmflxr, da, phi, mp)
824 guez 183 snow_con = 0.
825 guez 62 clwcon0 = qcondc
826 guez 71 mfu = upwd + dnwd
827 guez 3
828 guez 103 IF (thermcep) THEN
829     zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
830     zqsat = zqsat / (1. - retv * zqsat)
831     ELSE
832     zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
833     ENDIF
834 guez 3
835 guez 103 ! Properties of convective clouds
836 guez 71 clwcon0 = fact_cldcon * clwcon0
837 guez 62 call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
838     rnebcon0)
839 guez 72
840 guez 190 forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
841 guez 72 mfd = 0.
842     pen_u = 0.
843     pen_d = 0.
844     pde_d = 0.
845     pde_u = 0.
846 guez 182 else
847     conv_q = d_q_dyn + d_q_vdf / dtphys
848     conv_t = d_t_dyn + d_t_vdf / dtphys
849     z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
850     CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
851     q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
852     d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
853     mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
854     kdtop, pmflxr, pmflxs)
855     WHERE (rain_con < 0.) rain_con = 0.
856     WHERE (snow_con < 0.) snow_con = 0.
857     ibas_con = llm + 1 - kcbot
858     itop_con = llm + 1 - kctop
859 guez 69 END if
860 guez 3
861     DO k = 1, llm
862     DO i = 1, klon
863     t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
864     q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
865     u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
866     v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
867     ENDDO
868     ENDDO
869    
870 guez 190 IF (if_ebil >= 2) THEN
871 guez 62 tit = 'after convect'
872     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
873 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
874 guez 62 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
875 guez 98 zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec)
876 guez 3 END IF
877    
878     IF (check) THEN
879 guez 98 za = qcheck(paprs, q_seri, ql_seri)
880 guez 62 print *, "aprescon = ", za
881 guez 72 zx_t = 0.
882     za = 0.
883 guez 3 DO i = 1, klon
884     za = za + airephy(i)/REAL(klon)
885     zx_t = zx_t + (rain_con(i)+ &
886     snow_con(i))*airephy(i)/REAL(klon)
887     ENDDO
888 guez 47 zx_t = zx_t/za*dtphys
889 guez 62 print *, "Precip = ", zx_t
890 guez 3 ENDIF
891 guez 69
892 guez 182 IF (.not. conv_emanuel) THEN
893 guez 69 z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
894     z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
895 guez 3 DO k = 1, llm
896     DO i = 1, klon
897 guez 52 IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
898 guez 3 q_seri(i, k) = q_seri(i, k) * z_factor(i)
899     ENDIF
900     ENDDO
901     ENDDO
902     ENDIF
903    
904 guez 90 ! Convection s\`eche (thermiques ou ajustement)
905 guez 3
906 guez 51 d_t_ajs = 0.
907     d_u_ajs = 0.
908     d_v_ajs = 0.
909     d_q_ajs = 0.
910     fm_therm = 0.
911     entr_therm = 0.
912 guez 3
913 guez 47 if (iflag_thermals == 0) then
914     ! Ajustement sec
915     CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)
916 guez 13 t_seri = t_seri + d_t_ajs
917     q_seri = q_seri + d_q_ajs
918 guez 3 else
919 guez 47 ! Thermiques
920     call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
921     q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
922 guez 3 endif
923    
924 guez 190 IF (if_ebil >= 2) THEN
925 guez 62 tit = 'after dry_adjust'
926     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
927 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
928 guez 3 END IF
929    
930 guez 47 ! Caclul des ratqs
931 guez 3
932 guez 90 ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
933     ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
934 guez 3 if (iflag_cldcon == 1) then
935 guez 51 do k = 1, llm
936     do i = 1, klon
937 guez 3 if(ptconv(i, k)) then
938 guez 70 ratqsc(i, k) = ratqsbas + fact_cldcon &
939     * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
940 guez 3 else
941 guez 51 ratqsc(i, k) = 0.
942 guez 3 endif
943     enddo
944     enddo
945     endif
946    
947 guez 47 ! ratqs stables
948 guez 51 do k = 1, llm
949     do i = 1, klon
950 guez 70 ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
951 guez 190 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
952 guez 3 enddo
953     enddo
954    
955 guez 47 ! ratqs final
956 guez 69 if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
957 guez 47 ! les ratqs sont une conbinaison de ratqss et ratqsc
958     ! ratqs final
959     ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
960     ! relaxation des ratqs
961 guez 70 ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
962 guez 51 ratqs = max(ratqs, ratqsc)
963 guez 3 else
964 guez 47 ! on ne prend que le ratqs stable pour fisrtilp
965 guez 51 ratqs = ratqss
966 guez 3 endif
967    
968 guez 51 CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
969     d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
970     pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
971     psfl, rhcl)
972 guez 3
973     WHERE (rain_lsc < 0) rain_lsc = 0.
974     WHERE (snow_lsc < 0) snow_lsc = 0.
975     DO k = 1, llm
976     DO i = 1, klon
977     t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
978     q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
979     ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
980     cldfra(i, k) = rneb(i, k)
981     IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
982     ENDDO
983     ENDDO
984     IF (check) THEN
985 guez 98 za = qcheck(paprs, q_seri, ql_seri)
986 guez 62 print *, "apresilp = ", za
987 guez 72 zx_t = 0.
988     za = 0.
989 guez 3 DO i = 1, klon
990     za = za + airephy(i)/REAL(klon)
991     zx_t = zx_t + (rain_lsc(i) &
992     + snow_lsc(i))*airephy(i)/REAL(klon)
993     ENDDO
994 guez 47 zx_t = zx_t/za*dtphys
995 guez 62 print *, "Precip = ", zx_t
996 guez 3 ENDIF
997    
998 guez 190 IF (if_ebil >= 2) THEN
999 guez 62 tit = 'after fisrt'
1000     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1001 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1002 guez 62 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1003 guez 98 zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec)
1004 guez 3 END IF
1005    
1006 guez 47 ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1007 guez 3
1008     ! 1. NUAGES CONVECTIFS
1009    
1010 guez 174 IF (iflag_cldcon <= - 1) THEN
1011 guez 62 ! seulement pour Tiedtke
1012 guez 51 snow_tiedtke = 0.
1013 guez 174 if (iflag_cldcon == - 1) then
1014 guez 51 rain_tiedtke = rain_con
1015 guez 3 else
1016 guez 51 rain_tiedtke = 0.
1017     do k = 1, llm
1018     do i = 1, klon
1019 guez 7 if (d_q_con(i, k) < 0.) then
1020 guez 174 rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k)/dtphys &
1021 guez 17 *zmasse(i, k)
1022 guez 3 endif
1023     enddo
1024     enddo
1025     endif
1026    
1027     ! Nuages diagnostiques pour Tiedtke
1028 guez 69 CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1029     itop_con, diafra, dialiq)
1030 guez 3 DO k = 1, llm
1031     DO i = 1, klon
1032 guez 51 IF (diafra(i, k) > cldfra(i, k)) THEN
1033 guez 3 cldliq(i, k) = dialiq(i, k)
1034     cldfra(i, k) = diafra(i, k)
1035     ENDIF
1036     ENDDO
1037     ENDDO
1038     ELSE IF (iflag_cldcon == 3) THEN
1039 guez 72 ! On prend pour les nuages convectifs le maximum du calcul de
1040 guez 90 ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
1041 guez 72 ! d'un facteur facttemps.
1042     facteur = dtphys * facttemps
1043 guez 51 do k = 1, llm
1044     do i = 1, klon
1045 guez 70 rnebcon(i, k) = rnebcon(i, k) * facteur
1046 guez 72 if (rnebcon0(i, k) * clwcon0(i, k) &
1047     > rnebcon(i, k) * clwcon(i, k)) then
1048 guez 51 rnebcon(i, k) = rnebcon0(i, k)
1049     clwcon(i, k) = clwcon0(i, k)
1050 guez 3 endif
1051     enddo
1052     enddo
1053    
1054 guez 47 ! On prend la somme des fractions nuageuses et des contenus en eau
1055 guez 51 cldfra = min(max(cldfra, rnebcon), 1.)
1056     cldliq = cldliq + rnebcon*clwcon
1057 guez 3 ENDIF
1058    
1059 guez 51 ! 2. Nuages stratiformes
1060 guez 3
1061     IF (ok_stratus) THEN
1062 guez 47 CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
1063 guez 3 DO k = 1, llm
1064     DO i = 1, klon
1065 guez 51 IF (diafra(i, k) > cldfra(i, k)) THEN
1066 guez 3 cldliq(i, k) = dialiq(i, k)
1067     cldfra(i, k) = diafra(i, k)
1068     ENDIF
1069     ENDDO
1070     ENDDO
1071     ENDIF
1072    
1073     ! Precipitation totale
1074     DO i = 1, klon
1075     rain_fall(i) = rain_con(i) + rain_lsc(i)
1076     snow_fall(i) = snow_con(i) + snow_lsc(i)
1077     ENDDO
1078    
1079 guez 62 IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1080 guez 98 dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1081     d_qt, d_ec)
1082 guez 3
1083 guez 90 ! Humidit\'e relative pour diagnostic :
1084 guez 3 DO k = 1, llm
1085     DO i = 1, klon
1086     zx_t = t_seri(i, k)
1087     IF (thermcep) THEN
1088 guez 103 zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)
1089 guez 47 zx_qs = MIN(0.5, zx_qs)
1090 guez 174 zcor = 1./(1. - retv*zx_qs)
1091 guez 47 zx_qs = zx_qs*zcor
1092 guez 3 ELSE
1093 guez 7 IF (zx_t < t_coup) THEN
1094 guez 47 zx_qs = qsats(zx_t)/play(i, k)
1095 guez 3 ELSE
1096 guez 47 zx_qs = qsatl(zx_t)/play(i, k)
1097 guez 3 ENDIF
1098     ENDIF
1099     zx_rh(i, k) = q_seri(i, k)/zx_qs
1100 guez 51 zqsat(i, k) = zx_qs
1101 guez 3 ENDDO
1102     ENDDO
1103 guez 52
1104     ! Introduce the aerosol direct and first indirect radiative forcings:
1105 guez 192 tau_ae = 0.
1106     piz_ae = 0.
1107     cg_ae = 0.
1108 guez 3
1109 guez 97 ! Param\`etres optiques des nuages et quelques param\`etres pour
1110     ! diagnostics :
1111 guez 3 if (ok_newmicro) then
1112 guez 69 CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1113     cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1114     sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
1115 guez 3 else
1116 guez 52 CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1117     cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
1118     bl95_b1, cldtaupi, re, fl)
1119 guez 3 endif
1120    
1121 guez 154 IF (MOD(itap - 1, radpas) == 0) THEN
1122 guez 118 ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1123 guez 155 ! Calcul de l'abedo moyen par maille
1124     albsol = sum(falbe * pctsrf, dim = 2)
1125    
1126 guez 62 ! Rayonnement (compatible Arpege-IFS) :
1127 guez 155 CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &
1128     q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
1129     radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
1130     toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
1131     swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &
1132     solswad, cldtaupi, topswai, solswai)
1133 guez 3 ENDIF
1134 guez 118
1135 guez 3 ! Ajouter la tendance des rayonnements (tous les pas)
1136    
1137     DO k = 1, llm
1138     DO i = 1, klon
1139 guez 174 t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys/86400.
1140 guez 3 ENDDO
1141     ENDDO
1142    
1143 guez 190 IF (if_ebil >= 2) THEN
1144 guez 62 tit = 'after rad'
1145     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1146 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1147 guez 62 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1148 guez 98 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
1149 guez 3 END IF
1150    
1151     ! Calculer l'hydrologie de la surface
1152     DO i = 1, klon
1153 guez 72 zxqsurf(i) = 0.
1154     zxsnow(i) = 0.
1155 guez 3 ENDDO
1156     DO nsrf = 1, nbsrf
1157     DO i = 1, klon
1158     zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)
1159     zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)
1160     ENDDO
1161     ENDDO
1162    
1163 guez 90 ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1164 guez 3
1165     DO i = 1, klon
1166     bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1167     ENDDO
1168    
1169 guez 90 ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1170 guez 3
1171     IF (ok_orodr) THEN
1172 guez 174 ! S\'election des points pour lesquels le sch\'ema est actif :
1173 guez 51 igwd = 0
1174     DO i = 1, klon
1175     itest(i) = 0
1176 guez 174 IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
1177 guez 51 itest(i) = 1
1178     igwd = igwd + 1
1179 guez 3 ENDIF
1180     ENDDO
1181    
1182 guez 51 CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1183 guez 150 zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1184     zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1185 guez 3
1186 guez 47 ! ajout des tendances
1187 guez 3 DO k = 1, llm
1188     DO i = 1, klon
1189     t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
1190     u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
1191     v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
1192     ENDDO
1193     ENDDO
1194 guez 13 ENDIF
1195 guez 3
1196     IF (ok_orolf) THEN
1197 guez 90 ! S\'election des points pour lesquels le sch\'ema est actif :
1198 guez 51 igwd = 0
1199     DO i = 1, klon
1200     itest(i) = 0
1201 guez 174 IF (zpic(i) - zmea(i) > 100.) THEN
1202 guez 51 itest(i) = 1
1203     igwd = igwd + 1
1204 guez 3 ENDIF
1205     ENDDO
1206    
1207 guez 47 CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
1208     itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &
1209 guez 3 d_t_lif, d_u_lif, d_v_lif)
1210    
1211 guez 51 ! Ajout des tendances :
1212 guez 3 DO k = 1, llm
1213     DO i = 1, klon
1214     t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
1215     u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
1216     v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
1217     ENDDO
1218     ENDDO
1219 guez 49 ENDIF
1220 guez 3
1221 guez 90 ! Stress n\'ecessaires : toute la physique
1222 guez 3
1223     DO i = 1, klon
1224 guez 51 zustrph(i) = 0.
1225     zvstrph(i) = 0.
1226 guez 3 ENDDO
1227     DO k = 1, llm
1228     DO i = 1, klon
1229 guez 62 zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1230     * zmasse(i, k)
1231     zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1232     * zmasse(i, k)
1233 guez 3 ENDDO
1234     ENDDO
1235    
1236 guez 171 CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1237     zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1238 guez 3
1239 guez 62 IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1240 guez 98 2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1241     d_qt, d_ec)
1242 guez 3
1243 guez 47 ! Calcul des tendances traceurs
1244 guez 191 call phytrac(lmt_pas, julien, time, firstcal, lafin, dtphys, t, paprs, &
1245     play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, &
1246     yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
1247     tr_seri, zmasse, ncid_startphy)
1248 guez 3
1249 guez 190 IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &
1250     pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &
1251 guez 191 frac_impa, frac_nucl, pphis, airephy, dtphys)
1252 guez 3
1253     ! Calculer le transport de l'eau et de l'energie (diagnostique)
1254 guez 171 CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
1255 guez 3
1256 guez 31 ! diag. bilKP
1257 guez 3
1258 guez 178 CALL transp_lay(paprs, t_seri, q_seri, u_seri, v_seri, zphi, &
1259 guez 3 ve_lay, vq_lay, ue_lay, uq_lay)
1260    
1261     ! Accumuler les variables a stocker dans les fichiers histoire:
1262    
1263 guez 51 ! conversion Ec -> E thermique
1264 guez 3 DO k = 1, llm
1265     DO i = 1, klon
1266 guez 51 ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))
1267     d_t_ec(i, k) = 0.5 / ZRCPD &
1268     * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
1269     t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
1270     d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1271 guez 3 END DO
1272     END DO
1273 guez 51
1274 guez 190 IF (if_ebil >= 1) THEN
1275 guez 62 tit = 'after physic'
1276     CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1277 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1278 guez 190 ! Comme les tendances de la physique sont ajoute dans la dynamique,
1279 guez 47 ! on devrait avoir que la variation d'entalpie par la dynamique
1280     ! est egale a la variation de la physique au pas de temps precedent.
1281     ! Donc la somme de ces 2 variations devrait etre nulle.
1282 guez 62 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1283 guez 98 evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec)
1284 guez 51 d_h_vcol_phy = d_h_vcol
1285 guez 3 END IF
1286    
1287 guez 47 ! SORTIES
1288 guez 3
1289 guez 69 ! prw = eau precipitable
1290 guez 3 DO i = 1, klon
1291     prw(i) = 0.
1292     DO k = 1, llm
1293 guez 17 prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)
1294 guez 3 ENDDO
1295     ENDDO
1296    
1297     ! Convertir les incrementations en tendances
1298    
1299     DO k = 1, llm
1300     DO i = 1, klon
1301 guez 49 d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1302     d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1303     d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1304     d_qx(i, k, ivap) = (q_seri(i, k) - qx(i, k, ivap)) / dtphys
1305     d_qx(i, k, iliq) = (ql_seri(i, k) - qx(i, k, iliq)) / dtphys
1306 guez 3 ENDDO
1307     ENDDO
1308    
1309 guez 98 DO iq = 3, nqmx
1310     DO k = 1, llm
1311     DO i = 1, klon
1312 guez 174 d_qx(i, k, iq) = (tr_seri(i, k, iq - 2) - qx(i, k, iq)) / dtphys
1313 guez 3 ENDDO
1314     ENDDO
1315 guez 98 ENDDO
1316 guez 3
1317     ! Sauvegarder les valeurs de t et q a la fin de la physique:
1318     DO k = 1, llm
1319     DO i = 1, klon
1320     t_ancien(i, k) = t_seri(i, k)
1321     q_ancien(i, k) = q_seri(i, k)
1322     ENDDO
1323     ENDDO
1324    
1325 guez 191 CALL histwrite_phy("phis", pphis)
1326     CALL histwrite_phy("aire", airephy)
1327     CALL histwrite_phy("psol", paprs(:, 1))
1328     CALL histwrite_phy("precip", rain_fall + snow_fall)
1329     CALL histwrite_phy("plul", rain_lsc + snow_lsc)
1330     CALL histwrite_phy("pluc", rain_con + snow_con)
1331     CALL histwrite_phy("tsol", zxtsol)
1332     CALL histwrite_phy("t2m", zt2m)
1333     CALL histwrite_phy("q2m", zq2m)
1334     CALL histwrite_phy("u10m", zu10m)
1335     CALL histwrite_phy("v10m", zv10m)
1336     CALL histwrite_phy("snow", snow_fall)
1337     CALL histwrite_phy("cdrm", cdragm)
1338     CALL histwrite_phy("cdrh", cdragh)
1339     CALL histwrite_phy("topl", toplw)
1340     CALL histwrite_phy("evap", evap)
1341     CALL histwrite_phy("sols", solsw)
1342     CALL histwrite_phy("soll", sollw)
1343     CALL histwrite_phy("solldown", sollwdown)
1344     CALL histwrite_phy("bils", bils)
1345     CALL histwrite_phy("sens", - sens)
1346     CALL histwrite_phy("fder", fder)
1347     CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
1348     CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1349     CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1350     CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
1351 guez 3
1352 guez 191 DO nsrf = 1, nbsrf
1353     CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf)*100.)
1354     CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1355     CALL histwrite_phy("sens_"//clnsurf(nsrf), fluxt(:, 1, nsrf))
1356     CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1357     CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1358     CALL histwrite_phy("taux_"//clnsurf(nsrf), fluxu(:, 1, nsrf))
1359     CALL histwrite_phy("tauy_"//clnsurf(nsrf), fluxv(:, 1, nsrf))
1360     CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1361     CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1362     END DO
1363    
1364     CALL histwrite_phy("albs", albsol)
1365     CALL histwrite_phy("rugs", zxrugs)
1366     CALL histwrite_phy("s_pblh", s_pblh)
1367     CALL histwrite_phy("s_pblt", s_pblt)
1368     CALL histwrite_phy("s_lcl", s_lcl)
1369     CALL histwrite_phy("s_capCL", s_capCL)
1370     CALL histwrite_phy("s_oliqCL", s_oliqCL)
1371     CALL histwrite_phy("s_cteiCL", s_cteiCL)
1372     CALL histwrite_phy("s_therm", s_therm)
1373     CALL histwrite_phy("s_trmb1", s_trmb1)
1374     CALL histwrite_phy("s_trmb2", s_trmb2)
1375     CALL histwrite_phy("s_trmb3", s_trmb3)
1376     if (conv_emanuel) CALL histwrite_phy("ptop", ema_pct)
1377     CALL histwrite_phy("temp", t_seri)
1378     CALL histwrite_phy("vitu", u_seri)
1379     CALL histwrite_phy("vitv", v_seri)
1380     CALL histwrite_phy("geop", zphi)
1381     CALL histwrite_phy("pres", play)
1382     CALL histwrite_phy("dtvdf", d_t_vdf)
1383     CALL histwrite_phy("dqvdf", d_q_vdf)
1384     CALL histwrite_phy("rhum", zx_rh)
1385    
1386     if (ok_instan) call histsync(nid_ins)
1387    
1388 guez 157 IF (lafin) then
1389     call NF95_CLOSE(ncid_startphy)
1390 guez 175 CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1391 guez 157 fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1392     radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1393     t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1394     w01)
1395     end IF
1396 guez 3
1397 guez 35 firstcal = .FALSE.
1398    
1399 guez 3 END SUBROUTINE physiq
1400    
1401     end module physiq_m
  ViewVC Help
Powered by ViewVC 1.1.21