source: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90 @ 5609

Last change on this file since 5609 was 5609, checked in by aborella, 3 months ago
  • changed treatment of prognostic variables for prognostic clouds
  • adapted sedimentation and autoconversion for prognostic cirrus clouds
  • cloud mixing, ice sedimentation and ISSR diagnosis are now consistent with the water vapor PDF
  • simplified assumptions for ice crystals deposition / sublimation
  • first version of the coupling between prognostic cirrus clouds and deep convection
  • added persistent contrail cirrus clouds in radiative diagnostics
File size: 61.0 KB
Line 
1MODULE lmdz_lscp
2
3IMPLICIT NONE
4
5CONTAINS
6
7!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
9     paprs, pplay, omega, temp, qt, ql_seri, qi_seri,   &
10     ratqs, sigma_qtherm, ptconv, cfcon_old, qvcon_old, &
11     qccon_old, cfcon, qvcon, qccon,                    &
12     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
13     pfraclr, pfracld,                                  &
14     cldfraliq, sigma2_icefracturb,mean_icefracturb,    &
15     radocond, radicefrac, rain, snow,                  &
16     frac_impa, frac_nucl, beta,                        &
17     prfl, psfl, rhcl, qta, fraca,                      &
18     tv, pspsk, tla, thl, iflag_cld_th,                 &
19     iflag_ice_thermo, distcltop, temp_cltop,           &
20     tke, tke_dissip,                                   &
21     cell_area, stratomask,                             &
22     cf_seri, qvc_seri, u_seri, v_seri,                 &
23     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
24     dcf_sub, dcf_con, dcf_mix,                         &
25     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
26     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
27     cfa_seri, pcf_seri, qva_seri, qia_seri,flight_dist,&
28     flight_h2o, qice_perscont, Tcritcont,              &
29     qcritcont, potcontfraP, potcontfraNP,              &
30     cloudth_sth,                                       &
31     cloudth_senv, cloudth_sigmath, cloudth_sigmaenv,   &
32     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
33     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
34     dqsmelt, dqsfreez, dqised, dcfsed, dqvcsed)
35
36!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
37! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
38!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
39!--------------------------------------------------------------------------------
40! Date: 01/2021
41!--------------------------------------------------------------------------------
42! Aim: Large Scale Clouds and Precipitation (LSCP)
43!
44! This code is a new version of the fisrtilp.F90 routine, which is itself a
45! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
46! and 'ilp' (il pleut, L. Li)
47!
48! Compared to the original fisrtilp code, lscp
49! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
50! -> consider always precipitation thermalisation (fl_cor_ebil>0)
51! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
52! -> option "oldbug" by JYG has been removed
53! -> iflag_t_glace >0 always
54! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
55! -> rectangular distribution from L. Li no longer available
56! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
57!--------------------------------------------------------------------------------
58! References:
59!
60! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
61! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
62! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
63! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
64! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
65! - Touzze-Peifert Ludo, PhD thesis, p117-124
66! -------------------------------------------------------------------------------
67! Code structure:
68!
69! P0>     Thermalization of the precipitation coming from the overlying layer
70! P1>     Evaporation of the precipitation (falling from the k+1 level)
71! P2>     Cloud formation (at the k level)
72! P2.A.1> With the PDFs, calculation of cloud properties using the inital
73!         values of T and Q
74! P2.A.2> Coupling between condensed water and temperature
75! P2.A.3> Calculation of final quantities associated with cloud formation
76! P2.B>   Release of Latent heat after cloud formation
77! P3>     Autoconversion to precipitation (k-level)
78! P4>     Wet scavenging
79!------------------------------------------------------------------------------
80! Some preliminary comments (JBM) :
81!
82! The cloud water that the radiation scheme sees is not the same that the cloud
83! water used in the physics and the dynamics
84!
85! During the autoconversion to precipitation (P3 step), radocond (cloud water used
86! by the radiation scheme) is calculated as an average of the water that remains
87! in the cloud during the precipitation and not the water remaining at the end
88! of the time step. The latter is used in the rest of the physics and advected
89! by the dynamics.
90!
91! In summary:
92!
93! Radiation:
94! xflwc(newmicro)+xfiwc(newmicro) =
95!   radocond=lwcon(nc)+iwcon(nc)
96!
97! Notetheless, be aware of:
98!
99! radocond .NE. ocond(nc)
100! i.e.:
101! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
102! but oliq+(ocond-oliq) .EQ. ocond
103! (which is not trivial)
104!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
105!
106
107! USE de modules contenant des fonctions.
108USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
109USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
110USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
111USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
112USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
113USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
114
115! Use du module lmdz_lscp_ini contenant les constantes
116USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
117USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
118USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
119USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
120USE lmdz_lscp_ini, ONLY : t_glace_min, temp_nowater, min_frac_th_cld
121USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
122USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
123USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
124USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato
125USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_ice_sedim
126
127! Temporary call for Lamquin et al (2012) diagnostics
128USE phys_local_var_mod, ONLY : issrfra100to150, issrfra150to200, issrfra200to250
129USE phys_local_var_mod, ONLY : issrfra250to300, issrfra300to400, issrfra400to500
130USE phys_local_var_mod, ONLY : dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub
131USE phys_local_var_mod, ONLY : dcfa_cir, dqta_cir, dcfa_mix, dqia_mix, dqta_mix
132
133IMPLICIT NONE
134
135!===============================================================================
136! VARIABLES DECLARATION
137!===============================================================================
138
139! INPUT VARIABLES:
140!-----------------
141
142  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
143  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
144  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
145
146  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
147  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
148  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
149  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
150  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
151  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ql_seri         ! liquid specific content [kg/kg]
152  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qi_seri         ! ice specific content [kg/kg]
153  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
154  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
155                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
156  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
157  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon_old       ! cloud fraction from deep convection from previous timestep [-]
158  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon_old       ! in-cloud vapor specific humidity from deep convection from previous timestep [kg/kg]
159  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon_old       ! in-cloud condensed specific humidity from deep convection from previous timestep [kg/kg]
160  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon           ! cloud fraction from deep convection [-]
161  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon           ! in-cloud vapor specific humidity from deep convection [kg/kg]
162  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon           ! in-cloud condensed specific humidity from deep convection [kg/kg]
163
164  !Inputs associated with thermal plumes
165
166  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
167  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
168  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
169  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
170  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
171  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
172  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
173
174  ! INPUT/OUTPUT variables
175  !------------------------
176 
177  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
178  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs,sigma_qtherm        ! function of pressure that sets the large-scale
179
180
181  ! INPUT/OUTPUT condensation and ice supersaturation
182  !--------------------------------------------------
183  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
184  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qvc_seri         ! cloudy water vapor [kg/kg]
185  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
186  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
187  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
188  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: stratomask       ! fraction of stratosphere (0 or 1)
189
190  ! INPUT/OUTPUT aviation
191  !--------------------------------------------------
192  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfa_seri         ! linear contrails fraction [-]
193  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: pcf_seri         ! contrails induced cirrus fraction [-]
194  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qva_seri         ! linear contrails total specific humidity [kg/kg]
195  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qia_seri         ! linear contrails total specific humidity [kg/kg]
196  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
197  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
198 
199  ! OUTPUT variables
200  !-----------------
201
202  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
203  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
204  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
206  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
207  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
208  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
209  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
210  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
211  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
212  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
213  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
214  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
216  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
217  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
218  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
219  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
220  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
223
224  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
225 
226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
228 
229  ! for condensation and ice supersaturation
230
231  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
238  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
239  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
240  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
241  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
242  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
246  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
247  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
248  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
250
251  ! for contrails and aviation
252
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_perscont  !--condensed water in contrail induced cirrus used in the radiation scheme [kg/kg]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
256  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: potcontfraP    !--potential persistent contrail fraction [-]
257  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: potcontfraNP   !--potential non-persistent contrail fraction [-]
258
259
260  ! for POPRECIP
261
262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
263  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
264  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
265  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
266  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
267  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
268  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
269  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
270  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
271  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
272  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
273  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
274  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
275  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqised         !--ice water content tendency due to sedmentation of ice crystals [kg/kg/s]
276  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcfsed         !--cloud fraction tendency due to sedimentation of ice crystals [kg/kg/s]
277  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvcsed        !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
278
279  ! for thermals
280
281  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
282  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
283  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
284  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
285
286
287  ! LOCAL VARIABLES:
288  !----------------
289  REAL, DIMENSION(klon) :: qliq_in, qice_in, qvc_in, cldfra_in
290  REAL, DIMENSION(klon,klev) :: ctot
291  REAL, DIMENSION(klon,klev) :: ctot_vol
292  REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi
293  REAL :: zdelta
294  REAL, DIMENSION(klon) :: zdqsdT_raw
295  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
296  REAL, DIMENSION(klon) :: Tbef,qlbef,DT                          ! temperature, humidity and temp. variation during lognormal iteration
297  REAL :: num,denom
298  REAL :: cste
299  REAL, DIMENSION(klon) :: zfice_th
300  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
301  REAL, DIMENSION(klon) :: zrfl, zifl
302  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn, zqnl
303  REAL, DIMENSION(klon) :: zoliql, zoliqi
304  REAL, DIMENSION(klon) :: zt
305  REAL, DIMENSION(klon) :: zfice,zneb, znebl
306  REAL, DIMENSION(klon) :: dzfice
307  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
308  REAL, DIMENSION(klon) :: qtot, qzero
309  ! Variables precipitation energy conservation
310  REAL, DIMENSION(klon) :: zmqc
311  REAL :: zalpha_tr
312  REAL :: zfrac_lessi
313  REAL, DIMENSION(klon) :: zprec_cond
314  REAL, DIMENSION(klon) :: zlh_solid
315  REAL, DIMENSION(klon) :: ztupnew
316  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
317  REAL, DIMENSION(klon) :: cldfra_above, icesed_flux ! for sedimentation of ice crystals
318  REAL, DIMENSION(klon) :: zrflclr, zrflcld
319  REAL, DIMENSION(klon) :: ziflclr, ziflcld
320  REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
321  REAL, DIMENSION(klon) :: tot_zneb
322  REAL :: qlmpc, qimpc, rnebmpc
323  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
324  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
325
326  ! for quantity of condensates seen by radiation
327  REAL, DIMENSION(klon) :: zradocond, zradoice
328  REAL, DIMENSION(klon) :: zrho_up, zvelo_up
329 
330  ! for condensation and ice supersaturation
331  REAL, DIMENSION(klon) :: qvc, qvcl, shear
332  REAL :: delta_z
333  ! for contrails
334  REAL, DIMENSION(klon) :: contfra, perscontfra, qcont
335  LOGICAL, DIMENSION(klon) :: pt_pron_clds
336  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrail)
337  ! Constants used for calculating ratios that are advected (using a parent-child
338  ! formalism). This is not done in the dynamical core because at this moment,
339  ! only isotopes can use this parent-child formalism. Note that the two constants
340  ! are the same as the one use in the dynamical core, being also defined in
341  ! dyn3d_common/infotrac.F90
342  REAL :: min_qParent, min_ratio
343  !--for Lamquin et al 2012 diagnostics
344  REAL, DIMENSION(klon) :: issrfra100to150UP, issrfra150to200UP, issrfra200to250UP
345  REAL, DIMENSION(klon) :: issrfra250to300UP, issrfra300to400UP, issrfra400to500UP
346
347  INTEGER i, k, kk, iter
348  INTEGER, DIMENSION(klon) :: n_i
349  INTEGER ncoreczq
350  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
351  LOGICAL iftop
352
353  LOGICAL, DIMENSION(klon) :: lognormale
354  LOGICAL, DIMENSION(klon) :: keepgoing
355
356  CHARACTER (len = 20) :: modname = 'lscp'
357  CHARACTER (len = 80) :: abort_message
358 
359
360!===============================================================================
361! INITIALISATION
362!===============================================================================
363
364! Few initial checks
365
366
367IF (iflag_fisrtilp_qsat .LT. 0) THEN
368    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
369    CALL abort_physic(modname,abort_message,1)
370ENDIF
371
372! Few initialisations
373
374ctot_vol(1:klon,1:klev)=0.0
375rneblsvol(1:klon,1:klev)=0.0
376znebprecip(:)=0.0
377znebprecipclr(:)=0.0
378znebprecipcld(:)=0.0
379mpc_bl_points(:,:)=0
380
381IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
382
383! AA for 'safety' reasons
384zalpha_tr   = 0.
385zfrac_lessi = 0.
386beta(:,:)= 0.
387
388! Initialisation of variables:
389
390prfl(:,:) = 0.0
391psfl(:,:) = 0.0
392d_t(:,:) = 0.0
393d_q(:,:) = 0.0
394d_ql(:,:) = 0.0
395d_qi(:,:) = 0.0
396rneb(:,:) = 0.0
397pfraclr(:,:)=0.0
398pfracld(:,:)=0.0
399cldfraliq(:,:)=0.
400sigma2_icefracturb(:,:)=0.
401mean_icefracturb(:,:)=0.
402radocond(:,:) = 0.0
403radicefrac(:,:) = 0.0
404frac_nucl(:,:) = 1.0
405frac_impa(:,:) = 1.0
406rain(:) = 0.0
407snow(:) = 0.0
408zfice(:)=1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
409dzfice(:)=0.0
410zfice_turb(:)=0.0
411dzfice_turb(:)=0.0
412zrfl(:) = 0.0
413zifl(:) = 0.0
414zneb(:) = seuil_neb
415zrflclr(:) = 0.0
416ziflclr(:) = 0.0
417zrflcld(:) = 0.0
418ziflcld(:) = 0.0
419tot_zneb(:) = 0.0
420qzero(:) = 0.0
421zdistcltop(:)=0.0
422ztemp_cltop(:) = 0.0
423ztupnew(:)=0.0
424
425distcltop(:,:)=0.
426temp_cltop(:,:)=0.
427
428!--Ice supersaturation
429gamma_cond(:,:) = 1.
430qissr(:,:)      = 0.
431issrfra(:,:)    = 0.
432dcf_sub(:,:)    = 0.
433dcf_con(:,:)    = 0.
434dcf_mix(:,:)    = 0.
435dcfsed(:,:)     = 0.
436dqi_adj(:,:)    = 0.
437dqi_sub(:,:)    = 0.
438dqi_con(:,:)    = 0.
439dqi_mix(:,:)    = 0.
440dqised(:,:)     = 0.
441dqvc_adj(:,:)   = 0.
442dqvc_sub(:,:)   = 0.
443dqvc_con(:,:)   = 0.
444dqvc_mix(:,:)   = 0.
445dqvcsed(:,:)    = 0.
446qvc(:)          = 0.
447shear(:)        = 0.
448min_qParent     = 1.e-30
449min_ratio       = 1.e-16
450
451!--for Lamquin et al (2012) diagnostics
452issrfra100to150(:)   = 0.
453issrfra100to150UP(:) = 0.
454issrfra150to200(:)   = 0.
455issrfra150to200UP(:) = 0.
456issrfra200to250(:)   = 0.
457issrfra200to250UP(:) = 0.
458issrfra250to300(:)   = 0.
459issrfra250to300UP(:) = 0.
460issrfra300to400(:)   = 0.
461issrfra300to400UP(:) = 0.
462issrfra400to500(:)   = 0.
463issrfra400to500UP(:) = 0.
464
465!-- poprecip
466qraindiag(:,:)= 0.
467qsnowdiag(:,:)= 0.
468dqreva(:,:)   = 0.
469dqrauto(:,:)  = 0.
470dqrmelt(:,:)  = 0.
471dqrfreez(:,:) = 0.
472dqrcol(:,:)   = 0.
473dqssub(:,:)   = 0.
474dqsauto(:,:)  = 0.
475dqsrim(:,:)   = 0.
476dqsagg(:,:)   = 0.
477dqsfreez(:,:) = 0.
478dqsmelt(:,:)  = 0.
479zqupnew(:)    = 0.
480zqvapclr(:)   = 0.
481cldfra_above(:)= 0.
482icesed_flux(:)= 0.
483
484
485
486!c_iso: variable initialisation for iso
487
488!===============================================================================
489!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
490!===============================================================================
491
492ncoreczq=0
493
494DO k = klev, 1, -1
495
496    IF (k.LE.klev-1) THEN
497       iftop=.false.
498    ELSE
499       iftop=.true.
500    ENDIF
501
502    ! Initialisation temperature and specific humidity
503    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
504    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
505    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
506    ! d_t = temperature tendency due to lscp
507    ! The temperature of the overlying layer is updated here because needed for thermalization
508    DO i = 1, klon
509        zt(i)=temp(i,k)
510        zq(i)=qt(i,k)
511        qliq_in(i) = ql_seri(i,k)
512        qice_in(i) = qi_seri(i,k)
513        IF (.not. iftop) THEN
514           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
515           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
516           !--zqs(i) is the saturation specific humidity in the layer above
517           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
518        ENDIF
519        !c_iso init of iso
520    ENDDO
521    IF ( ok_ice_supersat ) THEN
522      cldfra_in(:) = cf_seri(:,k)
523      qvc_in(:) = qvc_seri(:,k)
524    ENDIF
525
526    ! --------------------------------------------------------------------
527    ! P1> Precipitation processes, before cloud formation:
528    !     Thermalization of precipitation falling from the overlying layer AND
529    !     Precipitation evaporation/sublimation/melting
530    !---------------------------------------------------------------------
531   
532    !================================================================
533    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
534    IF ( ok_poprecip ) THEN
535
536      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
537                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
538                        zqvapclr, zqupnew, icesed_flux, &
539                        cldfra_in, qvc_in, qliq_in, qice_in, &
540                        zrfl, zrflclr, zrflcld, &
541                        zifl, ziflclr, ziflcld, &
542                        dqreva(:,k), dqssub(:,k) &
543                        )
544
545    !================================================================
546    ELSE
547
548      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
549                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
550                        zrfl, zrflclr, zrflcld, &
551                        zifl, ziflclr, ziflcld, &
552                        dqreva(:,k), dqssub(:,k) &
553                        )
554
555    ENDIF ! (ok_poprecip)
556   
557    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
558    !-------------------------------------------------------
559
560     qtot(:)=zq(:)+zmqc(:)
561     CALL calc_qsat_ecmwf(klon,zt,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
562     DO i = 1, klon
563            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
564            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
565            IF (zq(i) .LT. 1.e-15) THEN
566                ncoreczq=ncoreczq+1
567                zq(i)=1.e-15
568            ENDIF
569            ! c_iso: do something similar for isotopes
570
571     ENDDO
572   
573    ! --------------------------------------------------------------------
574    ! P2> Cloud formation
575    !---------------------------------------------------------------------
576    !
577    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
578    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
579    ! is prescribed and depends on large scale variables and boundary layer
580    ! properties)
581    ! The decrease in condensed part due tu latent heating is taken into
582    ! account
583    ! -------------------------------------------------------------------
584
585        ! P2.1> With the PDFs (log-normal, bigaussian)
586        ! cloud properties calculation with the initial values of t and q
587        ! ----------------------------------------------------------------
588
589        ! initialise gammasat and qincloud_mpc
590        gammasat(:)=1.
591        qincloud_mpc(:)=0.
592
593        IF (iflag_cld_th.GE.5) THEN
594               ! Cloud cover and content in meshes affected by shallow convection,
595               ! are retrieved from a bi-gaussian distribution of the saturation deficit
596               ! following Jam et al. 2013
597
598               IF (iflag_cloudth_vert.LE.2) THEN
599                  ! Old version of Arnaud Jam
600
601                    CALL cloudth(klon,klev,k,tv,                  &
602                         zq,qta,fraca,                            &
603                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
604                         ratqs,zqs,temp,                              &
605                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
606
607
608                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
609                   ! Default version of Arnaud Jam
610                       
611                    CALL cloudth_v3(klon,klev,k,tv,                        &
612                         zq,qta,fraca,                                     &
613                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
614                         ratqs,sigma_qtherm,zqs,temp, &
615                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
616
617
618                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
619                   ! Jean Jouhaud's version, with specific separation between surface and volume
620                   ! cloud fraction Decembre 2018
621
622                    CALL cloudth_v6(klon,klev,k,tv,                        &
623                         zq,qta,fraca,                                     &
624                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
625                         ratqs,zqs,temp, &
626                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
627
628                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
629                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
630                   ! for boundary-layer mixed phase clouds
631                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
632                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
633                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
634                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
635
636                ENDIF
637
638
639                DO i=1,klon
640                    rneb(i,k)=ctot(i,k)
641                    rneblsvol(i,k)=ctot_vol(i,k)
642                    zqn(i)=qcloud(i)
643                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
644                    qvc(i) = rneb(i,k) * zqs(i)
645                ENDDO
646
647        ENDIF
648
649        IF (iflag_cld_th .LE. 4) THEN
650           
651                ! lognormal
652            lognormale(:) = .TRUE.
653
654        ELSEIF (iflag_cld_th .GE. 6) THEN
655           
656                ! lognormal distribution when no thermals
657            lognormale(:) = fraca(:,k) < min_frac_th_cld
658
659        ELSE
660                ! When iflag_cld_th=5, we always assume
661                ! bi-gaussian distribution
662            lognormale(:) = .FALSE.
663       
664        ENDIF
665
666
667        IF ( ok_ice_supersat ) THEN
668
669          !--Initialisation
670          IF ( ok_plane_contrail ) THEN
671            contfra(:)     = 0.
672            qcont(:)       = 0.
673            perscontfra(:) = 0.
674          ENDIF
675
676          DO i = 1, klon
677
678            pt_pron_clds(i) = ( ( ( ( zt(i) .LE. temp_nowater ) .OR. ok_weibull_warm_clouds ) &
679                .AND. ( .NOT. ok_no_issr_strato .OR. ( stratomask(i,k) .EQ. 0. ) ) ) &
680                .AND. ( cfcon(i,k) .LT. ( 1. - eps ) ) )
681
682            IF ( pt_pron_clds(i) ) THEN
683
684              !--If deep convection is activated, the condensation scheme activates
685              !--only in the environment. NB. the clear sky fraction will the be
686              !--maximised by 1. - cfcon(i,k)
687              IF ( ptconv(i,k) ) zq(i) = zq(i) - ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
688
689              IF ( ( cfcon(i,k) * qccon(i,k) ) .LT. ( cfcon_old(i,k) * qccon_old(i,k) ) ) THEN
690                !--If deep convection is weakening, we add the clouds that are not anymore
691                !--'in' deep convection to the advected clouds
692                cldfra_in(i) = cldfra_in(i) + MAX(0., cfcon_old(i,k) - cfcon(i,k))
693                qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * MAX(0., cfcon_old(i,k) - cfcon(i,k))
694                qice_in(i) = qice_in(i) + ( qccon_old(i,k) * cfcon_old(i,k) &
695                                          - qccon(i,k) * cfcon(i,k) )
696              ELSEIF ( cfcon(i,k) .GT. cfcon_old(i,k) ) THEN
697                !--Else if deep convection is strengthening, it consumes the existing cloud
698                !--fraction (which does not at this moment represent deep convection)
699                !--NB. if deep convection is strengthening while the fraction decreases,
700                !--clear sky water vapor will be transfered in priority
701                cldfra_in(i) = cldfra_in(i) * ( 1. &
702                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
703                qvc_in(i)    = qvc_in(i)    * ( 1. &
704                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
705                qice_in(i)   = qice_in(i)   * ( 1. &
706                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
707              ENDIF
708
709              !--Barriers
710              cldfra_in(i) = MAX(0., MIN(1. - cfcon(i,k), cldfra_in(i)))
711              qvc_in(i)    = MAX(0., MIN(zq(i), qvc_in(i)))
712              qice_in(i)   = MAX(0., MIN(zq(i) - qvc_in(i), qice_in(i)))
713
714              !--Calculate the shear value (input for condensation and ice supersat)
715              !--Cell thickness [m]
716              delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(i) * RD
717              IF ( iftop ) THEN
718                ! top
719                shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
720                               + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
721              ELSEIF ( k .EQ. 1 ) THEN
722                ! surface
723                shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
724                               + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
725              ELSE
726                ! other layers
727                shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
728                                   - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
729                               + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
730                                   - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
731              ENDIF
732            ENDIF
733          ENDDO
734        ENDIF
735
736        DT(:) = 0.
737        n_i(:)=0
738        Tbef(:)=zt(:)
739        qlbef(:)=0.
740
741        ! Treatment of non-boundary layer clouds (lognormale)
742        ! We iterate here to take into account the change in qsat(T) and ice fraction
743        ! during the condensation process
744        ! the increment in temperature is calculated using the first principle of
745        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
746        ! and a clear sky fraction)
747        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
748
749        DO iter=1,iflag_fisrtilp_qsat+1
750
751                keepgoing(:) = .FALSE.
752
753                DO i=1,klon
754
755                    ! keepgoing = .true. while convergence is not satisfied
756
757                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
758
759                        ! if not convergence:
760                        ! we calculate a new iteration
761                        keepgoing(i) = .TRUE.
762
763                        ! P2.2.1> cloud fraction and condensed water mass calculation
764                        ! Calculated variables:
765                        ! rneb : cloud fraction
766                        ! zqn : total water within the cloud
767                        ! zcond: mean condensed water within the mesh
768                        ! rhcl: clear-sky relative humidity
769                        !---------------------------------------------------------------
770
771                        ! new temperature that only serves in the iteration process:
772                        Tbef(i)=Tbef(i)+DT(i)
773
774                        ! Rneb, qzn and zcond for lognormal PDFs
775                        qtot(i)=zq(i)+zmqc(i)
776
777                      ENDIF
778
779                  ENDDO
780
781                  ! Calculation of saturation specific humidity and ice fraction
782                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
783                 
784                  IF (iflag_icefrac .GE. 3) THEN
785                  ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
786                  ! and cloud condensed water content. idea following Dietlitcher et al. 2018, GMD
787                  ! we make this option works only for the physically-based and tke-depdenent param (iflag_icefrac>=1)
788                      DO i=1,klon
789                           CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,1,.false.,zqsl,zdqsl)
790                           CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,2,.false.,zqsi,zdqsi)
791                           zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i)
792                           zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i)
793                      ENDDO
794                  ENDIF
795
796                  CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
797                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
798                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
799
800                  ! Cloud condensation based on subgrid pdf
801                  !--AB Activates a condensation scheme that allows for
802                  !--ice supersaturation and contrails evolution from aviation
803                  IF (ok_ice_supersat) THEN
804
805                    IF ( iftop ) THEN
806                      cldfra_above(:) = 0.
807                    ELSE
808                      cldfra_above(:) = rneb(:,k+1)
809                    ENDIF
810
811                    !---------------------------------------------
812                    !--   CONDENSATION AND ICE SUPERSATURATION  --
813                    !---------------------------------------------
814
815                    CALL condensation_ice_supersat( &
816                        klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), &
817                        cfcon(:,k), cldfra_in, qvc_in, qliq_in, qice_in, &
818                        shear, tke_dissip(:,k), cell_area, Tbef, zq, zqs, &
819                        gammasat, ratqs(:,k), keepgoing, pt_pron_clds, &
820                        cldfra_above, icesed_flux,&
821                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
822                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), dcfsed(:,k), &
823                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), dqised(:,k), &
824                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), dqvcsed(:,k), &
825                        cfa_seri(:,k), pcf_seri(:,k), qva_seri(:,k), qia_seri(:,k), &
826                        flight_dist(:,k), flight_h2o(:,k), contfra, perscontfra, qcont, &
827                        Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
828                        dcfa_ini(:,k), dqia_ini(:,k), dqta_ini(:,k), &
829                        dcfa_sub(:,k), dqia_sub(:,k), dqta_sub(:,k), &
830                        dcfa_cir(:,k), dqta_cir(:,k), &
831                        dcfa_mix(:,k), dqia_mix(:,k), dqta_mix(:,k))
832
833
834                  ELSE
835                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
836
837                   CALL condensation_lognormal( &
838                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
839                       keepgoing, rneb(:,k), zqn, qvc)
840
841
842                  ENDIF ! .NOT. ok_ice_supersat
843
844                  ! cloud phase determination
845                  IF (iflag_icefrac .GE. 2) THEN
846                     ! phase partitioning depending on temperature. activates here in the iteration process if iflag_icefrac > 2
847                     CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_in, ziflcld, zqn, &
848                     rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
849                  ELSE
850                     ! phase partitioning depending on temperature and eventually distance to cloud top
851                     IF (iflag_t_glace.GE.4) THEN
852                       ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
853                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
854                     ENDIF
855                     CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
856                  ENDIF
857
858
859                  DO i=1,klon
860                      IF (keepgoing(i)) THEN
861
862                        ! If vertical heterogeneity, change fraction by volume as well
863                        IF (iflag_cloudth_vert.GE.3) THEN
864                            ctot_vol(i,k)=rneb(i,k)
865                            rneblsvol(i,k)=ctot_vol(i,k)
866                        ENDIF
867
868
869                       ! P2.2.2> Approximative calculation of temperature variation DT
870                       !           due to condensation.
871                       ! Calculated variables:
872                       ! dT : temperature change due to condensation
873                       !---------------------------------------------------------------
874
875               
876                        IF (zfice(i).LT.1) THEN
877                            cste=RLVTT
878                        ELSE
879                            cste=RLSTT
880                        ENDIF
881                       
882                        IF ( ok_unadjusted_clouds ) THEN
883                          !--AB We relax the saturation adjustment assumption
884                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
885                          IF ( rneb(i,k) .GT. eps ) THEN
886                            qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
887                          ELSE
888                            qlbef(i) = 0.
889                          ENDIF
890                        ELSE
891                          qlbef(i)=max(0.,zqn(i)-zqs(i))
892                        ENDIF
893
894                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
895                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
896                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
897                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
898                              *qlbef(i)*dzfice(i)
899                        ! here we update a provisory temperature variable that only serves in the iteration
900                        ! process
901                        DT(i)=num/denom
902                        n_i(i)=n_i(i)+1
903
904                    ENDIF ! end keepgoing
905
906                ENDDO     ! end loop on i
907
908        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
909
910        ! P2.2> Final quantities calculation
911        ! Calculated variables:
912        !   rneb : cloud fraction
913        !   zcond: mean condensed water in the mesh
914        !   zqn  : mean water vapor in the mesh
915        !   zfice: ice fraction in clouds
916        !   zt   : temperature
917        !   rhcl : clear-sky relative humidity
918        ! ----------------------------------------------------------------
919
920
921        ! Cloud phase final determination
922        !--------------------------------
923        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
924        IF (iflag_t_glace.GE.4) THEN
925           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
926           distcltop(:,k)=zdistcltop(:)
927           temp_cltop(:,k)=ztemp_cltop(:)
928        ENDIF
929        ! Partition function depending on temperature for all clouds (shallow convective and stratiform)
930        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
931
932        ! Partition function depending on tke for non shallow-convective clouds, erase previous estimation
933        IF (iflag_icefrac .GE. 1) THEN
934           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_in, ziflcld, zqn, &
935           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
936        ENDIF
937
938        ! Water vapor update, subsequent latent heat exchange for each cloud type
939        !------------------------------------------------------------------------
940        DO i=1, klon
941            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
942            ! iflag_cloudth_vert=7 and specific param is activated
943            IF (mpc_bl_points(i,k) .GT. 0) THEN
944                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
945                ! following line is very strange and probably wrong
946                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
947                ! water vapor update and partition function if thermals
948                zq(i) = zq(i) - zcond(i)       
949                zfice(i)=zfice_th(i)
950            ELSE
951                ! Checks on rneb, rhcl and zqn
952                IF (rneb(i,k) .LE. 0.0) THEN
953                    zqn(i) = 0.0
954                    rneb(i,k) = 0.0
955                    zcond(i) = 0.0
956                    rhcl(i,k)=zq(i)/zqs(i)
957                ELSE IF (rneb(i,k) .GE. 1.0) THEN
958                    zqn(i) = zq(i)
959                    rneb(i,k) = 1.0
960                    IF ( ok_unadjusted_clouds ) THEN
961                      !--AB We relax the saturation adjustment assumption
962                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
963                      zcond(i) = MAX(0., zqn(i) - qvc(i))
964                    ELSE
965                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
966                    ENDIF
967                    rhcl(i,k)=1.0
968                ELSE
969                    IF ( ok_unadjusted_clouds ) THEN
970                      !--AB We relax the saturation adjustment assumption
971                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
972                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
973                    ELSE
974                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
975                    ENDIF
976                    ! following line is very strange and probably wrong:
977                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
978                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
979                    IF (iflag_icefrac .GE. 1) THEN
980                        IF (lognormale(i)) THEN 
981                           zcond(i)  = zqliq(i) + zqice(i)
982                           zfice(i)  = zfice_turb(i)
983                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
984                        ENDIF
985                    ENDIF
986                ENDIF
987
988                ! water vapor update
989                zq(i) = zq(i) - zcond(i)
990
991            ENDIF
992           
993           
994            ! temperature update due to phase change
995            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
996            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
997                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
998        ENDDO
999
1000        ! If vertical heterogeneity, change volume fraction
1001        IF (iflag_cloudth_vert .GE. 3) THEN
1002          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1003          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1004        ENDIF
1005
1006
1007    ! ----------------------------------------------------------------
1008    ! P3> Precipitation processes, after cloud formation
1009    !     - precipitation formation, melting/freezing
1010    ! ----------------------------------------------------------------
1011
1012    ! Initiate the quantity of liquid and solid condensates
1013    ! Note that in the following, zcond is the total amount of condensates
1014    ! including newly formed precipitations (i.e., condensates formed by the
1015    ! condensation process), while zoliq is the total amount of condensates in
1016    ! the cloud (i.e., on which precipitation processes have applied)
1017    DO i = 1, klon
1018      zoliq(i)  = zcond(i)
1019      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1020      zoliqi(i) = zcond(i) * zfice(i)
1021    ENDDO
1022
1023    IF (ok_plane_contrail) THEN
1024      !--Contrails precipitate as natural clouds. We save the partition of ice
1025      !--between natural clouds and contrails
1026      !--NB. we use qcont as a temporary variable to save this partition
1027      DO i = 1, klon
1028        IF ( zoliqi(i) .GT. 0. ) THEN
1029          qcont(i) = ( qcont(i) - zqs(i) * contfra(i) ) / zoliqi(i)
1030        ELSE
1031          qcont(i) = 0.
1032        ENDIF
1033      ENDDO
1034    ENDIF
1035
1036    !================================================================
1037    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
1038    IF (ok_poprecip) THEN
1039
1040      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1041                            ctot_vol(:,k), ptconv(:,k), &
1042                            zt, zq, zoliql, zoliqi, zfice, &
1043                            rneb(:,k), icesed_flux, znebprecipclr, znebprecipcld, &
1044                            zrfl, zrflclr, zrflcld, &
1045                            zifl, ziflclr, ziflcld, &
1046                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1047                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1048                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1049                            dqsmelt(:,k), dqsfreez(:,k), dqised(:,k) &
1050                            )
1051      DO i = 1, klon
1052          zoliq(i) = zoliql(i) + zoliqi(i)
1053      ENDDO
1054
1055    !================================================================
1056    ELSE
1057
1058      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1059                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
1060                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
1061                            rneb(:,k), znebprecipclr, znebprecipcld, &
1062                            zneb, tot_zneb, zrho_up, zvelo_up, &
1063                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
1064                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k), dqised(:,k) &
1065                            )
1066
1067    ENDIF ! ok_poprecip
1068   
1069    IF (ok_plane_contrail) THEN
1070      !--Contrails fraction is left unchanged, but contrails water has changed
1071      DO i = 1, klon
1072        IF ( zoliqi(i) .LE. 0. ) THEN
1073          contfra(i) = 0.
1074          qcont(i) = 0.
1075        ELSE
1076          qcont(i) = zqs(i) * contfra(i) + zoliqi(i) * qcont(i)
1077        ENDIF
1078      ENDDO
1079    ENDIF
1080
1081    ! End of precipitation processes after cloud formation
1082    ! ----------------------------------------------------
1083
1084    !----------------------------------------------------------------------
1085    ! P4> Calculation of cloud condensates amount seen by radiative scheme
1086    !----------------------------------------------------------------------
1087
1088    DO i=1,klon
1089
1090      IF (ok_poprecip) THEN
1091        IF (ok_radocond_snow) THEN
1092          radocond(i,k) = zoliq(i)
1093          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
1094        ELSE
1095          radocond(i,k) = zoliq(i)
1096          zradoice(i) = zoliqi(i)
1097        ENDIF
1098      ELSE
1099        radocond(i,k) = zradocond(i)
1100      ENDIF
1101
1102      ! calculate the percentage of ice in "radocond" so cloud+precip seen
1103      ! by the radiation scheme
1104      IF (radocond(i,k) .GT. 0.) THEN
1105        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
1106      ENDIF
1107    ENDDO
1108
1109    ! ----------------------------------------------------------------
1110    ! P5> Wet scavenging
1111    ! ----------------------------------------------------------------
1112
1113    !Scavenging through nucleation in the layer
1114
1115    DO i = 1,klon
1116       
1117        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1118            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1119        ELSE
1120            beta(i,k) = 0.
1121        ENDIF
1122
1123        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1124
1125        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1126
1127            IF (temp(i,k) .GE. t_glace_min) THEN
1128                zalpha_tr = a_tr_sca(3)
1129            ELSE
1130                zalpha_tr = a_tr_sca(4)
1131            ENDIF
1132
1133            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1134            frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1135
1136            ! Nucleation with a factor of -1 instead of -0.5
1137            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1138
1139        ENDIF
1140
1141    ENDDO
1142
1143    ! Scavenging through impaction in the underlying layer
1144
1145    DO kk = k-1, 1, -1
1146
1147        DO i = 1, klon
1148
1149            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1150
1151                IF (temp(i,kk) .GE. t_glace_min) THEN
1152                    zalpha_tr = a_tr_sca(1)
1153                ELSE
1154                    zalpha_tr = a_tr_sca(2)
1155                ENDIF
1156
1157              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1158              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1159
1160            ENDIF
1161
1162        ENDDO
1163
1164    ENDDO
1165   
1166    !------------------------------------------------------------
1167    ! P6 > write diagnostics and outputs
1168    !------------------------------------------------------------
1169
1170    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1171    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1172   
1173    !--AB Write diagnostics and tracers for ice supersaturation
1174    IF ( ok_plane_contrail ) THEN
1175      cfa_seri(:,k) = contfra(:)
1176      pcf_seri(:,k) = perscontfra(:)
1177      qva_seri(:,k) = zqs(:) * contfra(:)
1178      qia_seri(:,k) = qcont(:) - zqs(:) * contfra(:)
1179      DO i = 1, klon
1180        IF ( ( rneb(i,k) - cfa_seri(i,k) ) .GT. eps ) THEN
1181          qice_perscont(i,k) = ( radocond(i,k) - qia_seri(i,k) ) &
1182              * perscontfra(i) / ( rneb(i,k) - cfa_seri(i,k) )
1183        ELSE
1184          qice_perscont(i,k) = 0.
1185        ENDIF
1186      ENDDO
1187    ENDIF
1188
1189    IF ( ok_ice_supersat ) THEN
1190
1191      DO i = 1, klon
1192
1193        !--If prognostic clouds are activated, deep convection vapor is
1194        !--re-added to the total water vapor
1195        IF ( ptconv(i,k) .AND. pt_pron_clds(i) ) &
1196            zq(i) = zq(i) + ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
1197
1198        cf_seri(i,k) = rneb(i,k)
1199
1200        IF ( zoliq(i) .LE. 0. ) THEN
1201          !--If everything was precipitated, the remaining empty cloud is dissipated
1202          !--and everything is transfered to the subsaturated clear sky region
1203          !--NB. we do not change rneb, as it is a diagnostic only
1204          cf_seri(i,k) = 0.
1205          qvc(i) = 0.
1206        ENDIF
1207
1208        qvc_seri(i,k) = qvc(i)
1209
1210        !--Diagnostics
1211        gamma_cond(i,k) = gammasat(i)
1212        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1213        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1214        qcld(i,k) = qvc(i) + zoliq(i)
1215
1216        !--Calculation of the ice supersaturated fraction following Lamquin et al (2012)
1217        !--methodology: in each layer, we make a maximum random overlap assumption for
1218        !--ice supersaturation
1219        IF ( ( paprs(i,k) .GT. 10000. ) .AND. ( paprs(i,k) .LE. 15000. ) ) THEN
1220                IF ( issrfra100to150UP(i) .GT. ( 1. - eps ) ) THEN
1221                        issrfra100to150(i) = 1.
1222                ELSE
1223                        issrfra100to150(i) = 1. - ( 1. - issrfra100to150(i) ) * &
1224                                ( 1. - MAX( issrfra(i,k), issrfra100to150UP(i) ) ) &
1225                              / ( 1. - issrfra100to150UP(i) )
1226                        issrfra100to150UP(i) = issrfra(i,k)
1227                ENDIF
1228        ELSEIF ( ( paprs(i,k) .GT. 15000. ) .AND. ( paprs(i,k) .LE. 20000. ) ) THEN
1229                IF ( issrfra150to200UP(i) .GT. ( 1. - eps ) ) THEN
1230                        issrfra150to200(i) = 1.
1231                ELSE
1232                        issrfra150to200(i) = 1. - ( 1. - issrfra150to200(i) ) * &
1233                                ( 1. - MAX( issrfra(i,k), issrfra150to200UP(i) ) ) &
1234                              / ( 1. - issrfra150to200UP(i) )
1235                        issrfra150to200UP(i) = issrfra(i,k)
1236                ENDIF
1237        ELSEIF ( ( paprs(i,k) .GT. 20000. ) .AND. ( paprs(i,k) .LE. 25000. ) ) THEN
1238                IF ( issrfra200to250UP(i) .GT. ( 1. - eps ) ) THEN
1239                        issrfra200to250(i) = 1.
1240                ELSE
1241                        issrfra200to250(i) = 1. - ( 1. - issrfra200to250(i) ) * &
1242                                ( 1. - MAX( issrfra(i,k), issrfra200to250UP(i) ) ) &
1243                              / ( 1. - issrfra200to250UP(i) )
1244                        issrfra200to250UP(i) = issrfra(i,k)
1245                ENDIF
1246        ELSEIF ( ( paprs(i,k) .GT. 25000. ) .AND. ( paprs(i,k) .LE. 30000. ) ) THEN
1247                IF ( issrfra250to300UP(i) .GT. ( 1. - eps ) ) THEN
1248                        issrfra250to300(i) = 1.
1249                ELSE
1250                        issrfra250to300(i) = 1. - ( 1. - issrfra250to300(i) ) * &
1251                                ( 1. - MAX( issrfra(i,k), issrfra250to300UP(i) ) ) &
1252                              / ( 1. - issrfra250to300UP(i) )
1253                        issrfra250to300UP(i) = issrfra(i,k)
1254                ENDIF
1255        ELSEIF ( ( paprs(i,k) .GT. 30000. ) .AND. ( paprs(i,k) .LE. 40000. ) ) THEN
1256                IF ( issrfra300to400UP(i) .GT. ( 1. - eps ) ) THEN
1257                        issrfra300to400(i) = 1.
1258                ELSE
1259                        issrfra300to400(i) = 1. - ( 1. - issrfra300to400(i) ) * &
1260                                ( 1. - MAX( issrfra(i,k), issrfra300to400UP(i) ) ) &
1261                              / ( 1. - issrfra300to400UP(i) )
1262                        issrfra300to400UP(i) = issrfra(i,k)
1263                ENDIF
1264        ELSEIF ( ( paprs(i,k) .GT. 40000. ) .AND. ( paprs(i,k) .LE. 50000. ) ) THEN
1265                IF ( issrfra400to500UP(i) .GT. ( 1. - eps ) ) THEN
1266                        issrfra400to500(i) = 1.
1267                ELSE
1268                        issrfra400to500(i) = 1. - ( 1. - issrfra400to500(i) ) * &
1269                                ( 1. - MAX( issrfra(i,k), issrfra400to500UP(i) ) ) &
1270                              / ( 1. - issrfra400to500UP(i) )
1271                        issrfra400to500UP(i) = issrfra(i,k)
1272                ENDIF
1273        ENDIF
1274
1275      ENDDO
1276    ENDIF
1277
1278    ! Outputs:
1279    !-------------------------------
1280    ! Precipitation fluxes at layer interfaces
1281    ! + precipitation fractions +
1282    ! temperature and water species tendencies
1283    DO i = 1, klon
1284        psfl(i,k)=zifl(i)
1285        prfl(i,k)=zrfl(i)
1286        pfraclr(i,k)=znebprecipclr(i)
1287        pfracld(i,k)=znebprecipcld(i)
1288        d_q(i,k) = zq(i) - qt(i,k)
1289        d_t(i,k) = zt(i) - temp(i,k)
1290
1291        IF (ok_bug_phase_lscp) THEN
1292           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1293           d_qi(i,k) = zfice(i)*zoliq(i)
1294        ELSE
1295           d_ql(i,k) = zoliql(i)
1296           d_qi(i,k) = zoliqi(i)   
1297        ENDIF
1298
1299    ENDDO
1300
1301
1302ENDDO ! loop on k from top to bottom
1303
1304
1305  ! Rain or snow at the surface (depending on the first layer temperature)
1306  DO i = 1, klon
1307      snow(i) = zifl(i)
1308      rain(i) = zrfl(i)
1309      ! c_iso final output
1310  ENDDO
1311
1312  IF ( ok_ice_sedim ) THEN
1313    DO i = 1, klon
1314      snow(i) = snow(i) + icesed_flux(i)
1315    ENDDO
1316  ENDIF
1317
1318  IF (ncoreczq>0) THEN
1319      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1320  ENDIF
1321
1322
1323RETURN
1324
1325END SUBROUTINE lscp
1326!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1327
1328END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.