Index: /LMDZ6/trunk/DefLists/field_def_lmdz.xml
===================================================================
--- /LMDZ6/trunk/DefLists/field_def_lmdz.xml	(revision 4802)
+++ /LMDZ6/trunk/DefLists/field_def_lmdz.xml	(revision 4803)
@@ -643,5 +643,16 @@
 	<field id="rneblsvol" long_name="LS Cloud fraction by volume"    unit="-" />
         <field id="pfraclr" long_name="LS precip fraction clear-sky part"    unit="-" />		
-        <field id="pfracld" long_name="LS precip fraction cloudy part"    unit="-" />   
+	<field id="pfracld" long_name="LS precip fraction cloudy part"    unit="-" />  
+	<field id="dqreva" long_name="LS rain tendency evaporation"    unit="kg/m2/s" />
+        <field id="dqrauto" long_name="LS rain tendency autoconversion"    unit="kg/m2/s" />
+        <field id="dqrcol" long_name="LS rain tendency collection"    unit="kg/m2/s" />
+        <field id="dqrmelt" long_name="LS rain tendency melting"    unit="kg/m2/s" />
+        <field id="dqrfreez" long_name="LS rain tendency freezing"    unit="kg/m2/s" />
+        <field id="dqssub" long_name="LS snow tendency sublimation"    unit="kg/m2/s" />
+        <field id="dqsauto" long_name="LS snow tendency autoconversion"    unit="kg/m2/s" />
+        <field id="dqsagg" long_name="LS snow tendency aggregation"    unit="kg/m2/s" />
+        <field id="dqsmelt" long_name="LS snow tendency melting"    unit="kg/m2/s" />
+	<field id="dqsfreez" long_name="LS snow tendency freezing"    unit="kg/m2/s" />
+	<field id="dqsrim" long_name="LS snow tendency riming"    unit="kg/m2/s" />
         <field id="rhum"    long_name="Relative humidity"    unit="-" />
 	<field id="rhl" long_name="Relative humidity wrt liquid" unit="%" />
Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90	(revision 4802)
+++ /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)
Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90	(revision 4802)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90	(revision 4803)
@@ -9,5 +9,5 @@
   !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG)
 
-  REAL, SAVE :: seuil_neb=0.001                 ! cloud fraction threshold: a cloud really exists when exceeded
+  REAL, SAVE, PROTECTED :: seuil_neb=0.001      ! cloud fraction threshold: a cloud really exists when exceeded
   !$OMP THREADPRIVATE(seuil_neb)
 
@@ -15,8 +15,8 @@
   !$OMP THREADPRIVATE(lunout,prt_level)
 
-  INTEGER, SAVE :: niter_lscp=5                 ! number of iterations to calculate autoconversion to precipitation
+  INTEGER, SAVE, PROTECTED :: niter_lscp=5      ! number of iterations to calculate autoconversion to precipitation
   !$OMP THREADPRIVATE(niter_lscp)
 
-  INTEGER,SAVE :: iflag_evap_prec=1             ! precipitation evaporation flag. 0: nothing, 1: "old way", 
+  INTEGER, SAVE, PROTECTED :: iflag_evap_prec=1 ! precipitation evaporation flag. 0: nothing, 1: "old way", 
                                                 ! 2: Max cloud fraction above to calculate the max of reevaporation
                                                 ! >=4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only
@@ -25,111 +25,164 @@
   !$OMP THREADPRIVATE(iflag_evap_prec)
 
-  REAL t_coup                                   ! temperature threshold which determines the phase 
-  PARAMETER (t_coup=234.0)                      ! for which the saturation vapor pressure is calculated
-
-  REAL DDT0                                     ! iteration parameter 
-  PARAMETER (DDT0=.01)
-
-  REAL ztfondue                                 ! parameter to calculate melting fraction of precipitation
-  PARAMETER (ztfondue=278.15)
-
-  REAL temp_nowater                             ! temperature below which liquid water no longer exists
-  PARAMETER (temp_nowater=233.15)
-
-  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)
-
-  REAL, SAVE :: a_tr_sca(4)                     ! Variables for tracers temporary: alpha parameter for scavenging, 4 possible scavenging processes
+  REAL, SAVE, PROTECTED :: t_coup=234.0         ! temperature threshold which determines the phase 
+                                                ! for which the saturation vapor pressure is calculated
+  !$OMP THREADPRIVATE(t_coup)
+  REAL, SAVE, PROTECTED :: DDT0=0.01            ! iteration parameter 
+  !$OMP THREADPRIVATE(DDT0)
+
+  REAL, SAVE, PROTECTED :: ztfondue=278.15      ! parameter to calculate melting fraction of precipitation
+  !$OMP THREADPRIVATE(ztfondue)
+
+  REAL, SAVE, PROTECTED :: temp_nowater=233.15  ! temperature below which liquid water no longer exists
+  !$OMP THREADPRIVATE(temp_nowater)
+
+  REAL, SAVE, PROTECTED :: a_tr_sca(4)          ! Variables for tracers temporary: alpha parameter for scavenging, 4 possible scavenging processes
   !$OMP THREADPRIVATE(a_tr_sca)
   
-  INTEGER, SAVE ::  iflag_mpc_bl=0              ! flag to activate boundary layer mixed phase cloud param
+  INTEGER, SAVE, PROTECTED ::  iflag_mpc_bl=0   ! flag to activate boundary layer mixed phase cloud param
   !$OMP THREADPRIVATE(iflag_mpc_bl)
   
-  LOGICAL, SAVE :: ok_radocond_snow=.false.       ! take into account the mass of ice precip in the cloud ice content seen by radiation
+  LOGICAL, SAVE, PROTECTED :: ok_radocond_snow=.false. ! take into account the mass of ice precip in the cloud ice content seen by radiation
   !$OMP THREADPRIVATE(ok_radocond_snow)
 
-  REAL, SAVE :: t_glace_min=258.0                ! lower-bound temperature parameter for cloud phase determination
+  REAL, SAVE, PROTECTED :: t_glace_min=258.0    ! lower-bound temperature parameter for cloud phase determination
   !$OMP THREADPRIVATE(t_glace_min)
 
-  REAL, SAVE :: t_glace_max=273.15               ! upper-bound temperature parameter for cloud phase determination
+  REAL, SAVE, PROTECTED :: t_glace_max=273.15   ! upper-bound temperature parameter for cloud phase determination
   !$OMP THREADPRIVATE(t_glace_max)
 
-  REAL, SAVE :: exposant_glace=1.0               ! parameter for cloud phase determination
+  REAL, SAVE, PROTECTED :: exposant_glace=1.0   ! parameter for cloud phase determination
   !$OMP THREADPRIVATE(exposant_glace)
 
