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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 192 - (hide annotations)
Thu May 12 13:00:07 2016 UTC (8 years ago) by guez
File size: 48921 byte(s) 
Removed the possibility to read aerosol fields. This was not
operational. It required fields already regridded in the three
dimensions. It seems quite weird to me not to have online vertical
regridding, since the surface pressure varies. There was the
possibility of adding vertical regridding. But development is not in
the spirit of LMDZE. Furthermore, the treatment of aerosols that was
in LMDZE is completely obsolete in LMDZ. We could try importing the
up-to-date treatment of aerosols of LMDZ, but that carries LMDZE quite
far: there is the problem of the calendar and the problem of updated
radiative transfer required for updated aerosols.
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 97 CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &
822 guez 183 w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, &
823 guez 189 itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, qcondc, pmflxr, &
824     da, phi, mp)
825 guez 183 snow_con = 0.
826 guez 62 clwcon0 = qcondc
827 guez 71 mfu = upwd + dnwd
828 guez 3
829 guez 103 IF (thermcep) THEN
830     zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
831     zqsat = zqsat / (1. - retv * zqsat)
832     ELSE
833     zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play
834     ENDIF
835 guez 3
836 guez 103 ! Properties of convective clouds
837 guez 71 clwcon0 = fact_cldcon * clwcon0
838 guez 62 call clouds_gno(klon, llm, q_seri, zqsat, clwcon0, ptconv, ratqsc, &
839     rnebcon0)
840 guez 72
841 guez 190 forall (i = 1:klon) ema_pct(i) = paprs(i, itop_con(i) + 1)
842 guez 72 mfd = 0.
843     pen_u = 0.
844     pen_d = 0.
845     pde_d = 0.
846     pde_u = 0.
847 guez 182 else
848     conv_q = d_q_dyn + d_q_vdf / dtphys
849     conv_t = d_t_dyn + d_t_vdf / dtphys
850     z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
851     CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &
852     q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &
853     d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &
854     mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &
855     kdtop, pmflxr, pmflxs)
856     WHERE (rain_con < 0.) rain_con = 0.
857     WHERE (snow_con < 0.) snow_con = 0.
858     ibas_con = llm + 1 - kcbot
859     itop_con = llm + 1 - kctop
860 guez 69 END if
861 guez 3
862     DO k = 1, llm
863     DO i = 1, klon
864     t_seri(i, k) = t_seri(i, k) + d_t_con(i, k)
865     q_seri(i, k) = q_seri(i, k) + d_q_con(i, k)
866     u_seri(i, k) = u_seri(i, k) + d_u_con(i, k)
867     v_seri(i, k) = v_seri(i, k) + d_v_con(i, k)
868     ENDDO
869     ENDDO
870    
871 guez 190 IF (if_ebil >= 2) THEN
872 guez 62 tit = 'after convect'
873     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
874 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
875 guez 62 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
876 guez 98 zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec)
877 guez 3 END IF
878    
879     IF (check) THEN
880 guez 98 za = qcheck(paprs, q_seri, ql_seri)
881 guez 62 print *, "aprescon = ", za
882 guez 72 zx_t = 0.
883     za = 0.
884 guez 3 DO i = 1, klon
885     za = za + airephy(i)/REAL(klon)
886     zx_t = zx_t + (rain_con(i)+ &
887     snow_con(i))*airephy(i)/REAL(klon)
888     ENDDO
889 guez 47 zx_t = zx_t/za*dtphys
890 guez 62 print *, "Precip = ", zx_t
891 guez 3 ENDIF
892 guez 69
893 guez 182 IF (.not. conv_emanuel) THEN
894 guez 69 z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
895     z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
896 guez 3 DO k = 1, llm
897     DO i = 1, klon
898 guez 52 IF (z_factor(i) > 1. + 1E-8 .OR. z_factor(i) < 1. - 1E-8) THEN
899 guez 3 q_seri(i, k) = q_seri(i, k) * z_factor(i)
900     ENDIF
901     ENDDO
902     ENDDO
903     ENDIF
904    
905 guez 90 ! Convection s\`eche (thermiques ou ajustement)
906 guez 3
907 guez 51 d_t_ajs = 0.
908     d_u_ajs = 0.
909     d_v_ajs = 0.
910     d_q_ajs = 0.
911     fm_therm = 0.
912     entr_therm = 0.
913 guez 3
914 guez 47 if (iflag_thermals == 0) then
915     ! Ajustement sec
916     CALL ajsec(paprs, play, t_seri, q_seri, d_t_ajs, d_q_ajs)
917 guez 13 t_seri = t_seri + d_t_ajs
918     q_seri = q_seri + d_q_ajs
919 guez 3 else
920 guez 47 ! Thermiques
921     call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &
922     q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
923 guez 3 endif
924    
925 guez 190 IF (if_ebil >= 2) THEN
926 guez 62 tit = 'after dry_adjust'
927     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
928 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
929 guez 3 END IF
930    
931 guez 47 ! Caclul des ratqs
932 guez 3
933 guez 90 ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
934     ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
935 guez 3 if (iflag_cldcon == 1) then
936 guez 51 do k = 1, llm
937     do i = 1, klon
938 guez 3 if(ptconv(i, k)) then
939 guez 70 ratqsc(i, k) = ratqsbas + fact_cldcon &
940     * (q_seri(i, 1) - q_seri(i, k)) / q_seri(i, k)
941 guez 3 else
942 guez 51 ratqsc(i, k) = 0.
943 guez 3 endif
944     enddo
945     enddo
946     endif
947    
948 guez 47 ! ratqs stables
949 guez 51 do k = 1, llm
950     do i = 1, klon
951 guez 70 ratqss(i, k) = ratqsbas + (ratqshaut - ratqsbas) &
952 guez 190 * min((paprs(i, 1) - play(i, k)) / (paprs(i, 1) - 3e4), 1.)
953 guez 3 enddo
954     enddo
955    
956 guez 47 ! ratqs final
957 guez 69 if (iflag_cldcon == 1 .or. iflag_cldcon == 2) then
958 guez 47 ! les ratqs sont une conbinaison de ratqss et ratqsc
959     ! ratqs final
960     ! 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
961     ! relaxation des ratqs
962 guez 70 ratqs = max(ratqs * exp(- dtphys * facttemps), ratqss)
963 guez 51 ratqs = max(ratqs, ratqsc)
964 guez 3 else
965 guez 47 ! on ne prend que le ratqs stable pour fisrtilp
966 guez 51 ratqs = ratqss
967 guez 3 endif
968    
969 guez 51 CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &
970     d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &
971     pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &
972     psfl, rhcl)
973 guez 3
974     WHERE (rain_lsc < 0) rain_lsc = 0.
975     WHERE (snow_lsc < 0) snow_lsc = 0.
976     DO k = 1, llm
977     DO i = 1, klon
978     t_seri(i, k) = t_seri(i, k) + d_t_lsc(i, k)
979     q_seri(i, k) = q_seri(i, k) + d_q_lsc(i, k)
980     ql_seri(i, k) = ql_seri(i, k) + d_ql_lsc(i, k)
981     cldfra(i, k) = rneb(i, k)
982     IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
983     ENDDO
984     ENDDO
985     IF (check) THEN
986 guez 98 za = qcheck(paprs, q_seri, ql_seri)
987 guez 62 print *, "apresilp = ", za
988 guez 72 zx_t = 0.
989     za = 0.
990 guez 3 DO i = 1, klon
991     za = za + airephy(i)/REAL(klon)
992     zx_t = zx_t + (rain_lsc(i) &
993     + snow_lsc(i))*airephy(i)/REAL(klon)
994     ENDDO
995 guez 47 zx_t = zx_t/za*dtphys
996 guez 62 print *, "Precip = ", zx_t
997 guez 3 ENDIF
998    
999 guez 190 IF (if_ebil >= 2) THEN
1000 guez 62 tit = 'after fisrt'
1001     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1002 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1003 guez 62 call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &
1004 guez 98 zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec)
1005 guez 3 END IF
1006    
1007 guez 47 ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1008 guez 3
1009     ! 1. NUAGES CONVECTIFS
1010    
1011 guez 174 IF (iflag_cldcon <= - 1) THEN
1012 guez 62 ! seulement pour Tiedtke
1013 guez 51 snow_tiedtke = 0.
1014 guez 174 if (iflag_cldcon == - 1) then
1015 guez 51 rain_tiedtke = rain_con
1016 guez 3 else
1017 guez 51 rain_tiedtke = 0.
1018     do k = 1, llm
1019     do i = 1, klon
1020 guez 7 if (d_q_con(i, k) < 0.) then
1021 guez 174 rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k)/dtphys &
1022 guez 17 *zmasse(i, k)
1023 guez 3 endif
1024     enddo
1025     enddo
1026     endif
1027    
1028     ! Nuages diagnostiques pour Tiedtke
1029 guez 69 CALL diagcld1(paprs, play, rain_tiedtke, snow_tiedtke, ibas_con, &
1030     itop_con, diafra, dialiq)
1031 guez 3 DO k = 1, llm
1032     DO i = 1, klon
1033 guez 51 IF (diafra(i, k) > cldfra(i, k)) THEN
1034 guez 3 cldliq(i, k) = dialiq(i, k)
1035     cldfra(i, k) = diafra(i, k)
1036     ENDIF
1037     ENDDO
1038     ENDDO
1039     ELSE IF (iflag_cldcon == 3) THEN
1040 guez 72 ! On prend pour les nuages convectifs le maximum du calcul de
1041 guez 90 ! la convection et du calcul du pas de temps pr\'ec\'edent diminu\'e
1042 guez 72 ! d'un facteur facttemps.
1043     facteur = dtphys * facttemps
1044 guez 51 do k = 1, llm
1045     do i = 1, klon
1046 guez 70 rnebcon(i, k) = rnebcon(i, k) * facteur
1047 guez 72 if (rnebcon0(i, k) * clwcon0(i, k) &
1048     > rnebcon(i, k) * clwcon(i, k)) then
1049 guez 51 rnebcon(i, k) = rnebcon0(i, k)
1050     clwcon(i, k) = clwcon0(i, k)
1051 guez 3 endif
1052     enddo
1053     enddo
1054    
1055 guez 47 ! On prend la somme des fractions nuageuses et des contenus en eau
1056 guez 51 cldfra = min(max(cldfra, rnebcon), 1.)
1057     cldliq = cldliq + rnebcon*clwcon
1058 guez 3 ENDIF
1059    
1060 guez 51 ! 2. Nuages stratiformes
1061 guez 3
1062     IF (ok_stratus) THEN
1063 guez 47 CALL diagcld2(paprs, play, t_seri, q_seri, diafra, dialiq)
1064 guez 3 DO k = 1, llm
1065     DO i = 1, klon
1066 guez 51 IF (diafra(i, k) > cldfra(i, k)) THEN
1067 guez 3 cldliq(i, k) = dialiq(i, k)
1068     cldfra(i, k) = diafra(i, k)
1069     ENDIF
1070     ENDDO
1071     ENDDO
1072     ENDIF
1073    
1074     ! Precipitation totale
1075     DO i = 1, klon
1076     rain_fall(i) = rain_con(i) + rain_lsc(i)
1077     snow_fall(i) = snow_con(i) + snow_lsc(i)
1078     ENDDO
1079    
1080 guez 62 IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &
1081 guez 98 dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1082     d_qt, d_ec)
1083 guez 3
1084 guez 90 ! Humidit\'e relative pour diagnostic :
1085 guez 3 DO k = 1, llm
1086     DO i = 1, klon
1087     zx_t = t_seri(i, k)
1088     IF (thermcep) THEN
1089 guez 103 zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)
1090 guez 47 zx_qs = MIN(0.5, zx_qs)
1091 guez 174 zcor = 1./(1. - retv*zx_qs)
1092 guez 47 zx_qs = zx_qs*zcor
1093 guez 3 ELSE
1094 guez 7 IF (zx_t < t_coup) THEN
1095 guez 47 zx_qs = qsats(zx_t)/play(i, k)
1096 guez 3 ELSE
1097 guez 47 zx_qs = qsatl(zx_t)/play(i, k)
1098 guez 3 ENDIF
1099     ENDIF
1100     zx_rh(i, k) = q_seri(i, k)/zx_qs
1101 guez 51 zqsat(i, k) = zx_qs
1102 guez 3 ENDDO
1103     ENDDO
1104 guez 52
1105     ! Introduce the aerosol direct and first indirect radiative forcings:
1106 guez 192 tau_ae = 0.
1107     piz_ae = 0.
1108     cg_ae = 0.
1109 guez 3
1110 guez 97 ! Param\`etres optiques des nuages et quelques param\`etres pour
1111     ! diagnostics :
1112 guez 3 if (ok_newmicro) then
1113 guez 69 CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
1114     cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &
1115     sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)
1116 guez 3 else
1117 guez 52 CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
1118     cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &
1119     bl95_b1, cldtaupi, re, fl)
1120 guez 3 endif
1121    
1122 guez 154 IF (MOD(itap - 1, radpas) == 0) THEN
1123 guez 118 ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
1124 guez 155 ! Calcul de l'abedo moyen par maille
1125     albsol = sum(falbe * pctsrf, dim = 2)
1126    
1127 guez 62 ! Rayonnement (compatible Arpege-IFS) :
1128 guez 155 CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &
1129     q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
1130     radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
1131     toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
1132     swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &
1133     solswad, cldtaupi, topswai, solswai)
1134 guez 3 ENDIF
1135 guez 118
1136 guez 3 ! Ajouter la tendance des rayonnements (tous les pas)
1137    
1138     DO k = 1, llm
1139     DO i = 1, klon
1140 guez 174 t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys/86400.
1141 guez 3 ENDDO
1142     ENDDO
1143    
1144 guez 190 IF (if_ebil >= 2) THEN
1145 guez 62 tit = 'after rad'
1146     CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &
1147 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1148 guez 62 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &
1149 guez 98 zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)
1150 guez 3 END IF
1151    
1152     ! Calculer l'hydrologie de la surface
1153     DO i = 1, klon
1154 guez 72 zxqsurf(i) = 0.
1155     zxsnow(i) = 0.
1156 guez 3 ENDDO
1157     DO nsrf = 1, nbsrf
1158     DO i = 1, klon
1159     zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)
1160     zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)
1161     ENDDO
1162     ENDDO
1163    
1164 guez 90 ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
1165 guez 3
1166     DO i = 1, klon
1167     bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
1168     ENDDO
1169    
1170 guez 90 ! Param\'etrisation de l'orographie \`a l'\'echelle sous-maille :
1171 guez 3
1172     IF (ok_orodr) THEN
1173 guez 174 ! S\'election des points pour lesquels le sch\'ema est actif :
1174 guez 51 igwd = 0
1175     DO i = 1, klon
1176     itest(i) = 0
1177 guez 174 IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
1178 guez 51 itest(i) = 1
1179     igwd = igwd + 1
1180 guez 3 ENDIF
1181     ENDDO
1182    
1183 guez 51 CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &
1184 guez 150 zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &
1185     zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)
1186 guez 3
1187 guez 47 ! ajout des tendances
1188 guez 3 DO k = 1, llm
1189     DO i = 1, klon
1190     t_seri(i, k) = t_seri(i, k) + d_t_oro(i, k)
1191     u_seri(i, k) = u_seri(i, k) + d_u_oro(i, k)
1192     v_seri(i, k) = v_seri(i, k) + d_v_oro(i, k)
1193     ENDDO
1194     ENDDO
1195 guez 13 ENDIF
1196 guez 3
1197     IF (ok_orolf) THEN
1198 guez 90 ! S\'election des points pour lesquels le sch\'ema est actif :
1199 guez 51 igwd = 0
1200     DO i = 1, klon
1201     itest(i) = 0
1202 guez 174 IF (zpic(i) - zmea(i) > 100.) THEN
1203 guez 51 itest(i) = 1
1204     igwd = igwd + 1
1205 guez 3 ENDIF
1206     ENDDO
1207    
1208 guez 47 CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &
1209     itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &
1210 guez 3 d_t_lif, d_u_lif, d_v_lif)
1211    
1212 guez 51 ! Ajout des tendances :
1213 guez 3 DO k = 1, llm
1214     DO i = 1, klon
1215     t_seri(i, k) = t_seri(i, k) + d_t_lif(i, k)
1216     u_seri(i, k) = u_seri(i, k) + d_u_lif(i, k)
1217     v_seri(i, k) = v_seri(i, k) + d_v_lif(i, k)
1218     ENDDO
1219     ENDDO
1220 guez 49 ENDIF
1221 guez 3
1222 guez 90 ! Stress n\'ecessaires : toute la physique
1223 guez 3
1224     DO i = 1, klon
1225 guez 51 zustrph(i) = 0.
1226     zvstrph(i) = 0.
1227 guez 3 ENDDO
1228     DO k = 1, llm
1229     DO i = 1, klon
1230 guez 62 zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &
1231     * zmasse(i, k)
1232     zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &
1233     * zmasse(i, k)
1234 guez 3 ENDDO
1235     ENDDO
1236    
1237 guez 171 CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &
1238     zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)
1239 guez 3
1240 guez 62 IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &
1241 guez 98 2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &
1242     d_qt, d_ec)
1243 guez 3
1244 guez 47 ! Calcul des tendances traceurs
1245 guez 191 call phytrac(lmt_pas, julien, time, firstcal, lafin, dtphys, t, paprs, &
1246     play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, &
1247     yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
1248     tr_seri, zmasse, ncid_startphy)
1249 guez 3
1250 guez 190 IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &
1251     pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &
1252 guez 191 frac_impa, frac_nucl, pphis, airephy, dtphys)
1253 guez 3
1254     ! Calculer le transport de l'eau et de l'energie (diagnostique)
1255 guez 171 CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
1256 guez 3
1257 guez 31 ! diag. bilKP
1258 guez 3
1259 guez 178 CALL transp_lay(paprs, t_seri, q_seri, u_seri, v_seri, zphi, &
1260 guez 3 ve_lay, vq_lay, ue_lay, uq_lay)
1261    
1262     ! Accumuler les variables a stocker dans les fichiers histoire:
1263    
1264 guez 51 ! conversion Ec -> E thermique
1265 guez 3 DO k = 1, llm
1266     DO i = 1, klon
1267 guez 51 ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))
1268     d_t_ec(i, k) = 0.5 / ZRCPD &
1269     * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
1270     t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
1271     d_t_ec(i, k) = d_t_ec(i, k) / dtphys
1272 guez 3 END DO
1273     END DO
1274 guez 51
1275 guez 190 IF (if_ebil >= 1) THEN
1276 guez 62 tit = 'after physic'
1277     CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &
1278 guez 98 ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)
1279 guez 190 ! Comme les tendances de la physique sont ajoute dans la dynamique,
1280 guez 47 ! on devrait avoir que la variation d'entalpie par la dynamique
1281     ! est egale a la variation de la physique au pas de temps precedent.
1282     ! Donc la somme de ces 2 variations devrait etre nulle.
1283 guez 62 call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &
1284 guez 98 evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec)
1285 guez 51 d_h_vcol_phy = d_h_vcol
1286 guez 3 END IF
1287    
1288 guez 47 ! SORTIES
1289 guez 3
1290 guez 69 ! prw = eau precipitable
1291 guez 3 DO i = 1, klon
1292     prw(i) = 0.
1293     DO k = 1, llm
1294 guez 17 prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)
1295 guez 3 ENDDO
1296     ENDDO
1297    
1298     ! Convertir les incrementations en tendances
1299    
1300     DO k = 1, llm
1301     DO i = 1, klon
1302 guez 49 d_u(i, k) = (u_seri(i, k) - u(i, k)) / dtphys
1303     d_v(i, k) = (v_seri(i, k) - v(i, k)) / dtphys
1304     d_t(i, k) = (t_seri(i, k) - t(i, k)) / dtphys
1305     d_qx(i, k, ivap) = (q_seri(i, k) - qx(i, k, ivap)) / dtphys
1306     d_qx(i, k, iliq) = (ql_seri(i, k) - qx(i, k, iliq)) / dtphys
1307 guez 3 ENDDO
1308     ENDDO
1309    
1310 guez 98 DO iq = 3, nqmx
1311     DO k = 1, llm
1312     DO i = 1, klon
1313 guez 174 d_qx(i, k, iq) = (tr_seri(i, k, iq - 2) - qx(i, k, iq)) / dtphys
1314 guez 3 ENDDO
1315     ENDDO
1316 guez 98 ENDDO
1317 guez 3
1318     ! Sauvegarder les valeurs de t et q a la fin de la physique:
1319     DO k = 1, llm
1320     DO i = 1, klon
1321     t_ancien(i, k) = t_seri(i, k)
1322     q_ancien(i, k) = q_seri(i, k)
1323     ENDDO
1324     ENDDO
1325    
1326 guez 191 CALL histwrite_phy("phis", pphis)
1327     CALL histwrite_phy("aire", airephy)
1328     CALL histwrite_phy("psol", paprs(:, 1))
1329     CALL histwrite_phy("precip", rain_fall + snow_fall)
1330     CALL histwrite_phy("plul", rain_lsc + snow_lsc)
1331     CALL histwrite_phy("pluc", rain_con + snow_con)
1332     CALL histwrite_phy("tsol", zxtsol)
1333     CALL histwrite_phy("t2m", zt2m)
1334     CALL histwrite_phy("q2m", zq2m)
1335     CALL histwrite_phy("u10m", zu10m)
1336     CALL histwrite_phy("v10m", zv10m)
1337     CALL histwrite_phy("snow", snow_fall)
1338     CALL histwrite_phy("cdrm", cdragm)
1339     CALL histwrite_phy("cdrh", cdragh)
1340     CALL histwrite_phy("topl", toplw)
1341     CALL histwrite_phy("evap", evap)
1342     CALL histwrite_phy("sols", solsw)
1343     CALL histwrite_phy("soll", sollw)
1344     CALL histwrite_phy("solldown", sollwdown)
1345     CALL histwrite_phy("bils", bils)
1346     CALL histwrite_phy("sens", - sens)
1347     CALL histwrite_phy("fder", fder)
1348     CALL histwrite_phy("dtsvdfo", d_ts(:, is_oce))
1349     CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1350     CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1351     CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
1352 guez 3
1353 guez 191 DO nsrf = 1, nbsrf
1354     CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf)*100.)
1355     CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1356     CALL histwrite_phy("sens_"//clnsurf(nsrf), fluxt(:, 1, nsrf))
1357     CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1358     CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1359     CALL histwrite_phy("taux_"//clnsurf(nsrf), fluxu(:, 1, nsrf))
1360     CALL histwrite_phy("tauy_"//clnsurf(nsrf), fluxv(:, 1, nsrf))
1361     CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1362     CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1363     END DO
1364    
1365     CALL histwrite_phy("albs", albsol)
1366     CALL histwrite_phy("rugs", zxrugs)
1367     CALL histwrite_phy("s_pblh", s_pblh)
1368     CALL histwrite_phy("s_pblt", s_pblt)
1369     CALL histwrite_phy("s_lcl", s_lcl)
1370     CALL histwrite_phy("s_capCL", s_capCL)
1371     CALL histwrite_phy("s_oliqCL", s_oliqCL)
1372     CALL histwrite_phy("s_cteiCL", s_cteiCL)
1373     CALL histwrite_phy("s_therm", s_therm)
1374     CALL histwrite_phy("s_trmb1", s_trmb1)
1375     CALL histwrite_phy("s_trmb2", s_trmb2)
1376     CALL histwrite_phy("s_trmb3", s_trmb3)
1377     if (conv_emanuel) CALL histwrite_phy("ptop", ema_pct)
1378     CALL histwrite_phy("temp", t_seri)
1379     CALL histwrite_phy("vitu", u_seri)
1380     CALL histwrite_phy("vitv", v_seri)
1381     CALL histwrite_phy("geop", zphi)
1382     CALL histwrite_phy("pres", play)
1383     CALL histwrite_phy("dtvdf", d_t_vdf)
1384     CALL histwrite_phy("dqvdf", d_q_vdf)
1385     CALL histwrite_phy("rhum", zx_rh)
1386    
1387     if (ok_instan) call histsync(nid_ins)
1388    
1389 guez 157 IF (lafin) then
1390     call NF95_CLOSE(ncid_startphy)
1391 guez 175 CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &
1392 guez 157 fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &
1393     radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
1394     t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &
1395     w01)
1396     end IF
1397 guez 3
1398 guez 35 firstcal = .FALSE.
1399    
1400 guez 3 END SUBROUTINE physiq
1401    
1402     end module physiq_m
  ViewVC Help
Powered by ViewVC 1.1.21