Index: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90
===================================================================
--- LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90	(revision 5163)
+++ LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90	(revision 5202)
@@ -7,13 +7,16 @@
 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
-     paprs,pplay,temp,qt,ptconv,ratqs,                  &
+     paprs,pplay,temp,qt,qice_save,ptconv,ratqs,        &
      d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
-     pfraclr,pfracld,                                   &
+     pfraclr, pfracld,                                  &
+     cldfraliq, sigma2_icefracturb,mean_icefracturb,    &
      radocond, radicefrac, rain, snow,                  &
      frac_impa, frac_nucl, beta,                        &
-     prfl, psfl, rhcl, qta, fraca,                     &
-     tv, pspsk, tla, thl, iflag_cld_th,             &
-     iflag_ice_thermo, distcltop, temp_cltop, cell_area,&
-     cf_seri, rvc_seri, u_seri, v_seri, pbl_eps,        &
+     prfl, psfl, rhcl, qta, fraca,                      &
+     tv, pspsk, tla, thl, iflag_cld_th,                 &
+     iflag_ice_thermo, distcltop, temp_cltop,           &
+     tke, tke_dissip,                                   &
+     cell_area,                                         &
+     cf_seri, rvc_seri, u_seri, v_seri,                 &
      qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
      ratio_qi_qtot, dcf_sub, dcf_con, dcf_mix,          &
@@ -100,5 +103,6 @@
 ! USE de modules contenant des fonctions.
 USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc 
-USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat
+USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
+USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb
 USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
 USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
@@ -115,5 +119,5 @@
 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
 USE lmdz_lscp_ini, ONLY : ok_poprecip
-USE lmdz_lscp_ini, ONLY : ok_external_lognormal, ok_ice_supersat, ok_unadjusted_clouds
+USE lmdz_lscp_ini, ONLY : ok_external_lognormal, ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
 
 IMPLICIT NONE
@@ -134,4 +138,5 @@
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg] 
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qice_save       ! ice specific from previous time step [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
@@ -141,15 +146,17 @@
   !Inputs associated with thermal plumes
 
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv             ! virtual potential temperature [K]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta            ! specific humidity within thermals [kg/kg]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca          ! fraction of thermals within the mesh [-]
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk          ! exner potential (p/100000)**(R/cp) 
-  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla            ! liquid temperature within thermals [K]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp) 
+  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
+  REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
+  REAL, DIMENSION(klon,klev+1),      INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
 
   ! INPUT/OUTPUT variables
   !------------------------
   
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
-  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs        ! function of pressure that sets the large-scale
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl              ! liquid potential temperature [K]
+  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs            ! function of pressure that sets the large-scale
 
   ! INPUT/OUTPUT condensation and ice supersaturation
@@ -160,5 +167,4 @@
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
-  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: pbl_eps          ! TKE dissipation [?]
   REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
 
@@ -179,4 +185,7 @@
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
   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
@@ -190,8 +199,8 @@
   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
 
-  ! fraction of aerosol scavenging through impaction and nucleation (for on-line)
+  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
   
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa        ! scavenging fraction due tu impaction [-]
-  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl        ! scavenging fraction due tu nucleation [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
+  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
   
   ! for condensation and ice supersaturation
@@ -255,6 +264,5 @@
   ! LOCAL VARIABLES:
   !----------------
-
-  REAL,DIMENSION(klon) :: qsl, qsi
+  REAL,DIMENSION(klon) :: qsl, qsi                                ! saturation threshold at current vertical level
   REAL :: zct, zcl,zexpo
   REAL, DIMENSION(klon,klev) :: ctot
@@ -263,10 +271,10 @@
   REAL :: zdelta, zcor, zcvm5
   REAL, DIMENSION(klon) :: zdqsdT_raw
-  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                ! coefficient to make cold condensation at the correct RH and derivative wrt T
-  REAL, DIMENSION(klon) :: Tbef,qlbef,DT
+  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
+  REAL, DIMENSION(klon) :: Tbef,qlbef,DT                          ! temperature, humidity and temp. variation during lognormal iteration
   REAL :: num,denom
   REAL :: cste
-  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta
-  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2
+  REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta             ! lognormal parameters
+  REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2          ! lognormal intermediate variables
   REAL :: erf
   REAL, DIMENSION(klon) :: zfice_th
@@ -285,4 +293,5 @@
   REAL :: zmelt,zrain,zsnow,zprecip
   REAL, DIMENSION(klon) :: dzfice
+  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
   REAL :: zsolid
   REAL, DIMENSION(klon) :: qtot, qzero
@@ -315,5 +324,6 @@
   REAL, DIMENSION(klon,klev) :: radocondi, radocondl
   REAL :: effective_zneb
-  REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D
+  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
+  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
   
   ! for condensation and ice supersaturation
@@ -328,5 +338,4 @@
   REAL :: min_qParent, min_ratio
 
-
   INTEGER i, k, n, kk, iter
   INTEGER, DIMENSION(klon) :: n_i
@@ -382,4 +391,7 @@
 pfraclr(:,:)=0.0
 pfracld(:,:)=0.0
+cldfraliq(:,:)=0.
+sigma2_icefracturb(:,:)=0.
+mean_icefracturb(:,:)=0.
 radocond(:,:) = 0.0
 radicefrac(:,:) = 0.0
@@ -391,4 +403,6 @@
 zfice(:)=0.0
 dzfice(:)=0.0
+zfice_turb(:)=0.0
+dzfice_turb(:)=0.0
 zqprecl(:)=0.0
 zqpreci(:)=0.0
@@ -405,6 +419,6 @@
 d_tot_zneb(:) = 0.0
 qzero(:) = 0.0
-distcltop1D(:)=0.0
-temp_cltop1D(:) = 0.0
+zdistcltop(:)=0.0
+ztemp_cltop(:) = 0.0
 ztupnew(:)=0.0
 
@@ -459,5 +473,4 @@
 
 
-
 !c_iso: variable initialisation for iso 
 
@@ -478,4 +491,9 @@
 
     ! Initialisation temperature and specific humidity
+    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
+    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
+    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated 
+    ! d_t = temperature tendency due to lscp
+    ! The temperature of the overlying layer is updated here because needed for thermalization
     DO i = 1, klon
         zt(i)=temp(i,k)
@@ -812,5 +830,5 @@
                 ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
                    ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
-                   ! for boundary-layer mixed phase clouds following Vignon et al.  
+                   ! for boundary-layer mixed phase clouds 
                     CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
                                      pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
@@ -834,15 +852,15 @@
             
                 ! lognormal 
-            lognormale = .TRUE.
+            lognormale(:) = .TRUE.
 
         ELSEIF (iflag_cld_th .GE. 6) THEN
             
                 ! lognormal distribution when no thermals 
-            lognormale = fraca(:,k) < min_frac_th_cld
+            lognormale(:) = fraca(:,k) < min_frac_th_cld
 
         ELSE
                 ! When iflag_cld_th=5, we always assume
                 ! bi-gaussian distribution
-            lognormale = .FALSE.
+            lognormale(:) = .FALSE.
         
         ENDIF
@@ -900,8 +918,8 @@
                   IF (iflag_t_glace.GE.4) THEN
                   ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
-                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
+                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
                   ENDIF
-                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:))
-
+
+                  CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:),ztemp_cltop(:),zfice(:),dzfice(:))
 
                   !--AB Activates a condensation scheme that allows for
