Index: LMDZ6/trunk/libf/phylmd/lscp_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4397)
+++ LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4412)
@@ -9,5 +9,5 @@
      paprs,pplay,t,q,ptconv,ratqs,                      &
      d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri,  &
-     radliq, radicefrac, rain, snow,                    &
+     radocond, radicefrac, rain, snow,                    &
      frac_impa, frac_nucl, beta,                        &
      prfl, psfl, rhcl, zqta, fraca,                     &
@@ -41,7 +41,10 @@
 ! References: 
 !
+! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
 ! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
 ! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
+! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
 ! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
+! - Touzze-Peifert Ludo, PhD thesis, p117-124
 ! -------------------------------------------------------------------------------
 ! Code structure:
@@ -50,11 +53,9 @@
 ! P1>     Evaporation of the precipitation (falling from the k+1 level)
 ! P2>     Cloud formation (at the k level) 
-! P2.A.0> Cloud properties calculation from a rectangular pdf
-! P2.A.1> With the new PDFs, calculation of cloud properties using the inital
+! P2.A.1> With the PDFs, calculation of cloud properties using the inital
 !         values of T and Q
 ! P2.A.2> Coupling between condensed water and temperature
 ! P2.A.3> Calculation of final quantities associated with cloud formation
-! P2.B>   'All or nothing' cloud
-! P2.C>   Release of Latent heat after cloud formation
+! P2.B>   Release of Latent heat after cloud formation
 ! P3>     Autoconversion to precipitation (k-level)
 ! P4>     Wet scavenging
@@ -65,5 +66,5 @@
 ! water used in the physics and the dynamics 
 ! 
-! During the autoconversion to precipitation (P3 step), radliq (cloud water used
+! During the autoconversion to precipitation (P3 step), radocond (cloud water used
 ! by the radiation scheme) is calculated as an average of the water that remains
 ! in the cloud during the precipitation and not the water remaining at the end
@@ -75,11 +76,12 @@
 ! Radiation:
 ! xflwc(newmicro)+xfiwc(newmicro) =
-!   cldliq(physiq)=radliq(fisrt)=lwcon(nc)+iwcon(nc)
+!   radocond=lwcon(nc)+iwcon(nc)
 !
 ! Notetheless, be aware of:
 !
-! radliq(fisrt) .NE. ocond(nc)
+! radocond .NE. ocond(nc)
 ! i.e.:
 ! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
+! but oliq+(ocond-oliq) .EQ. ocond
 ! (which is not trivial)
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -91,5 +93,5 @@
 
 USE lscp_ini_mod, ONLY : seuil_neb, ninter, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
-USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radliq_snow, a_tr_sca
+USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca
 
 
@@ -115,5 +117,5 @@
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t               ! temperature (K)
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q               ! specific humidity [kg/kg] 
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q               ! total specific humidity (= vapor phase) [kg/kg] 
   INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
   INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
@@ -152,9 +154,9 @@
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-]  
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq           ! condensed water used in the radiation scheme [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
-  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! large-scale rainfall [kg/s/m2] 
-  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! large-scale snowfall [kg/s/m2]
+  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2] 
+  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl            ! saturation specific humidity wrt liquid [kg/kg]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsats            ! saturation specific humidity wrt ice [kg/kg]  
@@ -231,5 +233,5 @@
   REAL velo(klon,klev), vr
   REAL qlmpc, qimpc, rnebmpc
-  REAL radliqi(klon,klev), radliql(klon,klev)
+  REAL radocondi(klon,klev), radocondl(klon,klev)
 
   INTEGER i, k, n, kk
@@ -295,5 +297,5 @@
 d_qi(:,:) = 0.0
 rneb(:,:) = 0.0
-radliq(:,:) = 0.0
+radocond(:,:) = 0.0
 radicefrac(:,:) = 0.0
 frac_nucl(:,:) = 1.0 
@@ -350,11 +352,15 @@
     ! P0> Thermalization of precipitation falling from the overlying layer
     ! --------------------------------------------------------------------
-    ! Computes air temperature variation due to latent heat transported by
+    ! Computes air temperature variation due to enthalpy transported by
     ! precipitation. Precipitation is then thermalized with the air in the
     ! layer. 
     ! The precipitation should remain thermalized throughout the different
