Index: /LMDZ6/trunk/libf/phylmd/calcratqs.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/calcratqs.F90	(revision 4008)
+++ /LMDZ6/trunk/libf/phylmd/calcratqs.F90	(revision 4009)
@@ -3,7 +3,15 @@
            ratqsbas,ratqshaut,ratqsp0,ratqsdp, &
            tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
-           ptconv,ptconvth,clwcon0th, rnebcon0th,      &
-           paprs,pplay,q_seri,zqsat,fm_therm,          &
+           ptconv,ptconvth,clwcon0th, rnebcon0th,       &
+           paprs,pplay,t_seri,q_seri,                   &
+           qtc_cv, sigt_cv, zqsat,             &
+           tke,tke_dissip,lmix,wprime, &
+           t2m,q2m,fm_therm, &
            ratqs,ratqsc,ratqs_inter)
+
+
+USE indice_sol_mod
+USE phys_state_var_mod, ONLY: pctsrf
+USE calcratqs_multi_mod, ONLY: calcratqs_inter, calcratqs_oro, calcratqs_hetero, calcratqs_tke
 
 implicit none
@@ -23,11 +31,13 @@
 real,intent(in) :: pdtphys,ratqsbas,ratqshaut,fact_cldcon,tau_ratqs
 real,intent(in) :: ratqsp0, ratqsdp
-real, dimension(klon,klev+1),intent(in) :: paprs
-real, dimension(klon,klev),intent(in) :: pplay,q_seri,zqsat,fm_therm
+real, dimension(klon,klev+1),intent(in) :: paprs,tke,tke_dissip,lmix,wprime
+real, dimension(klon,klev),intent(in) :: pplay,t_seri,q_seri,zqsat,fm_therm, qtc_cv, sigt_cv
 logical, dimension(klon,klev),intent(in) :: ptconv
 real, dimension(klon,klev),intent(in) :: rnebcon0th,clwcon0th
 real, dimension(klon,klev),intent(in) :: wake_deltaq,wake_s
+real, dimension(klon,nbsrf),intent(in) :: t2m,q2m
 ! Output
 real, dimension(klon,klev),intent(inout) :: ratqs,ratqsc,ratqs_inter
+
 logical, dimension(klon,klev),intent(inout) :: ptconvth
 
@@ -36,4 +46,6 @@
 real, dimension(klon,klev) :: ratqss
 real facteur,zfratqs1,zfratqs2
+real, dimension(klon,klev) :: ratqs_hetero,ratqs_oro,ratqs_tke
+
 
 !-------------------------------------------------------------------------
@@ -131,13 +143,53 @@
          enddo
 
