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

Last change on this file since 5615 was 5615, checked in by aborella, 8 months ago

Comments + small bugfixes
File size: 60.5 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  !--for Lamquin et al 2012 diagnostics
337  REAL, DIMENSION(klon) :: issrfra100to150UP, issrfra150to200UP, issrfra200to250UP
338  REAL, DIMENSION(klon) :: issrfra250to300UP, issrfra300to400UP, issrfra400to500UP
339
340  INTEGER i, k, kk, iter
341  INTEGER, DIMENSION(klon) :: n_i
342  INTEGER ncoreczq
343  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
344  LOGICAL iftop
345
346  LOGICAL, DIMENSION(klon) :: lognormale
347  LOGICAL, DIMENSION(klon) :: keepgoing
348
349  CHARACTER (len = 20) :: modname = 'lscp'
350  CHARACTER (len = 80) :: abort_message
351 
352
353!===============================================================================
354! INITIALISATION
355!===============================================================================
356
357! Few initial checks
358
359
360IF (iflag_fisrtilp_qsat .LT. 0) THEN
361    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
362    CALL abort_physic(modname,abort_message,1)
363ENDIF
364
365! Few initialisations
366
367ctot_vol(1:klon,1:klev)=0.0
368rneblsvol(1:klon,1:klev)=0.0
369znebprecip(:)=0.0
370znebprecipclr(:)=0.0
371znebprecipcld(:)=0.0
372mpc_bl_points(:,:)=0
373
374IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
375
376! AA for 'safety' reasons
377zalpha_tr   = 0.
378zfrac_lessi = 0.
379beta(:,:)= 0.
380
381! Initialisation of variables:
382
383prfl(:,:) = 0.0
384psfl(:,:) = 0.0
385d_t(:,:) = 0.0
386d_q(:,:) = 0.0
387d_ql(:,:) = 0.0
388d_qi(:,:) = 0.0
389rneb(:,:) = 0.0
390pfraclr(:,:)=0.0
391pfracld(:,:)=0.0
392cldfraliq(:,:)=0.
393sigma2_icefracturb(:,:)=0.
394mean_icefracturb(:,:)=0.
395radocond(:,:) = 0.0
396radicefrac(:,:) = 0.0
397frac_nucl(:,:) = 1.0
398frac_impa(:,:) = 1.0
399rain(:) = 0.0
400snow(:) = 0.0
401zfice(:)=1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
402dzfice(:)=0.0
403zfice_turb(:)=0.0
404dzfice_turb(:)=0.0
405zrfl(:) = 0.0
406zifl(:) = 0.0
407zneb(:) = seuil_neb
408zrflclr(:) = 0.0
409ziflclr(:) = 0.0
410zrflcld(:) = 0.0
411ziflcld(:) = 0.0
412tot_zneb(:) = 0.0
413qzero(:) = 0.0
414zdistcltop(:)=0.0
415ztemp_cltop(:) = 0.0
416ztupnew(:)=0.0
417
418distcltop(:,:)=0.
419temp_cltop(:,:)=0.
420
421!--Ice supersaturation
422gamma_cond(:,:) = 1.
423qissr(:,:)      = 0.
424issrfra(:,:)    = 0.
425dcf_sub(:,:)    = 0.
426dcf_con(:,:)    = 0.
427dcf_mix(:,:)    = 0.
428dcfsed(:,:)     = 0.
429dqi_adj(:,:)    = 0.
430dqi_sub(:,:)    = 0.
431dqi_con(:,:)    = 0.
432dqi_mix(:,:)    = 0.
433dqised(:,:)     = 0.
434dqvc_adj(:,:)   = 0.
435dqvc_sub(:,:)   = 0.
436dqvc_con(:,:)   = 0.
437dqvc_mix(:,:)   = 0.
438dqvcsed(:,:)    = 0.
439qvc(:)          = 0.
440shear(:)        = 0.
441pt_pron_clds(:) = .FALSE.
442
443!--for Lamquin et al (2012) diagnostics
444issrfra100to150(:)   = 0.
445issrfra100to150UP(:) = 0.
446issrfra150to200(:)   = 0.
447issrfra150to200UP(:) = 0.
448issrfra200to250(:)   = 0.
449issrfra200to250UP(:) = 0.
450issrfra250to300(:)   = 0.
451issrfra250to300UP(:) = 0.
452issrfra300to400(:)   = 0.
453issrfra300to400UP(:) = 0.
454issrfra400to500(:)   = 0.
455issrfra400to500UP(:) = 0.
456
457!-- poprecip
458qraindiag(:,:)= 0.
459qsnowdiag(:,:)= 0.
460dqreva(:,:)   = 0.
461dqrauto(:,:)  = 0.
462dqrmelt(:,:)  = 0.
463dqrfreez(:,:) = 0.
464dqrcol(:,:)   = 0.
465dqssub(:,:)   = 0.
466dqsauto(:,:)  = 0.
467dqsrim(:,:)   = 0.
468dqsagg(:,:)   = 0.
469dqsfreez(:,:) = 0.
470dqsmelt(:,:)  = 0.
471zqupnew(:)    = 0.
472zqvapclr(:)   = 0.
473cldfra_above(:)= 0.
474icesed_flux(:)= 0.
475
476
477
478!c_iso: variable initialisation for iso
479
480!===============================================================================
481!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
482!===============================================================================
483
484ncoreczq=0
485
486DO k = klev, 1, -1
487
488    IF (k.LE.klev-1) THEN
489       iftop=.false.
490    ELSE
491       iftop=.true.
492    ENDIF
493
494    ! Initialisation temperature and specific humidity
495    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
496    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
497    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
498    ! d_t = temperature tendency due to lscp
499    ! The temperature of the overlying layer is updated here because needed for thermalization
500    DO i = 1, klon
501        zt(i)=temp(i,k)
502        zq(i)=qt(i,k)
503        qliq_in(i) = ql_seri(i,k)
504        qice_in(i) = qi_seri(i,k)
505        IF (.not. iftop) THEN
506           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
507           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
508           !--zqs(i) is the saturation specific humidity in the layer above
509           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
510        ENDIF
511        !c_iso init of iso
512    ENDDO
513    IF ( ok_ice_supersat ) THEN
514      cldfra_in(:) = cf_seri(:,k)
515      qvc_in(:) = qvc_seri(:,k)
516    ENDIF
517
518    ! --------------------------------------------------------------------
519    ! P1> Precipitation processes, before cloud formation:
520    !     Thermalization of precipitation falling from the overlying layer AND
521    !     Precipitation evaporation/sublimation/melting
522    !---------------------------------------------------------------------
523   
524    !================================================================
525    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
526    IF ( ok_poprecip ) THEN
527
528      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
529                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
530                        zqvapclr, zqupnew, icesed_flux, &
531                        cldfra_in, qvc_in, qliq_in, qice_in, &
532                        zrfl, zrflclr, zrflcld, &
533                        zifl, ziflclr, ziflcld, &
534                        dqreva(:,k), dqssub(:,k) &
535                        )
536
537    !================================================================
538    ELSE
539
540      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
541                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
542                        zrfl, zrflclr, zrflcld, &
543                        zifl, ziflclr, ziflcld, &
544                        dqreva(:,k), dqssub(:,k) &
545                        )
546
547    ENDIF ! (ok_poprecip)
548   
549    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
550    !-------------------------------------------------------
551
552     qtot(:)=zq(:)+zmqc(:)
553     CALL calc_qsat_ecmwf(klon,zt,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
554     DO i = 1, klon
555            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
556            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
557            IF (zq(i) .LT. 1.e-15) THEN
558                ncoreczq=ncoreczq+1
559                zq(i)=1.e-15
560            ENDIF
561            ! c_iso: do something similar for isotopes
562
563     ENDDO
564   
565    ! --------------------------------------------------------------------
566    ! P2> Cloud formation
567    !---------------------------------------------------------------------
568    !
569    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
570    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
571    ! is prescribed and depends on large scale variables and boundary layer
572    ! properties)
573    ! The decrease in condensed part due tu latent heating is taken into
574    ! account
575    ! -------------------------------------------------------------------
576
577        ! P2.1> With the PDFs (log-normal, bigaussian)
578        ! cloud properties calculation with the initial values of t and q
579        ! ----------------------------------------------------------------
580
581        ! initialise gammasat and qincloud_mpc
582        gammasat(:)=1.
583        qincloud_mpc(:)=0.
584
585        IF (iflag_cld_th.GE.5) THEN
586               ! Cloud cover and content in meshes affected by shallow convection,
587               ! are retrieved from a bi-gaussian distribution of the saturation deficit
588               ! following Jam et al. 2013
589
590               IF (iflag_cloudth_vert.LE.2) THEN
591                  ! Old version of Arnaud Jam
592
593                    CALL cloudth(klon,klev,k,tv,                  &
594                         zq,qta,fraca,                            &
595                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
596                         ratqs,zqs,temp,                              &
597                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
598
599
600                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
601                   ! Default version of Arnaud Jam
602                       
603                    CALL cloudth_v3(klon,klev,k,tv,                        &
604                         zq,qta,fraca,                                     &
605                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
606                         ratqs,sigma_qtherm,zqs,temp, &
607                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
608
609
610                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
611                   ! Jean Jouhaud's version, with specific separation between surface and volume
612                   ! cloud fraction Decembre 2018
613
614                    CALL cloudth_v6(klon,klev,k,tv,                        &
615                         zq,qta,fraca,                                     &
616                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
617                         ratqs,zqs,temp, &
618                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
619
620                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
621                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
622                   ! for boundary-layer mixed phase clouds
623                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
624                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
625                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
626                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
627
628                ENDIF
629
630
631                DO i=1,klon
632                    rneb(i,k)=ctot(i,k)
633                    rneblsvol(i,k)=ctot_vol(i,k)
634                    zqn(i)=qcloud(i)
635                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
636                    qvc(i) = rneb(i,k) * zqs(i)
637                ENDDO
638
639        ENDIF
640
641        IF (iflag_cld_th .LE. 4) THEN
642           
643                ! lognormal
644            lognormale(:) = .TRUE.
645
646        ELSEIF (iflag_cld_th .GE. 6) THEN
647           
648                ! lognormal distribution when no thermals
649            lognormale(:) = fraca(:,k) < min_frac_th_cld
650
651        ELSE
652                ! When iflag_cld_th=5, we always assume
653                ! bi-gaussian distribution
654            lognormale(:) = .FALSE.
655       
656        ENDIF
657
658
659        IF ( ok_ice_supersat ) THEN
660
661          !--Initialisation
662          IF ( ok_plane_contrail ) THEN
663            contfra(:)     = 0.
664            qcont(:)       = 0.
665            perscontfra(:) = 0.
666          ENDIF
667
668          DO i = 1, klon
669
670            pt_pron_clds(i) = ( ( ( ( zt(i) .LE. temp_nowater ) .OR. ok_weibull_warm_clouds ) &
671                .AND. ( .NOT. ok_no_issr_strato .OR. ( stratomask(i,k) .EQ. 0. ) ) ) &
672                .AND. ( cfcon(i,k) .LT. ( 1. - eps ) ) )
673
674            IF ( pt_pron_clds(i) ) THEN
675
676              !--If deep convection is activated, the condensation scheme activates
677              !--only in the environment. NB. the clear sky fraction will the be
678              !--maximised by 1. - cfcon(i,k)
679              IF ( ptconv(i,k) ) zq(i) = zq(i) - ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
680
681              IF ( ( cfcon(i,k) * qccon(i,k) ) .LT. ( cfcon_old(i,k) * qccon_old(i,k) ) ) THEN
682                !--If deep convection is weakening, we add the clouds that are not anymore
683                !--'in' deep convection to the advected clouds
684                cldfra_in(i) = cldfra_in(i) + MAX(0., cfcon_old(i,k) - cfcon(i,k))
685                qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * MAX(0., cfcon_old(i,k) - cfcon(i,k))
686                qice_in(i) = qice_in(i) + ( qccon_old(i,k) * cfcon_old(i,k) &
687                                          - qccon(i,k) * cfcon(i,k) )
688              ELSEIF ( cfcon(i,k) .GT. cfcon_old(i,k) ) THEN
689                !--Else if deep convection is strengthening, it consumes the existing cloud
690                !--fraction (which does not at this moment represent deep convection)
691                !--NB. if deep convection is strengthening while the fraction decreases,
692                !--clear sky water vapor will be transfered in priority
693                cldfra_in(i) = cldfra_in(i) * ( 1. &
694                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
695                qvc_in(i)    = qvc_in(i)    * ( 1. &
696                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
697                qice_in(i)   = qice_in(i)   * ( 1. &
698                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
699              ENDIF
700
701              !--Barriers
702              cldfra_in(i) = MAX(0., MIN(1. - cfcon(i,k), cldfra_in(i)))
703              qvc_in(i)    = MAX(0., MIN(zq(i), qvc_in(i)))
704              qice_in(i)   = MAX(0., MIN(zq(i) - qvc_in(i), qice_in(i)))
705
706              !--Calculate the shear value (input for condensation and ice supersat)
707              !--Cell thickness [m]
708              delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(i) * RD
709              IF ( iftop ) THEN
710                ! top
711                shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
712                               + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
713              ELSEIF ( k .EQ. 1 ) THEN
714                ! surface
715                shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
716                               + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
717              ELSE
718                ! other layers
719                shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
720                                   - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
721                               + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
722                                   - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
723              ENDIF
724            ENDIF
725          ENDDO
726        ENDIF
727
728        DT(:) = 0.
729        n_i(:)=0
730        Tbef(:)=zt(:)
731        qlbef(:)=0.
732
733        ! Treatment of non-boundary layer clouds (lognormale)
734        ! We iterate here to take into account the change in qsat(T) and ice fraction
735        ! during the condensation process
736        ! the increment in temperature is calculated using the first principle of
737        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
738        ! and a clear sky fraction)
739        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
740
741        DO iter=1,iflag_fisrtilp_qsat+1
742
743                keepgoing(:) = .FALSE.
744
745                DO i=1,klon
746
747                    ! keepgoing = .true. while convergence is not satisfied
748
749                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
750
751                        ! if not convergence:
752                        ! we calculate a new iteration
753                        keepgoing(i) = .TRUE.
754
755                        ! P2.2.1> cloud fraction and condensed water mass calculation
756                        ! Calculated variables:
757                        ! rneb : cloud fraction
758                        ! zqn : total water within the cloud
759                        ! zcond: mean condensed water within the mesh
760                        ! rhcl: clear-sky relative humidity
761                        !---------------------------------------------------------------
762
763                        ! new temperature that only serves in the iteration process:
764                        Tbef(i)=Tbef(i)+DT(i)
765
766                        ! Rneb, qzn and zcond for lognormal PDFs
767                        qtot(i)=zq(i)+zmqc(i)
768
769                      ENDIF
770
771                  ENDDO
772
773                  ! Calculation of saturation specific humidity and ice fraction
774                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
775                 
776                  IF (iflag_icefrac .GE. 3) THEN
777                  ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
778                  ! and cloud condensed water content. idea following Dietlitcher et al. 2018, GMD
779                  ! we make this option works only for the physically-based and tke-depdenent param (iflag_icefrac>=1)
780                      DO i=1,klon
781                           CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,1,.false.,zqsl,zdqsl)
782                           CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,2,.false.,zqsi,zdqsi)
783                           zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i)
784                           zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i)
785                      ENDDO
786                  ENDIF
787
788                  CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
789                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
790                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
791
792                  ! Cloud condensation based on subgrid pdf
793                  !--AB Activates a condensation scheme that allows for
794                  !--ice supersaturation and contrails evolution from aviation
795                  IF (ok_ice_supersat) THEN
796
797                    IF ( iftop ) THEN
798                      cldfra_above(:) = 0.
799                    ELSE
800                      cldfra_above(:) = rneb(:,k+1)
801                    ENDIF
802
803                    !---------------------------------------------
804                    !--   CONDENSATION AND ICE SUPERSATURATION  --
805                    !---------------------------------------------
806
807                    CALL condensation_ice_supersat( &
808                        klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), &
809                        cfcon(:,k), cldfra_in, qvc_in, qliq_in, qice_in, &
810                        shear, tke_dissip(:,k), cell_area, Tbef, zq, zqs, &
811                        gammasat, ratqs(:,k), keepgoing, pt_pron_clds, &
812                        cldfra_above, icesed_flux,&
813                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
814                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), dcfsed(:,k), &
815                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), dqised(:,k), &
816                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), dqvcsed(:,k), &
817                        cfa_seri(:,k), pcf_seri(:,k), qva_seri(:,k), qia_seri(:,k), &
818                        flight_dist(:,k), flight_h2o(:,k), contfra, perscontfra, qcont, &
819                        Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
820                        dcfa_ini(:,k), dqia_ini(:,k), dqta_ini(:,k), &
821                        dcfa_sub(:,k), dqia_sub(:,k), dqta_sub(:,k), &
822                        dcfa_cir(:,k), dqta_cir(:,k), &
823                        dcfa_mix(:,k), dqia_mix(:,k), dqta_mix(:,k))
824
825
826                  ELSE
827                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
828
829                   CALL condensation_lognormal( &
830                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
831                       keepgoing, rneb(:,k), zqn, qvc)
832
833
834                  ENDIF ! .NOT. ok_ice_supersat
835
836                  ! cloud phase determination
837                  IF (iflag_icefrac .GE. 2) THEN
838                     ! phase partitioning depending on temperature. activates here in the iteration process if iflag_icefrac > 2
839                     CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_in, ziflcld, zqn, &
840                     rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
841                  ELSE
842                     ! phase partitioning depending on temperature and eventually distance to cloud top
843                     IF (iflag_t_glace.GE.4) THEN
844                       ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
845                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
846                     ENDIF
847                     CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
848                  ENDIF
849
850
851                  DO i=1,klon
852                      IF (keepgoing(i)) THEN
853
854                        ! If vertical heterogeneity, change fraction by volume as well
855                        IF (iflag_cloudth_vert.GE.3) THEN
856                            ctot_vol(i,k)=rneb(i,k)
857                            rneblsvol(i,k)=ctot_vol(i,k)
858                        ENDIF
859
860
861                       ! P2.2.2> Approximative calculation of temperature variation DT
862                       !           due to condensation.
863                       ! Calculated variables:
864                       ! dT : temperature change due to condensation
865                       !---------------------------------------------------------------
866
867               
868                        IF (zfice(i).LT.1) THEN
869                            cste=RLVTT
870                        ELSE
871                            cste=RLSTT
872                        ENDIF
873                       
874                        IF ( ok_unadjusted_clouds ) THEN
875                          !--AB We relax the saturation adjustment assumption
876                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
877                          IF ( rneb(i,k) .GT. eps ) THEN
878                            qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
879                          ELSE
880                            qlbef(i) = 0.
881                          ENDIF
882                        ELSE
883                          qlbef(i)=max(0.,zqn(i)-zqs(i))
884                        ENDIF
885
886                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
887                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
888                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
889                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
890                              *qlbef(i)*dzfice(i)
891                        ! here we update a provisory temperature variable that only serves in the iteration
892                        ! process
893                        DT(i)=num/denom
894                        n_i(i)=n_i(i)+1
895
896                    ENDIF ! end keepgoing
897
898                ENDDO     ! end loop on i
899
900        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
901
902        ! P2.2> Final quantities calculation
903        ! Calculated variables:
904        !   rneb : cloud fraction
905        !   zcond: mean condensed water in the mesh
906        !   zqn  : mean water vapor in the mesh
907        !   zfice: ice fraction in clouds
908        !   zt   : temperature
909        !   rhcl : clear-sky relative humidity
910        ! ----------------------------------------------------------------
911
912
913        ! Cloud phase final determination
914        !--------------------------------
915        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
916        IF (iflag_t_glace.GE.4) THEN
917           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
918           distcltop(:,k)=zdistcltop(:)
919           temp_cltop(:,k)=ztemp_cltop(:)
920        ENDIF
921        ! Partition function depending on temperature for all clouds (shallow convective and stratiform)
922        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
923
924        ! Partition function depending on tke for non shallow-convective clouds, erase previous estimation
925        IF (iflag_icefrac .GE. 1) THEN
926           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_in, ziflcld, zqn, &
927           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
928        ENDIF
929
930        ! Water vapor update, subsequent latent heat exchange for each cloud type
931        !------------------------------------------------------------------------
932        DO i=1, klon
933            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
934            ! iflag_cloudth_vert=7 and specific param is activated
935            IF (mpc_bl_points(i,k) .GT. 0) THEN
936                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
937                ! following line is very strange and probably wrong
938                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
939                ! water vapor update and partition function if thermals
940                zq(i) = zq(i) - zcond(i)       
941                zfice(i)=zfice_th(i)
942            ELSE
943                ! Checks on rneb, rhcl and zqn
944                IF (rneb(i,k) .LE. 0.0) THEN
945                    zqn(i) = 0.0
946                    rneb(i,k) = 0.0
947                    zcond(i) = 0.0
948                    rhcl(i,k)=zq(i)/zqs(i)
949                ELSE IF (rneb(i,k) .GE. 1.0) THEN
950                    zqn(i) = zq(i)
951                    rneb(i,k) = 1.0
952                    IF ( ok_unadjusted_clouds ) THEN
953                      !--AB We relax the saturation adjustment assumption
954                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
955                      zcond(i) = MAX(0., zqn(i) - qvc(i))
956                    ELSE
957                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
958                    ENDIF
959                    rhcl(i,k)=1.0
960                ELSE
961                    IF ( ok_unadjusted_clouds ) THEN
962                      !--AB We relax the saturation adjustment assumption
963                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
964                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
965                    ELSE
966                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
967                    ENDIF
968                    ! following line is very strange and probably wrong:
969                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
970                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
971                    IF (iflag_icefrac .GE. 1) THEN
972                        IF (lognormale(i)) THEN 
973                           zcond(i)  = zqliq(i) + zqice(i)
974                           zfice(i)  = zfice_turb(i)
975                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
976                        ENDIF
977                    ENDIF
978                ENDIF
979
980                ! water vapor update
981                zq(i) = zq(i) - zcond(i)
982
983            ENDIF
984           
985           
986            ! temperature update due to phase change
987            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
988            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
989                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
990        ENDDO
991
992        ! If vertical heterogeneity, change volume fraction
993        IF (iflag_cloudth_vert .GE. 3) THEN
994          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
995          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
996        ENDIF
997
998
999    ! ----------------------------------------------------------------
1000    ! P3> Precipitation processes, after cloud formation
1001    !     - precipitation formation, melting/freezing
1002    ! ----------------------------------------------------------------
1003
1004    ! Initiate the quantity of liquid and solid condensates
1005    ! Note that in the following, zcond is the total amount of condensates
1006    ! including newly formed precipitations (i.e., condensates formed by the
1007    ! condensation process), while zoliq is the total amount of condensates in
1008    ! the cloud (i.e., on which precipitation processes have applied)
1009    DO i = 1, klon
1010      zoliq(i)  = zcond(i)
1011      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1012      zoliqi(i) = zcond(i) * zfice(i)
1013    ENDDO
1014
1015    IF (ok_plane_contrail) THEN
1016      !--Contrails precipitate as natural clouds. We save the partition of ice
1017      !--between natural clouds and contrails
1018      !--NB. we use qcont as a temporary variable to save this partition
1019      DO i = 1, klon
1020        IF ( zoliqi(i) .GT. 0. ) THEN
1021          qcont(i) = ( qcont(i) - zqs(i) * contfra(i) ) / zoliqi(i)
1022        ELSE
1023          qcont(i) = 0.
1024        ENDIF
1025      ENDDO
1026    ENDIF
1027
1028    !================================================================
1029    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
1030    IF (ok_poprecip) THEN
1031
1032      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1033                            ctot_vol(:,k), ptconv(:,k), &
1034                            zt, zq, zoliql, zoliqi, zfice, &
1035                            rneb(:,k), icesed_flux, znebprecipclr, znebprecipcld, &
1036                            zrfl, zrflclr, zrflcld, &
1037                            zifl, ziflclr, ziflcld, &
1038                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1039                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1040                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1041                            dqsmelt(:,k), dqsfreez(:,k), dqised(:,k) &
1042                            )
1043      DO i = 1, klon
1044          zoliq(i) = zoliql(i) + zoliqi(i)
1045      ENDDO
1046
1047    !================================================================
1048    ELSE
1049
1050      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1051                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
1052                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
1053                            rneb(:,k), znebprecipclr, znebprecipcld, &
1054                            zneb, tot_zneb, zrho_up, zvelo_up, &
1055                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
1056                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k), dqised(:,k) &
1057                            )
1058
1059    ENDIF ! ok_poprecip
1060   
1061    IF (ok_plane_contrail) THEN
1062      !--Contrails fraction is left unchanged, but contrails water has changed
1063      DO i = 1, klon
1064        IF ( zoliqi(i) .LE. 0. ) THEN
1065          contfra(i) = 0.
1066          qcont(i) = 0.
1067        ELSE
1068          qcont(i) = zqs(i) * contfra(i) + zoliqi(i) * qcont(i)
1069        ENDIF
1070      ENDDO
1071    ENDIF
1072
1073    ! End of precipitation processes after cloud formation
1074    ! ----------------------------------------------------
1075
1076    !----------------------------------------------------------------------
1077    ! P4> Calculation of cloud condensates amount seen by radiative scheme
1078    !----------------------------------------------------------------------
1079
1080    DO i=1,klon
1081
1082      IF (ok_poprecip) THEN
1083        IF (ok_radocond_snow) THEN
1084          radocond(i,k) = zoliq(i)
1085          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
1086        ELSE
1087          radocond(i,k) = zoliq(i)
1088          zradoice(i) = zoliqi(i)
1089        ENDIF
1090      ELSE
1091        radocond(i,k) = zradocond(i)
1092      ENDIF
1093
1094      ! calculate the percentage of ice in "radocond" so cloud+precip seen
1095      ! by the radiation scheme
1096      IF (radocond(i,k) .GT. 0.) THEN
1097        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
1098      ENDIF
1099    ENDDO
1100
1101    ! ----------------------------------------------------------------
1102    ! P5> Wet scavenging
1103    ! ----------------------------------------------------------------
1104
1105    !Scavenging through nucleation in the layer
1106
1107    DO i = 1,klon
1108       
1109        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1110            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1111        ELSE
1112            beta(i,k) = 0.
1113        ENDIF
1114
1115        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1116
1117        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1118
1119            IF (temp(i,k) .GE. t_glace_min) THEN
1120                zalpha_tr = a_tr_sca(3)
1121            ELSE
1122                zalpha_tr = a_tr_sca(4)
1123            ENDIF
1124
1125            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1126            frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1127
1128            ! Nucleation with a factor of -1 instead of -0.5
1129            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1130
1131        ENDIF
1132
1133    ENDDO
1134
1135    ! Scavenging through impaction in the underlying layer
1136
1137    DO kk = k-1, 1, -1
1138
1139        DO i = 1, klon
1140
1141            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1142
1143                IF (temp(i,kk) .GE. t_glace_min) THEN
1144                    zalpha_tr = a_tr_sca(1)
1145                ELSE
1146                    zalpha_tr = a_tr_sca(2)
1147                ENDIF
1148
1149              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1150              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1151
1152            ENDIF
1153
1154        ENDDO
1155
1156    ENDDO
1157   
1158    !------------------------------------------------------------
1159    ! P6 > write diagnostics and outputs
1160    !------------------------------------------------------------
1161
1162    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1163    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1164   
1165    !--AB Write diagnostics and tracers for ice supersaturation
1166    IF ( ok_plane_contrail ) THEN
1167      cfa_seri(:,k) = contfra(:)
1168      pcf_seri(:,k) = perscontfra(:)
1169      qva_seri(:,k) = zqs(:) * contfra(:)
1170      qia_seri(:,k) = qcont(:) - zqs(:) * contfra(:)
1171      DO i = 1, klon
1172        IF ( ( rneb(i,k) - cfa_seri(i,k) ) .GT. eps ) THEN
1173          qice_perscont(i,k) = ( radocond(i,k) - qia_seri(i,k) ) &
1174              * perscontfra(i) / ( rneb(i,k) - cfa_seri(i,k) )
1175        ELSE
1176          qice_perscont(i,k) = 0.
1177        ENDIF
1178      ENDDO
1179    ENDIF
1180
1181    IF ( ok_ice_supersat ) THEN
1182
1183      DO i = 1, klon
1184
1185        !--If prognostic clouds are activated, deep convection vapor is
1186        !--re-added to the total water vapor
1187        IF ( ptconv(i,k) .AND. pt_pron_clds(i) ) &
1188            zq(i) = zq(i) + ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
1189
1190        cf_seri(i,k) = rneb(i,k)
1191
1192        IF ( zoliq(i) .LE. 0. ) THEN
1193          !--If everything was precipitated, the remaining empty cloud is dissipated
1194          !--and everything is transfered to the subsaturated clear sky region
1195          !--NB. we do not change rneb, as it is a diagnostic only
1196          cf_seri(i,k) = 0.
1197          qvc(i) = 0.
1198        ENDIF
1199
1200        qvc_seri(i,k) = qvc(i)
1201
1202        !--Diagnostics
1203        gamma_cond(i,k) = gammasat(i)
1204        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1205        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1206        qcld(i,k) = qvc(i) + zoliq(i)
1207
1208        !--Calculation of the ice supersaturated fraction following Lamquin et al (2012)
1209        !--methodology: in each layer, we make a maximum random overlap assumption for
1210        !--ice supersaturation
1211        IF ( ( paprs(i,k) .GT. 10000. ) .AND. ( paprs(i,k) .LE. 15000. ) ) THEN
1212                IF ( issrfra100to150UP(i) .GT. ( 1. - eps ) ) THEN
1213                        issrfra100to150(i) = 1.
1214                ELSE
1215                        issrfra100to150(i) = 1. - ( 1. - issrfra100to150(i) ) * &
1216                                ( 1. - MAX( issrfra(i,k), issrfra100to150UP(i) ) ) &
1217                              / ( 1. - issrfra100to150UP(i) )
1218                        issrfra100to150UP(i) = issrfra(i,k)
1219                ENDIF
1220        ELSEIF ( ( paprs(i,k) .GT. 15000. ) .AND. ( paprs(i,k) .LE. 20000. ) ) THEN
1221                IF ( issrfra150to200UP(i) .GT. ( 1. - eps ) ) THEN
1222                        issrfra150to200(i) = 1.
1223                ELSE
1224                        issrfra150to200(i) = 1. - ( 1. - issrfra150to200(i) ) * &
1225                                ( 1. - MAX( issrfra(i,k), issrfra150to200UP(i) ) ) &
1226                              / ( 1. - issrfra150to200UP(i) )
1227                        issrfra150to200UP(i) = issrfra(i,k)
1228                ENDIF
1229        ELSEIF ( ( paprs(i,k) .GT. 20000. ) .AND. ( paprs(i,k) .LE. 25000. ) ) THEN
1230                IF ( issrfra200to250UP(i) .GT. ( 1. - eps ) ) THEN
1231                        issrfra200to250(i) = 1.
1232                ELSE
1233                        issrfra200to250(i) = 1. - ( 1. - issrfra200to250(i) ) * &
1234                                ( 1. - MAX( issrfra(i,k), issrfra200to250UP(i) ) ) &
1235                              / ( 1. - issrfra200to250UP(i) )
1236                        issrfra200to250UP(i) = issrfra(i,k)
1237                ENDIF
1238        ELSEIF ( ( paprs(i,k) .GT. 25000. ) .AND. ( paprs(i,k) .LE. 30000. ) ) THEN
1239                IF ( issrfra250to300UP(i) .GT. ( 1. - eps ) ) THEN
1240                        issrfra250to300(i) = 1.
1241                ELSE
1242                        issrfra250to300(i) = 1. - ( 1. - issrfra250to300(i) ) * &
1243                                ( 1. - MAX( issrfra(i,k), issrfra250to300UP(i) ) ) &
1244                              / ( 1. - issrfra250to300UP(i) )
1245                        issrfra250to300UP(i) = issrfra(i,k)
1246                ENDIF
1247        ELSEIF ( ( paprs(i,k) .GT. 30000. ) .AND. ( paprs(i,k) .LE. 40000. ) ) THEN
1248                IF ( issrfra300to400UP(i) .GT. ( 1. - eps ) ) THEN
1249                        issrfra300to400(i) = 1.
1250                ELSE
1251                        issrfra300to400(i) = 1. - ( 1. - issrfra300to400(i) ) * &
1252                                ( 1. - MAX( issrfra(i,k), issrfra300to400UP(i) ) ) &
1253                              / ( 1. - issrfra300to400UP(i) )
1254                        issrfra300to400UP(i) = issrfra(i,k)
1255                ENDIF
1256        ELSEIF ( ( paprs(i,k) .GT. 40000. ) .AND. ( paprs(i,k) .LE. 50000. ) ) THEN
1257                IF ( issrfra400to500UP(i) .GT. ( 1. - eps ) ) THEN
1258                        issrfra400to500(i) = 1.
1259                ELSE
1260                        issrfra400to500(i) = 1. - ( 1. - issrfra400to500(i) ) * &
1261                                ( 1. - MAX( issrfra(i,k), issrfra400to500UP(i) ) ) &
1262                              / ( 1. - issrfra400to500UP(i) )
1263                        issrfra400to500UP(i) = issrfra(i,k)
1264                ENDIF
1265        ENDIF
1266
1267      ENDDO
1268    ENDIF
1269
1270    ! Outputs:
1271    !-------------------------------
1272    ! Precipitation fluxes at layer interfaces
1273    ! + precipitation fractions +
1274    ! temperature and water species tendencies
1275    DO i = 1, klon
1276        psfl(i,k)=zifl(i)
1277        prfl(i,k)=zrfl(i)
1278        pfraclr(i,k)=znebprecipclr(i)
1279        pfracld(i,k)=znebprecipcld(i)
1280        d_q(i,k) = zq(i) - qt(i,k)
1281        d_t(i,k) = zt(i) - temp(i,k)
1282
1283        IF (ok_bug_phase_lscp) THEN
1284           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1285           d_qi(i,k) = zfice(i)*zoliq(i)
1286        ELSE
1287           d_ql(i,k) = zoliql(i)
1288           d_qi(i,k) = zoliqi(i)   
1289        ENDIF
1290
1291    ENDDO
1292
1293
1294ENDDO ! loop on k from top to bottom
1295
1296
1297  ! Rain or snow at the surface (depending on the first layer temperature)
1298  DO i = 1, klon
1299      snow(i) = zifl(i)
1300      rain(i) = zrfl(i)
1301      ! c_iso final output
1302  ENDDO
1303
1304  IF ( ok_ice_sedim ) THEN
1305    DO i = 1, klon
1306      snow(i) = snow(i) + icesed_flux(i)
1307    ENDDO
1308  ENDIF
1309
1310  IF (ncoreczq>0) THEN
1311      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1312  ENDIF
1313
1314
1315RETURN
1316
1317END SUBROUTINE lscp
1318!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1319
1320END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.