-    ! thermodynamical transformations. The corresponding water mass should
+    ! thermodynamical transformations. 
+    ! The corresponding water mass should
     ! be added when calculating the layer's enthalpy change with 
     ! temperature
+    ! See lmdzpedia page todoan
+    ! todoan: check consistency with ice phase
+    ! todoan: understand why several steps
     ! ---------------------------------------------------------------------
     
@@ -409,5 +415,5 @@
             IF (zrfl(i)+zifl(i).GT.0.) THEN
 
-            ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
+            ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
             ! c_iso: likely important to distinguish cs from neb isotope precipitation
 
@@ -475,4 +481,5 @@
                 ! 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)
@@ -494,4 +501,5 @@
                 zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
                 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+
                 ! precip thermalization
                 zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
@@ -525,16 +533,25 @@
 
                 ! Melting:
-                zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! JYG
+                zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
                 ! precip fraction that is melted
                 zmelt = MIN(MAX(zmelt,0.),1.)
 
-                ! Ice melting
+                ! update of rainfall and snowfall due to melting
                 IF (iflag_evap_prec.EQ.4) THEN
                     zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
                     zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
                     zrfl(i)=zrflclr(i)+zrflcld(i)
+
+                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
+                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
+                    zifl(i)=ziflclr(i)+ziflcld(i)
+
                 ELSE
                     zrfl(i)=zrfl(i)+zmelt*zifl(i)
-                ENDIF
+
+                    zifl(i)=zifl(i)*(1.-zmelt)
+                ENDIF
+
+
                 ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
 
@@ -543,12 +560,4 @@
                 *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
                 
-                IF (iflag_evap_prec.EQ.4) THEN
-                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
-                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
-                    zifl(i)=ziflclr(i)+ziflcld(i)
-                ELSE
-                    zifl(i)=zifl(i)*(1.-zmelt)
-                ENDIF
-                !c_iso: same for isotopic precip.
 
             ELSE
@@ -908,5 +917,5 @@
     ENDIF
 
-    ! Initialisation of zoliq and radliq variables
+    ! Initialisation of zoliq and radocond variables
 
     DO i = 1, klon
@@ -919,8 +928,8 @@
             iwc(i)   = 0. 
             zneb(i)  = MAX(rneb(i,k),seuil_neb)
-            radliq(i,k) = zoliq(i)/REAL(ninter+1)
+            radocond(i,k) = zoliq(i)/REAL(ninter+1)
             radicefrac(i,k) = zfice(i)
-            radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
-            radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
+            radocondi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
+            radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
     ENDDO
 
@@ -992,7 +1001,7 @@
                 ! c_iso: call isotope_conversion (for readibility) or duplicate 
 
-                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)
+                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(ninter+1)
+                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(ninter+1)
+                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(ninter+1)
 
             ENDIF ! rneb >0
@@ -1109,11 +1118,11 @@
 
     ! Calculation of the concentration of condensates seen by the radiation scheme
-    ! for liquid phase, we take radliql
-    ! for ice phase, we take radliqi if we neglect snowfall, otherwise (ok_radliq_snow=true) 
-    ! we recaulate radliqi to account for contributions from the precipitation flux
-
-    IF ((ok_radliq_snow) .AND. (k .LT. klev)) THEN
+    ! 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
         ! for the solid phase (crystals + snowflakes)
-        ! we recalculate radliqi by summing
+        ! we recalculate radocondi by summing
         ! the ice content calculated in the mesh
         ! + the contribution of the non-evaporated snowfall 
@@ -1121,16 +1130,16 @@
         DO i=1,klon
            IF (ziflprev(i) .NE. 0.0) THEN
-              radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
+              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
            ELSE
-              radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
+              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
            ENDIF
-           radliq(i,k)=radliql(i,k)+radliqi(i,k)
+           radocond(i,k)=radocondl(i,k)+radocondi(i,k)
         ENDDO
     ENDIF
 
-    ! caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
+    ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
     DO i=1,klon
-        IF (radliq(i,k) .GT. 0.) THEN
-            radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
+        IF (radocond(i,k) .GT. 0.) THEN
+            radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
         ENDIF
     ENDDO