-      else if (iflag_ratqs==10) then ! ratqs interactif dépendant de la présence de poches froides
-         call calcratqs_inter(klon,klev,pdtphys,ratqsbas,wake_deltaq,wake_s,q_seri,ratqs_inter)
-         do k=1,klev
-            do i=1, klon
-                ratqss(i,k)=ratqs_inter(i,k)+0.5*(ratqshaut-ratqs_inter(i,k)) &
-                *( tanh( (ratqsp0-pplay(i,k))/ratqsdp) + 1.)
-            enddo
-         enddo
-
+       else if (iflag_ratqs .GT. 9) then
+ 
+       ! interactive ratqs calculations that depend on cold pools, orography, surface heterogeneity and small-scale turbulence
+       ! This should help getting a more realistic ratqs in the low and mid troposphere
+       ! We however need a "background" ratqs to account for subgrid distribution of qt (or qt/qs)
+       ! in the high troposphere
+       
+       ! background ratqs and initialisations
+          do k=1,klev
+             do i=1,klon
+              ratqss(i,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
+              *( tanh( (ratqsp0-pplay(i,k))/ratqsdp) + 1.)
+              ratqss(i,k)=max(ratqss(i,k),0.0)
+
+              ratqs_hetero(i,k)=0.
+              ratqs_oro(i,k)=0.
+              ratqs_tke(i,k)=0.
+              ratqs_inter(i,k)=0
+             enddo
+          enddo
+      
+          if (iflag_ratqs .EQ. 10) then 
+             ! interactive ratqs in presence of cold pools     
+             call calcratqs_inter(klon,klev,iflag_ratqs,pdtphys,ratqsbas,wake_deltaq,wake_s,q_seri,qtc_cv, sigt_cv,ratqs_inter)
+             do k=1,klev
+                do i=1,klon
+                    ratqs_inter(i,k)=ratqs_inter(i,k)-0.5*ratqs_inter(i,k)*(tanh((ratqsp0-pplay(i,k))/ratqsdp)+1.)
+                enddo
+             enddo
+             ratqss=ratqss+ratqs_inter
+          else if (iflag_ratqs .EQ. 11) then
+            ! interactive ratqs with several sources
+            call calcratqs_inter(klon,klev,iflag_ratqs,pdtphys,ratqsbas,wake_deltaq,wake_s,q_seri,qtc_cv, sigt_cv,ratqs_inter)
+             ratqss=ratqss+ratqs_inter
+          else if (iflag_ratqs .EQ. 12) then
+             ! contribution of surface heterogeneities to ratqs
+             call calcratqs_hetero(klon,klev,t2m,q2m,t_seri,q_seri,pplay,paprs,ratqs_hetero)
+             ratqss=ratqss+ratqs_hetero
+          else if (iflag_ratqs .EQ. 13) then
+             ! contribution of ubgrid orography to ratqs
+             call calcratqs_oro(klon,klev,zqsat,t_seri,pplay,paprs,ratqs_oro)
+             ratqss=ratqss+ratqs_oro
+          else if (iflag_ratqs .EQ. 14) then
+             ! effect of subgrid-scale TKE on ratqs (in development)
+             call calcratqs_tke(klon,klev,pdtphys,t_seri,q_seri,zqsat,pplay,paprs,tke,tke_dissip,lmix,wprime,ratqs_tke)     
+             ratqss=ratqss+ratqs_tke
+          endif
+          
+      
       endif
 
Index: DZ6/trunk/libf/phylmd/calcratqs_inter.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/calcratqs_inter.F90	(revision 4008)
+++ 	(revision )
@@ -1,69 +1,0 @@
-SUBROUTINE calcratqs_inter(klon,klev,pdtphys, &
-           ratqsbas, wake_deltaq, wake_s, q_seri,     &
-           ratqs_inter)
-USE ioipsl_getin_p_mod, ONLY : getin_p
-implicit none
-
-!========================================================================
-! Subroutine écrite par L. d'Alençon le 25/02/2021
-! Cette subroutine calcule une valeur de ratqsbas interactive dépendant de la présence de poches froides dans l'environnement. 
-! Elle est appelée par la subroutine calcratqs lorsque iflag_ratqs = 10. 
-!========================================================================
-
-! Declarations
-
-
-LOGICAL, SAVE :: first = .TRUE.
-!$OMP THREADPRIVATE(first)
-! Input
-integer,intent(in) :: klon,klev
-real,intent(in) :: pdtphys,ratqsbas
-real, dimension(klon,klev),intent(in) :: wake_deltaq, q_seri
-real, dimension(klon),intent(in) :: wake_s
-
-! Output
-real, dimension(klon,klev),intent(inout) :: ratqs_inter
-
-! local
-integer i,k
-real, dimension(klon,klev) :: wake_dq
-REAL, SAVE             :: tau_ratqs_wake
-!$OMP THREADPRIVATE(tau_ratqs_wake) 
-REAL, SAVE             :: a_ratqs_wake
-!$OMP THREADPRIVATE(a_ratqs_wake)
-real, dimension(klon) :: max_wake_dq, max_dqconv
-!-------------------------------------------------------------------------
-!  Caclul de ratqs_inter
-!-------------------------------------------------------------------------
-
-!  
-      if (first) then
-         tau_ratqs_wake = 3600. ! temps de relaxation de la variabilité
-         a_ratqs_wake = 1.85    ! paramètre pilotant l'importance du terme dépendant des poches froides
-         CALL getin_p('tau_ratqs_wake', tau_ratqs_wake)
-         CALL getin_p('a_ratqs_wake', a_ratqs_wake)         
-         first=.false.
-      endif       
-      max_wake_dq(:) = 0.
-      max_dqconv (:) = 0
-      
-      do k=1,klev
-        do i=1,klon
-         max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i))
-        enddo
-      enddo
-      
-      do k=1,klev
-        do i=1,klon
-         ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) +   &
-         a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1)   
-         if (ratqs_inter(i,k)<ratqsbas) then
-            ratqs_inter(i,k) = ratqsbas
-         endif       
-        enddo
-      enddo
-
-
-
-return
-end
Index: /LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90	(revision 4009)
+++ /LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90	(revision 4009)
@@ -0,0 +1,367 @@
+MODULE calcratqs_multi_mod
+
+!=============================================
+! module containing subroutines that take 
+! into account the effect of convection, orography,
+! surface heterogeneities and subgrid-scale
+! turbulence on ratqs, i.e. on the width of the
+! total water subgrid distribution.
+!=============================================
+
+IMPLICIT NONE
+
+! Include
+!=============================================
+    INCLUDE "YOETHF.h"
+    INCLUDE "YOMCST.h"
+
+
+      CONTAINS
+
+
+!========================================================================    
+SUBROUTINE calcratqs_inter(klon,klev,iflag_ratqs,pdtphys, &
+           ratqsbas, wake_deltaq, wake_s, q_seri,qtc_cv, sigt_cv,     &
+           ratqs_inter)
+USE ioipsl_getin_p_mod, ONLY : getin_p
+implicit none
+
+!========================================================================
+! L. d'Alençon, 25/02/2021
+! Cette subroutine calcule une valeur de ratqsbas interactive dépendant de la présence de poches froides dans l'environnement. 
+! Elle est appelée par la subroutine calcratqs lorsque iflag_ratqs = 10. 
+!========================================================================
+
+! Declarations
+
+
+LOGICAL, SAVE :: first = .TRUE.
+!$OMP THREADPRIVATE(first)
+! Input
+integer,intent(in) :: klon,klev,iflag_ratqs
+real,intent(in) :: pdtphys,ratqsbas
+real, dimension(klon,klev),intent(in) :: wake_deltaq, q_seri,qtc_cv, sigt_cv
+real, dimension(klon),intent(in) :: wake_s
+
+! Output
+real, dimension(klon,klev),intent(inout) :: ratqs_inter
+
+! local
+integer i,k
+real, dimension(klon,klev) :: wake_dq
+REAL, SAVE             :: a_ratqs_cv
+!$OMP THREADPRIVATE(tau_ratqs_wake)
+REAL, SAVE             :: tau_ratqs_wake
+!$OMP THREADPRIVATE(tau_ratqs_wake) 
+REAL, SAVE             :: a_ratqs_wake
+!$OMP THREADPRIVATE(a_ratqs_wake)
+real, dimension(klon) :: max_wake_dq, max_dqconv,max_sigt
+!-------------------------------------------------------------------------
+!  Caclul de ratqs_inter
+!-------------------------------------------------------------------------
+
+!  
+      if (first) then
+         tau_ratqs_wake = 3600. ! temps de relaxation de la variabilité
+         a_ratqs_wake = 3.    ! paramètre pilotant l'importance du terme dépendant des poches froides
+         a_ratqs_cv = 0.5
+         CALL getin_p('tau_ratqs_wake', tau_ratqs_wake)
+         CALL getin_p('a_ratqs_wake', a_ratqs_wake) 
+         CALL getin_p('a_ratqs_cv', a_ratqs_cv)         
+         first=.false.
+      endif       
+      max_wake_dq(:) = 0.
+      max_dqconv (:) = 0
+      max_sigt(:)    = 0. 
+      if (iflag_ratqs.eq.10) then
+        do k=1,klev
+          do i=1,klon
+           max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i))
+           max_sigt(i) = max(abs(sigt_cv(i,k)),max_sigt(i))
+           max_dqconv(i) = max(abs(q_seri(i,k) - qtc_cv(i,k)),max_dqconv(i))
+          enddo
+        enddo
+      
+        do k=1,klev
+          do i=1,klon
+           ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) +   &
+           a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1)   
+           if (ratqs_inter(i,k)<ratqsbas) then
+              ratqs_inter(i,k) = ratqsbas
+           endif       
+          enddo
+        enddo
+      endif
+
+      if (iflag_ratqs.eq.11) then
+        do k=1,klev
+          do i=1,klon
+           max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i))
+           max_sigt(i) = max(abs(sigt_cv(i,k)),max_sigt(i))
+           max_dqconv(i) = max(abs(q_seri(i,k) - qtc_cv(i,k)),max_dqconv(i))
+          enddo
+        enddo
+      
+        do k=1,klev
+          do i=1,klon
+           ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) +   &
+           a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1)  +   &
+           a_ratqs_cv*max_dqconv(i)*max_sigt(i)*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1)
+!           if (ratqs_inter(i,k)>0) then
+!              ratqs_inter(i,k) = abs(q_seri(i,k) - qtc_cv(i,k))
+!           endif       
+          enddo
+        enddo
+      endif
+return
+end          
+
+
+!------------------------------------------------------------------
+SUBROUTINE calcratqs_oro(klon,klev,qsat,temp,pplay,paprs,ratqs_oro)
+
+! Etienne Vignon, November 2021: effect of subgrid orography on ratqs
+
+USE phys_state_var_mod, ONLY: zstd
+USE phys_state_var_mod, ONLY: pctsrf
+USE indice_sol_mod, only: nbsrf, is_lic, is_ter
+
+IMPLICIT NONE
+
+! Declarations
+!--------------
+
+! INPUTS
+
+INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
+INTEGER, INTENT(IN) :: klev                       ! number of vertical layers
+REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat    ! saturation specific humidity [kg/kg]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: temp    ! air temperature [K]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay    ! air pressure, layer's center [Pa]
+REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs    ! air pressure, lower inteface [Pa]
+
+! OUTPUTS
+
+REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_oro ! ratqs profile due to subgrid orography
+
+
+! LOCAL
+
+INTEGER :: i,k
+REAL, DIMENSION(klon) :: orogradT,xsi0
+REAL, DIMENSION (klon,klev) :: zlay
+REAL :: Lvs, temp0
+
+
+! Calculation of the near-surface temperature gradient along the topography
+!--------------------------------------------------------------------------
+
+! at the moment, we fix it at a constant value (moist adiab. lapse rate)
+
+orogradT(:)=-6.5/1000. ! K/m
+
+! Calculation of near-surface surface ratqs
+!-------------------------------------------
+
+DO i=1,klon
+    temp0=temp(i,1)
+    IF (temp0 .LT. RTT) THEN
+        Lvs=RLSTT
+    ELSE
+        Lvs=RLVTT
+    ENDIF
+    xsi0(i)=zstd(i)*ABS(orogradT(i))*Lvs/temp0/temp0/RV
+    ratqs_oro(i,1)=xsi0(i)
+END DO
+
+! Vertical profile of ratqs assuming an exponential decrease with height
+!------------------------------------------------------------------------
+     
+! calculation of geop. height AGL        
+zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) &
+           *(paprs(:,1)-pplay(:,1))/RG
+
+DO k=2,klev
+   DO i = 1, klon
+      zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) &
+               /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG
+               
+      ratqs_oro(i,k)=MAX(0.0,pctsrf(i,is_ter)*xsi0(i)*exp(-zlay(i,k)/MAX(zstd(i),1.)))   
+    END DO
+END DO
+
+
+
+
+END SUBROUTINE calcratqs_oro
+
+!=============================================
+
+SUBROUTINE calcratqs_hetero(klon,klev,t2m,q2m,temp,q,pplay,paprs,ratqs_hetero)
+
+! Etienne Vignon, November 2021
+! Effect of subgrid surface heterogeneities on ratqs
+
+USE phys_local_var_mod, ONLY: s_pblh
+USE phys_state_var_mod, ONLY: pctsrf
+USE indice_sol_mod
+USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF
+
+IMPLICIT NONE
+
+include "YOMCST.h"
+
+! INPUTS
+
+
+INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
+INTEGER, INTENT(IN) :: klev                       ! number of vertical layers
+REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: t2m    ! 2m temperature for each tile [K]
+REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: q2m    ! 2m specific humidity for each tile [kg/kg]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: temp    ! air temperature [K]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: q       ! specific humidity [kg/kg]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay   ! air pressure, layer's center [Pa]
+REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa]
+
+! OUTPUTS
+
+REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_hetero ! ratsq profile due to surface heterogeneities
+
+
+INTEGER :: i,k,nsrf
+REAL :: ratiom, qsat2m, dqsatdT
+REAL, DIMENSION(klon) :: xsi0
+REAL, DIMENSION (klon,klev) :: zlay
+
+
+
+! Calculation of near-surface surface ratqs
+!-------------------------------------------
+
+DO i=1,klon
+    
+    ratiom=0.
+    xsi0(i)=0.
+    
+    DO nsrf=1,nbsrf
+    CALL CALC_QSAT_ECMWF(t2m(i,nsrf),q2m(i,nsrf),paprs(i,1),RTT,0,.false.,qsat2m,dqsatdT)
+    ratiom=ratiom+pctsrf(i,nsrf)*(q2m(i,nsrf)/qsat2m)
+    xsi0(i)=xsi0(i)+pctsrf(i,nsrf)*((q2m(i,nsrf)/qsat2m-ratiom)**2)
+    END DO
+    
+    xsi0(i)=sqrt(xsi0(i))/(ratiom+1E-6)
+END DO
+
+
+
+! Vertical profile of ratqs assuming an exponential decrease with height
+!------------------------------------------------------------------------
+        
+! calculation of geop. height AGL
+
+zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) &
+           *(paprs(:,1)-pplay(:,1))/RG
+ratqs_hetero(:,1)=xsi0(:)
+
+DO k=2,klev
+   DO i = 1, klon
+      zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) &
+               /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG
+               
+      ratqs_hetero(i,k)=MAX(xsi0(i)*exp(-zlay(i,k)/(s_pblh(i)+1.0)),0.0)   
+    END DO
+END DO
+
+END SUBROUTINE calcratqs_hetero
+
+!=============================================
+
+SUBROUTINE calcratqs_tke(klon,klev,pdtphys,temp,q,qsat,pplay,paprs,tke,tke_dissip,lmix,wprime,ratqs_tke)
+
+! References:
+!
+! Etienne Vignon: effect of subgrid turbulence on ratqs
+!
+! Field, P.R., Hill, A., Furtado, K., Korolev, A., 2014b. Mixed-phase clouds in a turbulent environment. Part
+! 2: analytic treatment. Q. J. R. Meteorol. Soc. 21, 2651–2663. https://doi.org/10.1002/qj.2175.
+! 
+! Furtado, K., Field, P.R., Boutle, I.A., Morcrette, C.R., Wilkinson, J., 2016. A physically-based, subgrid
+! parametrization for the production and maintenance of mixed-phase clouds in a general circulation
+! model. J. Atmos. Sci. 73, 279–291. https://doi.org/10.1175/JAS-D-15-0021.
+
+USE phys_local_var_mod, ONLY: omega
+
+IMPLICIT NONE
+
+! INPUTS
+
+INTEGER, INTENT(IN) :: klon                             ! number of horizontal grid points
+INTEGER, INTENT(IN) :: klev                             ! number of vertical layers
+REAL, INTENT(IN) :: pdtphys                             ! physics time step [s]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: temp          ! air temperature [K]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: q             ! specific humidity [kg/kg]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat          ! saturation specific humidity [kg/kg]
+REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay         ! air pressure, layer's center [Pa]
+REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs       ! air pressure, lower inteface [Pa]
+REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke         ! Turbulent Kinetic Energy [m2/s2]
+REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip  ! Turbulent Kinetic Energy Dissipation rate [m2/s3]
+REAL, DIMENSION(klon,klev+1), INTENT(IN) :: lmix  ! Turbulent mixing length
+REAL, DIMENSION(klon,klev+1), INTENT(IN) :: wprime      ! Turbulent vertical velocity scale [m/s]
+
+! OUTPUTS
+
+REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_tke  ! ratsq profile due to subgrid TKE
+
+! LOCAL
+INTEGER :: i, k
+REAL :: AA, DD, NW, AAprime, VARLOG,rho,Lvs,taue,lhomo,dissmin,maxvarlog
+REAL, DIMENSION(klon,klev) :: sigmaw,w
+REAL, PARAMETER :: C0=10.0
+REAL, PARAMETER :: lmin=0.001
+REAL, PARAMETER :: ratqsmin=1E-6
+REAL, PARAMETER :: ratqsmax=0.5
+
+
+! Calculation of large scale and turbulent vertical velocities
+!---------------------------------------------------------------
+
+DO k=1,klev
+    DO i=1,klon
+        rho=pplay(i,k)/temp(i,k)/RD
+        w(i,k)=-rho*RG*omega(i,k)
+        sigmaw(i,k)=0.5*(wprime(i,k+1)+wprime(i,k)) ! turbulent vertical velocity at the middle of model layers. 
+    END DO
+END DO
+
+! Calculation of ratqs
+!---------------------------------------------------------------
+ratqs_tke(:,1)=ratqsmin ! set to a very low value to avoid division by 0 in order parts
+                        ! of the code
+DO k=2,klev ! we start from second model level since TKE is not defined at k=1
+    DO i=1,klon
+
+       IF (temp(i,k) .LT. RTT) THEN
+           Lvs=RLSTT
+       ELSE
+           Lvs=RLVTT
+       ENDIF
+       dissmin=0.01*(0.5*(tke(i,k)+tke(i,k+1))/pdtphys)
+       maxvarlog=LOG(1.0+ratqsmax**2)! to prevent ratqs from exceeding an arbitrary threshold value 
+       AA=RG*(Lvs/(RCPD*temp(i,k)*temp(i,k)*RV) - 1./(RD*temp(i,k)))
+       lhomo=MAX(0.5*(lmix(i,k)+lmix(i,k+1)),lmin)
+       taue=(lhomo*lhomo/MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin))**(1./3) ! Fields et al. 2014
+       DD=1.0/taue
+       NW=(sigmaw(i,k)**2)*SQRT(2./(C0*MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin)))
+       AAprime=AA*NW
+       VARLOG=AAprime/2./DD
+       VARLOG=MIN(VARLOG,maxvarlog)
+       ratqs_tke(i,k)=SQRT(MAX(EXP(VARLOG)-1.0,ratqsmin))
+       END DO
+END DO
+END SUBROUTINE calcratqs_tke
+
+
+
+
+
+END MODULE calcratqs_multi_mod
Index: /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4008)
+++ /LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 4009)
@@ -16,6 +16,6 @@
       REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
       !$OMP THREADPRIVATE(u_seri, v_seri)
