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

Last change on this file since 5006 was 4915, checked in by evignon, 8 months ago

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