-  INTEGER, SAVE :: iflag_vice=0                  ! which expression for ice crystall fall velocity
+  INTEGER, SAVE, PROTECTED :: iflag_vice=0      ! which expression for ice crystall fall velocity
   !$OMP THREADPRIVATE(iflag_vice)
 
-  INTEGER, SAVE :: iflag_t_glace=0               ! which expression for cloud phase partitioning
+  INTEGER, SAVE, PROTECTED :: iflag_t_glace=0   ! which expression for cloud phase partitioning
   !$OMP THREADPRIVATE(iflag_t_glace)
 
-  INTEGER, SAVE :: iflag_cloudth_vert=0          ! option for determining cloud fraction and content in convective boundary layers
+  INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert=0         ! option for determining cloud fraction and content in convective boundary layers
   !$OMP THREADPRIVATE(iflag_cloudth_vert)
 
-  INTEGER, SAVE :: iflag_gammasat=0              ! which threshold for homogeneous nucleation below -40oC
+  INTEGER, SAVE, PROTECTED :: iflag_gammasat=0             ! which threshold for homogeneous nucleation below -40oC
   !$OMP THREADPRIVATE(iflag_gammasat)
 
-  INTEGER, SAVE :: iflag_rain_incloud_vol=0      ! use of volume cloud fraction for rain autoconversion 
+  INTEGER, SAVE, PROTECTED :: iflag_rain_incloud_vol=0     ! use of volume cloud fraction for rain autoconversion 
   !$OMP THREADPRIVATE(iflag_rain_incloud_vol)
 
-  INTEGER, SAVE :: iflag_bergeron=0              ! bergeron effect for liquid precipitation treatment  
+  INTEGER, SAVE, PROTECTED :: iflag_bergeron=0             ! bergeron effect for liquid precipitation treatment  
   !$OMP THREADPRIVATE(iflag_bergeron)
 
-  INTEGER, SAVE :: iflag_fisrtilp_qsat=0         ! qsat adjustment (iterative) during autoconversion
+  INTEGER, SAVE, PROTECTED :: iflag_fisrtilp_qsat=0        ! qsat adjustment (iterative) during autoconversion
   !$OMP THREADPRIVATE(iflag_fisrtilp_qsat)
 
-  INTEGER, SAVE :: iflag_pdf=0                   ! type of subgrid scale qtot pdf
+  INTEGER, SAVE, PROTECTED :: iflag_pdf=0                  ! type of subgrid scale qtot pdf
   !$OMP THREADPRIVATE(iflag_pdf)
 
-  INTEGER, SAVE :: iflag_autoconversion=0        ! autoconversion option
+  INTEGER, SAVE, PROTECTED :: iflag_autoconversion=0       ! autoconversion option
   !$OMP THREADPRIVATE(iflag_autoconversion)
 
-  LOGICAL, SAVE :: reevap_ice=.false.            ! no liquid precip for T< threshold
+  LOGICAL, SAVE, PROTECTED :: reevap_ice=.false.           ! no liquid precip for T< threshold
   !$OMP THREADPRIVATE(reevap_ice)
 
-  REAL, SAVE :: cld_lc_lsc=2.6e-4                ! liquid autoconversion coefficient, stratiform rain
+  REAL, SAVE, PROTECTED :: cld_lc_lsc=2.6e-4               ! liquid autoconversion coefficient, stratiform rain
   !$OMP THREADPRIVATE(cld_lc_lsc)
 
-  REAL, SAVE :: cld_lc_con=2.6e-4                ! liquid autoconversion coefficient, convective rain
+  REAL, SAVE, PROTECTED :: cld_lc_con=2.6e-4                ! liquid autoconversion coefficient, convective rain
   !$OMP THREADPRIVATE(cld_lc_con)
 
-  REAL, SAVE :: cld_tau_lsc=3600.                ! liquid autoconversion timescale, stratiform rain
+  REAL, SAVE, PROTECTED :: cld_tau_lsc=3600.                ! liquid autoconversion timescale, stratiform rain
   !$OMP THREADPRIVATE(cld_tau_lsc)
 
-  REAL, SAVE :: cld_tau_con=3600.                ! liquid autoconversion timescale, convective rain
+  REAL, SAVE, PROTECTED :: cld_tau_con=3600.                ! liquid autoconversion timescale, convective rain
   !$OMP THREADPRIVATE(cld_tau_con)
 
-  REAL, SAVE :: cld_expo_lsc=2.                  ! liquid autoconversion threshold exponent, stratiform rain
+  REAL, SAVE, PROTECTED :: cld_expo_lsc=2.                  ! liquid autoconversion threshold exponent, stratiform rain
   !$OMP THREADPRIVATE(cld_expo_lsc)
 
-  REAL, SAVE :: cld_expo_con=2.                  ! liquid autoconversion threshold exponent, convective rain
+  REAL, SAVE, PROTECTED :: cld_expo_con=2.                  ! liquid autoconversion threshold exponent, convective rain
   !$OMP THREADPRIVATE(cld_expo_con)
 
-  REAL, SAVE :: ffallv_lsc=1.                    ! tuning coefficient crystal fall velocity, stratiform 
+  REAL, SAVE, PROTECTED :: ffallv_lsc=1.                    ! tuning coefficient crystal fall velocity, stratiform 
   !$OMP THREADPRIVATE(ffallv_lsc)
 
-  REAL, SAVE :: ffallv_con=1.                    ! tuning coefficient crystal fall velocity, convective
+  REAL, SAVE, PROTECTED :: ffallv_con=1.                    ! tuning coefficient crystal fall velocity, convective
   !$OMP THREADPRIVATE(ffallv_con)
 
-  REAL, SAVE :: coef_eva=2e-5                    ! tuning coefficient liquid precip evaporation
+  REAL, SAVE, PROTECTED :: coef_eva=2e-5                    ! tuning coefficient liquid precip evaporation
   !$OMP THREADPRIVATE(coef_eva)
 
-  REAL, SAVE :: coef_eva_i                       ! tuning coefficient ice precip sublimation
+  REAL, SAVE, PROTECTED :: coef_eva_i                       ! tuning coefficient ice precip sublimation
   !$OMP THREADPRIVATE(coef_eva_i)
 