-      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:), l_mix(:,:,:), tke_dissip(:,:,:)
-      !$OMP THREADPRIVATE(l_mixmin, l_mix, tke_dissip)
+      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),tke_dissip(:,:,:),wprime(:,:,:)
+      !$OMP THREADPRIVATE(l_mixmin, l_mix, tke_dissip,wprime)
       REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
       !$OMP THREADPRIVATE(tr_seri)
@@ -559,6 +559,6 @@
       ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev))
       ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
-      ALLOCATE(l_mixmin(klon,klev+1,nbsrf), l_mix(klon,klev+1,nbsrf), tke_dissip(klon,klev+1,nbsrf))
-      l_mix(:,:,:)=0. ; l_mixmin(:,:,:)=0. ; tke_dissip(:,:,:)=0. ! doit etre initialse car pas toujours remplis
+      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),tke_dissip(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
+      l_mix(:,:,:)=0.;l_mixmin(:,:,:)=0.;tke_dissip(:,:,:)=0.;wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
 
       ALLOCATE(tr_seri(klon,klev,nbtr))
@@ -878,5 +878,5 @@
       DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri)
       DEALLOCATE(u_seri,v_seri)
-      DEALLOCATE(l_mixmin,l_mix, tke_dissip)
+      DEALLOCATE(l_mixmin,l_mix, tke_dissip,wprime)
 
       DEALLOCATE(tr_seri)
