source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90 @ 5408

Last change on this file since 5408 was 5408, checked in by evignon, 5 hours ago

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