@@ -938,5 +956,5 @@
                         pplay(:,k), paprs(:,k), paprs(:,k+1), &
                         cf_seri(:,k), rvc_seri(:,k), ratio_qi_qtot(:,k), &
-                        shear(:), pbl_eps(:,k), cell_area(:), &
+                        shear(:), tke_dissip(:,k), cell_area(:), &
                         Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:,k), keepgoing(:), &
                         rneb(:,k), zqn(:), qvc(:), issrfra(:,k), qissr(:,k), &
@@ -1017,5 +1035,6 @@
                             cste=RLSTT
                         ENDIF
-
+                        
+                        ! LEA_R : check formule
                         IF ( ok_unadjusted_clouds ) THEN
                           !--AB We relax the saturation adjustment assumption
@@ -1059,17 +1078,25 @@
         ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
         IF (iflag_t_glace.GE.4) THEN
-            CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D)
-            distcltop(:,k)=distcltop1D(:)
-            temp_cltop(:,k)=temp_cltop1D(:)
-        ENDIF   
-        ! 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 distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
+           distcltop(:,k)=zdistcltop(:)
+           temp_cltop(:,k)=ztemp_cltop(:)
+        ENDIF
+
+        ! Partition function depending on temperature
+        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
+
+        ! Partition function depending on tke for non shallow-convective clouds 
+        IF (iflag_icefrac .GE. 1) THEN
+
+           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), qice_save(:,k), ziflcld, zqn, &
+           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
+
+        ENDIF
 
         ! Water vapor update, Phase determination and subsequent latent heat exchange
         DO i=1, klon
-
+            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the 
+            ! iflag_cloudth_vert=7 and specific param is activated
             IF (mpc_bl_points(i,k) .GT. 0) THEN
-               
                 zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
                 ! following line is very strange and probably wrong
@@ -1078,7 +1105,5 @@
                 zq(i) = zq(i) - zcond(i)        
                 zfice(i)=zfice_th(i)
-
             ELSE
-
                 ! Checks on rneb, rhcl and zqn
                 IF (rneb(i,k) .LE. 0.0) THEN
@@ -1108,4 +1133,12 @@
                     ! following line is very strange and probably wrong:
                     rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
+                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param) 
+                    IF (iflag_icefrac .GE. 1) THEN
+                        IF (lognormale(i)) THEN  
+                           zcond(i)  = zqliq(i) + zqice(i)
+                           zfice(i)=zfice_turb(i)
+                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
+                        ENDIF
+                    ENDIF
                 ENDIF
 
@@ -1493,5 +1526,8 @@
                 znebprecipcld(i)=0.0
             ENDIF
-
+        !IF ( ((1-zfice(i))*zoliq(i) .GT. 0.) .AND. (zt(i) .LE. 233.15) ) THEN
+        !print*,'WARNING LEA OLIQ A <-40°C '
+        !print*,'zt,Tbef,oliq,oice,cldfraliq,icefrac,rneb',zt(i),Tbef(i),(1-zfice(i))*zoliq(i),zfice(i)*zoliq(i),cldfraliq(i,k),zfice(i),rneb(i,k)
+        !ENDIF 
         ENDDO
 