Index: /LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4008)
+++ /LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4009)
@@ -220,5 +220,5 @@
        d_t_vdf_x, d_t_vdf_w, &
        d_q_vdf_x, d_q_vdf_w, &
-       pbl_tke_input, &
+       pbl_tke_input, tke_dissip, l_mix, wprime,&
        t_therm, q_therm, u_therm, v_therm, &
        cdragh_x, cdragh_w, &
@@ -1186,4 +1186,5 @@
     integer iostat
 
+    REAL, dimension(klon,klev+1) :: tke_dissip_ave, l_mix_ave, wprime_ave
     REAL zzz
     !albedo SB >>>
@@ -3498,4 +3499,19 @@
     ! Computation of ratqs, the width (normalized) of the subrid scale 
     ! water distribution
+
+    tke_dissip_ave(:,:)=0.
+    l_mix_ave(:,:)=0.
+    wprime_ave(:,:)=0.
+
+
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
+          l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
+          wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
+       ENDDO
+    ENDDO
+
+
     CALL  calcratqs(klon,klev,prt_level,lunout,        &
          iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
@@ -3503,5 +3519,6 @@
          tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
          ptconv,ptconvth,clwcon0th, rnebcon0th,     &
-         paprs,pplay,q_seri,zqsat,fm_therm, &
+         paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
+         pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
          ratqs,ratqsc,ratqs_inter)
 
