Index: LMDZ6/trunk/libf/phylmd/lscp_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4062)
+++ LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4072)
@@ -6,7 +6,7 @@
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-SUBROUTINE LSCP(dtime,missing_val,                      & 
+SUBROUTINE LSCP(dtime,missing_val,                      &
      paprs,pplay,t,q,ptconv,ratqs,                      &
-     d_t, d_q, d_ql, d_qi, rneb, rneb_seri,             & 
+     d_t, d_q, d_ql, d_qi, rneb, rneb_seri,             &
      radliq, radicefrac, rain, snow,                    &
      pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
@@ -25,5 +25,5 @@
 !      
 ! This code is a new version of the fisrtilp.F90 routine, which is itself a 
-! fusion of 'first' (superrsaturation physics, P. LeVan K. Laval)
+! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
 ! and 'ilp' (il pleut, L. Li)
 !
@@ -126,5 +126,5 @@
   INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
                                                                    ! CR: if iflag_ice_thermo=2, only convection is active    
-  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated 
+  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated
 
   LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
@@ -138,5 +138,5 @@
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk          ! exner potential (p/100000)**(R/cp) 
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztla            ! liquid temperature within thermals [K]
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl            ! liquid potential temperature [K]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl         ! liquid potential temperature [K]
 
   ! INPUT/OUTPUT variables
@@ -145,4 +145,5 @@
   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs            ! function of pressure that sets the large-scale
                                                                     ! cloud PDF (sigma=ratqs*qt)
+  
 
   ! Input sursaturation en glace
@@ -204,18 +205,21 @@
   PARAMETER (ztfondue=278.15)
 
-  REAL, SAVE    :: rain_int_min=0.001                               ! Minimum local rain intensity [mm/s] before the decrease in  associates precipitation fraction 
+  REAL, SAVE    :: rain_int_min=0.001                               ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction 
   !$OMP THREADPRIVATE(rain_int_min)
 
-
-  
+  LOGICAL, SAVE :: ok_radliq_precip=.false.                         ! Inclusion of all precipitation in radliq (cloud water seen by radiation)
+  !$OMP THREADPRIVATE(ok_radliq_precip)
+
+
+
+
   ! LOCAL VARIABLES:
   !----------------
 
 
-  REAL qsl, qsi
+  REAL qsl(klon), qsi(klon)
   REAL zct, zcl
   REAL ctot(klon,klev)
   REAL ctot_vol(klon,klev)
-  INTEGER mpc_bl_points(klon,klev)
   REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5  
   REAL zdqsdT_raw(klon)
@@ -230,11 +234,13 @@
   REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
   REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
-  REAL zoliqp(klon), zoliqi(klon)
+  REAL zoliql(klon), zoliqi(klon)
   REAL zt(klon)
-  REAL zdz(klon),zrho(klon),ztot, zrhol(klon)
+  REAL zdz(klon),zrho(klon),iwc(klon)
   REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
-  REAL zmelt,zpluie,zice
+  REAL zmelt,zrain,zsnow,zprecip
   REAL dzfice(klon)
-  REAL zsolid,qtot,dqsl,dqsi
+  REAL zsolid
+  REAL qtot(klon), qzero(klon)
+  REAL dqsl(klon),dqsi(klon)
   REAL smallestreal
   !  Variables for Bergeron process
@@ -262,11 +268,12 @@
   REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)                                              
   REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
-  REAL velo(klon)
+  REAL velo(klon,klev), vr
   REAL qlmpc, qimpc, rnebmpc
-  REAL radliqi(klon,klev)
+  REAL radliqi(klon,klev), radliql(klon,klev)
 
   INTEGER i, k, n, kk
   INTEGER n_i(klon), iter
   INTEGER ncoreczq  
+  INTEGER mpc_bl_points(klon,klev)
   INTEGER,SAVE :: itap=0
   !$OMP THREADPRIVATE(itap)
@@ -334,8 +341,10 @@
     CALL getin_p('seuil_neb',seuil_neb)
     CALL getin_p('rain_int_min',rain_int_min) 
+    CALL getin_p('ok_radliq_precip',ok_radliq_precip)
     WRITE(lunout,*) 'lscp, ninter:', ninter
     WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec
     WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb
     WRITE(lunout,*) 'lscp, rain_int_min:', rain_int_min
+    WRITE(lunout,*) 'lscp, ok_radliq_precip:', ok_radliq_precip
 
     ! check for precipitation sub-time steps
@@ -397,4 +406,5 @@
 tot_znebn(:) = 0.0                                                                               
 d_tot_zneb(:) = 0.0   
+qzero(:)=0.0
 
 !--ice sursaturation
@@ -468,5 +478,4 @@
     ! --------------------------------------------------------------------
     ! A part of the precipitation coming from above is evaporated/sublimated.
-    !
     ! --------------------------------------------------------------------
 
@@ -474,4 +483,11 @@
     IF (iflag_evap_prec.GE.1) THEN
 
+        ! Calculation of saturation specific humidity
+        ! depending on temperature:
+        CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
+        ! wrt liquid water
+        CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
+        ! wrt ice
+        CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
 
         DO i = 1, klon
@@ -480,6 +496,4 @@
             IF (zrfl(i)+zifl(i).GT.0.) THEN
                     
-                CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
-
             ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
                 IF (iflag_evap_prec.EQ.4) THEN                                                          
@@ -503,10 +517,5 @@
                 ENDIF   
 
-                ! A. JAM
-                ! We consider separately qsatice and qstal
-                ! qsat wrt liquid phase according to thermcep
-                CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,1,.false.,qsl,dqsl)
-
-         
+
                 ! Evaporation of liquid precipitation coming from above
                 ! dP/dz=beta*(1-q/qsat)*sqrt(P)
