Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90	(revision 4686)
+++ LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90	(revision 4803)
@@ -18,5 +18,7 @@
      qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
      Tcontr, qcontr, qcontr2, fcontrN, fcontrP,         &
-      cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
+     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
+     dqreva,dqssub,dqrauto,dqrcol,dqrmelt,dqrfreez,dqsauto, &
+     dqsagg,dqsrim,dqsmelt,dqsfreez)
 
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -96,4 +98,5 @@
 USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
 USE ice_sursat_mod, ONLY : ice_sursat
+USE lmdz_lscp_poprecip, ONLY : poprecip_evapsub, poprecip_postcld
 
 ! Use du module lmdz_lscp_ini contenant les constantes
@@ -106,4 +109,5 @@
 USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
+USE lmdz_lscp_ini, ONLY : ok_poprecip
 
 IMPLICIT NONE
@@ -197,4 +201,17 @@
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv ! std of saturation deficit in environment
 
+  ! for POPRECIP
+
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !-- rain tendendy due to evaporation [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !-- snow tendency due to sublimation [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !-- rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !-- snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !-- rain tendency due to autoconversion of cloud liquid [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !-- snow tendency due to autoconversion of cloud ice [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !-- snow tendency due to riming [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !-- snow tendency due to melting [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !-- rain tendency due to melting [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !-- snow tendency due to freezing [kg/kg/s]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !-- rain tendency due to freezing [kg/kg/s]
 
 
@@ -243,7 +260,8 @@
   REAL :: zfrac_lessi
   REAL, DIMENSION(klon) :: zprec_cond
-  REAL, DIMENSION(klon) :: zmair
+  REAL :: zmair
   REAL :: zcpair, zcpeau
   REAL, DIMENSION(klon) :: zlh_solid
+  REAL, DIMENSION(klon) :: ztupnew
   REAL :: zm_solid         ! for liquid -> solid conversion
   REAL, DIMENSION(klon) :: zrflclr, zrflcld
@@ -266,4 +284,5 @@
   INTEGER ncoreczq
   INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
+  LOGICAL iftop
 
   LOGICAL, DIMENSION(klon) :: lognormale
@@ -339,4 +358,5 @@
 distcltop1D(:)=0.0
 temp_cltop1D(:) = 0.0
+ztupnew(:)=0.0
 !--ice supersaturation
 gamma_ss(:,:) = 1.0
@@ -350,4 +370,19 @@
 distcltop(:,:)=0.
 temp_cltop(:,:)=0.
+!-- poprecip
+dqreva(:,:)=0.0
+dqrauto(:,:)=0.0
+dqrmelt(:,:)=0.0
+dqrfreez(:,:)=0.0
+dqrcol(:,:)=0.0
+dqssub(:,:)=0.0
+dqsauto(:,:)=0.0
+dqsrim(:,:)=0.0
+dqsagg(:,:)=0.0
+dqsfreez(:,:)=0.0
+dqsmelt(:,:)=0.0
+
+
+
 
 !c_iso: variable initialisation for iso 
@@ -361,4 +396,10 @@
 
 DO k = klev, 1, -1
+
+    IF (k.LE.klev-1) THEN
+       iftop=.false. 
+    ELSE 
+       iftop=.true.
+    ENDIF
 
     ! Initialisation temperature and specific humidity
@@ -366,9 +407,26 @@
         zt(i)=temp(i,k)
         zq(i)=qt(i,k)
+        IF (.not. iftop) THEN
+           ztupnew(i)=temp(i,k+1)+d_t(i,k+1)
+        ENDIF
         !c_iso init of iso
     ENDDO
-    
+   
+    !================================================================
+    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
+    IF (ok_poprecip) THEN
+
+            CALL poprecip_evapsub(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
+                              zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
+                              zrfl, zrflclr, zrflcld, &
+                              zifl, ziflclr, ziflcld, &
+                              dqreva(:,k),dqssub(:,k) &
+                             )
+
+    !================================================================
+    ELSE
+
     ! --------------------------------------------------------------------
-    ! P0> Thermalization of precipitation falling from the overlying layer
+    ! P1> Thermalization of precipitation falling from the overlying layer
     ! --------------------------------------------------------------------
     ! Computes air temperature variation due to enthalpy transported by
@@ -385,9 +443,15 @@
     ! ---------------------------------------------------------------------
     
-    IF (k.LE.klev-1) THEN 
+    IF (iftop) THEN
+
+        DO i = 1, klon
+            zmqc(i) = 0.
+        ENDDO
+
+    ELSE
 
         DO i = 1, klon
  
-            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+            zmair=(paprs(i,k)-paprs(i,k+1))/RG
             ! no condensed water so cp=cp(vapor+dry air)
             ! RVTMP2=rcpv/rcpd-1
@@ -398,16 +462,9 @@
             ! layer's air so that precipitation at the ground has the
             ! same temperature as the lowermost layer
-            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
+            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
             ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
-            zt(i) = ( (temp(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
+            zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) &
              / (zcpair + zmqc(i)*zcpeau)
 
-        ENDDO
- 
-    ELSE
- 
-        DO i = 1, klon 
-            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG 
-            zmqc(i) = 0.
         ENDDO
  
@@ -415,9 +472,9 @@
 
     ! --------------------------------------------------------------------
-    ! P1> Precipitation evaporation/sublimation
+    ! P2> Precipitation evaporation/sublimation/melting
     ! --------------------------------------------------------------------
-    ! A part of the precipitation coming from above is evaporated/sublimated.
+    ! A part of the precipitation coming from above is evaporated/sublimated/melted.
     ! --------------------------------------------------------------------
-
+    
     IF (iflag_evap_prec.GE.1) THEN 
 
@@ -506,5 +563,4 @@
                 ! redistribute zqev0 conserving the ratio liquid/ice 
 
-                ! todoan: check the consistency of this formula
                 IF (zqevt+zqevti.GT.zqev0) THEN
                     zqev=zqev0*zqevt/(zqevt+zqevti)
@@ -524,8 +580,12 @@
 
                 ! vapor, temperature, precip fluxes update
+                ! vapor is updated after evaporation/sublimation (it is increased)
                 zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
                 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+                ! zmqc is the total condensed water in the precip flux (it is decreased)
                 zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
                 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+                ! air and precip temperature (i.e., gridbox temperature)
+                ! is updated due to latent heat cooling
                 zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
                 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
@@ -579,5 +639,6 @@
                 ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
 
-                ! Latent heat of melting with precipitation thermalization
+                ! Latent heat of melting because of precipitation melting
+                ! NB: the air + precip temperature is simultaneously updated
                 zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
                 *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
@@ -590,7 +651,9 @@
             ENDIF ! (zrfl(i)+zifl(i).GT.0.)
 
-        ENDDO
+        ENDDO ! loop on klon
 
     ENDIF ! (iflag_evap_prec>=1)
+
+    ENDIF ! (ok_poprecip)
 
     ! --------------------------------------------------------------------
@@ -647,7 +710,4 @@
 
 
-
-
-
                 ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
 
@@ -861,5 +921,5 @@
 
         ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
-        CALL icefrac_lscp(klon,zt(:),iflag_ice_thermo,distcltop1D(:),temp_cltop1D(:),zfice(:), dzfice(:))
+        CALL icefrac_lscp(klon,zt,iflag_ice_thermo,distcltop1D,temp_cltop1D,zfice, dzfice)
 
 
@@ -918,4 +978,38 @@
     ! ----------------------------------------------------------------
 
+    !================================================================
+    IF (ok_poprecip) THEN
+
+      DO i = 1, klon
+        zoliql(i) = zcond(i) * ( 1. - zfice(i) )
+        zoliqi(i) = zcond(i) * zfice(i)
+      ENDDO
+
+      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
+                            ctot_vol(:,k), ptconv(:,k), &
+                            zt, zq, zoliql, zoliqi, zfice, &
+                            rneb(:,k), znebprecipclr, znebprecipcld, &
+                            zrfl, zrflclr, zrflcld, &
+                            zifl, ziflclr, ziflcld, &
+                            dqrauto(:,k),dqrcol(:,k),dqrmelt(:,k),dqrfreez(:,k), &
+                            dqsauto(:,k),dqsagg(:,k),dqsrim(:,k),dqsmelt(:,k),dqsfreez(:,k) &
+                            )
+
+      DO i = 1, klon
+        zoliq(i) = zoliql(i) + zoliqi(i)
+        IF ( zoliq(i) .GT. 0. ) THEN 
+                zfice(i) = zoliqi(i)/zoliq(i) 
+        ELSE  
+                zfice(i) = 0.0
+        ENDIF
+        ! when poprecip activated, radiation does not see any precipitation content
+        radocond(i,k) = zoliq(i)
+        radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
+        radocondi(i,k)= radocond(i,k)*zfice(i)
+      ENDDO
+
+    !================================================================
+    ELSE
+
     ! LTP:
     IF (iflag_evap_prec .GE. 4) THEN
@@ -982,5 +1076,4 @@
             zneb(i)  = MAX(rneb(i,k),seuil_neb)
             radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
-            radicefrac(i,k) = zfice(i)
             radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
             radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
@@ -1195,32 +1288,22 @@
         ENDDO
 
-
     ENDIF
+
+    ENDIF ! ok_poprecip
 
     ! End of precipitation formation
     ! --------------------------------
 
-    ! Outputs:
-    ! Precipitation fluxes at layer interfaces
-    ! + precipitation fractions +
-    ! temperature and water species tendencies
-    DO i = 1, klon
-        psfl(i,k)=zifl(i) 
-        prfl(i,k)=zrfl(i)
-        pfraclr(i,k)=znebprecipclr(i)
-        pfracld(i,k)=znebprecipcld(i)
-        d_ql(i,k) = (1-zfice(i))*zoliq(i)
-        d_qi(i,k) = zfice(i)*zoliq(i)
-        d_q(i,k) = zq(i) - qt(i,k)
-        ! c_iso: same for isotopes
-        d_t(i,k) = zt(i) - temp(i,k)
-    ENDDO
+
+    ! Calculation of cloud condensates amount seen by radiative scheme
+    !-----------------------------------------------------------------
 
     ! Calculation of the concentration of condensates seen by the radiation scheme
     ! for liquid phase, we take radocondl
     ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true) 
-    ! we recaulate radocondi to account for contributions from the precipitation flux
-
-    IF ((ok_radocond_snow) .AND. (k .LT. klev)) THEN
+    ! we recalculate radocondi to account for contributions from the precipitation flux
+    ! TODO: for the moment, we deactivate this option when ok_poprecip=.true.
+
+    IF ((ok_radocond_snow) .AND. (k .LT. klev) .AND. (.NOT. ok_poprecip)) THEN
         ! for the solid phase (crystals + snowflakes)
         ! we recalculate radocondi by summing
@@ -1319,9 +1402,26 @@
    CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
 
+
+
+    ! Outputs:
+    !-------------------------------
+    ! Precipitation fluxes at layer interfaces
+    ! + precipitation fractions +
+    ! temperature and water species tendencies
+    DO i = 1, klon
+        psfl(i,k)=zifl(i)
+        prfl(i,k)=zrfl(i)
+        pfraclr(i,k)=znebprecipclr(i)
+        pfracld(i,k)=znebprecipcld(i)
+        d_ql(i,k) = (1-zfice(i))*zoliq(i)
+        d_qi(i,k) = zfice(i)*zoliq(i)
+        d_q(i,k) = zq(i) - qt(i,k)
+        ! c_iso: same for isotopes
+        d_t(i,k) = zt(i) - temp(i,k)
+    ENDDO
+
+
 ENDDO
 
-!======================================================================
-!                      END OF VERTICAL LOOP
-!======================================================================
 
   ! Rain or snow at the surface (depending on the first layer temperature)