Index: /LMDZ6/trunk/libf/phylmd/yamada4.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/yamada4.F90	(revision 4008)
+++ /LMDZ6/trunk/libf/phylmd/yamada4.F90	(revision 4009)
@@ -6,5 +6,5 @@
   USE dimphy
   USE ioipsl_getin_p_mod, ONLY : getin_p
-  USE phys_local_var_mod, only: tke_dissip
+  USE phys_local_var_mod, only: tke_dissip,wprime
   
   IMPLICIT NONE
@@ -726,5 +726,5 @@
       lyam(1:ngrid, 2:klev)*5.17*kn(1:ngrid, 2:klev)*n2(1:ngrid, 2:klev)/ &
       sqrt(q2(1:ngrid,2:klev))
-
+ 
     t2yam(1:ngrid, 2:klev) = 9.1*kn(1:ngrid, 2:klev)* &
       dtetadz(1:ngrid, 2:klev)**2/sqrt(q2(1:ngrid,2:klev))* &
@@ -750,15 +750,17 @@
 
 !============================================================================
-! Diagnostique de la dissipation
+! Diagnostique de la dissipation et vitesse verticale
 !============================================================================
 
 ! Diagnostics
  tke_dissip(1:ngrid,:,nsrf)=0.
-! DO k=2,klev
-!    DO ig=1,ngrid
-!       jg=ni(ig)
-!       tke_dissip(jg,k,nsrf)=dissip(ig,k)
-!    ENDDO
-! ENDDO
+ wprime(1:ngrid,:,nsrf)=0.
+ DO k=2,klev
+    DO ig=1,ngrid
+       jg=ni(ig)
+       wprime(jg,k,nsrf)=sqrt(MAX(1./3*q2(ig,k),0.))
+       tke_dissip(jg,k,nsrf)=dissip(ig,k)
+    ENDDO
+ ENDDO
  
 !=============================================================================