@@ -515,13 +524,13 @@
 
                 IF (iflag_evap_prec.EQ.3) THEN
-                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) &
+                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
                     *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 ELSE IF (iflag_evap_prec.EQ.4) THEN
-                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl) &
+                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
                     *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 ELSE
-                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
+                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 ENDIF
@@ -531,18 +540,15 @@
                 *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 
          
-                ! qsat wrt ice according to thermcep
-                CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,2,.false.,qsi,dqsi)
-
                 ! sublimation of the solid precipitation coming from above
                 IF (iflag_evap_prec.EQ.3) THEN
-                    zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) &
+                    zqevti = znebprecip(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) &
                     *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 ELSE IF (iflag_evap_prec.EQ.4) THEN
-                     zqevti = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsi) &
+                     zqevti = znebprecipclr(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) &
                     *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 ELSE
-                    zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
+                    zqevti = 1.*coef_eva_i*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
                     *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
                 ENDIF
@@ -566,19 +572,4 @@
 
 
-                ! EV: agrees with JL's comments below so correct and comment
-                ! JLD: I do not understand the lines below. We distribute the liquid 
-                ! and solid parts of the precipitation even if the layer does not get
-                ! saturated. To my opinion, we should all replace with:
-                ! zqev=zqevt
-                ! zqevi=zqevti
-                !    IF (zqevt+zqevti.GT.0.) THEN
-                !        zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
-                !        zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
-                !    ELSE
-                !        zqev=0.
-                !        zqevi=0.
-                !    ENDIF
-                ! ENDIF
-
                 ! New solid and liquid precipitation fluxes
                 zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
@@ -659,11 +650,10 @@
     ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
     !-------------------------------------------------------
