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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 178 - (hide annotations)
Fri Mar 11 18:47:26 2016 UTC (8 years, 1 month ago) by guez
File size: 57235 byte(s) 
Moved variables date0, deltat, datasz_max, ncvar_ids, point, buff_pos,
buffer, regular from module histcom_var to modules where they are
defined.

Removed procedure ioipslmpp, useless for a sequential program.

Added argument datasz_max to histwrite_real (to avoid circular
dependency with histwrite).

Removed useless variables and computations everywhere.

Changed real litteral constants from default kind to double precision
in lwb, lwu, lwvn, sw1s, swtt, swtt1, swu.

Removed unused arguments: paer of sw, sw1s, sw2s, swclr; pcldsw of
sw1s, sw2s; pdsig, prayl of swr; co2_ppm of clmain, clqh; tsol of
transp_lay; nsrf of screenp; kcrit and kknu of gwstress; pstd of
orosetup.

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