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

Diff of /trunk/phylmd/physiq.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/Sources/phylmd/physiq.f revision 191 by guez, Mon May 9 19:56:28 2016 UTC trunk/phylmd/physiq.f revision 305 by guez, Tue Sep 11 11:08:38 2018 UTC
# Line 16  contains Line 16  contains
16    
17      use aaam_bud_m, only: aaam_bud      use aaam_bud_m, only: aaam_bud
18      USE abort_gcm_m, ONLY: abort_gcm      USE abort_gcm_m, ONLY: abort_gcm
     use aeropt_m, only: aeropt  
19      use ajsec_m, only: ajsec      use ajsec_m, only: ajsec
20      use calltherm_m, only: calltherm      use calltherm_m, only: calltherm
21      USE clesphys, ONLY: cdhmax, cdmmax, ecrit_hf, ecrit_ins, ecrit_mth, &      USE clesphys, ONLY: cdhmax, cdmmax, ecrit_ins, ok_instan
22           ecrit_reg, ecrit_tra, ksta, ksta_ter, ok_kzmin, ok_instan      USE clesphys2, ONLY: conv_emanuel, nbapp_rad, new_oliq, ok_orodr, ok_orolf
23      USE clesphys2, ONLY: cycle_diurne, conv_emanuel, nbapp_rad, new_oliq, &      USE conf_interface_m, ONLY: conf_interface
24           ok_orodr, ok_orolf      USE pbl_surface_m, ONLY: pbl_surface
     USE clmain_m, ONLY: clmain  
25      use clouds_gno_m, only: clouds_gno      use clouds_gno_m, only: clouds_gno
26      use comconst, only: dtphys      use comconst, only: dtphys
27      USE comgeomphy, ONLY: airephy      USE comgeomphy, ONLY: airephy
28      USE concvl_m, ONLY: concvl      USE concvl_m, ONLY: concvl
29      USE conf_gcm_m, ONLY: offline, day_step, iphysiq      USE conf_gcm_m, ONLY: lmt_pas
30      USE conf_phys_m, ONLY: conf_phys      USE conf_phys_m, ONLY: conf_phys
31      use conflx_m, only: conflx      use conflx_m, only: conflx
32      USE ctherm, ONLY: iflag_thermals, nsplit_thermals      USE ctherm, ONLY: iflag_thermals, nsplit_thermals
33      use diagcld2_m, only: diagcld2      use diagcld2_m, only: diagcld2
34      use diagetpq_m, only: diagetpq      USE dimensions, ONLY: llm, nqmx
     use diagphy_m, only: diagphy  
     USE dimens_m, ONLY: llm, nqmx  
35      USE dimphy, ONLY: klon      USE dimphy, ONLY: klon
36      USE dimsoil, ONLY: nsoilmx      USE dimsoil, ONLY: nsoilmx
37      use drag_noro_m, only: drag_noro      use drag_noro_m, only: drag_noro
38      use dynetat0_m, only: day_ref, annee_ref      use dynetat0_m, only: day_ref, annee_ref
39      USE fcttre, ONLY: foeew, qsatl, qsats, thermcep      USE fcttre, ONLY: foeew
40      use fisrtilp_m, only: fisrtilp      use fisrtilp_m, only: fisrtilp
41      USE hgardfou_m, ONLY: hgardfou      USE hgardfou_m, ONLY: hgardfou
42      USE histsync_m, ONLY: histsync      USE histsync_m, ONLY: histsync
# Line 48  contains Line 44  contains
44      USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &      USE indicesol, ONLY: clnsurf, epsfra, is_lic, is_oce, is_sic, is_ter, &
45           nbsrf           nbsrf
46      USE ini_histins_m, ONLY: ini_histins, nid_ins      USE ini_histins_m, ONLY: ini_histins, nid_ins
47        use lift_noro_m, only: lift_noro
48      use netcdf95, only: NF95_CLOSE      use netcdf95, only: NF95_CLOSE
49      use newmicro_m, only: newmicro      use newmicro_m, only: newmicro
50      use nr_util, only: assert      use nr_util, only: assert
51      use nuage_m, only: nuage      use nuage_m, only: nuage
52      USE orbite_m, ONLY: orbite      USE orbite_m, ONLY: orbite
53      USE ozonecm_m, ONLY: ozonecm      USE ozonecm_m, ONLY: ozonecm
54      USE phyetat0_m, ONLY: phyetat0, rlat, rlon      USE phyetat0_m, ONLY: phyetat0
55      USE phyredem_m, ONLY: phyredem      USE phyredem_m, ONLY: phyredem
56      USE phyredem0_m, ONLY: phyredem0      USE phyredem0_m, ONLY: phyredem0
     USE phystokenc_m, ONLY: phystokenc  
57      USE phytrac_m, ONLY: phytrac      USE phytrac_m, ONLY: phytrac
     USE qcheck_m, ONLY: qcheck  
58      use radlwsw_m, only: radlwsw      use radlwsw_m, only: radlwsw
     use readsulfate_m, only: readsulfate  
     use readsulfate_preind_m, only: readsulfate_preind  
59      use yoegwd, only: sugwd      use yoegwd, only: sugwd
60      USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt      USE suphec_m, ONLY: rcpd, retv, rg, rlvtt, romega, rsigma, rtt, rmo3, md
61      use time_phylmdz, only: itap, increment_itap      use time_phylmdz, only: itap, increment_itap
62      use transp_m, only: transp      use transp_m, only: transp
63      use transp_lay_m, only: transp_lay      use transp_lay_m, only: transp_lay
# Line 92  contains Line 85  contains
85      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol      REAL, intent(in):: pphis(:) ! (klon) géopotentiel du sol
86    
87      REAL, intent(in):: u(:, :) ! (klon, llm)      REAL, intent(in):: u(:, :) ! (klon, llm)
88      ! vitesse dans la direction X (de O a E) en m/s      ! vitesse dans la direction X (de O a E) en m / s
89    
90      REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m/s      REAL, intent(in):: v(:, :) ! (klon, llm) vitesse Y (de S a N) en m / s
91      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)      REAL, intent(in):: t(:, :) ! (klon, llm) temperature (K)
92    
93      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)      REAL, intent(in):: qx(:, :, :) ! (klon, llm, nqmx)
94      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)      ! (humidit\'e sp\'ecifique et fractions massiques des autres traceurs)
95    
96      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa/s      REAL, intent(in):: omega(:, :) ! (klon, llm) vitesse verticale en Pa / s
97      REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)      REAL, intent(out):: d_u(:, :) ! (klon, llm) tendance physique de "u" (m s-2)
98      REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)      REAL, intent(out):: d_v(:, :) ! (klon, llm) tendance physique de "v" (m s-2)
99      REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K/s)      REAL, intent(out):: d_t(:, :) ! (klon, llm) tendance physique de "t" (K / s)
100    
101      REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)      REAL, intent(out):: d_qx(:, :, :) ! (klon, llm, nqmx)
102      ! tendance physique de "qx" (s-1)      ! tendance physique de "qx" (s-1)
# Line 112  contains Line 105  contains
105    
106      LOGICAL:: firstcal = .true.      LOGICAL:: firstcal = .true.
107    
     LOGICAL, PARAMETER:: check = .FALSE.  
     ! Verifier la conservation du modele en eau  
   
108      LOGICAL, PARAMETER:: ok_stratus = .FALSE.      LOGICAL, PARAMETER:: ok_stratus = .FALSE.
109      ! Ajouter artificiellement les stratus      ! Ajouter artificiellement les stratus
110    
111      ! pour phsystoke avec thermiques      ! pour phystoke avec thermiques
112      REAL fm_therm(klon, llm + 1)      REAL fm_therm(klon, llm + 1)
113      REAL entr_therm(klon, llm)      REAL entr_therm(klon, llm)
114      real, save:: q2(klon, llm + 1, nbsrf)      real, save:: q2(klon, llm + 1, nbsrf)
# Line 129  contains Line 119  contains
119      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)      REAL, save:: t_ancien(klon, llm), q_ancien(klon, llm)
120      LOGICAL, save:: ancien_ok      LOGICAL, save:: ancien_ok
121    
122      REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K/s)      REAL d_t_dyn(klon, llm) ! tendance dynamique pour "t" (K / s)
123      REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg/kg/s)      REAL d_q_dyn(klon, llm) ! tendance dynamique pour "q" (kg / kg / s)
124    
125      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)      real da(klon, llm), phi(klon, llm, llm), mp(klon, llm)
126    
127      REAL swdn0(klon, llm + 1), swdn(klon, llm + 1)      REAL, save:: swdn0(klon, llm + 1), swdn(klon, llm + 1)
128      REAL swup0(klon, llm + 1), swup(klon, llm + 1)      REAL, save:: swup0(klon, llm + 1), swup(klon, llm + 1)
129      SAVE swdn0, swdn, swup0, swup  
130        REAL, save:: lwdn0(klon, llm + 1), lwdn(klon, llm + 1)
131      REAL lwdn0(klon, llm + 1), lwdn(klon, llm + 1)      REAL, save:: lwup0(klon, llm + 1), lwup(klon, llm + 1)
     REAL lwup0(klon, llm + 1), lwup(klon, llm + 1)  
     SAVE lwdn0, lwdn, lwup0, lwup  
132    
133      ! prw: precipitable water      ! prw: precipitable water
134      real prw(klon)      real prw(klon)
135    
136      ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)      ! flwp, fiwp = Liquid Water Path & Ice Water Path (kg / m2)
137      ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)      ! flwc, fiwc = Liquid Water Content & Ice Water Content (kg / kg)
138      REAL flwp(klon), fiwp(klon)      REAL flwp(klon), fiwp(klon)
139      REAL flwc(klon, llm), fiwc(klon, llm)      REAL flwc(klon, llm), fiwc(klon, llm)
140    
# Line 156  contains Line 144  contains
144      ! Radiative transfer computations are made every "radpas" call to      ! Radiative transfer computations are made every "radpas" call to
145      ! "physiq".      ! "physiq".
146    
147      REAL radsol(klon)      REAL, save:: radsol(klon) ! bilan radiatif au sol calcule par code radiatif
     SAVE radsol ! bilan radiatif au sol calcule par code radiatif  
   
