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

Last change on this file since 5631 was 5631, checked in by aborella, 2 months ago

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