Index: /LMDZ6/trunk/libf/phylmdiso/calcratqs_multi_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/calcratqs_multi_mod.F90	(revision 4009)
+++ /LMDZ6/trunk/libf/phylmdiso/calcratqs_multi_mod.F90	(revision 4009)
@@ -0,0 +1,1 @@
+link ../phylmd/calcratqs_multi_mod.F90
Index: /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4008)
+++ /LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90	(revision 4009)
@@ -16,6 +16,6 @@
       REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
       !$OMP THREADPRIVATE(u_seri, v_seri)
-      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:), l_mix(:,:,:), tke_dissip(:,:,:)
-      !$OMP THREADPRIVATE(l_mixmin, l_mix, tke_dissip)
+      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:), l_mix(:,:,:), tke_dissip(:,:,:), wprime(:,:,:)
+      !$OMP THREADPRIVATE(l_mixmin, l_mix, tke_dissip, wprime)
       REAL, SAVE, ALLOCATABLE :: tr_seri(:,:,:)
       !$OMP THREADPRIVATE(tr_seri)
@@ -693,6 +693,6 @@
       ALLOCATE(t_seri(klon,klev),q_seri(klon,klev),ql_seri(klon,klev),qs_seri(klon,klev))
       ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
-      ALLOCATE(l_mixmin(klon,klev+1,nbsrf), l_mix(klon,klev+1,nbsrf), tke_dissip(klon,klev+1,nbsrf))
-      l_mix(:,:,:)=0. ; l_mixmin(:,:,:)=0. ; tke_dissip(:,:,:)=0. ! doit etre initialse car pas toujours remplis
+      ALLOCATE(l_mixmin(klon,klev+1,nbsrf), l_mix(klon,klev+1,nbsrf), tke_dissip(klon,klev+1,nbsrf), wprime(klon,klev+1,nbsrf))
+      l_mix(:,:,:)=0. ; l_mixmin(:,:,:)=0. ; tke_dissip(:,:,:)=0. ; wprime(:,:,:)=0. ! doit etre initialse car pas toujours remplis
 
       ALLOCATE(tr_seri(klon,klev,nbtr))
