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

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

Merge with trunk testing r5597. We have convergence in prod and debug in NPv7.0.1c
File size: 61.9 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                ! Correct calculation of clear-sky relative humidity. To activate once
940                ! issues related to zqn>zq in bi-gaussian clouds are corrected
941                !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
942                !--This is slighly approximated, the actual formula is
943                !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
944                !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
945                !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
946                ! rhcl(i,k) = ( zq(i) - qincloud_mpc(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
947                ! water vapor update and partition function if thermals
948                zq(i) = zq(i) - zcond(i)       
949                zfice(i)=zfice_th(i)
950            ELSE
951                ! Checks on rneb, rhcl and zqn
952                IF (rneb(i,k) .LE. 0.0) THEN
953                    zqn(i) = 0.0
954                    rneb(i,k) = 0.0
955                    zcond(i) = 0.0
956                    rhcl(i,k)=zq(i)/zqs(i)
957                ELSE IF (rneb(i,k) .GE. 1.0) THEN
958                    zqn(i) = zq(i)
959                    rneb(i,k) = 1.0
960                    IF ( ok_unadjusted_clouds ) THEN
961                      !--AB We relax the saturation adjustment assumption
962                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
963                      zcond(i) = MAX(0., zqn(i) - qvc(i))
964                    ELSE
965                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
966                    ENDIF
967                    rhcl(i,k)=1.0
968                ELSE
969                    IF ( ok_unadjusted_clouds ) THEN
970                      !--AB We relax the saturation adjustment assumption
971                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
972                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
973                    ELSE
974                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
975                    ENDIF
976                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
977                    IF (iflag_icefrac .GE. 1) THEN
978                        IF (lognormale(i)) THEN 
979                           zcond(i)  = zqliq(i) + zqice(i)
980                           zfice(i)  = zfice_turb(i)
981                        ENDIF
982                    ENDIF
983
984                    ! following line is very strange and probably wrong
985                    rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
986                    ! Correct calculation of clear-sky relative humidity. To activate once
987                    ! issues related to zqn>zq in bi-gaussian clouds are corrected
988                    !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
989                    !--This is slighly approximated, the actual formula is
990                    !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
991                    !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
992                    !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
993                    ! rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
994                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
995
996                ENDIF
997
998                ! water vapor update
999                zq(i) = zq(i) - zcond(i)
1000
1001            ENDIF
1002           
1003           
1004            ! temperature update due to phase change
1005            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1006            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1007                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1008        ENDDO
1009
1010        ! If vertical heterogeneity, change volume fraction
1011        IF (iflag_cloudth_vert .GE. 3) THEN
1012          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
1013          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
1014        ENDIF
1015
1016
1017    ! ----------------------------------------------------------------
1018    ! P3> Precipitation processes, after cloud formation
1019    !     - precipitation formation, melting/freezing
1020    ! ----------------------------------------------------------------
1021
1022    ! Initiate the quantity of liquid and solid condensates
1023    ! Note that in the following, zcond is the total amount of condensates
1024    ! including newly formed precipitations (i.e., condensates formed by the
1025    ! condensation process), while zoliq is the total amount of condensates in
1026    ! the cloud (i.e., on which precipitation processes have applied)
1027    DO i = 1, klon
1028      zoliq(i)  = zcond(i)
1029      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1030      zoliqi(i) = zcond(i) * zfice(i)
1031    ENDDO
1032
1033    IF (ok_plane_contrail) THEN
1034      !--Contrails precipitate as natural clouds. We save the partition of ice
1035      !--between natural clouds and contrails
1036      !--NB. we use qcont as a temporary variable to save this partition
1037      DO i = 1, klon
1038        IF ( zoliqi(i) .GT. 0. ) THEN
1039          qcont(i) = ( qcont(i) - zqs(i) * contfra(i) ) / zoliqi(i)
1040        ELSE
1041          qcont(i) = 0.
1042        ENDIF
1043      ENDDO
1044    ENDIF
1045
1046    !================================================================
1047    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
1048    IF (ok_poprecip) THEN
1049
1050      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1051                            ctot_vol(:,k), ptconv(:,k), &
1052                            zt, zq, zoliql, zoliqi, zfice, &
1053                            rneb(:,k), icesed_flux, znebprecipclr, znebprecipcld, &
1054                            zrfl, zrflclr, zrflcld, &
1055                            zifl, ziflclr, ziflcld, &
1056                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1057                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1058                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1059                            dqsmelt(:,k), dqsfreez(:,k), dqised(:,k) &
1060                            )
1061      DO i = 1, klon
1062          zoliq(i) = zoliql(i) + zoliqi(i)
1063      ENDDO
1064
1065    !================================================================
1066    ELSE
1067
1068      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
1069                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
1070                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
1071                            rneb(:,k), znebprecipclr, znebprecipcld, &
1072                            zneb, tot_zneb, zrho_up, zvelo_up, &
1073                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
1074                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k), dqised(:,k) &
1075                            )
1076
1077    ENDIF ! ok_poprecip
1078   
1079    IF (ok_plane_contrail) THEN
1080      !--Contrails fraction is left unchanged, but contrails water has changed
1081      DO i = 1, klon
1082        IF ( zoliqi(i) .LE. 0. ) THEN
1083          contfra(i) = 0.
1084          qcont(i) = 0.
1085        ELSE
1086          qcont(i) = zqs(i) * contfra(i) + zoliqi(i) * qcont(i)
1087        ENDIF
1088      ENDDO
1089    ENDIF
1090
1091    ! End of precipitation processes after cloud formation
1092    ! ----------------------------------------------------
1093
1094    !----------------------------------------------------------------------
1095    ! P4> Calculation of cloud condensates amount seen by radiative scheme
1096    !----------------------------------------------------------------------
1097
1098    DO i=1,klon
1099
1100      IF (ok_poprecip) THEN
1101        IF (ok_radocond_snow) THEN
1102          radocond(i,k) = zoliq(i)
1103          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
1104        ELSE
1105          radocond(i,k) = zoliq(i)
1106          zradoice(i) = zoliqi(i)
1107        ENDIF
1108      ELSE
1109        radocond(i,k) = zradocond(i)
1110      ENDIF
1111
1112      ! calculate the percentage of ice in "radocond" so cloud+precip seen
1113      ! by the radiation scheme
1114      IF (radocond(i,k) .GT. 0.) THEN
1115        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
1116      ENDIF
1117    ENDDO
1118
1119    ! ----------------------------------------------------------------
1120    ! P5> Wet scavenging
1121    ! ----------------------------------------------------------------
1122
1123    !Scavenging through nucleation in the layer
1124
1125    DO i = 1,klon
1126       
1127        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1128            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1129        ELSE
1130            beta(i,k) = 0.
1131        ENDIF
1132
1133        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1134
1135        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1136
1137            IF (temp(i,k) .GE. t_glace_min) THEN
1138                zalpha_tr = a_tr_sca(3)
1139            ELSE
1140                zalpha_tr = a_tr_sca(4)
1141            ENDIF
1142
1143            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1144            frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1145
1146            ! Nucleation with a factor of -1 instead of -0.5
1147            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1148
1149        ENDIF
1150
1151    ENDDO
1152
1153    ! Scavenging through impaction in the underlying layer
1154
1155    DO kk = k-1, 1, -1
1156
1157        DO i = 1, klon
1158
1159            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1160
1161                IF (temp(i,kk) .GE. t_glace_min) THEN
1162                    zalpha_tr = a_tr_sca(1)
1163                ELSE
1164                    zalpha_tr = a_tr_sca(2)
1165                ENDIF
1166
1167              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1168              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1169
1170            ENDIF
1171
1172        ENDDO
1173
1174    ENDDO
1175   
1176    !------------------------------------------------------------
1177    ! P6 > write diagnostics and outputs
1178    !------------------------------------------------------------
1179
1180    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1181    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1182   
1183    !--AB Write diagnostics and tracers for ice supersaturation
1184    IF ( ok_plane_contrail ) THEN
1185      cfa_seri(:,k) = contfra(:)
1186      pcf_seri(:,k) = perscontfra(:)
1187      qva_seri(:,k) = zqs(:) * contfra(:)
1188      qia_seri(:,k) = qcont(:) - zqs(:) * contfra(:)
1189      DO i = 1, klon
1190        IF ( ( rneb(i,k) - cfa_seri(i,k) ) .GT. eps ) THEN
1191          qice_perscont(i,k) = ( radocond(i,k) - qia_seri(i,k) ) &
1192              * perscontfra(i) / ( rneb(i,k) - cfa_seri(i,k) )
1193        ELSE
1194          qice_perscont(i,k) = 0.
1195        ENDIF
1196      ENDDO
1197    ENDIF
1198
1199    IF ( ok_ice_supersat ) THEN
1200
1201      DO i = 1, klon
1202
1203        !--If prognostic clouds are activated, deep convection vapor is
1204        !--re-added to the total water vapor
1205        IF ( ptconv(i,k) .AND. pt_pron_clds(i) ) &
1206            zq(i) = zq(i) + ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
1207
1208        cf_seri(i,k) = rneb(i,k)
1209
1210        IF ( zoliq(i) .LE. 0. ) THEN
1211          !--If everything was precipitated, the remaining empty cloud is dissipated
1212          !--and everything is transfered to the subsaturated clear sky region
1213          !--NB. we do not change rneb, as it is a diagnostic only
1214          cf_seri(i,k) = 0.
1215          qvc(i) = 0.
1216        ENDIF
1217
1218        qvc_seri(i,k) = qvc(i)
1219
1220        !--Diagnostics
1221        gamma_cond(i,k) = gammasat(i)
1222        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1223        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1224        qcld(i,k) = qvc(i) + zoliq(i)
1225
1226        !--Calculation of the ice supersaturated fraction following Lamquin et al (2012)
1227        !--methodology: in each layer, we make a maximum random overlap assumption for
1228        !--ice supersaturation
1229        IF ( ( paprs(i,k) .GT. 10000. ) .AND. ( paprs(i,k) .LE. 15000. ) ) THEN
1230                IF ( issrfra100to150UP(i) .GT. ( 1. - eps ) ) THEN
1231                        issrfra100to150(i) = 1.
1232                ELSE
1233                        issrfra100to150(i) = 1. - ( 1. - issrfra100to150(i) ) * &
1234                                ( 1. - MAX( issrfra(i,k), issrfra100to150UP(i) ) ) &
1235                              / ( 1. - issrfra100to150UP(i) )
1236                        issrfra100to150UP(i) = issrfra(i,k)
1237                ENDIF
1238        ELSEIF ( ( paprs(i,k) .GT. 15000. ) .AND. ( paprs(i,k) .LE. 20000. ) ) THEN
1239                IF ( issrfra150to200UP(i) .GT. ( 1. - eps ) ) THEN
1240                        issrfra150to200(i) = 1.
1241                ELSE
1242                        issrfra150to200(i) = 1. - ( 1. - issrfra150to200(i) ) * &
1243                                ( 1. - MAX( issrfra(i,k), issrfra150to200UP(i) ) ) &
1244                              / ( 1. - issrfra150to200UP(i) )
1245                        issrfra150to200UP(i) = issrfra(i,k)
1246                ENDIF
1247        ELSEIF ( ( paprs(i,k) .GT. 20000. ) .AND. ( paprs(i,k) .LE. 25000. ) ) THEN
1248                IF ( issrfra200to250UP(i) .GT. ( 1. - eps ) ) THEN
1249                        issrfra200to250(i) = 1.
1250                ELSE
1251                        issrfra200to250(i) = 1. - ( 1. - issrfra200to250(i) ) * &
1252                                ( 1. - MAX( issrfra(i,k), issrfra200to250UP(i) ) ) &
1253                              / ( 1. - issrfra200to250UP(i) )
1254                        issrfra200to250UP(i) = issrfra(i,k)
1255                ENDIF
1256        ELSEIF ( ( paprs(i,k) .GT. 25000. ) .AND. ( paprs(i,k) .LE. 30000. ) ) THEN
1257                IF ( issrfra250to300UP(i) .GT. ( 1. - eps ) ) THEN
1258                        issrfra250to300(i) = 1.
1259                ELSE
1260                        issrfra250to300(i) = 1. - ( 1. - issrfra250to300(i) ) * &
1261                                ( 1. - MAX( issrfra(i,k), issrfra250to300UP(i) ) ) &
1262                              / ( 1. - issrfra250to300UP(i) )
1263                        issrfra250to300UP(i) = issrfra(i,k)
1264                ENDIF
1265        ELSEIF ( ( paprs(i,k) .GT. 30000. ) .AND. ( paprs(i,k) .LE. 40000. ) ) THEN
1266                IF ( issrfra300to400UP(i) .GT. ( 1. - eps ) ) THEN
1267                        issrfra300to400(i) = 1.
1268                ELSE
1269                        issrfra300to400(i) = 1. - ( 1. - issrfra300to400(i) ) * &
1270                                ( 1. - MAX( issrfra(i,k), issrfra300to400UP(i) ) ) &
1271                              / ( 1. - issrfra300to400UP(i) )
1272                        issrfra300to400UP(i) = issrfra(i,k)
1273                ENDIF
1274        ELSEIF ( ( paprs(i,k) .GT. 40000. ) .AND. ( paprs(i,k) .LE. 50000. ) ) THEN
1275                IF ( issrfra400to500UP(i) .GT. ( 1. - eps ) ) THEN
1276                        issrfra400to500(i) = 1.
1277                ELSE
1278                        issrfra400to500(i) = 1. - ( 1. - issrfra400to500(i) ) * &
1279                                ( 1. - MAX( issrfra(i,k), issrfra400to500UP(i) ) ) &
1280                              / ( 1. - issrfra400to500UP(i) )
1281                        issrfra400to500UP(i) = issrfra(i,k)
1282                ENDIF
1283        ENDIF
1284
1285      ENDDO
1286    ENDIF
1287
1288    ! Outputs:
1289    !-------------------------------
1290    ! Precipitation fluxes at layer interfaces
1291    ! + precipitation fractions +
1292    ! temperature and water species tendencies
1293    DO i = 1, klon
1294        psfl(i,k)=zifl(i)
1295        prfl(i,k)=zrfl(i)
1296        pfraclr(i,k)=znebprecipclr(i)
1297        pfracld(i,k)=znebprecipcld(i)
1298        d_q(i,k) = zq(i) - qt(i,k)
1299        d_t(i,k) = zt(i) - temp(i,k)
1300
1301        IF (ok_bug_phase_lscp) THEN
1302           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1303           d_qi(i,k) = zfice(i)*zoliq(i)
1304        ELSE
1305           d_ql(i,k) = zoliql(i)
1306           d_qi(i,k) = zoliqi(i)   
1307        ENDIF
1308
1309    ENDDO
1310
1311
1312ENDDO ! loop on k from top to bottom
1313
1314
1315  ! Rain or snow at the surface (depending on the first layer temperature)
1316  DO i = 1, klon
1317      snow(i) = zifl(i)
1318      rain(i) = zrfl(i)
1319      ! c_iso final output
1320  ENDDO
1321
1322  IF ( ok_ice_sedim ) THEN
1323    DO i = 1, klon
1324      snow(i) = snow(i) + icesed_flux(i)
1325    ENDDO
1326  ENDIF
1327
1328  IF (ncoreczq>0) THEN
1329      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1330  ENDIF
1331
1332
1333RETURN
1334
1335END SUBROUTINE lscp
1336!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1337
1338END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.