-  REAL cice_velo                                 ! factor in the ice fall velocity formulation
-  PARAMETER (cice_velo=1.645)
-
-  REAL dice_velo                                 ! exponent in the ice fall velocity formulation
-  PARAMETER (dice_velo=0.16)
-
-  REAL, SAVE :: dist_liq=300.                    ! typical deph of cloud-top liquid layer in mpcs
+  REAL, SAVE, PROTECTED :: expo_eva=0.5                     ! tuning coefficient liquid precip evaporation
+  !$OMP THREADPRIVATE(expo_eva)
+
+  REAL, SAVE, PROTECTED :: expo_eva_i                       ! tuning coefficient ice precip sublimation
+  !$OMP THREADPRIVATE(expo_eva_i)
+
+  REAL, SAVE, PROTECTED :: cice_velo=1.645                  ! factor in the ice fall velocity formulation
+  !$OMP THREADPRIVATE(cice_velo)
+
+  REAL, SAVE, PROTECTED ::  dice_velo=0.16                  ! exponent in the ice fall velocity formulation
+  !$OMP THREADPRIVATE(dice_velo)
+
+  REAL, SAVE, PROTECTED :: dist_liq=300.                    ! typical deph of cloud-top liquid layer in mpcs
   !$OMP THREADPRIVATE(dist_liq)
 
-  REAL, SAVE    :: tresh_cl=0.0                  ! cloud fraction threshold for cloud top search
+  REAL, SAVE, PROTECTED    :: tresh_cl=0.0                  ! cloud fraction threshold for cloud top search
   !$OMP THREADPRIVATE(tresh_cl)
 
+  !--Parameters for poprecip
+  LOGICAL, SAVE, PROTECTED :: ok_poprecip=.FALSE.           ! use the processes-oriented formulation of precipitations
+  !$OMP THREADPRIVATE(ok_poprecip)
+
+  REAL, SAVE, PROTECTED :: rain_int_min=0.001               ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction
+  !$OMP THREADPRIVATE(rain_int_min)
+
+  REAL, SAVE, PROTECTED :: thresh_precip_frac=1.E-6         ! precipitation fraction threshold TODO [-]
+  !$OMP THREADPRIVATE(thresh_precip_frac)
+
+  REAL, SAVE, PROTECTED :: gamma_col=1.                     ! A COMMENTER TODO [-]
+  !$OMP THREADPRIVATE(gamma_col)
+
+  REAL, SAVE, PROTECTED :: gamma_agg=1.                     ! A COMMENTER TODO [-]
+  !$OMP THREADPRIVATE(gamma_agg)
+
+  REAL, SAVE, PROTECTED :: gamma_rim=1.                     ! A COMMENTER TODO [-]
+  !$OMP THREADPRIVATE(gamma_rim)
+
+  REAL, SAVE, PROTECTED :: rho_rain=1000.                    ! A COMMENTER TODO [kg/m3]
+  !$OMP THREADPRIVATE(rho_rain)
+
+  REAL, SAVE, PROTECTED :: rho_snow=500.                   ! A COMMENTER TODO [kg/m3]
+  !$OMP THREADPRIVATE(rho_snow)
+
+  REAL, SAVE, PROTECTED :: r_rain=100.E-6                   ! A COMMENTER TODO [m]
+  !$OMP THREADPRIVATE(r_rain)
+
+  REAL, SAVE, PROTECTED :: r_snow=100.E-6                    ! A COMMENTER TODO [m]
+  !$OMP THREADPRIVATE(r_snow)
+
+  REAL, SAVE, PROTECTED :: Eff_rain_liq=1.0                  ! A COMMENTER TODO [-]
+  !$OMP THREADPRIVATE(Eff_rain_liq)
+
+  REAL, SAVE, PROTECTED :: Eff_snow_ice=0.5                ! A COMMENTER TODO [-]
+  !$OMP THREADPRIVATE(Eff_snow_ice)
+
+  REAL, SAVE, PROTECTED :: Eff_snow_liq=1.0              ! A COMMENTER TODO [-]
+  !$OMP THREADPRIVATE(Eff_snow_liq)
+
+  REAL, SAVE, PROTECTED :: tau_auto_snow_min=1800.          ! A COMMENTER TODO [s]
+  !$OMP THREADPRIVATE(tau_auto_snow_min)
+
+  REAL, SAVE, PROTECTED :: tau_auto_snow_max=7200.          ! A COMMENTER TODO [s]
+  !$OMP THREADPRIVATE(tau_auto_snow_max)
+
+  REAL, SAVE, PROTECTED :: eps=1.E-10                       ! A COMMENTER TODO [-]
+  !$OMP THREADPRIVATE(eps)
+  !--End of the parameters for poprecip
+
 ! Two parameters used for lmdz_lscp_old only
-  INTEGER, SAVE :: iflag_oldbug_fisrtilp=0, fl_cor_ebil 
+  INTEGER, SAVE, PROTECTED :: iflag_oldbug_fisrtilp=0, fl_cor_ebil 
   !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp,fl_cor_ebil)
 
@@ -173,5 +226,4 @@
     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
     CALL getin_p('seuil_neb',seuil_neb)
-    CALL getin_p('rain_int_min',rain_int_min) 
     CALL getin_p('iflag_mpc_bl',iflag_mpc_bl)
     CALL getin_p('ok_radocond_snow',ok_radocond_snow) 
@@ -199,9 +251,15 @@
     coef_eva_i=coef_eva
     CALL getin_p('coef_eva_i',coef_eva_i)
+    CALL getin_p('expo_eva',expo_eva)
+    expo_eva_i=expo_eva
+    CALL getin_p('expo_eva_i',expo_eva_i)
     CALL getin_p('iflag_autoconversion',iflag_autoconversion)
     CALL getin_p('dist_liq',dist_liq)
     CALL getin_p('tresh_cl',tresh_cl)
     CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
-
+    CALL getin_p('ok_poprecip',ok_poprecip)
+    CALL getin_p('rain_int_min',rain_int_min) 
+    CALL getin_p('gamma_agg',gamma_agg)
+    CALL getin_p('gamma_col',gamma_col)
 
 
@@ -210,5 +268,4 @@
     WRITE(lunout,*) 'lscp_ini, iflag_evap_prec:', iflag_evap_prec
     WRITE(lunout,*) 'lscp_ini, seuil_neb:', seuil_neb
-    WRITE(lunout,*) 'lscp_ini, rain_int_min:', rain_int_min
     WRITE(lunout,*) 'lscp_ini, iflag_mpc_bl:', iflag_mpc_bl
     WRITE(lunout,*) 'lscp_ini, ok_radocond_snow:', ok_radocond_snow
@@ -235,4 +292,6 @@
     WRITE(lunout,*) 'lscp_ini, coef_eva', coef_eva
     WRITE(lunout,*) 'lscp_ini, coef_eva_i', coef_eva_i
+    WRITE(lunout,*) 'lscp_ini, expo_eva', expo_eva
+    WRITE(lunout,*) 'lscp_ini, expo_eva_i', expo_eva_i
     WRITE(lunout,*) 'lscp_ini, iflag_autoconversion', iflag_autoconversion
     WRITE(lunout,*) 'lscp_ini, dist_liq', dist_liq
