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

Last change on this file since 4803 was 4803, checked in by evignon, 4 months ago

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