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

Annotation of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log

Revision 175 - (hide annotations)
Fri Feb 5 16:02:34 2016 UTC (8 years, 3 months ago) by guez
Original Path: trunk/Sources/phylmd/physiq.f
File size: 57159 byte(s) 
Added argument itau_phy to ini_histins, phyetat0, phytrac and
phyredem0. Removed variable itau_phy of module temps. Avoiding side
effect in etat0 and phyetat0. The procedures ini_histins, phyetat0,
phytrac and phyredem0 are all called by physiq so there is no
cascading variable penalty.

In procedure inifilr, made the condition on colat0 weaker to allow for
rounding error.

Removed arguments flux_o, flux_g and t_slab of clmain, flux_o and
flux_g of clqh and interfsurf_hq, tslab and seaice of phyetat0 and
phyredem. NetCDF variables TSLAB and SEAICE no longer in
restartphy.nc. All these variables were related to the not-implemented
slab ocean. seaice and tslab were just set to 0 in phyetat0 and never
used nor changed. flux_o and flux_g were computed in clmain but never
used in physiq.

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