@@ -240,4 +299,8 @@
     WRITE(lunout,*) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp
     WRITE(lunout,*) 'lscp_ini, fl_cor_ebil', fl_cor_ebil
+    WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
+    WRITE(lunout,*) 'lscp_ini, rain_int_min:', rain_int_min
+    WRITE(lunout,*) 'lscp_ini, gamma_agg:', gamma_agg
+    WRITE(lunout,*) 'lscp_ini, gamma_col:', gamma_col
 
 
@@ -262,4 +325,15 @@
 
 
+    ! check consistency between the use of the processes-oriented precipitation formulations
+    ! and other options
+   
+    IF (ok_poprecip) THEN
+        IF ((iflag_evap_prec .LT. 4) .OR. (niter_lscp .GT. 1)) THEN
+           abort_message = 'in lscp, ok_poprecip=y requires iflag_evap_prec >= 4 and niter_lscp=1'
+           CALL abort_physic (modname,abort_message,1) 
+        ENDIF 
+    ENDIF
+
+
     !AA Temporary initialisation
     a_tr_sca(1) = -0.5
Index: /LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90	(revision 4803)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90	(revision 4803)
@@ -0,0 +1,552 @@
+MODULE lmdz_lscp_poprecip
+!----------------------------------------------------------------
+! Module for the process-oriented treament of precipitation
+! that are called in LSCP
+! Authors: Atelier Nuage (G. Riviere, L. Raillard, M. Wimmer,
+! N. Dutrievoz, E. Vignon, A. Borella, et al.)
+! Jan. 2024
+
+
+IMPLICIT NONE
+
+CONTAINS
+
+!----------------------------------------------------------------
+! Computes the processes-oriented precipitation formulations for
+! evaporation and sublimation
+! 
+SUBROUTINE poprecip_evapsub( &
+           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
+           qprecip, precipfracclr, precipfraccld, &
+           rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
+           )
+
+USE lmdz_lscp_ini, ONLY : prt_level, lunout
+USE lmdz_lscp_ini, ONLY : coef_eva, coef_eva_i, expo_eva, expo_eva_i, thresh_precip_frac
+USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
+USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
+
+IMPLICIT NONE
+
+
+INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
+REAL,    INTENT(IN)                     :: dtime          !--time step [s]
+LOGICAL, INTENT(IN)                     :: iftop          !--if top of the column
+
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: tempupnew      !--updated temperature of the overlying layer [K]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qprecip        !--specific humidity in the precipitation falling from the upper layer [kg/kg]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain in clear sky coming from the layer above [kg/s/m2]
+REAL,    INTENT(IN),    DIMENSION(klon) :: raincld        !--flux of rain in cloudy air coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow in clear sky coming from the layer above [kg/s/m2]
+REAL,    INTENT(IN),    DIMENSION(klon) :: snowcld        !--flux of snow in cloudy air coming from the layer above [kg/s/m2]
+
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
+
+
+
+
+! integer for interating over klon
+INTEGER :: i
+
+! saturation values
+REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
+! fluxes tendencies because of evaporation
+REAL :: flevapmax, flevapl, flevapi, flevaptot
+! specific humidity tendencies because of evaporation
+REAL :: dqevapl, dqevapi
+! specific heat constant
+REAL :: cpair, cpw
+
+qzero(:)  = 0.0
+dqreva(:) = 0.0
+dqssub(:) = 0.0
+dqevapl=0.0
+dqevapi=0.0
+
+! Calculation of saturation specific humidity
+! depending on temperature:
+CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:))
+! wrt liquid water
+CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))
+! wrt ice
+CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))
+
+
+
+! First step consists in "thermalizing" the layer:
+! as the flux of precip from layer above "advects" some heat (as the precip is at the temperature 
+! of the overlying layer) we recalculate a mean temperature that both the air and the precip in the 
+! layer have.
+
+IF (iftop) THEN
+
+   DO i = 1, klon
+      qprecip(i) = 0.
+   ENDDO
+
+ELSE
+
+   DO i = 1, klon
+       ! no condensed water so cp=cp(vapor+dry air)
+       ! RVTMP2=rcpv/rcpd-1
+       cpair=RCPD*(1.0+RVTMP2*qvap(i))
+       cpw=RCPD*RVTMP2
+       ! qprecip has to be thermalized with 
+       ! layer's air so that precipitation at the ground has the
+       ! same temperature as the lowermost layer
+       qprecip(i) = (rain(i)+snow(i))*dtime/((paprsdn(i)-paprsup(i))/RG)
+       ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
+       temp(i) = ( (tempupnew(i))*qprecip(i)*cpw + cpair*temp(i) ) &
+             / (cpair + qprecip(i)*cpw)
+   ENDDO
+
+ENDIF
+
+
+DO i = 1, klon
+
+  ! if precipitation from the layer above
+  IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
+
+    ! Evaporation of liquid precipitation coming from above
+    ! dP/dz=beta*(1-q/qsat)*(P**expo_eva) (lines 1-2)
+    ! multiplying by dz = - dP / g / rho (line 3-4)
+    ! formula from Sundqvist 1988, Klemp & Wilhemson 1978
+    ! LTP: evaporation only in the clear sky part
+
+    flevapl = precipfracclr(i) * coef_eva * (1.0 - qvap(i) / qsatl(i)) &
+            * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva &
+            * temp(i) * RD / pplay(i) &
+            * ( paprsdn(i) - paprsup(i) ) / RG
+
+    ! evaporation is limited by 0 and by the total water amount in
+    ! the precipitation
+    flevapl = MAX(0.0, MIN(flevapl, rainclr(i)))
+
+
+    ! sublimation of the solid precipitation coming from above
+    ! (same formula as for liquid precip)
+    flevapi = precipfracclr(i) * coef_eva_i * (1.0 - qvap(i) / qsati(i)) &
+            * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva_i &
+            * temp(i) * RD / pplay(i) &
+            * ( paprsdn(i) - paprsup(i) ) / RG
+
+    ! sublimation is limited by 0 and by the total water amount in
+    ! the precipitation
+    ! TODO: change max when we will allow for vapor deposition in supersaturated regions
+    flevapi = MAX(0.0, MIN(flevapi, snowclr(i)))
+
+    ! Evaporation limit: we ensure that the layer's fraction below
+    ! the clear sky does not reach saturation. In this case, we 
+    ! redistribute the maximum flux flevapmax conserving the ratio liquid/ice 
+    ! Max evaporation is computed not to saturate the clear sky precip fraction
+    ! (i.e., the fraction where evaporation occurs)
+    ! It is expressed as a max flux flevapmax
+    !
+    flevapmax = MAX(0.0, ( qsat(i) - qvap(i) ) * precipfracclr(i)) &
+                * ( paprsdn(i) - paprsup(i) ) / RG / dtime
+    flevaptot = flevapl + flevapi
+
+    IF ( flevaptot .GT. flevapmax ) THEN
+        flevapl = flevapmax * flevapl / flevaptot
+        flevapi = flevapmax * flevapi / flevaptot
+    ENDIF
+
+
+    ! New solid and liquid precipitation fluxes after evap and sublimation
+    dqevapl = flevapl / ( paprsdn(i) - paprsup(i) ) * RG * dtime
+    dqevapi = flevapi / ( paprsdn(i) - paprsup(i) ) * RG * dtime
+
+
+    ! vapor is updated after evaporation/sublimation (it is increased)
+    qvap(i) = qvap(i) + dqevapl + dqevapi
+    ! qprecip is the total condensed water in the precip flux (it is decreased)
+    qprecip(i) = qprecip(i) - dqevapl - dqevapi
+    ! air and precip temperature (i.e., gridbox temperature)
+    ! is updated due to latent heat cooling
+    temp(i) = temp(i) &
+            - dqevapl * RLVTT / RCPD &
+            / ( 1.0 + RVTMP2 * ( qvap(i) + qprecip(i) ) ) &
+            - dqevapi * RLSTT / RCPD &
+            / ( 1.0 + RVTMP2 * ( qvap(i) + qprecip(i) ) )
+
+    ! New values of liquid and solid precipitation 
+    rainclr(i) = rainclr(i) - flevapl
+    snowclr(i) = snowclr(i) - flevapi
+
+    ! if there is no more precip fluxes, the precipitation fraction in clear
+    ! sky is set to 0
+    IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0.
+
+    ! calculation of the total fluxes
+    rain(i) = rainclr(i) + raincld(i)
+    snow(i) = snowclr(i) + snowcld(i)
+
+  ELSE
+    ! if no precip, we reinitialize the cloud fraction used for the precip to 0 
+    precipfraccld(i) = 0.
+    precipfracclr(i) = 0.
+
+  ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. )
+
+
+
+  ! write output tendencies for rain and snow
+
+  dqssub(i)  = -dqevapi/dtime
+  dqreva(i)  = -dqevapl/dtime
+
+ENDDO ! loop on klon
+
+
+END SUBROUTINE poprecip_evapsub
+
+!----------------------------------------------------------------
+! Computes the processes-oriented precipitation formulations for
+! - autoconversion (auto) via a deposition process
+! - aggregation (agg)
+! - riming (rim)
+! - collection (coll)
+! - melting (melt)
+! - freezing (free)
+! 
+SUBROUTINE poprecip_postcld( &
+           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
+           temp, qvap, qliq, qice, icefrac, cldfra, &
+           precipfracclr, precipfraccld, &
+           rain, rainclr, raincld, snow, snowclr, snowcld, &
+           dqrauto,dqrcol,dqrmelt,dqrfreez,dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez)
+
+USE lmdz_lscp_ini, ONLY : prt_level, lunout
+USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
+USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
+
+USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb,    &
+                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & 
+                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
+                          rho_rain, rho_snow, r_rain, r_snow, Eff_rain_liq,    &
+                          Eff_snow_ice, Eff_snow_liq, tau_auto_snow_min,       &
+                          tau_auto_snow_max, thresh_precip_frac, eps,          &
+                          iflag_cloudth_vert, iflag_rain_incloud_vol
+
+IMPLICIT NONE
+
+INTEGER, INTENT(IN)                     :: klon           !--number of horizontal grid points [-]
+REAL,    INTENT(IN)                     :: dtime          !--time step [s]
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsdn        !--pressure at the bottom interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: paprsup        !--pressure at the top interface of the layer [Pa]
+REAL,    INTENT(IN),    DIMENSION(klon) :: pplay          !--pressure in the middle of the layer [Pa]
+
+REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--
+LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: temp           !--current temperature [K]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qvap           !--current water vapor specific humidity [kg/kg]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--current liquid water specific humidity [kg/kg]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--current ice water specific humidity [kg/kg]
+REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--
+REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
+                                                          !--NB. at the end of the routine, becomes the fraction of precip
+                                                          !--in the current layer
+
+REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: rainclr        !--flux of rain in clear sky coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: raincld        !--flux of rain in cloudy air coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: snow           !--flux of snow coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: snowclr        !--flux of snow in clear sky coming from the layer above [kg/s/m2]
+REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow in cloudy air coming from the layer above [kg/s/m2]
+
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrcol         !-- rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsagg         !-- snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !-- rain tendency due to autoconversion of cloud liquid [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !-- snow tendency due to autoconversion of cloud ice [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsrim         !-- snow tendency due to riming [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsmelt        !-- snow tendency due to melting [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrmelt        !-- rain tendency due to melting [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !-- snow tendency due to freezing [kg/kg/s]
+REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !-- rain tendency due to freezing [kg/kg/s]
+
+
+
+!--Local variables
+
+INTEGER :: i
+
+REAL :: hum_to_flux
+REAL :: dcldfra
+REAL :: precipfractot
+REAL :: dprecipfracclr, dprecipfraccld
+REAL :: drainclr, dsnowclr
+REAL :: draincld, dsnowcld
+REAL :: eff_cldfra
+REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain
+REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow
+REAL ::  dqlcol           ! loss of liquid cloud content due to collection by rain [kg/kg/s]
+REAL ::  dqiagg           ! loss of ice cloud content due to collection by aggregation [kg/kg/s]
+REAL ::  dqlauto          ! loss of liquid cloud content due to autoconversion to rain [kg/kg/s]
+REAL ::  dqiauto          ! loss of ice cloud content due to autoconversion to snow [kg/kg/s]
+REAL ::  dqlrim           ! loss of liquid cloud content due to riming on snow[kg/kg/s]
+
+
+!--Initialisation of variables
+
+
+dqlrim=0.0
+dqlcol=0.0
+dqiagg=0.0
+dqiauto=0.0
+dqlauto=0.0
+
+
+dqrcol(:) = 0.
+dqsagg(:) = 0.
+dqsauto(:) = 0.
+dqrauto(:) = 0.
+dqsrim(:) = 0.
+dqrmelt(:) = 0.
+dqsmelt(:) = 0.
+dqrfreez(:) = 0.
+dqsfreez(:) = 0.
+
+
+
+!Partitionning between precipitation coming from clouds and that coming from CS
+
+!0) Calculate tot_zneb, total cloud fraction above the cloud
+!assuming a maximum-random overlap (voir Jakob and Klein, 2000)
+
+DO i = 1, klon
+
+  !--Initialisation
+  drainclr = 0.
+  dsnowclr = 0.
+  hum_to_flux = ( paprsdn(i) - paprsup(i) ) / RG / dtime
+  precipfractot = precipfracclr(i) + precipfraccld(i)
+
+  ! formule maximum random
+  precipfractot = 1. - ( 1. - precipfractot ) * &
+                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
+               / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
+
+
+  ! precipfraccld est égal a rneb au niveau d'au dessus
+  dcldfra = cldfra(i) - precipfraccld(i)
+  ! add max to guaranttee precip fraction >=0
+  dprecipfracclr = MAX(0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i)
+  dprecipfraccld = MAX( dcldfra , -precipfraccld(i) ) 
+
+
+  IF ( dprecipfraccld .GT. 0. ) THEN
+    IF ( precipfracclr(i) .GT. 0. ) THEN
+      drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i) 
+      dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i) 
+    ENDIF
+  ELSE
+    IF  ( precipfraccld(i) .GT. 0. ) THEN
+      draincld = dprecipfraccld / precipfraccld(i) * raincld(i) 
+      dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i) 
+      drainclr = - draincld
+      dsnowclr = - dsnowcld
+    ENDIF
+  ENDIF
+
+
+  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
+  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
+  rainclr(i) = rainclr(i) + drainclr
+  snowclr(i) = snowclr(i) + dsnowclr
+  raincld(i) = raincld(i) - drainclr
+  snowcld(i) = snowcld(i) - dsnowclr
+  
+  
+  IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN
+    eff_cldfra = ctot_vol(i)
+  ELSE
+    eff_cldfra = cldfra(i)
+  ENDIF 
+
+  IF ( cldfra(i) .GE. seuil_neb ) THEN
+
+    ! if vertical heterogeneity is taken into account, we use
+    ! the "true" volume fraction instead of a modified 
+    ! surface fraction (which is larger and artificially
+    ! reduces the in-cloud water).
+
+    ! agg et col avant autoconv car on ne veut pas agg et col la nouvelle precip
+    !---------------------------------------------------------
+    !--           COLLECTION AND AGGREGATION
+    !---------------------------------------------------------
+
+    ! remplacer la premiere ligne par "coef_col" ?
+    dqlcol = - gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq &
+              * qliq(i) / eff_cldfra * raincld(i) * dtime
+
+    ! remplacer la premiere ligne par "coef_agg" ?
+    dqiagg = - gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice &
+              * qice(i) / eff_cldfra * snowcld(i) * dtime
+
+    dqiagg = MAX( - qice(i), dqiagg )
+    dqlcol = MAX( - qliq(i), dqlcol )
+
+    qliq(i) = qliq(i) + dqlcol
+    qice(i) = qice(i) + dqiagg
+
+    raincld(i) = raincld(i) - dqlcol * hum_to_flux
+    snowcld(i) = snowcld(i) - dqiagg * hum_to_flux
+
+
+    !---------------------------------------------------------
+    !--                  AUTOCONVERSION
+    !---------------------------------------------------------
+
+    ! TODO
+    IF ( ptconv(i) ) THEN ! if convective point
+        qthresh_auto_rain = cld_lc_con
+        qthresh_auto_snow = cld_lc_con
+
+        tau_auto_rain = cld_tau_con
+        tau_auto_snow = tau_auto_snow_max &
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+
+        expo_auto_rain = cld_expo_con
+        expo_auto_snow = cld_expo_con
+    ELSE
+        qthresh_auto_rain = cld_lc_lsc
+        qthresh_auto_snow = cld_lc_lsc
+
+        tau_auto_rain = cld_tau_lsc
+        tau_auto_snow = tau_auto_snow_max &
+                      + ( tau_auto_snow_min - tau_auto_snow_max ) * ( 1. - icefrac(i) )
+
+        expo_auto_rain = cld_expo_lsc
+        expo_auto_snow = cld_expo_lsc
+    ENDIF
+
+
+    ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
+    ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
+    !.........................................................
+    ! we first treat the second term (with exponential) in an explicit way
+    ! and then treat the first term (-q/tau) in an exact way
+
+    dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( &
+               - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) )
+
+    dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( &
+               - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) )
+
+
+    dqlauto = MAX( - qliq(i), dqlauto )
+    dqiauto = MAX( - qice(i), dqiauto )
+
+    qliq(i) = qliq(i) + dqlauto
+    qice(i) = qice(i) + dqiauto
+
+    raincld(i) = raincld(i) - dqlauto * hum_to_flux
+    snowcld(i) = snowcld(i) - dqiauto * hum_to_flux
+
+
+    ! FOLLOWING PROCESSES IMPLY A PHASE CHANGE SO A TEMPERATURE 
+    ! ADJUSTMENT 
+
+    !---------------------------------------------------------
+    !--                  RIMING
+    !---------------------------------------------------------
+
+    ! remplacer la premiere ligne par "coef_rim" ?
+    dqlrim = - gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq &
+              * qliq(i) / eff_cldfra * snowcld(i) * dtime
+
+    dqlrim = MAX( - qliq(i), dqlrim )
+
+    qliq(i) = qliq(i) + dqlrim
+
+    snowcld(i) = snowcld(i) - dqlrim * hum_to_flux
+
+
+  ENDIF ! rneb .GE. seuil_neb
+
+
+
+  !---------------------------------------------------------
+  !--                  FREEZING
+  !---------------------------------------------------------
+
+  !dqrfree_max = MINOUMAX(0., ( RTT - temp(i) ) / RLMLT * RCPD / ( 1. + RVTMP2 * ( qtot(i) + qprecip(i) ) ))
+
+  !---------------------------------------------------------
+  !--                  MELTING
+  !---------------------------------------------------------
+
+  !flux = velocity * N_0 * 4. / 3. * PI * r_snow**3. * rho_snow
+
+  !IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN
+  !  dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD / ( 1. + RVTMP2 * ( qtot(i) + qprecip(i) ) ))
+  !  dsnowtotmelt_max = dqsmelt_max * hum_to_flux(i)
+  !  
+  !  dsnowtotmelt = - nb_snowflake * 4. * PI * mol_diff_vap * snowflake_capa / RLMLT * coef_ventil &
+  !              * MAX(0., ticebulb - RTT) &
+  !              * ( paprsdn(i) - paprsup(i) ) / RG
+  !  ! max bec. negative values
+  !  dsnowtotmelt = MAX(dsnowtotmelt, dsnowmelt_max)
+  !  dsnowclrmelt = dsnowtotmelt * snowclr(i) / ( snowclr(i) + snowcld(i) )
+  !  dsnowcldmelt = dsnowtotmelt - dsnowclrmelt
+
+
+  !  ! update of rainfall and snowfall due to melting
+  !  rainclr(i) = rainclr(i) - dsnowclrmelt(i)
+  !  raincld(i) = raincld(i) - dsnowcldmelt(i)
+  !  snowclr(i) = snowclr(i) + dsnowclrmelt(i)
+  !  snowcld(i) = snowcld(i) + dsnowcldmelt(i)
+  !ENDIF
+  !
+  !! Latent heat of melting with precipitation thermalization
+  !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)))
+
+
+  !! MISE A JOUR DES FRACTIONS PRECIP CLD et CS
+  ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
+
+  precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / rain_int_min )
+  precipfraccld(i) = MIN( precipfraccld(i), ( raincld(i) + snowcld(i) ) / rain_int_min )
+
+  rain(i) = rainclr(i) + raincld(i)
+  snow(i) = snowclr(i) + snowcld(i)
+
+! write output tendencies for rain and snow
+
+  dqsrim(i)  = -dqlrim/dtime
+  dqrcol(i)  = -dqlcol/dtime
+  dqsagg(i)  = -dqiagg/dtime
+  dqsauto(i) = -dqiauto/dtime
+  dqrauto(i) = -dqlauto/dtime
+
+
+
+ENDDO
+
+
+
+END SUBROUTINE poprecip_postcld
+
+END MODULE lmdz_lscp_poprecip
Index: /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4802)
+++ /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4803)
@@ -538,4 +538,32 @@
 
 
+!--POPRECIP variables
+      REAL, SAVE, ALLOCATABLE :: dqreva(:,:)
+      !$OMP THREADPRIVATE(dqreva)
+      REAL, SAVE, ALLOCATABLE :: dqrauto(:,:)
+      !$OMP THREADPRIVATE(dqrauto)
+      REAL, SAVE, ALLOCATABLE :: dqrcol(:,:)
+      !$OMP THREADPRIVATE(dqrcol)
+      REAL, SAVE, ALLOCATABLE :: dqrmelt(:,:)
+      !$OMP THREADPRIVATE(dqrmelt)
+       REAL, SAVE, ALLOCATABLE :: dqrfreez(:,:)
+      !$OMP THREADPRIVATE(dqrfreez)
+      REAL, SAVE, ALLOCATABLE :: dqssub(:,:)
+      !$OMP THREADPRIVATE(dqssub)
+      REAL, SAVE, ALLOCATABLE :: dqsauto(:,:)
+      !$OMP THREADPRIVATE(dqsauto)
+      REAL, SAVE, ALLOCATABLE :: dqsagg(:,:)
+      !$OMP THREADPRIVATE(dqsagg)
+      REAL, SAVE, ALLOCATABLE :: dqsrim(:,:)
+      !$OMP THREADPRIVATE(dqsrim)
+      REAL, SAVE, ALLOCATABLE :: dqsmelt(:,:)
+      !$OMP THREADPRIVATE(dqsmelt)
+      REAL, SAVE, ALLOCATABLE :: dqsfreez(:,:)
+      !$OMP THREADPRIVATE(dqsfreez)
+
+
+      
+
+
 #ifdef CPP_StratAer
 !