@@ -1090,5 +1090,5 @@
       DEALLOCATE(t_seri,q_seri,ql_seri,qs_seri)
       DEALLOCATE(u_seri,v_seri)
-      DEALLOCATE(l_mixmin,l_mix, tke_dissip)
+      DEALLOCATE(l_mixmin,l_mix, tke_dissip, wprime)
 
       DEALLOCATE(tr_seri)
Index: /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4008)
+++ /LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90	(revision 4009)
@@ -264,5 +264,5 @@
        d_t_vdf_x, d_t_vdf_w, &
        d_q_vdf_x, d_q_vdf_w, &
-       pbl_tke_input, &
+       pbl_tke_input, tke_dissip, l_mix, wprime, &
        t_therm, q_therm, u_therm, v_therm, &
        cdragh_x, cdragh_w, &
@@ -1288,4 +1288,5 @@
     integer iostat
 
+    REAL, dimension(klon,klev+1) :: tke_dissip_ave, l_mix_ave, wprime_ave
     REAL zzz
     !albedo SB >>>
@@ -4468,11 +4469,26 @@
     ! Computation of ratqs, the width (normalized) of the subrid scale 
     ! water distribution
-    CALL  calcratqs(klon,klev,prt_level,lunout,        &
+
+    tke_dissip_ave(:,:)=0.
+    l_mix_ave(:,:)=0.
+    wprime_ave(:,:)=0.
+
+    DO nsrf = 1, nbsrf
+       DO i = 1, klon
+          tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
+          l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
+          wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
+       ENDDO
+    ENDDO
+
+   CALL  calcratqs(klon,klev,prt_level,lunout,        &
          iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
          ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
          tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
          ptconv,ptconvth,clwcon0th, rnebcon0th,     &
-         paprs,pplay,q_seri,zqsat,fm_therm, &
+         paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
+         pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
          ratqs,ratqsc,ratqs_inter)
+
 
     !
