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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 191 - (hide annotations)
Mon May 9 19:56:28 2016 UTC (8 years ago) by guez
File size: 49430 byte(s)
Extracted the call to read_comdissnew out of conf_gcm.

Made ok_instan a variable of module clesphys, itau_phy a variable of
module phyetat0_m, nid_ins a variable of module ini_histins_m, itap a
variable of new module time_phylmdz, so that histwrite_phy can be
called from any procedure without the need to cascade those variables
into that procedure. Made itau_w a variable of module time_phylmdz so
that it is computed only once per time step of physics.

Extracted variables of module clesphys which were in namelist
conf_phys_nml into their own namelist, clesphys_nml, and created
procedure read_clesphys reading clesphys_nml, to avoid side effect.

No need for double precision in procedure getso4fromfile. Assume there
is a single variable for the whole year in the NetCDF file instead of
one variable per month.

Created generic procedure histwrite_phy and removed procedure
write_histins, following LMDZ. histwrite_phy has only two arguments,
can be called from anywhere, and should manage the logic of writing or
not writing into various history files with various operations. So the
test on ok_instan goes inside histwrite_phy.

Test for raz_date in phyetat0 instead of physiq to avoid side effect.

Created procedure increment_itap to avoid side effect.

Removed unnecessary differences between procedures readsulfate and
readsulfate_pi.

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

  ViewVC Help
Powered by ViewVC 1.1.21