@@ -928,4 +956,9 @@
       ALLOCATE(zqsatl(klon,klev), zqsats(klon,klev))
       ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev), fcontrN(klon,klev), fcontrP(klon,klev))
+
+!--POPRECIP variables
+      ALLOCATE(dqreva(klon,klev),dqssub(klon,klev))
+      ALLOCATE(dqrauto(klon,klev), dqrcol(klon,klev), dqrmelt(klon,klev), dqrfreez(klon,klev))
+      ALLOCATE(dqsauto(klon,klev), dqsagg(klon,klev), dqsrim(klon,klev), dqsmelt(klon,klev), dqsfreez(klon,klev))
 
 #ifdef CPP_StratAer
@@ -1239,4 +1272,9 @@
       DEALLOCATE(Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
 
+!--POPRECIP variables
+      DEALLOCATE(dqreva,dqssub)
+      DEALLOCATE(dqrauto,dqrcol,dqrmelt,dqrfreez)
+      DEALLOCATE(dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez)
+
 #ifdef CPP_StratAer
 ! variables for strat. aerosol CK
Index: /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 4802)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 4803)
@@ -1556,4 +1556,26 @@
   TYPE(ctrl_out), SAVE :: o_pfracld = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'pfracld', 'LS precipitation fraction cloudy part', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqreva = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqreva', 'LS rain tendency due to evaporation', 'kg/m2/s', (/ ('', i=1, 10) /))