-       
-        DO i = 1, klon
+      
+     qtot(:)=zq(:)+zmqc(:)
+     CALL CALC_QSAT_ECMWF(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
+     DO i = 1, klon
             zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
-            qtot=zq(i)+zmqc(i)
-            CALL CALC_QSAT_ECMWF(zt(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
             zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
-         
             IF (zq(i) .LT. 1.e-15) THEN
                 ncoreczq=ncoreczq+1
@@ -671,5 +661,5 @@
             ENDIF
 
-        ENDDO
+     ENDDO
     
 
@@ -777,5 +767,5 @@
                 DO i=1,klon
 
-                    ! convergence = .true. since convergence is not satisfied
+                    ! convergence = .true. until when convergence is not satisfied
                     convergence(i)=ABS(DT(i)).GT.DDT0
                    
@@ -796,10 +786,22 @@
 
                         ! Rneb, qzn and zcond for lognormal PDFs
-                        qtot=zq(i)+zmqc(i)
-                        CALL CALC_QSAT_ECMWF(Tbef(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))
-                        CALL CALC_GAMMASAT(Tbef(i),qtot,pplay(i,k),gammasat(i),dgammasatdt(i))
-                        
-                        ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
-                        zdqs(i) = gammasat(i)*zdqs(i)+zqs(i)*dgammasatdt(i) 
+                        qtot(i)=zq(i)+zmqc(i)
+                    
+                      ENDIF
+
+                  ENDDO
+
+                  ! Calculation of saturation specific humidity and ce fraction
+                  CALL CALC_QSAT_ECMWF(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
+                  CALL CALC_GAMMASAT(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
+                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
+                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
+                  CALL ICEFRAC_LSCP(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
+
+
+
+                  DO i=1,klon
+
+                      IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
 
                         zpdf_sig(i)=ratqs(i,k)*zq(i)
@@ -826,17 +828,18 @@
                           ENDIF
 
-                          rnebss(i,k)=0.0   !--ajout OB (necessaire car boucle de convergence sur le temps)
+                          rnebss(i,k)=0.0   !--added by OB (needed because of the convergence loop on time)
                           fcontrN(i,k)=0.0  !--idem
                           fcontrP(i,k)=0.0  !--idem
                           qss(i,k)=0.0      !--idem
-
+                       
                         ELSE
+               
                         !------------------------------------
-                        ! SURSATURATION EN GLACE
+                        ! ICE SUPERSATURATION
                         !------------------------------------
 
                         CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
-                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),               &
-                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                 &
+                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),               &  
+                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                 & 
                              Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
 
@@ -850,12 +853,10 @@
 
 
-            ! P2.2.2> Approximative calculation of temperature variation DT 
-            !           due to condensation.
-            ! Calculated variables: 
-            ! dT : temperature change due to condensation
-            ! ---------------------------------------------------------------
-
-                       ! EV: calculation of icefrac in one sole function
-                        CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
+                       ! P2.2.2> Approximative calculation of temperature variation DT 
+                       !           due to condensation.
+                       ! Calculated variables: 
+                       ! dT : temperature change due to condensation
+                       !---------------------------------------------------------------
+
                 
                         IF (zfice(i).LT.1) THEN
@@ -896,5 +897,5 @@
 
         ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
-        CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
+        CALL ICEFRAC_LSCP(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
 
         DO i=1, klon
@@ -1002,14 +1003,15 @@
 
     DO i = 1, klon
-        IF (rneb(i,k).GT.0.0) THEN
             zoliq(i) = zcond(i)
-            zrho(i) = pplay(i,k) / zt(i) / RD
-            zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
-    
-            zneb(i) = MAX(rneb(i,k), seuil_neb)
+            zoliqi(i)= zoliq(i)*zfice(i)
+            zoliql(i)= zoliq(i)*(1.-zfice(i))
+            zrho(i)  = pplay(i,k) / zt(i) / RD
+            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
+            iwc(i)   = 0. 
+            zneb(i)  = MAX(rneb(i,k),seuil_neb)
             radliq(i,k) = zoliq(i)/REAL(ninter+1)
             radicefrac(i,k) = zfice(i)
             radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
-        ENDIF
+            radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
     ENDDO
 
@@ -1024,9 +1026,9 @@
         DO i=1,klon
             IF (rneb(i,k).GT.0.0) THEN
-                zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+                iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
             ENDIF
         ENDDO
 
-        CALL FALLICE_VELOCITY(klon,zrhol(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:))
+        CALL FALLICE_VELOCITY(klon,iwc(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:,k))
 
         DO i = 1, klon
@@ -1034,13 +1036,13 @@
             IF (rneb(i,k).GT.0.0) THEN
 
-                ! Initialization of zpluie, zice and ztot:
-                zpluie=0.
-                zice=0.
-                ztot=0.
+                ! Initialization of zrain, zsnow and zprecip:
+                zrain=0.
+                zsnow=0.
+                zprecip=0.
 
                 IF (zneb(i).EQ.seuil_neb) THEN
-                    ztot = 0.0
-                    zice = 0.0
-                    zpluie= 0.0
+                    zprecip = 0.0
+                    zsnow = 0.0
+                    zrain= 0.0
                 ELSE
                     ! water quantity to remove: zchau (Sundqvist, 1978)
@@ -1049,10 +1051,10 @@
                         zcl=cld_lc_con
                         zct=1./cld_tau_con
-                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i)*zfice(i)
+                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*zfice(i)
                     ELSE
                         zcl=cld_lc_lsc
                         zct=1./cld_tau_lsc
                         zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
-                            *velo(i) * zfice(i)
+                            *velo(i,k) * zfice(i)
                     ENDIF
 
@@ -1061,4 +1063,5 @@
                     ! surface fraction (which is larger and artificially
                     ! reduces the in-cloud water).
+
                     IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
                         zchau = zct   *dtime/REAL(ninter) * zoliq(i) &
@@ -1069,22 +1072,19 @@
                     ENDIF
 
-                    zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
-                    zice   = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 
-                    ztot   = MAX(zpluie + zice,0.0)
+                    zrain   = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
+                    zsnow   = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 
+                    zprecip = MAX(zrain + zsnow,0.0)
 
                 ENDIF
             
-                ztot    = MIN(ztot,zoliq(i))
-                zice    = MIN(zice,ztot)
-                zpluie  = MIN(zpluie,ztot)
-
-                zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie  , 0.0)
-                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice  , 0.0)
-                zoliq(i)  = MAX(zoliq(i)-ztot , 0.0)
+                zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
+                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
+                zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
 
                 radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
+                radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1)
                 radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
            
-            ENDIF
+            ENDIF ! rneb >0
 
         ENDDO  ! i = 1,klon
@@ -1093,7 +1093,17 @@
 
     
-    ! Caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
+    ! Include the contribution to non-evaporated precipitation to radliq 
+    IF ((ok_radliq_precip) .AND. (k .LT. klev)) THEN
+        ! for rain drops, we assume a fallspeed of 5m/s
+        vr=5.0
+        radliql(:,k)=radliql(:,k)+zrfl(:)/vr/zrho(:)
+        ! for ice crystals, we take the fallspeed from level above 
+        radliqi(:,k)=radliqi(:,k)+zifl(:)/zrho(:)/velo(:,k+1)
+        radliq(:,k)=radliql(:,k)+radliqi(:,k)
+    ENDIF
+
+    ! caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
     DO i=1,klon
-        IF (radliq(i,k) .GT. 0) THEN
+        IF (radliq(i,k) .GT. 0.) THEN
             radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
         ENDIF
@@ -1101,45 +1111,44 @@
 
 
-    ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
-    ! If T<0C, liquid precip are converted into ice, which leads to an increase in
-    ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
-    ! taken into account through a linearization of the equations and by approximating
-    ! the liquid precipitation process with a "threshold" process. We assume that 
-    ! condensates are not modified during this operation. Liquid precipitation is
-    ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. 
-    ! Water vapor increases as well
-    ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
-
-
+    ! Precipitation flux calculation
 
     DO i = 1, klon
 
             IF (rneb(i,k) .GT. 0.0) THEN
-
-                zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
-                zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
-                zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
-                coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
-                ! Computation of DT if all the liquid precip freezes
-                DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
-                ! T should not exceed the freezing point
-                ! that is Delta > RTT-zt(i)
-                DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
-                zt(i)  = zt(i) + DeltaT
-                ! water vaporization due to temp. increase
-                Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
-                ! we add this vaporized water to the total vapor and we remove it from the precip
-                zq(i) = zq(i) +  Deltaq
-                ! The three "max" lines herebelow protect from rounding errors
-                zcond(i) = max( zcond(i)- Deltaq, 0. )
-                ! liquid precipitation freezes oe evaporates
-                Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
-                zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
-                ! iced water budget
-                zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
+                
+
+               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
+               ! If T<0C, liquid precip are converted into ice, which leads to an increase in
+               ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
+               ! taken into account through a linearization of the equations and by approximating
+               ! the liquid precipitation process with a "threshold" process. We assume that
+               ! condensates are not modified during this operation. Liquid precipitation is
+               ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
+               ! Water vapor increases as well
+               ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
+
+                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
+                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
+                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
+                    coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
+                    ! Computation of DT if all the liquid precip freezes
+                    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
+                    ! T should not exceed the freezing point
+                    ! that is Delta > RTT-zt(i)
+                    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
+                    zt(i)  = zt(i) + DeltaT
+                    ! water vaporization due to temp. increase
+                    Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
+                    ! we add this vaporized water to the total vapor and we remove it from the precip
+                    zq(i) = zq(i) +  Deltaq
+                    ! The three "max" lines herebelow protect from rounding errors
+                    zcond(i) = max( zcond(i)- Deltaq, 0. )
+                    ! liquid precipitation converted to ice precip 
+                    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
+                    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
+                    ! iced water budget
+                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
+
            
-                d_ql(i,k) = (1-zfice(i))*zoliq(i)
-                d_qi(i,k) = zfice(i)*zoliq(i)
-        
                 IF (iflag_evap_prec.EQ.4) THEN                                                            
                     zrflcld(i) = zrflcld(i)+zqprecl(i) &                                                   
@@ -1165,5 +1174,5 @@
     IF (iflag_evap_prec.EQ.4) THEN  
 
-        DO i=1, klon                                        
+        DO i=1,klon                                        
                 
             IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN     
@@ -1192,8 +1201,10 @@
     ! Outputs:
     ! Precipitation fluxes at layer interfaces
-    ! and temperature and vapor tendencies
+    ! and temperature and water species tendencies
     DO i = 1, klon
         psfl(i,k)=zifl(i) 
         prfl(i,k)=zrfl(i)
+        d_ql(i,k) = (1-zfice(i))*zoliq(i)
+        d_qi(i,k) = zfice(i)*zoliq(i)
         d_q(i,k) = zq(i) - q(i,k)
         d_t(i,k) = zt(i) - t(i,k)
@@ -1264,24 +1275,22 @@
     ENDDO
      
-    !--save some variables for ice sursaturation
+    !--save some variables for ice supersaturation
     !
     DO i = 1, klon
-        ! pour la mémoire
+        ! for memory 
         rneb_seri(i,k) = rneb(i,k)
 
-        ! pour les diagnostics
+        ! for diagnostics
         rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
 
         qvc(i,k) = zqs(i) * rneb(i,k)
-        qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--ajout OB a cause de cas pathologiques avec lognormale=F
+        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
         qcld(i,k) = qvc(i,k) + zcond(i)
-
-        !q_sat
-        CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,1,.false.,zqsatl(i,k),zdqs(i))
-        CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,2,.false.,zqsats(i,k),zdqs(i))
-
-     ENDDO
-
-ENDDO
+   END DO
+   !q_sat
+   CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:))
+   CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:))     
+
+END DO
 
 !======================================================================
