source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90 @ 5051

Last change on this file since 5051 was 5007, checked in by evignon, 3 months ago

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