+   TYPE(ctrl_out), SAVE :: o_dqrauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqrauto', 'LS rain tendency due to autoconversion', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqrcol = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqrcol', 'LS rain tendency due to collection', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqrmelt = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqrmelt', 'LS rain tendency due to melting', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqrfreez = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqrfreez', 'LS rain tendency due to freezing', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqssub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqssub', 'LS snow tendency due to sublimation', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsauto', 'LS snow tendency due to autoconversion', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsagg = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsagg', 'LS snow tendency due to aggragation', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsrim = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsrim', 'LS snow tendency due to riming', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsmelt = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsmelt', 'LS snow tendency due to melting', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsfreez = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsfreez', 'LS snow tendency due to freezing', 'kg/m2/s', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /))
Index: /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4802)
+++ /LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 4803)
@@ -141,4 +141,6 @@
          o_rnebls, o_rneblsvol, o_rhum, o_rhl, o_rhi, o_ozone, o_ozone_light, &
          o_pfraclr, o_pfracld, &
+         o_dqreva, o_dqrauto, o_dqrcol, o_dqrmelt, o_dqrfreez, &
+         o_dqssub, o_dqsauto, o_dqsagg, o_dqsrim, o_dqsmelt, o_dqsfreez, &
          o_duphy, o_dtphy, o_dqphy, o_dqphy2d, o_dqlphy, o_dqlphy2d, &
          o_dqsphy, o_dqsphy2d, o_dqbsphy, o_dqbsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