148      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction      REAL, save:: ftsol(klon, nbsrf) ! skin temperature of surface fraction
149    
150      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)      REAL, save:: ftsoil(klon, nsoilmx, nbsrf)
151      ! soil temperature of surface fraction      ! soil temperature of surface fraction
152    
     REAL, save:: fevap(klon, nbsrf) ! evaporation  
153      REAL fluxlat(klon, nbsrf)      REAL fluxlat(klon, nbsrf)
     SAVE fluxlat  
154    
155      REAL, save:: fqsurf(klon, nbsrf)      REAL, save:: fqsurf(klon, nbsrf)
156      ! humidite de l'air au contact de la surface      ! humidite de l'air au contact de la surface
157    
158      REAL, save:: qsol(klon)      REAL, save:: qsol(klon) ! column-density of water in soil, in kg m-2
159      ! column-density of water in soil, in kg m-2      REAL, save:: fsnow(klon, nbsrf) ! \'epaisseur neigeuse
   
     REAL, save:: fsnow(klon, nbsrf) ! epaisseur neigeuse  
160      REAL, save:: falbe(klon, nbsrf) ! albedo visible par type de surface      REAL, save:: falbe(klon, nbsrf) ! albedo visible par type de surface
161    
162      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :      ! Param\`etres de l'orographie \`a l'\'echelle sous-maille (OESM) :
# Line 187  contains Line 169  contains
169      REAL, save:: zval(klon) ! Minimum de l'OESM      REAL, save:: zval(klon) ! Minimum de l'OESM
170      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM      REAL, save:: rugoro(klon) ! longueur de rugosite de l'OESM
171      REAL zulow(klon), zvlow(klon)      REAL zulow(klon), zvlow(klon)
172      INTEGER igwd, itest(klon)      INTEGER ktest(klon)
173    
174      REAL, save:: agesno(klon, nbsrf) ! age de la neige      REAL, save:: agesno(klon, nbsrf) ! age de la neige
175      REAL, save:: run_off_lic_0(klon)      REAL, save:: run_off_lic_0(klon)
176    
177      ! Variables li\'ees \`a la convection d'Emanuel :      ! Variables li\'ees \`a la convection d'Emanuel :
178      REAL, save:: Ma(klon, llm) ! undilute upward mass flux      REAL, save:: Ma(klon, llm) ! undilute upward mass flux
     REAL, save:: qcondc(klon, llm) ! in-cld water content from convect  
179      REAL, save:: sig1(klon, llm), w01(klon, llm)      REAL, save:: sig1(klon, llm), w01(klon, llm)
180    
181      ! Variables pour la couche limite (Alain Lahellec) :      ! Variables pour la couche limite (Alain Lahellec) :
182      REAL cdragh(klon) ! drag coefficient pour T and Q      REAL cdragh(klon) ! drag coefficient pour T and Q
183      REAL cdragm(klon) ! drag coefficient pour vent      REAL cdragm(klon) ! drag coefficient pour vent
184    
185      ! Pour phytrac :      REAL coefh(klon, 2:llm) ! coef d'echange pour phytrac
186      REAL ycoefh(klon, llm) ! coef d'echange pour phytrac  
187      REAL yu1(klon) ! vents dans la premiere couche U      REAL, save:: ffonte(klon, nbsrf)
188      REAL yv1(klon) ! vents dans la premiere couche V      ! flux thermique utilise pour fondre la neige
     REAL ffonte(klon, nbsrf) ! flux thermique utilise pour fondre la neige  
189    
190      REAL fqcalving(klon, nbsrf)      REAL fqcalving(klon, nbsrf)
191      ! flux d'eau "perdue" par la surface et necessaire pour limiter la      ! flux d'eau "perdue" par la surface et n\'ecessaire pour limiter
192      ! hauteur de neige, en kg/m2/s      ! la hauteur de neige, en kg / m2 / s
193    
194        REAL zxffonte(klon)
195    
196        REAL, save:: pfrac_impa(klon, llm)! Produits des coefs lessivage impaction
197        REAL, save:: pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation
198    
199      REAL zxffonte(klon), zxfqcalving(klon)      REAL, save:: pfrac_1nucl(klon, llm)
200        ! Produits des coefs lessi nucl (alpha = 1)
201    
202      REAL pfrac_impa(klon, llm)! Produits des coefs lessivage impaction      REAL frac_impa(klon, llm) ! fraction d'a\'erosols lessiv\'es (impaction)
     save pfrac_impa  
     REAL pfrac_nucl(klon, llm)! Produits des coefs lessivage nucleation  
     save pfrac_nucl  
     REAL pfrac_1nucl(klon, llm)! Produits des coefs lessi nucl (alpha = 1)  
     save pfrac_1nucl  
     REAL frac_impa(klon, llm) ! fractions d'aerosols lessivees (impaction)  
203      REAL frac_nucl(klon, llm) ! idem (nucleation)      REAL frac_nucl(klon, llm) ! idem (nucleation)
204    
205      REAL, save:: rain_fall(klon)      REAL, save:: rain_fall(klon)
206      ! liquid water mass flux (kg/m2/s), positive down      ! liquid water mass flux (kg / m2 / s), positive down
207    
208      REAL, save:: snow_fall(klon)      REAL, save:: snow_fall(klon)
209      ! solid water mass flux (kg/m2/s), positive down      ! solid water mass flux (kg / m2 / s), positive down
210    
211      REAL rain_tiedtke(klon), snow_tiedtke(klon)      REAL rain_tiedtke(klon), snow_tiedtke(klon)
212    
213      REAL evap(klon), devap(klon) ! evaporation and its derivative      REAL evap(klon) ! flux d'\'evaporation au sol
214      REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee      real dflux_q(klon) ! derivative of the evaporation flux at the surface
215      REAL dlw(klon) ! derivee infra rouge      REAL sens(klon) ! flux de chaleur sensible au sol
216      SAVE dlw      real dflux_t(klon) ! derivee du flux de chaleur sensible au sol
217        REAL, save:: dlw(klon) ! derivative of infra-red flux
218      REAL bils(klon) ! bilan de chaleur au sol      REAL bils(klon) ! bilan de chaleur au sol
219      REAL, save:: fder(klon) ! Derive de flux (sensible et latente)      REAL fder(klon) ! Derive de flux (sensible et latente)
220      REAL ve(klon) ! integr. verticale du transport meri. de l'energie      REAL ve(klon) ! integr. verticale du transport meri. de l'energie
221      REAL vq(klon) ! integr. verticale du transport meri. de l'eau      REAL vq(klon) ! integr. verticale du transport meri. de l'eau
222      REAL ue(klon) ! integr. verticale du transport zonal de l'energie      REAL ue(klon) ! integr. verticale du transport zonal de l'energie
# Line 247  contains Line 228  contains
228      ! Conditions aux limites      ! Conditions aux limites
229    
230      INTEGER julien      INTEGER julien
     INTEGER, SAVE:: lmt_pas ! number of time steps of "physics" per day  
231      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface      REAL, save:: pctsrf(klon, nbsrf) ! percentage of surface
232      REAL pctsrf_new(klon, nbsrf) ! pourcentage surfaces issus d'ORCHIDEE      REAL, save:: albsol(klon) ! albedo du sol total, visible, moyen par maille
     REAL, save:: albsol(klon) ! albedo du sol total visible  
233      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU      REAL, SAVE:: wo(klon, llm) ! column density of ozone in a cell, in kDU
234        real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
235    
236      real, save:: clwcon(klon, llm), rnebcon(klon, llm)      real, save:: clwcon(klon, llm), rnebcon(klon, llm)
237      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)      real, save:: clwcon0(klon, llm), rnebcon0(klon, llm)
238    
239      REAL rhcl(klon, llm) ! humiditi relative ciel clair      REAL rhcl(klon, llm) ! humidit\'e relative ciel clair
240      REAL dialiq(klon, llm) ! eau liquide nuageuse      REAL dialiq(klon, llm) ! eau liquide nuageuse
241      REAL diafra(klon, llm) ! fraction nuageuse      REAL diafra(klon, llm) ! fraction nuageuse
242      REAL cldliq(klon, llm) ! eau liquide nuageuse      REAL cldliq(klon, llm) ! eau liquide nuageuse
# Line 264  contains Line 244  contains
244      REAL cldtau(klon, llm) ! epaisseur optique      REAL cldtau(klon, llm) ! epaisseur optique
245      REAL cldemi(klon, llm) ! emissivite infrarouge      REAL cldemi(klon, llm) ! emissivite infrarouge
246    
247      REAL fluxq(klon, llm, nbsrf) ! flux turbulent d'humidite      REAL flux_q(klon, nbsrf) ! flux turbulent d'humidite à la surface
248      REAL fluxt(klon, llm, nbsrf) ! flux turbulent de chaleur      REAL flux_t(klon, nbsrf) ! flux turbulent de chaleur à la surface
249      REAL fluxu(klon, llm, nbsrf) ! flux turbulent de vitesse u  
250      REAL fluxv(klon, llm, nbsrf) ! flux turbulent de vitesse v      REAL flux_u(klon, nbsrf), flux_v(klon, nbsrf)
251        ! tension du vent (flux turbulent de vent) à la surface, en Pa
     REAL zxfluxt(klon, llm)  
     REAL zxfluxq(klon, llm)  
     REAL zxfluxu(klon, llm)  
     REAL zxfluxv(klon, llm)  
252    
253      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que      ! Le rayonnement n'est pas calcul\'e tous les pas, il faut donc que
254      ! les variables soient r\'emanentes.      ! les variables soient r\'emanentes.
# Line 288  contains Line 264  contains
264      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface      REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous-surface
265      REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface      REAL fsolsw(klon, nbsrf) ! flux solaire absorb\'e pour chaque sous-surface
266    
267      REAL conv_q(klon, llm) ! convergence de l'humidite (kg/kg/s)      REAL conv_q(klon, llm) ! convergence de l'humidite (kg / kg / s)
268      REAL conv_t(klon, llm) ! convergence of temperature (K/s)      REAL conv_t(klon, llm) ! convergence of temperature (K / s)
269    
270      REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut      REAL cldl(klon), cldm(klon), cldh(klon) ! nuages bas, moyen et haut
271      REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree      REAL cldt(klon), cldq(klon) ! nuage total, eau liquide integree
272    
273      REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)      REAL zxfluxlat(klon)
   
274      REAL dist, mu0(klon), fract(klon)      REAL dist, mu0(klon), fract(klon)
275      real longi      real longi
276      REAL z_avant(klon), z_apres(klon), z_factor(klon)      REAL z_avant(klon), z_apres(klon), z_factor(klon)
277      REAL za, zb      REAL zb
278      REAL zx_t, zx_qs, zcor      REAL zx_t, zx_qs, zcor
279      real zqsat(klon, llm)      real zqsat(klon, llm)
280      INTEGER i, k, iq, nsrf      INTEGER i, k, iq, nsrf
     REAL, PARAMETER:: t_coup = 234.  
281      REAL zphi(klon, llm)      REAL zphi(klon, llm)
282    
283      ! cf. Anne Mathieu variables pour la couche limite atmosphérique (hbtm)      ! cf. Anne Mathieu, variables pour la couche limite atmosphérique (hbtm)
284    
285      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite      REAL, SAVE:: pblh(klon, nbsrf) ! Hauteur de couche limite
286      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA      REAL, SAVE:: plcl(klon, nbsrf) ! Niveau de condensation de la CLA
287      REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite      REAL, SAVE:: capCL(klon, nbsrf) ! CAPE de couche limite
288      REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite      REAL, SAVE:: oliqCL(klon, nbsrf) ! eau_liqu integree de couche limite
289      REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite      REAL, SAVE:: cteiCL(klon, nbsrf) ! cloud top instab. crit. couche limite
290      REAL, SAVE:: pblt(klon, nbsrf) ! T a la Hauteur de couche limite      REAL, SAVE:: pblt(klon, nbsrf) ! T \`a la hauteur de couche limite
291      REAL, SAVE:: therm(klon, nbsrf)      REAL, SAVE:: therm(klon, nbsrf)
     REAL, SAVE:: trmb1(klon, nbsrf) ! deep_cape  
     REAL, SAVE:: trmb2(klon, nbsrf) ! inhibition  
     REAL, SAVE:: trmb3(klon, nbsrf) ! Point Omega  
292      ! Grandeurs de sorties      ! Grandeurs de sorties
293      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)      REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
294      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)      REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
295      REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)      REAL s_therm(klon)
     REAL s_trmb3(klon)  
296    
297      ! Variables pour la convection de K. Emanuel :      ! Variables pour la convection de K. Emanuel :
298    
299      REAL upwd(klon, llm) ! saturated updraft mass flux      REAL upwd(klon, llm) ! saturated updraft mass flux
300      REAL dnwd(klon, llm) ! saturated downdraft mass flux      REAL dnwd(klon, llm) ! saturated downdraft mass flux
301      REAL dnwd0(klon, llm) ! unsaturated downdraft mass flux      REAL, save:: cape(klon)
     REAL cape(klon) ! CAPE  
     SAVE cape  
302    
303      INTEGER iflagctrl(klon) ! flag fonctionnement de convect      INTEGER iflagctrl(klon) ! flag fonctionnement de convect
304    
# Line 342  contains Line 310  contains
310      ! eva: \'evaporation de l'eau liquide nuageuse      ! eva: \'evaporation de l'eau liquide nuageuse
311      ! vdf: vertical diffusion in boundary layer      ! vdf: vertical diffusion in boundary layer
312      REAL d_t_con(klon, llm), d_q_con(klon, llm)      REAL d_t_con(klon, llm), d_q_con(klon, llm)
313      REAL d_u_con(klon, llm), d_v_con(klon, llm)      REAL, save:: d_u_con(klon, llm), d_v_con(klon, llm)
314      REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)      REAL d_t_lsc(klon, llm), d_q_lsc(klon, llm), d_ql_lsc(klon, llm)
315      REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)      REAL d_t_ajs(klon, llm), d_q_ajs(klon, llm)
316      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)      REAL d_u_ajs(klon, llm), d_v_ajs(klon, llm)
# Line 358  contains Line 326  contains
326      INTEGER, save:: ibas_con(klon), itop_con(klon)      INTEGER, save:: ibas_con(klon), itop_con(klon)
327      real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa      real ema_pct(klon) ! Emanuel pressure at cloud top, in Pa
328    
329      REAL rain_con(klon), rain_lsc(klon)      REAL rain_con(klon)
330      REAL, save:: snow_con(klon) ! neige (mm / s)      real rain_lsc(klon)
331        REAL snow_con(klon) ! neige (mm / s)
332      real snow_lsc(klon)      real snow_lsc(klon)
333      REAL d_ts(klon, nbsrf)      REAL d_ts(klon, nbsrf) ! variation of ftsol
334    
335      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)      REAL d_u_vdf(klon, llm), d_v_vdf(klon, llm)
336      REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)      REAL d_t_vdf(klon, llm), d_q_vdf(klon, llm)
# Line 395  contains Line 364  contains
364    
365      REAL zustrdr(klon), zvstrdr(klon)      REAL zustrdr(klon), zvstrdr(klon)
366      REAL zustrli(klon), zvstrli(klon)      REAL zustrli(klon), zvstrli(klon)
     REAL zustrph(klon), zvstrph(klon)  
367      REAL aam, torsfc      REAL aam, torsfc
368    
369      REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.      REAL ve_lay(klon, llm) ! transport meri. de l'energie a chaque niveau vert.
# Line 403  contains Line 371  contains
371      REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.      REAL ue_lay(klon, llm) ! transport zonal de l'energie a chaque niveau vert.
372      REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.      REAL uq_lay(klon, llm) ! transport zonal de l'eau a chaque niveau vert.
373    
374      real date0      REAL tsol(klon)
375    
376      ! Variables li\'ees au bilan d'\'energie et d'enthalpie :      REAL d_t_ec(klon, llm)
377      REAL ztsol(klon)      ! tendance due \`a la conversion d'\'energie cin\'etique en
378      REAL d_h_vcol, d_qt, d_ec      ! énergie thermique
379      REAL, SAVE:: d_h_vcol_phy  
380      REAL zero_v(klon)      REAL, save:: t2m(klon, nbsrf), q2m(klon, nbsrf)
381      CHARACTER(LEN = 20) tit      ! temperature and humidity at 2 m
382      INTEGER:: ip_ebil = 0 ! print level for energy conservation diagnostics  
383      INTEGER:: if_ebil = 0 ! verbosity for diagnostics of energy conservation      REAL, save:: u10m_srf(klon, nbsrf), v10m_srf(klon, nbsrf)
384        ! composantes du vent \`a 10 m
385      REAL d_t_ec(klon, llm) ! tendance due \`a la conversion Ec -> E thermique      
386      REAL ZRCPD      REAL zt2m(klon), zq2m(klon) ! température, humidité 2 m moyenne sur 1 maille
387        REAL u10m(klon), v10m(klon) ! vent \`a 10 m moyenn\' sur les sous-surfaces
     REAL t2m(klon, nbsrf), q2m(klon, nbsrf) ! temperature and humidity at 2 m  
     REAL u10m(klon, nbsrf), v10m(klon, nbsrf) ! vents a 10 m  
     REAL zt2m(klon), zq2m(klon) ! temp., hum. 2 m moyenne s/ 1 maille  
     REAL zu10m(klon), zv10m(klon) ! vents a 10 m moyennes s/1 maille  
388    
389      ! Aerosol effects:      ! Aerosol effects:
390    
391      REAL sulfate(klon, llm) ! SO4 aerosol concentration (micro g/m3)      REAL, save:: topswad(klon), solswad(klon) ! aerosol direct effect
   
     REAL, save:: sulfate_pi(klon, llm)  
     ! SO4 aerosol concentration, in \mu g/m3, pre-industrial value  
   
     REAL cldtaupi(klon, llm)  
     ! cloud optical thickness for pre-industrial aerosols  
   
     REAL re(klon, llm) ! Cloud droplet effective radius  
     REAL fl(klon, llm) ! denominator of re  
   
     ! Aerosol optical properties  
     REAL, save:: tau_ae(klon, llm, 2), piz_ae(klon, llm, 2)  
     REAL, save:: cg_ae(klon, llm, 2)  
   
     REAL topswad(klon), solswad(klon) ! aerosol direct effect  
     REAL topswai(klon), solswai(klon) ! aerosol indirect effect  
   
     REAL aerindex(klon) ! POLDER aerosol index  
   
392      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect      LOGICAL:: ok_ade = .false. ! apply aerosol direct effect
     LOGICAL:: ok_aie = .false. ! apply aerosol indirect effect  
393    
394      REAL:: bl95_b0 = 2., bl95_b1 = 0.2      REAL:: bl95_b0 = 2., bl95_b1 = 0.2
395      ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus      ! Parameters in equation (D) of Boucher and Lohmann (1995, Tellus
396      ! B). They link cloud droplet number concentration to aerosol mass      ! B). They link cloud droplet number concentration to aerosol mass
397      ! concentration.      ! concentration.
398    
     SAVE u10m  
     SAVE v10m  
     SAVE t2m  
     SAVE q2m  
     SAVE ffonte  
     SAVE fqcalving  
     SAVE rain_con  
     SAVE topswai  
     SAVE topswad  
     SAVE solswai  
     SAVE solswad  
     SAVE d_u_con  
     SAVE d_v_con  
   
399      real zmasse(klon, llm)      real zmasse(klon, llm)
400      ! (column-density of mass of air in a cell, in kg m-2)      ! (column-density of mass of air in a cell, in kg m-2)
401    
402      integer, save:: ncid_startphy      integer, save:: ncid_startphy
403    
404      namelist /physiq_nml/ fact_cldcon, facttemps, ok_newmicro, &      namelist /physiq_nml/ fact_cldcon, facttemps, ok_newmicro, iflag_cldcon, &
405           iflag_cldcon, ratqsbas, ratqshaut, if_ebil, ok_ade, ok_aie, bl95_b0, &           ratqsbas, ratqshaut, ok_ade, bl95_b0, bl95_b1, iflag_thermals, &
406           bl95_b1, iflag_thermals, nsplit_thermals           nsplit_thermals
407    
408      !----------------------------------------------------------------      !----------------------------------------------------------------
409    
     IF (if_ebil >= 1) zero_v = 0.  
410      IF (nqmx < 2) CALL abort_gcm('physiq', &      IF (nqmx < 2) CALL abort_gcm('physiq', &
411           'eaux vapeur et liquide sont indispensables')           'eaux vapeur et liquide sont indispensables')
412    
413      test_firstcal: IF (firstcal) THEN      test_firstcal: IF (firstcal) THEN
414         ! initialiser         ! initialiser
415         u10m = 0.         u10m_srf = 0.
416         v10m = 0.         v10m_srf = 0.
417         t2m = 0.         t2m = 0.
418         q2m = 0.         q2m = 0.
419         ffonte = 0.         ffonte = 0.
        fqcalving = 0.  
        piz_ae = 0.  
        tau_ae = 0.  
        cg_ae = 0.  
        rain_con = 0.  
        snow_con = 0.  
        topswai = 0.  
        topswad = 0.  
        solswai = 0.  
        solswad = 0.  
   
420         d_u_con = 0.         d_u_con = 0.
421         d_v_con = 0.         d_v_con = 0.
422         rnebcon0 = 0.         rnebcon0 = 0.
423         clwcon0 = 0.         clwcon0 = 0.
424         rnebcon = 0.         rnebcon = 0.
425         clwcon = 0.         clwcon = 0.
   
426         pblh =0. ! Hauteur de couche limite         pblh =0. ! Hauteur de couche limite
427         plcl =0. ! Niveau de condensation de la CLA         plcl =0. ! Niveau de condensation de la CLA
428         capCL =0. ! CAPE de couche limite         capCL =0. ! CAPE de couche limite
429         oliqCL =0. ! eau_liqu integree de couche limite         oliqCL =0. ! eau_liqu integree de couche limite
430         cteiCL =0. ! cloud top instab. crit. couche limite         cteiCL =0. ! cloud top instab. crit. couche limite
431         pblt =0. ! T a la Hauteur de couche limite         pblt =0.
432         therm =0.         therm =0.
        trmb1 =0. ! deep_cape  
        trmb2 =0. ! inhibition  
        trmb3 =0. ! Point Omega  
   
        IF (if_ebil >= 1) d_h_vcol_phy = 0.  
433    
434         iflag_thermals = 0         iflag_thermals = 0
435         nsplit_thermals = 1         nsplit_thermals = 1
# Line 531  contains Line 443  contains
443    
444         frugs = 0.         frugs = 0.
445         CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &         CALL phyetat0(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
446              fevap, rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, &              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, &
447              agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &              zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &
448              q_ancien, ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &              ancien_ok, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01, &
449              w01, ncid_startphy)              ncid_startphy)
450    
451         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial         ! ATTENTION : il faudra a terme relire q2 dans l'etat initial
452         q2 = 1e-8         q2 = 1e-8
453    
        lmt_pas = day_step / iphysiq  
        print *, 'Number of time steps of "physics" per day: ', lmt_pas  
   
454         radpas = lmt_pas / nbapp_rad         radpas = lmt_pas / nbapp_rad
455         print *, "radpas = ", radpas         print *, "radpas = ", radpas
456    
# Line 558  contains Line 467  contains
467            rugoro = 0.            rugoro = 0.
468         ENDIF         ENDIF
469    
        ecrit_ins = NINT(ecrit_ins/dtphys)  
        ecrit_hf = NINT(ecrit_hf/dtphys)  
        ecrit_mth = NINT(ecrit_mth/dtphys)  
        ecrit_tra = NINT(86400.*ecrit_tra/dtphys)  
        ecrit_reg = NINT(ecrit_reg/dtphys)  
   
470         ! Initialisation des sorties         ! Initialisation des sorties
471           call ini_histins(ok_newmicro)
472         call ini_histins(dtphys)         CALL phyredem0
473         CALL ymds2ju(annee_ref, 1, day_ref, 0., date0)         call conf_interface
        ! Positionner date0 pour initialisation de ORCHIDEE  
        print *, 'physiq date0: ', date0  
        CALL phyredem0(lmt_pas)  
474      ENDIF test_firstcal      ENDIF test_firstcal
475    
476      ! We will modify variables *_seri and we will not touch variables      ! We will modify variables *_seri and we will not touch variables
# Line 582  contains Line 482  contains
482      ql_seri = qx(:, :, iliq)      ql_seri = qx(:, :, iliq)
483      tr_seri = qx(:, :, 3:nqmx)      tr_seri = qx(:, :, 3:nqmx)
484    
485      ztsol = sum(ftsol * pctsrf, dim = 2)      tsol = sum(ftsol * pctsrf, dim = 2)
   
     IF (if_ebil >= 1) THEN  
        tit = 'after dynamics'  
        CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &  
             ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)  
        ! Comme les tendances de la physique sont ajout\'es dans la  
        ! dynamique, la variation d'enthalpie par la dynamique devrait  
        ! \^etre \'egale \`a la variation de la physique au pas de temps  
        ! pr\'ec\'edent. Donc la somme de ces 2 variations devrait \^etre  
        ! nulle.  
        call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &  
             zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol + d_h_vcol_phy, &  
             d_qt, 0.)  
     END IF  
486    
487      ! Diagnostic de la tendance dynamique :      ! Diagnostic de la tendance dynamique :
488      IF (ancien_ok) THEN      IF (ancien_ok) THEN
# Line 632  contains Line 518  contains
518    
519      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg      forall (k = 1: llm) zmasse(:, k) = (paprs(:, k) - paprs(:, k + 1)) / rg
520    
     ! Prescrire l'ozone :  
     wo = ozonecm(REAL(julien), paprs)  
   
521      ! \'Evaporation de l'eau liquide nuageuse :      ! \'Evaporation de l'eau liquide nuageuse :
522      DO k = 1, llm      DO k = 1, llm
523         DO i = 1, klon         DO i = 1, klon
# Line 646  contains Line 529  contains
529      ENDDO      ENDDO
530      ql_seri = 0.      ql_seri = 0.
531    
     IF (if_ebil >= 2) THEN  
        tit = 'after reevap'  
        CALL diagetpq(airephy, tit, ip_ebil, 2, 1, dtphys, t_seri, q_seri, &  
             ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)  
        call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &  
             zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)  
     END IF  
   
532      frugs = MAX(frugs, 0.000015)      frugs = MAX(frugs, 0.000015)
533      zxrugs = sum(frugs * pctsrf, dim = 2)      zxrugs = sum(frugs * pctsrf, dim = 2)
534    
# Line 661  contains Line 536  contains
536      ! la surface.      ! la surface.
537    
538      CALL orbite(REAL(julien), longi, dist)      CALL orbite(REAL(julien), longi, dist)
539      IF (cycle_diurne) THEN      CALL zenang(longi, time, dtphys * radpas, mu0, fract)
        CALL zenang(longi, time, dtphys * radpas, mu0, fract)  
     ELSE  
        mu0 = - 999.999  
     ENDIF  
   
     ! Calcul de l'abedo moyen par maille  
540      albsol = sum(falbe * pctsrf, dim = 2)      albsol = sum(falbe * pctsrf, dim = 2)
541    
542      ! R\'epartition sous maille des flux longwave et shortwave      ! R\'epartition sous maille des flux longwave et shortwave
543      ! R\'epartition du longwave par sous-surface lin\'earis\'ee      ! R\'epartition du longwave par sous-surface lin\'earis\'ee
544    
545      forall (nsrf = 1: nbsrf)      forall (nsrf = 1: nbsrf)
546         fsollw(:, nsrf) = sollw + 4. * RSIGMA * ztsol**3 &         fsollw(:, nsrf) = sollw + 4. * RSIGMA * tsol**3 &
547              * (ztsol - ftsol(:, nsrf))              * (tsol - ftsol(:, nsrf))
548         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)         fsolsw(:, nsrf) = solsw * (1. - falbe(:, nsrf)) / (1. - albsol)
549      END forall      END forall
550    
551      fder = dlw      CALL pbl_surface(pctsrf, t_seri, q_seri, u_seri, v_seri, julien, mu0, &
552             ftsol, cdmmax, cdhmax, ftsoil, qsol, paprs, play, fsnow, fqsurf, &
553      ! Couche limite:           falbe, fluxlat, rain_fall, snow_fall, fsolsw, fsollw, frugs, agesno, &
554             rugoro, d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, flux_t, flux_q, &
555      CALL clmain(dtphys, pctsrf, pctsrf_new, t_seri, q_seri, u_seri, v_seri, &           flux_u, flux_v, cdragh, cdragm, q2, dflux_t, dflux_q, coefh, t2m, &
556           julien, mu0, ftsol, cdmmax, cdhmax, ksta, ksta_ter, ok_kzmin, &           q2m, u10m_srf, v10m_srf, pblh, capCL, oliqCL, cteiCL, pblT, therm, &
557           ftsoil, qsol, paprs, play, fsnow, fqsurf, fevap, falbe, fluxlat, &           plcl, fqcalving, ffonte, run_off_lic_0)
          rain_fall, snow_fall, fsolsw, fsollw, fder, rlat, frugs, firstcal, &  
          agesno, rugoro, d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_ts, fluxt, &  
          fluxq, fluxu, fluxv, cdragh, cdragm, q2, dsens, devap, ycoefh, yu1, &  
          yv1, t2m, q2m, u10m, v10m, pblh, capCL, oliqCL, cteiCL, pblT, therm, &  
          trmb1, trmb2, trmb3, plcl, fqcalving, ffonte, run_off_lic_0)  
558    
559      ! Incr\'ementation des flux      ! Incr\'ementation des flux
560    
561      zxfluxt = 0.      sens = - sum(flux_t * pctsrf, dim = 2)
562      zxfluxq = 0.      evap = - sum(flux_q * pctsrf, dim = 2)
563      zxfluxu = 0.      fder = dlw + dflux_t + dflux_q
     zxfluxv = 0.  
     DO nsrf = 1, nbsrf  
        DO k = 1, llm  
           DO i = 1, klon  
              zxfluxt(i, k) = zxfluxt(i, k) + fluxt(i, k, nsrf) * pctsrf(i, nsrf)  
              zxfluxq(i, k) = zxfluxq(i, k) + fluxq(i, k, nsrf) * pctsrf(i, nsrf)  
              zxfluxu(i, k) = zxfluxu(i, k) + fluxu(i, k, nsrf) * pctsrf(i, nsrf)  
              zxfluxv(i, k) = zxfluxv(i, k) + fluxv(i, k, nsrf) * pctsrf(i, nsrf)  
           END DO  
        END DO  
     END DO  
     DO i = 1, klon  
        sens(i) = - zxfluxt(i, 1) ! flux de chaleur sensible au sol  
        evap(i) = - zxfluxq(i, 1) ! flux d'\'evaporation au sol  
        fder(i) = dlw(i) + dsens(i) + devap(i)  
     ENDDO  
564    
565      DO k = 1, llm      DO k = 1, llm
566         DO i = 1, klon         DO i = 1, klon
# Line 723  contains Line 571  contains
571         ENDDO         ENDDO
572      ENDDO      ENDDO
573    
     IF (if_ebil >= 2) THEN  
        tit = 'after clmain'  
        CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)  
        call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &  
             sens, evap, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)  
     END IF  
   
     ! Update surface temperature:  
   
     DO i = 1, klon  
        zxtsol(i) = 0.  
        zxfluxlat(i) = 0.  
   
        zt2m(i) = 0.  
        zq2m(i) = 0.  
        zu10m(i) = 0.  
        zv10m(i) = 0.  
        zxffonte(i) = 0.  
        zxfqcalving(i) = 0.  
   
        s_pblh(i) = 0.  
        s_lcl(i) = 0.  
        s_capCL(i) = 0.  
        s_oliqCL(i) = 0.  
        s_cteiCL(i) = 0.  
        s_pblT(i) = 0.  
        s_therm(i) = 0.  
        s_trmb1(i) = 0.  
        s_trmb2(i) = 0.  
        s_trmb3(i) = 0.  
     ENDDO  
   
574      call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')      call assert(abs(sum(pctsrf, dim = 2) - 1.) <= EPSFRA, 'physiq: pctsrf')
575        ftsol = ftsol + d_ts ! update surface temperature
576        tsol = sum(ftsol * pctsrf, dim = 2)
577        zxfluxlat = sum(fluxlat * pctsrf, dim = 2)
578        zt2m = sum(t2m * pctsrf, dim = 2)
579        zq2m = sum(q2m * pctsrf, dim = 2)
580        u10m = sum(u10m_srf * pctsrf, dim = 2)
581        v10m = sum(v10m_srf * pctsrf, dim = 2)
582        zxffonte = sum(ffonte * pctsrf, dim = 2)
583        s_pblh = sum(pblh * pctsrf, dim = 2)
584        s_lcl = sum(plcl * pctsrf, dim = 2)
585        s_capCL = sum(capCL * pctsrf, dim = 2)
586        s_oliqCL = sum(oliqCL * pctsrf, dim = 2)
587        s_cteiCL = sum(cteiCL * pctsrf, dim = 2)
588        s_pblT = sum(pblT * pctsrf, dim = 2)
589        s_therm = sum(therm * pctsrf, dim = 2)
590    
591        ! Si une sous-fraction n'existe pas, elle prend la valeur moyenne :
592      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
593         DO i = 1, klon         DO i = 1, klon
594            ftsol(i, nsrf) = ftsol(i, nsrf) + d_ts(i, nsrf)            IF (pctsrf(i, nsrf) < epsfra) then
595            zxtsol(i) = zxtsol(i) + ftsol(i, nsrf)*pctsrf(i, nsrf)               ftsol(i, nsrf) = tsol(i)
596            zxfluxlat(i) = zxfluxlat(i) + fluxlat(i, nsrf)*pctsrf(i, nsrf)               t2m(i, nsrf) = zt2m(i)
597                 q2m(i, nsrf) = zq2m(i)
598            zt2m(i) = zt2m(i) + t2m(i, nsrf)*pctsrf(i, nsrf)               u10m_srf(i, nsrf) = u10m(i)
599            zq2m(i) = zq2m(i) + q2m(i, nsrf)*pctsrf(i, nsrf)               v10m_srf(i, nsrf) = v10m(i)
600            zu10m(i) = zu10m(i) + u10m(i, nsrf)*pctsrf(i, nsrf)               ffonte(i, nsrf) = zxffonte(i)
601            zv10m(i) = zv10m(i) + v10m(i, nsrf)*pctsrf(i, nsrf)               pblh(i, nsrf) = s_pblh(i)
602            zxffonte(i) = zxffonte(i) + ffonte(i, nsrf)*pctsrf(i, nsrf)               plcl(i, nsrf) = s_lcl(i)
603            zxfqcalving(i) = zxfqcalving(i) + &               capCL(i, nsrf) = s_capCL(i)
604                 fqcalving(i, nsrf)*pctsrf(i, nsrf)               oliqCL(i, nsrf) = s_oliqCL(i)
605            s_pblh(i) = s_pblh(i) + pblh(i, nsrf)*pctsrf(i, nsrf)               cteiCL(i, nsrf) = s_cteiCL(i)
606            s_lcl(i) = s_lcl(i) + plcl(i, nsrf)*pctsrf(i, nsrf)               pblT(i, nsrf) = s_pblT(i)
607            s_capCL(i) = s_capCL(i) + capCL(i, nsrf) *pctsrf(i, nsrf)               therm(i, nsrf) = s_therm(i)
608            s_oliqCL(i) = s_oliqCL(i) + oliqCL(i, nsrf) *pctsrf(i, nsrf)            end IF
           s_cteiCL(i) = s_cteiCL(i) + cteiCL(i, nsrf) *pctsrf(i, nsrf)  
           s_pblT(i) = s_pblT(i) + pblT(i, nsrf) *pctsrf(i, nsrf)  
           s_therm(i) = s_therm(i) + therm(i, nsrf) *pctsrf(i, nsrf)  
           s_trmb1(i) = s_trmb1(i) + trmb1(i, nsrf) *pctsrf(i, nsrf)  
           s_trmb2(i) = s_trmb2(i) + trmb2(i, nsrf) *pctsrf(i, nsrf)  
           s_trmb3(i) = s_trmb3(i) + trmb3(i, nsrf) *pctsrf(i, nsrf)  
609         ENDDO         ENDDO
610      ENDDO      ENDDO
611    
612      ! Si une sous-fraction n'existe pas, elle prend la température moyenne :      dlw = - 4. * RSIGMA * tsol**3
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           IF (pctsrf(i, nsrf) < epsfra) ftsol(i, nsrf) = zxtsol(i)  
   
           IF (pctsrf(i, nsrf) < epsfra) t2m(i, nsrf) = zt2m(i)  
           IF (pctsrf(i, nsrf) < epsfra) q2m(i, nsrf) = zq2m(i)  
           IF (pctsrf(i, nsrf) < epsfra) u10m(i, nsrf) = zu10m(i)  
           IF (pctsrf(i, nsrf) < epsfra) v10m(i, nsrf) = zv10m(i)  
           IF (pctsrf(i, nsrf) < epsfra) ffonte(i, nsrf) = zxffonte(i)  
           IF (pctsrf(i, nsrf) < epsfra) &  
                fqcalving(i, nsrf) = zxfqcalving(i)  
           IF (pctsrf(i, nsrf) < epsfra) pblh(i, nsrf) = s_pblh(i)  
           IF (pctsrf(i, nsrf) < epsfra) plcl(i, nsrf) = s_lcl(i)  
           IF (pctsrf(i, nsrf) < epsfra) capCL(i, nsrf) = s_capCL(i)  
           IF (pctsrf(i, nsrf) < epsfra) oliqCL(i, nsrf) = s_oliqCL(i)  
           IF (pctsrf(i, nsrf) < epsfra) cteiCL(i, nsrf) = s_cteiCL(i)  
           IF (pctsrf(i, nsrf) < epsfra) pblT(i, nsrf) = s_pblT(i)  
           IF (pctsrf(i, nsrf) < epsfra) therm(i, nsrf) = s_therm(i)  
           IF (pctsrf(i, nsrf) < epsfra) trmb1(i, nsrf) = s_trmb1(i)  
           IF (pctsrf(i, nsrf) < epsfra) trmb2(i, nsrf) = s_trmb2(i)  
           IF (pctsrf(i, nsrf) < epsfra) trmb3(i, nsrf) = s_trmb3(i)  
        ENDDO  
     ENDDO  
   
     ! Calculer la dérive du flux infrarouge  
   
     DO i = 1, klon  
        dlw(i) = - 4. * RSIGMA * zxtsol(i)**3  
     ENDDO  
   
     IF (check) print *, "avantcon = ", qcheck(paprs, q_seri, ql_seri)  
613    
614      ! Appeler la convection      ! Appeler la convection
615    
616      if (conv_emanuel) then      if (conv_emanuel) then
617         da = 0.         CALL concvl(paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, w01, &
618         mp = 0.              d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, itop_con, &
619         phi = 0.              upwd, dnwd, Ma, cape, iflagctrl, clwcon0, pmflxr, da, phi, mp)
        CALL concvl(dtphys, paprs, play, t_seri, q_seri, u_seri, v_seri, sig1, &  
             w01, d_t_con, d_q_con, d_u_con, d_v_con, rain_con, ibas_con, &  
             itop_con, upwd, dnwd, dnwd0, Ma, cape, iflagctrl, qcondc, pmflxr, &  
             da, phi, mp)  
620         snow_con = 0.         snow_con = 0.
        clwcon0 = qcondc  
621         mfu = upwd + dnwd         mfu = upwd + dnwd
622    
623         IF (thermcep) THEN         zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)
624            zqsat = MIN(0.5, r2es * FOEEW(t_seri, rtt >= t_seri) / play)         zqsat = zqsat / (1. - retv * zqsat)
           zqsat = zqsat / (1. - retv * zqsat)  
        ELSE  
           zqsat = merge(qsats(t_seri), qsatl(t_seri), t_seri < t_coup) / play  
        ENDIF  
625    
626         ! Properties of convective clouds         ! Properties of convective clouds
627         clwcon0 = fact_cldcon * clwcon0         clwcon0 = fact_cldcon * clwcon0
# Line 853  contains Line 638  contains
638         conv_q = d_q_dyn + d_q_vdf / dtphys         conv_q = d_q_dyn + d_q_vdf / dtphys
639         conv_t = d_t_dyn + d_t_vdf / dtphys         conv_t = d_t_dyn + d_t_vdf / dtphys
640         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)         z_avant = sum((q_seri + ql_seri) * zmasse, dim=2)
641         CALL conflx(dtphys, paprs, play, t_seri(:, llm:1:- 1), &         CALL conflx(paprs, play, t_seri(:, llm:1:- 1), q_seri(:, llm:1:- 1), &
642              q_seri(:, llm:1:- 1), conv_t, conv_q, zxfluxq(:, 1), omega, &              conv_t, conv_q, - evap, omega, d_t_con, d_q_con, rain_con, &
643              d_t_con, d_q_con, rain_con, snow_con, mfu(:, llm:1:- 1), &              snow_con, mfu(:, llm:1:- 1), mfd(:, llm:1:- 1), pen_u, pde_u, &
644              mfd(:, llm:1:- 1), pen_u, pde_u, pen_d, pde_d, kcbot, kctop, &              pen_d, pde_d, kcbot, kctop, kdtop, pmflxr, pmflxs)
             kdtop, pmflxr, pmflxs)  
645         WHERE (rain_con < 0.) rain_con = 0.         WHERE (rain_con < 0.) rain_con = 0.
646         WHERE (snow_con < 0.) snow_con = 0.         WHERE (snow_con < 0.) snow_con = 0.
647         ibas_con = llm + 1 - kcbot         ibas_con = llm + 1 - kcbot
# Line 873  contains Line 657  contains
657         ENDDO         ENDDO
658      ENDDO      ENDDO
659    
     IF (if_ebil >= 2) THEN  
        tit = 'after convect'  
        CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)  
        call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &  
             zero_v, zero_v, rain_con, snow_con, ztsol, d_h_vcol, d_qt, d_ec)  
     END IF  
   
     IF (check) THEN  
        za = qcheck(paprs, q_seri, ql_seri)  
        print *, "aprescon = ", za  
        zx_t = 0.  
        za = 0.  
        DO i = 1, klon  
           za = za + airephy(i)/REAL(klon)  
           zx_t = zx_t + (rain_con(i)+ &  
                snow_con(i))*airephy(i)/REAL(klon)  
        ENDDO  
        zx_t = zx_t/za*dtphys  
        print *, "Precip = ", zx_t  
     ENDIF  
   
660      IF (.not. conv_emanuel) THEN      IF (.not. conv_emanuel) THEN
661         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)         z_apres = sum((q_seri + ql_seri) * zmasse, dim=2)
662         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres         z_factor = (z_avant - (rain_con + snow_con) * dtphys) / z_apres
# Line 922  contains Line 684  contains
684         t_seri = t_seri + d_t_ajs         t_seri = t_seri + d_t_ajs
685         q_seri = q_seri + d_q_ajs         q_seri = q_seri + d_q_ajs
686      else      else
687         ! Thermiques         call calltherm(play, paprs, pphi, u_seri, v_seri, t_seri, q_seri, &
688         call calltherm(dtphys, play, paprs, pphi, u_seri, v_seri, t_seri, &              d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)
             q_seri, d_u_ajs, d_v_ajs, d_t_ajs, d_q_ajs, fm_therm, entr_therm)  
689      endif      endif
690    
     IF (if_ebil >= 2) THEN  
        tit = 'after dry_adjust'  
        CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)  
     END IF  
   
691      ! Caclul des ratqs      ! Caclul des ratqs
692    
     ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q  
     ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno  
693      if (iflag_cldcon == 1) then      if (iflag_cldcon == 1) then
694           ! ratqs convectifs \`a l'ancienne en fonction de (q(z = 0) - q) / q
695           ! on \'ecrase le tableau ratqsc calcul\'e par clouds_gno
696         do k = 1, llm         do k = 1, llm
697            do i = 1, klon            do i = 1, klon
698               if(ptconv(i, k)) then               if(ptconv(i, k)) then
# Line 971  contains Line 726  contains
726         ratqs = ratqss         ratqs = ratqss
727      endif      endif
728    
729      CALL fisrtilp(dtphys, paprs, play, t_seri, q_seri, ptconv, ratqs, &      CALL fisrtilp(paprs, play, t_seri, q_seri, ptconv, ratqs, d_t_lsc, &
730           d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, &           d_q_lsc, d_ql_lsc, rneb, cldliq, rain_lsc, snow_lsc, pfrac_impa, &
731           pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, &           pfrac_nucl, pfrac_1nucl, frac_impa, frac_nucl, prfl, psfl, rhcl)
          psfl, rhcl)  
732    
733      WHERE (rain_lsc < 0) rain_lsc = 0.      WHERE (rain_lsc < 0) rain_lsc = 0.
734      WHERE (snow_lsc < 0) snow_lsc = 0.      WHERE (snow_lsc < 0) snow_lsc = 0.
# Line 987  contains Line 741  contains
741            IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)            IF (.NOT.new_oliq) cldliq(i, k) = ql_seri(i, k)
742         ENDDO         ENDDO
743      ENDDO      ENDDO
     IF (check) THEN  
        za = qcheck(paprs, q_seri, ql_seri)  
        print *, "apresilp = ", za  
        zx_t = 0.  
        za = 0.  
        DO i = 1, klon  
           za = za + airephy(i)/REAL(klon)  
           zx_t = zx_t + (rain_lsc(i) &  
                + snow_lsc(i))*airephy(i)/REAL(klon)  
        ENDDO  
        zx_t = zx_t/za*dtphys  
        print *, "Precip = ", zx_t  
     ENDIF  
   
     IF (if_ebil >= 2) THEN  
        tit = 'after fisrt'  
        CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)  
        call diagphy(airephy, tit, ip_ebil, zero_v, zero_v, zero_v, zero_v, &  
             zero_v, zero_v, rain_lsc, snow_lsc, ztsol, d_h_vcol, d_qt, d_ec)  
     END IF  
744    
745      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT      ! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
746    
# Line 1023  contains Line 756  contains
756            do k = 1, llm            do k = 1, llm
757               do i = 1, klon               do i = 1, klon
758                  if (d_q_con(i, k) < 0.) then                  if (d_q_con(i, k) < 0.) then
759                     rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k)/dtphys &                     rain_tiedtke(i) = rain_tiedtke(i) - d_q_con(i, k) / dtphys &
760                          *zmasse(i, k)                          * zmasse(i, k)
761                  endif                  endif
762               enddo               enddo
763            enddo            enddo
# Line 1059  contains Line 792  contains
792    
793         ! On prend la somme des fractions nuageuses et des contenus en eau         ! On prend la somme des fractions nuageuses et des contenus en eau
794         cldfra = min(max(cldfra, rnebcon), 1.)         cldfra = min(max(cldfra, rnebcon), 1.)
795         cldliq = cldliq + rnebcon*clwcon         cldliq = cldliq + rnebcon * clwcon
796      ENDIF      ENDIF
797    
798      ! 2. Nuages stratiformes      ! 2. Nuages stratiformes
# Line 1082  contains Line 815  contains
815         snow_fall(i) = snow_con(i) + snow_lsc(i)         snow_fall(i) = snow_con(i) + snow_lsc(i)
816      ENDDO      ENDDO
817    
     IF (if_ebil >= 2) CALL diagetpq(airephy, "after diagcld", ip_ebil, 2, 2, &  
          dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &  
          d_qt, d_ec)  
   
818      ! Humidit\'e relative pour diagnostic :      ! Humidit\'e relative pour diagnostic :
819      DO k = 1, llm      DO k = 1, llm
820         DO i = 1, klon         DO i = 1, klon
821            zx_t = t_seri(i, k)            zx_t = t_seri(i, k)
822            IF (thermcep) THEN            zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t) / play(i, k)
823               zx_qs = r2es * FOEEW(zx_t, rtt >= zx_t)/play(i, k)            zx_qs = MIN(0.5, zx_qs)
824               zx_qs = MIN(0.5, zx_qs)            zcor = 1. / (1. - retv * zx_qs)
825               zcor = 1./(1. - retv*zx_qs)            zx_qs = zx_qs * zcor
826               zx_qs = zx_qs*zcor            zx_rh(i, k) = q_seri(i, k) / zx_qs
           ELSE  
              IF (zx_t < t_coup) THEN  
                 zx_qs = qsats(zx_t)/play(i, k)  
              ELSE  
                 zx_qs = qsatl(zx_t)/play(i, k)  
              ENDIF  
           ENDIF  
           zx_rh(i, k) = q_seri(i, k)/zx_qs  
827            zqsat(i, k) = zx_qs            zqsat(i, k) = zx_qs
828         ENDDO         ENDDO
829      ENDDO      ENDDO
830    
     ! Introduce the aerosol direct and first indirect radiative forcings:  
     IF (ok_ade .OR. ok_aie) THEN  
        ! Get sulfate aerosol distribution :  
        CALL readsulfate(dayvrai, time, firstcal, sulfate)  
        CALL readsulfate_preind(dayvrai, time, firstcal, sulfate_pi)  
   
        CALL aeropt(play, paprs, t_seri, sulfate, rhcl, tau_ae, piz_ae, cg_ae, &  
             aerindex)  
     ELSE  
        tau_ae = 0.  
        piz_ae = 0.  
        cg_ae = 0.  
     ENDIF  
   
831      ! Param\`etres optiques des nuages et quelques param\`etres pour      ! Param\`etres optiques des nuages et quelques param\`etres pour
832      ! diagnostics :      ! diagnostics :
833      if (ok_newmicro) then      if (ok_newmicro) then
834         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &         CALL newmicro(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, &
835              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc, ok_aie, &              cldh, cldl, cldm, cldt, cldq, flwp, fiwp, flwc, fiwc)
             sulfate, sulfate_pi, bl95_b0, bl95_b1, cldtaupi, re, fl)  
836      else      else
837         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &         CALL nuage(paprs, play, t_seri, cldliq, cldfra, cldtau, cldemi, cldh, &
838              cldl, cldm, cldt, cldq, ok_aie, sulfate, sulfate_pi, bl95_b0, &              cldl, cldm, cldt, cldq)
             bl95_b1, cldtaupi, re, fl)  
839      endif      endif
840    
841      IF (MOD(itap - 1, radpas) == 0) THEN      IF (MOD(itap - 1, radpas) == 0) THEN
842         ! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.         wo = ozonecm(REAL(julien), paprs)
        ! Calcul de l'abedo moyen par maille  
843         albsol = sum(falbe * pctsrf, dim = 2)         albsol = sum(falbe * pctsrf, dim = 2)
844           CALL radlwsw(dist, mu0, fract, paprs, play, tsol, albsol, t_seri, &
        ! Rayonnement (compatible Arpege-IFS) :  
        CALL radlwsw(dist, mu0, fract, paprs, play, zxtsol, albsol, t_seri, &  
845              q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &              q_seri, wo, cldfra, cldemi, cldtau, heat, heat0, cool, cool0, &
846              radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &              radsol, albpla, topsw, toplw, solsw, sollw, sollwdown, topsw0, &
847              toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &              toplw0, solsw0, sollw0, lwdn0, lwdn, lwup0, lwup, swdn0, swdn, &
848              swup0, swup, ok_ade, ok_aie, tau_ae, piz_ae, cg_ae, topswad, &              swup0, swup, ok_ade, topswad, solswad)
             solswad, cldtaupi, topswai, solswai)  
849      ENDIF      ENDIF
850    
851      ! Ajouter la tendance des rayonnements (tous les pas)      ! Ajouter la tendance des rayonnements (tous les pas)
   
852      DO k = 1, llm      DO k = 1, llm
853         DO i = 1, klon         DO i = 1, klon
854            t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys/86400.            t_seri(i, k) = t_seri(i, k) + (heat(i, k) - cool(i, k)) * dtphys &
855         ENDDO                 / 86400.
     ENDDO  
   
     IF (if_ebil >= 2) THEN  
        tit = 'after rad'  
        CALL diagetpq(airephy, tit, ip_ebil, 2, 2, dtphys, t_seri, q_seri, &  
             ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)  
        call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, &  
             zero_v, zero_v, zero_v, zero_v, ztsol, d_h_vcol, d_qt, d_ec)  
     END IF  
   
     ! Calculer l'hydrologie de la surface  
     DO i = 1, klon  
        zxqsurf(i) = 0.  
        zxsnow(i) = 0.  
     ENDDO  
     DO nsrf = 1, nbsrf  
        DO i = 1, klon  
           zxqsurf(i) = zxqsurf(i) + fqsurf(i, nsrf)*pctsrf(i, nsrf)  
           zxsnow(i) = zxsnow(i) + fsnow(i, nsrf)*pctsrf(i, nsrf)  
856         ENDDO         ENDDO
857      ENDDO      ENDDO
858    
859      ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)      ! Calculer le bilan du sol et la d\'erive de temp\'erature (couplage)
   
860      DO i = 1, klon      DO i = 1, klon
861         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)         bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
862      ENDDO      ENDDO
# Line 1185  contains Line 865  contains
865    
866      IF (ok_orodr) THEN      IF (ok_orodr) THEN
867         ! S\'election des points pour lesquels le sch\'ema est actif :         ! S\'election des points pour lesquels le sch\'ema est actif :
        igwd = 0  
868         DO i = 1, klon         DO i = 1, klon
869            itest(i) = 0            ktest(i) = 0
870            IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN            IF (zpic(i) - zmea(i) > 100. .AND. zstd(i) > 10.) THEN
871               itest(i) = 1               ktest(i) = 1
              igwd = igwd + 1  
872            ENDIF            ENDIF
873         ENDDO         ENDDO
874    
875         CALL drag_noro(klon, llm, dtphys, paprs, play, zmea, zstd, zsig, zgam, &         CALL drag_noro(paprs, play, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
876              zthe, zpic, zval, itest, t_seri, u_seri, v_seri, zulow, zvlow, &              ktest, t_seri, u_seri, v_seri, zulow, zvlow, zustrdr, zvstrdr, &
877              zustrdr, zvstrdr, d_t_oro, d_u_oro, d_v_oro)              d_t_oro, d_u_oro, d_v_oro)
878    
879         ! ajout des tendances         ! ajout des tendances
880         DO k = 1, llm         DO k = 1, llm
# Line 1210  contains Line 888  contains
888    
889      IF (ok_orolf) THEN      IF (ok_orolf) THEN
890         ! S\'election des points pour lesquels le sch\'ema est actif :         ! S\'election des points pour lesquels le sch\'ema est actif :
        igwd = 0  
891         DO i = 1, klon         DO i = 1, klon
892            itest(i) = 0            ktest(i) = 0
893            IF (zpic(i) - zmea(i) > 100.) THEN            IF (zpic(i) - zmea(i) > 100.) THEN
894               itest(i) = 1               ktest(i) = 1
              igwd = igwd + 1  
895            ENDIF            ENDIF
896         ENDDO         ENDDO
897    
898         CALL lift_noro(klon, llm, dtphys, paprs, play, rlat, zmea, zstd, zpic, &         CALL lift_noro(paprs, play, zmea, zstd, zpic, ktest, t_seri, u_seri, &
899              itest, t_seri, u_seri, v_seri, zulow, zvlow, zustrli, zvstrli, &              v_seri, zulow, zvlow, zustrli, zvstrli, d_t_lif, d_u_lif, d_v_lif)
             d_t_lif, d_u_lif, d_v_lif)  
900    
901         ! Ajout des tendances :         ! Ajout des tendances :
902         DO k = 1, llm         DO k = 1, llm
# Line 1233  contains Line 908  contains
908         ENDDO         ENDDO
909      ENDIF      ENDIF
910    
911      ! Stress n\'ecessaires : toute la physique      CALL aaam_bud(rg, romega, pphis, zustrdr, zustrli, &
912             sum((u_seri - u) / dtphys * zmasse, dim = 2), zvstrdr, &
913      DO i = 1, klon           zvstrli, sum((v_seri - v) / dtphys * zmasse, dim = 2), paprs, u, v, &
914         zustrph(i) = 0.           aam, torsfc)
        zvstrph(i) = 0.  
     ENDDO  
     DO k = 1, llm  
        DO i = 1, klon  
           zustrph(i) = zustrph(i) + (u_seri(i, k) - u(i, k)) / dtphys &  
                * zmasse(i, k)  
           zvstrph(i) = zvstrph(i) + (v_seri(i, k) - v(i, k)) / dtphys &  
                * zmasse(i, k)  
        ENDDO  
     ENDDO  
   
     CALL aaam_bud(rg, romega, rlat, rlon, pphis, zustrdr, zustrli, zustrph, &  
          zvstrdr, zvstrli, zvstrph, paprs, u, v, aam, torsfc)  
   
     IF (if_ebil >= 2) CALL diagetpq(airephy, 'after orography', ip_ebil, 2, &  
          2, dtphys, t_seri, q_seri, ql_seri, u_seri, v_seri, paprs, d_h_vcol, &  
          d_qt, d_ec)  
915    
916      ! Calcul des tendances traceurs      ! Calcul des tendances traceurs
917      call phytrac(lmt_pas, julien, time, firstcal, lafin, dtphys, t, paprs, &      call phytrac(julien, time, firstcal, lafin, t, paprs, play, mfu, mfd, &
918           play, mfu, mfd, pde_u, pen_d, ycoefh, fm_therm, entr_therm, yu1, &           pde_u, pen_d, coefh, cdragh, fm_therm, entr_therm, u(:, 1), v(:, 1), &
919           yv1, ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &           ftsol, pctsrf, frac_impa, frac_nucl, da, phi, mp, upwd, dnwd, &
920           tr_seri, zmasse, ncid_startphy)           tr_seri, zmasse, ncid_startphy)
921    
     IF (offline) call phystokenc(dtphys, t, mfu, mfd, pen_u, pde_u, pen_d, &  
          pde_d, fm_therm, entr_therm, ycoefh, yu1, yv1, ftsol, pctsrf, &  
          frac_impa, frac_nucl, pphis, airephy, dtphys)  
   
922      ! Calculer le transport de l'eau et de l'energie (diagnostique)      ! Calculer le transport de l'eau et de l'energie (diagnostique)
923      CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)      CALL transp(paprs, t_seri, q_seri, u_seri, v_seri, zphi, ve, vq, ue, uq)
924    
# Line 1275  contains Line 929  contains
929    
930      ! Accumuler les variables a stocker dans les fichiers histoire:      ! Accumuler les variables a stocker dans les fichiers histoire:
931    
932      ! conversion Ec -> E thermique      ! conversion Ec en énergie thermique
933      DO k = 1, llm      DO k = 1, llm
934         DO i = 1, klon         DO i = 1, klon
935            ZRCPD = RCPD * (1. + RVTMP2 * q_seri(i, k))            d_t_ec(i, k) = 0.5 / (RCPD * (1. + RVTMP2 * q_seri(i, k))) &
           d_t_ec(i, k) = 0.5 / ZRCPD &  
936                 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)                 * (u(i, k)**2 + v(i, k)**2 - u_seri(i, k)**2 - v_seri(i, k)**2)
937            t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)            t_seri(i, k) = t_seri(i, k) + d_t_ec(i, k)
938            d_t_ec(i, k) = d_t_ec(i, k) / dtphys            d_t_ec(i, k) = d_t_ec(i, k) / dtphys
939         END DO         END DO
940      END DO      END DO
941    
     IF (if_ebil >= 1) THEN  
        tit = 'after physic'  
        CALL diagetpq(airephy, tit, ip_ebil, 1, 1, dtphys, t_seri, q_seri, &  
             ql_seri, u_seri, v_seri, paprs, d_h_vcol, d_qt, d_ec)  
        ! Comme les tendances de la physique sont ajoute dans la dynamique,  
        ! on devrait avoir que la variation d'entalpie par la dynamique  
        ! est egale a la variation de la physique au pas de temps precedent.  
        ! Donc la somme de ces 2 variations devrait etre nulle.  
        call diagphy(airephy, tit, ip_ebil, topsw, toplw, solsw, sollw, sens, &  
             evap, rain_fall, snow_fall, ztsol, d_h_vcol, d_qt, d_ec)  
        d_h_vcol_phy = d_h_vcol  
     END IF  
   
942      ! SORTIES      ! SORTIES
943    
944      ! prw = eau precipitable      ! prw = eau precipitable
945      DO i = 1, klon      DO i = 1, klon
946         prw(i) = 0.         prw(i) = 0.
947         DO k = 1, llm         DO k = 1, llm
948            prw(i) = prw(i) + q_seri(i, k)*zmasse(i, k)            prw(i) = prw(i) + q_seri(i, k) * zmasse(i, k)
949         ENDDO         ENDDO
950      ENDDO      ENDDO
951    
# Line 1343  contains Line 983  contains
983      CALL histwrite_phy("precip", rain_fall + snow_fall)      CALL histwrite_phy("precip", rain_fall + snow_fall)
984      CALL histwrite_phy("plul", rain_lsc + snow_lsc)      CALL histwrite_phy("plul", rain_lsc + snow_lsc)
985      CALL histwrite_phy("pluc", rain_con + snow_con)      CALL histwrite_phy("pluc", rain_con + snow_con)
986      CALL histwrite_phy("tsol", zxtsol)      CALL histwrite_phy("tsol", tsol)
987      CALL histwrite_phy("t2m", zt2m)      CALL histwrite_phy("t2m", zt2m)
988      CALL histwrite_phy("q2m", zq2m)      CALL histwrite_phy("q2m", zq2m)
989      CALL histwrite_phy("u10m", zu10m)      CALL histwrite_phy("u10m", u10m)
990      CALL histwrite_phy("v10m", zv10m)      CALL histwrite_phy("v10m", v10m)
991      CALL histwrite_phy("snow", snow_fall)      CALL histwrite_phy("snow", snow_fall)
992      CALL histwrite_phy("cdrm", cdragm)      CALL histwrite_phy("cdrm", cdragm)
993      CALL histwrite_phy("cdrh", cdragh)      CALL histwrite_phy("cdrh", cdragh)
# Line 1363  contains Line 1003  contains
1003      CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))      CALL histwrite_phy("dtsvdft", d_ts(:, is_ter))
1004      CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))      CALL histwrite_phy("dtsvdfg", d_ts(:, is_lic))
1005      CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))      CALL histwrite_phy("dtsvdfi", d_ts(:, is_sic))
1006        CALL histwrite_phy("zxfqcalving", sum(fqcalving * pctsrf, dim = 2))
1007    
1008      DO nsrf = 1, nbsrf      DO nsrf = 1, nbsrf
1009         CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf)*100.)         CALL histwrite_phy("pourc_"//clnsurf(nsrf), pctsrf(:, nsrf) * 100.)
1010         CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))         CALL histwrite_phy("fract_"//clnsurf(nsrf), pctsrf(:, nsrf))
1011         CALL histwrite_phy("sens_"//clnsurf(nsrf), fluxt(:, 1, nsrf))         CALL histwrite_phy("sens_"//clnsurf(nsrf), flux_t(:, nsrf))
1012         CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))         CALL histwrite_phy("lat_"//clnsurf(nsrf), fluxlat(:, nsrf))
1013         CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))         CALL histwrite_phy("tsol_"//clnsurf(nsrf), ftsol(:, nsrf))
1014         CALL histwrite_phy("taux_"//clnsurf(nsrf), fluxu(:, 1, nsrf))         CALL histwrite_phy("taux_"//clnsurf(nsrf), flux_u(:, nsrf))
1015         CALL histwrite_phy("tauy_"//clnsurf(nsrf), fluxv(:, 1, nsrf))         CALL histwrite_phy("tauy_"//clnsurf(nsrf), flux_v(:, nsrf))
1016         CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))         CALL histwrite_phy("rugs_"//clnsurf(nsrf), frugs(:, nsrf))
1017         CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))         CALL histwrite_phy("albe_"//clnsurf(nsrf), falbe(:, nsrf))
1018           CALL histwrite_phy("u10m_"//clnsurf(nsrf), u10m_srf(:, nsrf))
1019           CALL histwrite_phy("v10m_"//clnsurf(nsrf), v10m_srf(:, nsrf))
1020      END DO      END DO
1021    
1022      CALL histwrite_phy("albs", albsol)      CALL histwrite_phy("albs", albsol)
1023        CALL histwrite_phy("tro3", wo * dobson_u * 1e3 / zmasse / rmo3 * md)
1024      CALL histwrite_phy("rugs", zxrugs)      CALL histwrite_phy("rugs", zxrugs)
1025      CALL histwrite_phy("s_pblh", s_pblh)      CALL histwrite_phy("s_pblh", s_pblh)
1026      CALL histwrite_phy("s_pblt", s_pblt)      CALL histwrite_phy("s_pblt", s_pblt)
# Line 1385  contains Line 1029  contains
1029      CALL histwrite_phy("s_oliqCL", s_oliqCL)      CALL histwrite_phy("s_oliqCL", s_oliqCL)
1030      CALL histwrite_phy("s_cteiCL", s_cteiCL)      CALL histwrite_phy("s_cteiCL", s_cteiCL)
1031      CALL histwrite_phy("s_therm", s_therm)      CALL histwrite_phy("s_therm", s_therm)
1032      CALL histwrite_phy("s_trmb1", s_trmb1)  
1033      CALL histwrite_phy("s_trmb2", s_trmb2)      if (conv_emanuel) then
1034      CALL histwrite_phy("s_trmb3", s_trmb3)         CALL histwrite_phy("ptop", ema_pct)
1035      if (conv_emanuel) CALL histwrite_phy("ptop", ema_pct)         CALL histwrite_phy("dnwd0", - mp)
1036        end if
1037    
1038      CALL histwrite_phy("temp", t_seri)      CALL histwrite_phy("temp", t_seri)
1039      CALL histwrite_phy("vitu", u_seri)      CALL histwrite_phy("vitu", u_seri)
1040      CALL histwrite_phy("vitv", v_seri)      CALL histwrite_phy("vitv", v_seri)
# Line 1397  contains Line 1043  contains
1043      CALL histwrite_phy("dtvdf", d_t_vdf)      CALL histwrite_phy("dtvdf", d_t_vdf)
1044      CALL histwrite_phy("dqvdf", d_q_vdf)      CALL histwrite_phy("dqvdf", d_q_vdf)
1045      CALL histwrite_phy("rhum", zx_rh)      CALL histwrite_phy("rhum", zx_rh)
1046        CALL histwrite_phy("d_t_ec", d_t_ec)
1047        CALL histwrite_phy("dtsw0", heat0 / 86400.)
1048        CALL histwrite_phy("dtlw0", - cool0 / 86400.)
1049        CALL histwrite_phy("msnow", sum(fsnow * pctsrf, dim = 2))
1050        call histwrite_phy("qsurf", sum(fqsurf * pctsrf, dim = 2))
1051    
1052      if (ok_instan) call histsync(nid_ins)      if (ok_instan) call histsync(nid_ins)
1053    
1054      IF (lafin) then      IF (lafin) then
1055         call NF95_CLOSE(ncid_startphy)         call NF95_CLOSE(ncid_startphy)
1056         CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, &         CALL phyredem(pctsrf, ftsol, ftsoil, fqsurf, qsol, fsnow, falbe, &
1057              fsnow, falbe, fevap, rain_fall, snow_fall, solsw, sollw, dlw, &              rain_fall, snow_fall, solsw, sollw, dlw, radsol, frugs, agesno, &
1058              radsol, frugs, agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &              zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, q_ancien, &
1059              t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0, sig1, &              rnebcon, ratqs, clwcon, run_off_lic_0, sig1, w01)
             w01)  
1060      end IF      end IF
1061    
1062      firstcal = .FALSE.      firstcal = .FALSE.
Legend:
Removed from v.191  
changed lines
  Added in v.305
  ViewVC Help
Powered by ViewVC 1.1.21