@@ -255,5 +257,6 @@
 
     USE ice_sursat_mod, ONLY: flight_m, flight_h2o
-    
+    USE lmdz_lscp_ini, ONLY: ok_poprecip
+
     USE phys_output_ctrlout_mod, ONLY: o_heat_volc, o_cool_volc !NL
     USE phys_state_var_mod, ONLY: heat_volc, cool_volc !NL
@@ -357,5 +360,9 @@
          zphi, u_seri, v_seri, omega, cldfra, &
          rneb, rnebjn, rneblsvol, zx_rh, zx_rhl, zx_rhi, &
-         pfraclr, pfracld, d_t_dyn,  & 
+         pfraclr, pfracld,  &
+         dqreva, dqssub, &
+         dqrauto,dqrcol,dqrmelt,dqrfreez, &
+         dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
+         d_t_dyn,  & 
          d_q_dyn,  d_ql_dyn, d_qs_dyn, d_qbs_dyn,  &
          d_q_dyn2d,  d_ql_dyn2d, d_qs_dyn2d, d_qbs_dyn2d, &
@@ -1942,4 +1949,17 @@
            CALL histwrite_phy(o_pfraclr, pfraclr)
            CALL histwrite_phy(o_pfracld, pfracld)
+           IF (ok_poprecip) THEN
+           CALL histwrite_phy(o_dqreva, dqreva)
+           CALL histwrite_phy(o_dqrauto, dqrauto)
+           CALL histwrite_phy(o_dqrcol, dqrcol)
+           CALL histwrite_phy(o_dqrmelt, dqrmelt)
+           CALL histwrite_phy(o_dqrfreez, dqrfreez)
+           CALL histwrite_phy(o_dqssub, dqssub)
+           CALL histwrite_phy(o_dqsauto, dqsauto)
+           CALL histwrite_phy(o_dqsagg, dqsagg)
+           CALL histwrite_phy(o_dqsmelt, dqsmelt)
+           CALL histwrite_phy(o_dqsfreez, dqsfreez)
+           CALL histwrite_phy(o_dqsrim, dqsrim)
+           ENDIF
        ENDIF
 
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4802)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4803)
@@ -195,4 +195,8 @@
        east_gwstress,west_gwstress, &
        d_q_ch4, &
+       ! proprecip
+       dqreva, dqssub, &
+       dqrauto,dqrcol,dqrmelt,dqrfreez, &
+       dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
        !  Special RRTM
        ZLWFT0_i,ZSWFT0_i,ZFLDN0,  &
@@ -3887,5 +3891,6 @@
          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)
 
 
Index: /LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_poprecip.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_poprecip.F90	(revision 4803)
+++ /LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_poprecip.F90	(revision 4803)
@@ -0,0 +1,1 @@
+link ../phylmd/lmdz_lscp_poprecip.F90
Index: /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4802)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4803)
@@ -639,4 +639,32 @@
       !$OMP THREADPRIVATE(fcontrP)
 
+!--POPRECIP variables
+      REAL, SAVE, ALLOCATABLE :: dqreva(:,:)
+      !$OMP THREADPRIVATE(dqreva)
+      REAL, SAVE, ALLOCATABLE :: dqrauto(:,:)
+      !$OMP THREADPRIVATE(dqrauto)
+      REAL, SAVE, ALLOCATABLE :: dqrcol(:,:)
+      !$OMP THREADPRIVATE(dqrcol)
+      REAL, SAVE, ALLOCATABLE :: dqrmelt(:,:)
+      !$OMP THREADPRIVATE(dqrmelt)
+       REAL, SAVE, ALLOCATABLE :: dqrfreez(:,:)
+      !$OMP THREADPRIVATE(dqrfreez)
+      REAL, SAVE, ALLOCATABLE :: dqssub(:,:)
+      !$OMP THREADPRIVATE(dqssub)
+      REAL, SAVE, ALLOCATABLE :: dqsauto(:,:)
+      !$OMP THREADPRIVATE(dqsauto)
+      REAL, SAVE, ALLOCATABLE :: dqsagg(:,:)
+      !$OMP THREADPRIVATE(dqsagg)
+      REAL, SAVE, ALLOCATABLE :: dqsrim(:,:)
+      !$OMP THREADPRIVATE(dqsrim)
+      REAL, SAVE, ALLOCATABLE :: dqsmelt(:,:)
+      !$OMP THREADPRIVATE(dqsmelt)
+      REAL, SAVE, ALLOCATABLE :: dqsfreez(:,:)
+      !$OMP THREADPRIVATE(dqsfreez)
+
+
+
+
+
 #ifdef CPP_StratAer
 !
@@ -1085,4 +1113,9 @@
       ALLOCATE(zqsatl(klon,klev), zqsats(klon,klev))
       ALLOCATE(Tcontr(klon,klev), qcontr(klon,klev), qcontr2(klon,klev), fcontrN(klon,klev), fcontrP(klon,klev))
+
+!--POPRECIP variables
+      ALLOCATE(dqreva(klon,klev),dqssub(klon,klev))
+      ALLOCATE(dqrauto(klon,klev), dqrcol(klon,klev), dqrmelt(klon,klev), dqrfreez(klon,klev))
+      ALLOCATE(dqsauto(klon,klev), dqsagg(klon,klev), dqsrim(klon,klev), dqsmelt(klon,klev), dqsfreez(klon,klev))
 
 #ifdef CPP_StratAer
@@ -1444,4 +1477,9 @@
       DEALLOCATE(Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
 
+!--POPRECIP variables
+      DEALLOCATE(dqreva,dqssub)
+      DEALLOCATE(dqrauto,dqrcol,dqrmelt,dqrfreez)
+      DEALLOCATE(dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez)
+
 #ifdef CPP_StratAer
 ! variables for strat. aerosol CK
Index: /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 4802)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90	(revision 4803)
@@ -1477,4 +1477,26 @@
   TYPE(ctrl_out), SAVE :: o_pfracld = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     'pfracld', 'LS precipitation fraction cloudy part', '-', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqreva = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqreva', 'LS rain tendency due to evaporation', 'kg/m2/s', (/ ('', i=1, 10) /))
+   TYPE(ctrl_out), SAVE :: o_dqrauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqrauto', 'LS rain tendency due to autoconversion', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqrcol = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqrcol', 'LS rain tendency due to collection', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqrmelt = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqrmelt', 'LS rain tendency due to melting', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqrfreez = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqrfreez', 'LS rain tendency due to freezing', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqssub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqssub', 'LS snow tendency due to sublimation', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsauto', 'LS snow tendency due to autoconversion', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsagg = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsagg', 'LS snow tendency due to aggragation', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsrim = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsrim', 'LS snow tendency due to riming', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsmelt = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsmelt', 'LS snow tendency due to melting', 'kg/m2/s', (/ ('', i=1, 10) /))
+  TYPE(ctrl_out), SAVE :: o_dqsfreez = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
+    'dqsfreez', 'LS snow tendency due to freezing', 'kg/m2/s', (/ ('', i=1, 10) /))
   TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
     'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /))
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4802)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4803)
@@ -229,4 +229,8 @@
        east_gwstress,west_gwstress, &
        d_q_ch4, &
+       ! proprecip
+       dqreva, dqssub, &
+       dqrauto,dqrcol,dqrmelt,dqrfreez, &
+       dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
        !  Special RRTM
        ZLWFT0_i,ZSWFT0_i,ZFLDN0,  &
@@ -4876,5 +4880,5 @@
     CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
          t_seri, q_seri,ptconv,ratqs, &
-         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
+         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, & 
          pfraclr,pfracld, &
          radocond, picefra, rain_lsc, snow_lsc, &
@@ -4885,5 +4889,6 @@
          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)
 
 
