Index: LMDZ6/trunk/libf/phylmd/calcratqs.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/calcratqs.F90	(revision 4608)
+++ 	(revision )
@@ -1,264 +1,0 @@
-SUBROUTINE 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,t_seri,q_seri,                   &
-           qtc_cv, sigt_cv, zqsat,             &
-           tke,tke_dissip,lmix,wprime, &
-           t2m,q2m,fm_therm,cell_area,&
-           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
-
-!========================================================================
-! Computation of ratqs, the width of the subrid scale water distribution
-! (normalized by the mean value)
-! Various options controled by flags iflag_con and iflag_ratqs
-! F Hourdin 2012/12/06
-!========================================================================
-
-! Declarations
-
-! Input
-integer,intent(in) :: klon,klev,prt_level,lunout
-integer,intent(in) :: iflag_con,iflag_cld_th,iflag_ratqs
-real,intent(in) :: pdtphys,ratqsbas,ratqshaut,fact_cldcon,tau_ratqs
-real,intent(in) :: ratqsp0, ratqsdp
-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
-real, dimension(klon), intent(in) :: cell_area
-! Output
-real, dimension(klon,klev),intent(inout) :: ratqs,ratqsc,ratqs_inter
-
-logical, dimension(klon,klev),intent(inout) :: ptconvth
-
-! local
-integer i,k
-real, dimension(klon,klev) :: ratqss
-real facteur,zfratqs1,zfratqs2
-real, dimension(klon,klev) :: ratqs_hetero,ratqs_oro,ratqs_tke
-real resol,resolmax,fact
-
-!-------------------------------------------------------------------------
-!  Caclul des ratqs
-!-------------------------------------------------------------------------
-
-!      print*,'calcul des ratqs'
-!   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
-!   ----------------
-!   on ecrase le tableau ratqsc calcule par clouds_gno
-      if (iflag_cld_th.eq.1) then
-         do k=1,klev
-         do i=1,klon
-            if(ptconv(i,k)) then
-              ratqsc(i,k)=ratqsbas &
-              +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
-            else
-               ratqsc(i,k)=0.
-            endif
-         enddo
-         enddo
-
-!-----------------------------------------------------------------------
-!  par nversion de la fonction log normale
-!-----------------------------------------------------------------------
-      else if (iflag_cld_th.eq.4) then
-         ptconvth(:,:)=.false.
-         ratqsc(:,:)=0.
-         if(prt_level.ge.9) print*,'avant clouds_gno thermique'
-         call clouds_gno &
-         (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
-         if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
-       
-       endif
-
-!   ratqs stables
-!   -------------
-
-      if (iflag_ratqs.eq.0) then
-
-! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
-         do k=1,klev
-            do i=1, klon
-               ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)* &
-               min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.) 
-            enddo 
-         enddo
-
-! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de 
-! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
-! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
-! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
-! Il s'agit de differents tests dans la phase de reglage du modele
-! avec thermiques.
-
-      else if (iflag_ratqs.eq.1) then
-
-         do k=1,klev
-            do i=1, klon
-               if (pplay(i,k).ge.60000.) then
-                  ratqss(i,k)=ratqsbas
-               else if ((pplay(i,k).ge.30000.).and.(pplay(i,k).lt.60000.)) then
-                  ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*(60000.-pplay(i,k))/(60000.-30000.)
-               else
-                  ratqss(i,k)=ratqshaut
-               endif
-            enddo
-         enddo
-
-      else if (iflag_ratqs.eq.2) then
-
-         do k=1,klev
-            do i=1, klon
-               if (pplay(i,k).ge.60000.) then
-                  ratqss(i,k)=ratqsbas*(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
-               else if ((pplay(i,k).ge.30000.).and.(pplay(i,k).lt.60000.)) then
-                    ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*(60000.-pplay(i,k))/(60000.-30000.)
-               else
-                    ratqss(i,k)=ratqshaut
-               endif
-            enddo
-         enddo
-
-      else if (iflag_ratqs==3) then
-         do k=1,klev
-           ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas) &
-           *min( ((paprs(:,1)-pplay(:,k))/70000.)**2 , 1. )
-         enddo
-
-      else if (iflag_ratqs==4) then 
-         do k=1,klev
-           ratqss(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
-!          *( tanh( (50000.-pplay(:,k))/20000.) + 1.)
-           *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
-         enddo
-
-
-      else if (iflag_ratqs==5) then
-! Dependency of ratqs on model resolution
-! Audran, Meryl, Lea, Gwendal and Etienne
-! April 2023 
-        resolmax=sqrt(maxval(cell_area))
-         do k=1,klev
-            do i=1,klon
-              resol=sqrt(cell_area(i))
-              fact=sqrt(resol/resolmax)
-              ratqss(i,k)=ratqsbas*fact+0.5*(ratqshaut-ratqsbas)*fact &
-              *( 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
-
-
-!  ratqs final
-!  -----------
-
-      if (iflag_cld_th.eq.1 .or.iflag_cld_th.eq.2.or.iflag_cld_th.eq.4) then
-
-! On ajoute une constante au ratqsc*2 pour tenir compte de 
-! fluctuations turbulentes de petite echelle
-
-         do k=1,klev
-            do i=1,klon
-               if ((fm_therm(i,k).gt.1.e-10)) then
-                  ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
-               endif
-            enddo
-         enddo
-
-!   les ratqs sont une combinaison de ratqss et ratqsc
-       if(prt_level.ge.9) write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
-
-         if (tau_ratqs>1.e-10) then
-            facteur=exp(-pdtphys/tau_ratqs)
-         else
-            facteur=0.
-         endif
-         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! FH 22/09/2009
-! La ligne ci-dessous faisait osciller le modele et donnait une solution
-! assymptotique bidon et dépendant fortement du pas de temps.
-!        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-         ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
-      else if (iflag_cld_th<=6) then
-!   on ne prend que le ratqs stable pour fisrtilp
-         ratqs(:,:)=ratqss(:,:)
-      else
-          zfratqs1=exp(-pdtphys/10800.)
-          zfratqs2=exp(-pdtphys/10800.)
-          do k=1,klev
-             do i=1,klon
-                if (ratqsc(i,k).gt.1.e-10) then
-                   ratqs(i,k)=ratqs(i,k)*zfratqs2+(iflag_cld_th/100.)*ratqsc(i,k)*(1.-zfratqs2)
-                endif
-                ratqs(i,k)=min(ratqs(i,k)*zfratqs1+ratqss(i,k)*(1.-zfratqs1),0.5)
-             enddo
-          enddo
-      endif
-
-
-return
-end
Index: LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/calcratqs_multi_mod.F90	(revision 4608)
+++ 	(revision )
@@ -1,364 +1,0 @@
-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(a_ratqs_cv)
-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, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT
-REAL, DIMENSION (klon,klev) :: zlay
-
-
-
-! Calculation of near-surface surface ratqs
-!-------------------------------------------
-
-    
-    ratiom(:)=0.
-    xsi0(:)=0.
-    
-    DO nsrf=1,nbsrf
-    CALL CALC_QSAT_ECMWF(klon,t2m(:,nsrf),q2m(:,nsrf),paprs(:,1),RTT,0,.false.,qsat2m,dqsatdT)
-    ratiom(:)=ratiom(:)+pctsrf(:,nsrf)*(q2m(:,nsrf)/qsat2m(:))
-    xsi0(:)=xsi0(:)+pctsrf(:,nsrf)*((q2m(:,nsrf)/qsat2m(:)-ratiom(:))**2)
-    END DO
-    
-    xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6)
-
-
-
-! 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/cloudth_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/cloudth_mod.F90	(revision 4613)
@@ -30,4 +30,5 @@
       INTEGER itap,ind1,ind2
       INTEGER ngrid,klev,klon,l,ig 
+      INTEGER iflag_ratqs
       
       REAL ztv(ngrid,klev)
@@ -274,4 +275,5 @@
       INTEGER itap,ind1,ind2
       INTEGER ngrid,klev,klon,l,ig 
+      INTEGER iflag_ratqs
       
       REAL ztv(ngrid,klev)
@@ -613,4 +615,5 @@
       INTEGER itap,ind1,ind2
       INTEGER ngrid,klev,klon,l,ig 
+      INTEGER iflag_ratqs
       
       REAL ztv(ngrid,klev)
@@ -835,4 +838,5 @@
       INTEGER itap,ind1,ind2
       INTEGER ngrid,klev,klon,l,ig 
+      INTEGER iflag_ratqs
       
       REAL ztv(ngrid,klev)
@@ -1022,4 +1026,7 @@
       ENDIF
       sigma1s = sigma1s_fraca + sigma1s_ratqs
+           IF (iflag_ratqs.eq.11) then
+              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
+           ENDIF
       sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
 !      tests
@@ -1600,4 +1607,5 @@
 
       INTEGER itap,ind1,l,ig,iter,k 
+      INTEGER iflag_ratqs
       INTEGER iflag_topthermals, niter
       LOGICAL falseklon(klon)
@@ -1894,4 +1902,7 @@
  
            sigma1s = sigma1s_fraca + sigma1s_ratqs
+           IF (iflag_ratqs.eq.11) then
+              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
+           ENDIF
            sigma2s=(sigma2s_factor*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**sigma2s_power))+0.002*zqta(ind1,ind2)
 
@@ -2136,4 +2147,10 @@
 
             sigma1s = sigma1s_fraca + sigma1s_ratqs
+           IF (iflag_ratqs.eq.11) then
+              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenv
+           ENDIF
+           IF (iflag_ratqs.eq.11) then
+              sigma1s = ratqs(ind1,ind2)*po(ind1)*aenvl
+           ENDIF
             deltasenv=aenvl*vert_alpha*sigma1s
             xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
Index: LMDZ6/trunk/libf/phylmd/concvl.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/concvl.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/concvl.F90	(revision 4613)
@@ -3,5 +3,5 @@
                   t, q, t_wake, q_wake, s_wake, u, v, tra, ntra, &
                   Ale, Alp, sig1, w01, &
-                  d_t, d_q, d_u, d_v, d_tra, &
+                  d_t, d_q, d_qcomp, d_u, d_v, d_tra, &
                   rain, snow, kbas, ktop, sigd, &
                   cbmf, plcl, plfc, wbeff, convoccur, &
@@ -16,5 +16,5 @@
                   dd_t, dd_q, lalim_conv, wght_th, &                 ! RomP
                   evap, ep, epmlmMm, eplaMm, &                       ! RomP
-                  wdtrainA, wdtrainS, wdtrainM, wght, qtc, sigt, &
+                  wdtrainA, wdtrainS, wdtrainM, wght, qtc, sigt, detrain, &
                   tau_cld_cv, coefw_cld_cv, &                           ! RomP+RL, AJ
 !RomP <<<
@@ -105,5 +105,5 @@
   REAL, DIMENSION(klon,klev),   INTENT(INOUT)   :: sig1, w01
 
-  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: d_t, d_q, d_u, d_v
+  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: d_t, d_q, d_qcomp, d_u, d_v
   REAL, DIMENSION(klon,klev, nbtr),INTENT(OUT)  :: d_tra
   REAL, DIMENSION(klon),        INTENT(OUT)     :: rain, snow
@@ -145,5 +145,5 @@
   REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: wght                       !RL
   REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: qtc
-  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: sigt
+  REAL, DIMENSION(klon,klev),   INTENT(OUT)     :: sigt, detrain
   REAL,                         INTENT(OUT)     :: tau_cld_cv, coefw_cld_cv
   REAL, DIMENSION(klon),        INTENT(OUT)     :: epmax_diag                ! epmax_cape
@@ -429,5 +429,5 @@
                     Ale, Alp, omega, &
                     em_sig1feed, em_sig2feed, em_wght, &
-                    iflag, d_t, d_q, d_u, d_v, d_tra, rain, kbas, ktop, &
+                    iflag, d_t, d_q, d_qcomp, d_u, d_v, d_tra, rain, kbas, ktop, &
                     cbmf, plcl, plfc, wbeff, sig1, w01, ptop2, sigd, &
                     Ma, mip, Vprecip, Vprecipi, upwd, dnwd, dnwdbis, qcondc, wd, &
@@ -440,5 +440,5 @@
                     da, phi, mp, phii, d1a, dam, sij, wght, &           ! RomP+RL
                     qta, clw, elij, evap, ep, epmlmMm, eplaMm, &        ! RomP+RL
-                    wdtrainA, wdtrainS, wdtrainM, qtc, sigt, &
+                    wdtrainA, wdtrainS, wdtrainM, qtc, sigt, detrain, &
                     tau_cld_cv, coefw_cld_cv, &                         ! RomP,AJ
 !AC!+!RomP+jyg
Index: LMDZ6/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/conf_phys_m.F90	(revision 4613)
@@ -16,5 +16,5 @@
        fact_cldcon, facttemps,ok_newmicro,iflag_radia,&
        iflag_cld_th, &
-       iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
+       ratqsbas,ratqshaut,tau_ratqs, &
        ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, & 
        chemistry_couple, flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
@@ -85,5 +85,5 @@
     REAL                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut,tau_ratqs
     INTEGER              :: iflag_cld_th
-    INTEGER              :: iflag_ratqs
+     
 
     CHARACTER (len = 6), SAVE  :: type_ocean_omp, version_ocean_omp, ocean_omp
Index: LMDZ6/trunk/libf/phylmd/cv3_routines.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cv3_routines.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/cv3_routines.F90	(revision 4613)
@@ -3448,9 +3448,9 @@
                      tv, tvp, wghti, &
                      iflag, precip, Vprecip, Vprecipi, &     ! jyg: Vprecipi
-                     ft, fr, fu, fv, ftra, &                 ! jyg
+                     ft, fr, fr_comp, fu, fv, ftra, &                 ! jyg
                      cbmf, upwd, dnwd, dnwd0, ma, mip, &
 !!                     tls, tps,                             ! useless . jyg
                      qcondc, wd, &
-                     ftd, fqd, qta, qtc, sigt, tau_cld_cv, coefw_cld_cv)
+                     ftd, fqd, qta, qtc, sigt, detrain, tau_cld_cv, coefw_cld_cv)
 
     USE print_control_mod, ONLY: lunout, prt_level
@@ -3505,5 +3505,5 @@
 !outputs:
       REAL, DIMENSION (nloc), INTENT (OUT)               :: precip
-      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ft, fr, fu, fv
+      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ft, fr, fu, fv , fr_comp
       REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: ftd, fqd
       REAL, DIMENSION (nloc, nd, ntra), INTENT (OUT)     :: ftra
@@ -3515,4 +3515,5 @@
       REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: qcondc                      ! cld
       REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: qtc, sigt                   ! cld
+      REAL, DIMENSION (nloc, nd), INTENT (OUT)           :: detrain                     ! Louis : pour le calcul de Klein du terme de variance qui dÃ©traine dans lenvironnement
       REAL, DIMENSION (nloc), INTENT (OUT)               :: wd                          ! gust
       REAL, DIMENSION (nloc), INTENT (OUT)               :: cbmf
@@ -3541,4 +3542,5 @@
       REAL, DIMENSION (nloc)                             :: sument
       REAL, DIMENSION (nloc, nd)                         :: sigment, qtment             ! cld
+      REAL, DIMENSION (nloc, nd, nd)                     :: qdet
       REAL sumdq !jyg
 ! 
@@ -3570,4 +3572,5 @@
       ft(il, i) = 0.0
       fr(il, i) = 0.0
+      fr_comp(il,i) = 0.0
       fu(il, i) = 0.0
       fv(il, i) = 0.0
@@ -3580,4 +3583,6 @@
       sigment(il, i) = 0.0 ! cld
       sigt(il, i) = 0.0 ! cld
+      qdet(il,i,:) = 0.0 ! cld
+      detrain(il, i) = 0.0 ! cld
       nqcond(il, i) = 0.0 ! cld
     END DO
@@ -3811,4 +3816,5 @@
       IF (j<=inb(il) .AND. iflag(il)<=1) THEN
         fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il,j,1)-rr(il,1))
+        fr_comp(il,1) = fr_comp(il,1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il,j,1)-rr(il,1))
         fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il,j,1)-u(il,1))
         fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il,j,1)-v(il,1))
@@ -4213,4 +4219,5 @@
           fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i) * &
                                                        (qent(il,k,i)-awat(il)-rr(il,i))
+          fr_comp(il,i) = fr_comp(il,i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-awat(il)-rr(il,i))
           fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k,i)-u(il,i))
           fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k,i)-v(il,i))
@@ -4218,4 +4225,5 @@
 ! (saturated updrafts resulting from mixing)                                   ! cld
           qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il))                ! cld
+          qdet(il,k,i) = (qent(il,k,i)-awat(il))                               ! cld Louis : specific humidity in detraining water
           qtment(il, i) = qtment(il, i) + qent(il,k,i)                         ! cld
           nqcond(il, i) = nqcond(il, i) + 1.                                   ! cld
@@ -4298,4 +4306,5 @@
 ! (saturated downdrafts resulting from mixing)                                 ! cld
           qcond(il, i) = qcond(il, i) + elij(il, k, i)                         ! cld
+          qdet(il,k,i) = qent(il,k,i)                                          ! cld Louis : specific humidity in detraining water
           qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
           nqcond(il, i) = nqcond(il, i) + 1.                                   ! cld
@@ -4316,4 +4325,5 @@
 !      PROBLEM: Should not qent(il,i,i) be taken into account even if nent(il,i)/=0?
 !!        qtment(il, i) = qent(il,k,i) + qtment(il,i)                            ! cld
+        qdet(il,i,i) = qent(il,i,i)                                            ! cld Louis : specific humidity in detraining water
         qtment(il, i) = qent(il,i,i) + qtment(il,i)                            ! cld
 !>jyg
@@ -4887,4 +4897,5 @@
         IF  (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN   ! cld
           sument(il) =sument(il) + abs(ment(il,k,i))
+          detrain(il,i) = detrain(il,i) + abs(ment(il,k,i))*(qdet(il,k,i) - rr(il,i))*(qdet(il,k,i) - rr(il,i)) ! Louis terme de dÃ©trainement dans le bilan de variance
         ENDIF
       ENDDO     ! il
Index: LMDZ6/trunk/libf/phylmd/cv3a_uncompress.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cv3a_uncompress.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/cv3a_uncompress.F90	(revision 4613)
@@ -2,5 +2,5 @@
                            iflag, kbas, ktop, &
                            precip, cbmf, plcl, plfc, wbeff, sig, w0, ptop2, &
-                           ft, fq, fu, fv, ftra,  &
+                           ft, fq, fqcomp, fu, fv, ftra,  &
                            sigd, ma, mip, vprecip, vprecipi, upwd, dnwd, dnwd0, &
                            qcondc, wd, cape, cin, &
@@ -12,9 +12,9 @@
                            qta, clw, elij, evap, ep, epmlmMm, eplaMm, &         ! RomP+jyg
                            wdtrainA, wdtrainS, wdtrainM, &                      ! RomP
-                           qtc, sigt,          &
+                           qtc, sigt, detrain,         &
                            epmax_diag, & ! epmax_cape
                            iflag1, kbas1, ktop1, &
                            precip1, cbmf1, plcl1, plfc1, wbeff1, sig1, w01, ptop21, &
-                           ft1, fq1, fu1, fv1, ftra1, &
+                           ft1, fq1, fqcomp1, fu1, fv1, ftra1, &
                            sigd1, ma1, mip1, vprecip1, vprecipi1, upwd1, dnwd1, dnwd01, &
                            qcondc1, wd1, cape1, cin1, &
@@ -26,5 +26,5 @@
                            qta1, clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, &  ! RomP+jyg
                            wdtrainA1, wdtrainS1, wdtrainM1, &                   ! RomP
-                           qtc1, sigt1, &
+                           qtc1, sigt1, detrain1, &
                            epmax_diag1) ! epmax_cape
 
@@ -54,5 +54,5 @@
   REAL, DIMENSION (nloc), INTENT (IN)                :: ptop2
   REAL, DIMENSION (nloc), INTENT (IN)                :: epmax_diag
-  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: ft, fq, fu, fv
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: ft, fq, fqcomp, fu, fv
   REAL, DIMENSION (nloc, nd, ntra), INTENT (IN)      :: ftra
   REAL, DIMENSION (nloc), INTENT (IN)                :: sigd
@@ -81,5 +81,5 @@
   REAL, DIMENSION (nloc, nd, nd), INTENT (IN)        :: epmlmMm                !RomP+jyg
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: eplamM                 !RomP+jyg
-  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qtc, sigt              !RomP
+  REAL, DIMENSION (nloc, nd), INTENT (IN)            :: qtc, sigt, detrain              !RomP
   REAL, DIMENSION (nloc, nd), INTENT (IN)            :: wdtrainA, wdtrainS, wdtrainM     !RomP
 
@@ -91,5 +91,5 @@
   REAL, DIMENSION (len), INTENT (OUT)                :: epmax_diag1 ! epmax_cape
   REAL, DIMENSION (len), INTENT (OUT)                :: ptop21
-  REAL, DIMENSION (len, nd), INTENT (OUT)            :: ft1, fq1, fu1, fv1
+  REAL, DIMENSION (len, nd), INTENT (OUT)            :: ft1, fq1, fqcomp1, fu1, fv1
   REAL, DIMENSION (len, nd, ntra), INTENT (OUT)      :: ftra1
   REAL, DIMENSION (len), INTENT (OUT)                :: sigd1
@@ -118,5 +118,5 @@
   REAL, DIMENSION (len, nd, nd), INTENT (OUT)        :: epmlmMm1               !RomP+jyg
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: eplamM1                !RomP+jyg
-  REAL, DIMENSION (len, nd), INTENT (OUT)            :: qtc1, sigt1            !RomP
+  REAL, DIMENSION (len, nd), INTENT (OUT)            :: qtc1, sigt1, detrain1            !RomP
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: wdtrainA1, wdtrainS1, wdtrainM1   !RomP
 
@@ -158,4 +158,5 @@
         ft1(idcum(i), k) = ft(i, k)
         fq1(idcum(i), k) = fq(i, k)
+        fqcomp1(idcum(i), k) = fqcomp(i, k)
         fu1(idcum(i), k) = fu(i, k)
         fv1(idcum(i), k) = fv(i, k)
@@ -187,4 +188,5 @@
         qtc1(idcum(i), k) = qtc(i, k)
         sigt1(idcum(i), k) = sigt(i, k)
+        detrain1(idcum(i), k) = detrain(i, k)
    
       END DO
@@ -285,4 +287,5 @@
       ft1(:, 1:nl) = ft(:, 1:nl)
       fq1(:, 1:nl) = fq(:, 1:nl)
+      fqcomp1(:, 1:nl) = fqcomp(:, 1:nl)
       fu1(:, 1:nl) = fu(:, 1:nl)
       fv1(:, 1:nl) = fv(:, 1:nl)
@@ -314,4 +317,5 @@
       qtc1(:, 1:nl) = qtc(:, 1:nl)
       sigt1(:, 1:nl) = sigt(:, 1:nl)
+      detrain1(:, 1:nl) = detrain(:, 1:nl)
 !
       ma1(:, nlp) = 0.
Index: LMDZ6/trunk/libf/phylmd/cva_driver.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/cva_driver.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/cva_driver.F90	(revision 4613)
@@ -11,5 +11,5 @@
                       Ale1, Alp1, omega1, &
                       sig1feed1, sig2feed1, wght1, &
-                      iflag1, ft1, fq1, fu1, fv1, ftra1, &
+                      iflag1, ft1, fq1, fqcomp1, fu1, fv1, ftra1, &
                       precip1, kbas1, ktop1, &
                       cbmf1, plcl1, plfc1, wbeff1, &
@@ -26,5 +26,5 @@
                       da1, phi1, mp1, phi21, d1a1, dam1, sigij1, wghti1, & ! RomP, RL
                       qta1, clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, &  ! RomP, RL
-                      wdtrainA1, wdtrainS1, wdtrainM1, qtc1, sigt1, tau_cld_cv, &     !!jygprl
+                      wdtrainA1, wdtrainS1, wdtrainM1, qtc1, sigt1, detrain1, tau_cld_cv, &     !!jygprl
                       coefw_cld_cv, &                                      ! RomP, AJ
                       epmax_diag1)  ! epmax_cape
@@ -88,4 +88,5 @@
 ! ft1           Real           Output       temp tend
 ! fq1           Real           Output       spec hum tend
+! fqcomp1       Real           Output       spec hum tend (only mixed draughts)
 ! fu1           Real           Output       u-wind tend
 ! fv1           Real           Output       v-wind tend
@@ -135,4 +136,5 @@
 ! qtc1          Real           Output     specific humidity in convection
 ! sigt1         Real           Output     surface fraction in adiabatic updrafts                                         
+! detrain1      Real           Output     detrainment terme klein
 ! phi21         Real           Output     used in tracer transport (cvltr)
                                          
@@ -209,4 +211,5 @@
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: ft1
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: fq1
+  REAL, DIMENSION (len, nd), INTENT (OUT)            :: fqcomp1
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: fu1
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: fv1
@@ -248,4 +251,5 @@
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: qtc1    ! in cloud water content (intensive)   ! cld
   REAL, DIMENSION (len, nd), INTENT (OUT)            :: sigt1   ! fract. cloud area (intensive)        ! cld
+  REAL, DIMENSION (len, nd), INTENT (OUT)            :: detrain1   ! detrainement term of mixed draughts in environment
 
 ! RomP >>>
@@ -502,5 +506,5 @@
   REAL, DIMENSION(len,nd)     :: ice, fondue, b
   REAL, DIMENSION(len,nd)     :: frac_a, frac_s, faci               !!jygprl
-  REAL ft(nloc, nd), fq(nloc, nd)
+  REAL ft(nloc, nd), fq(nloc, nd), fqcomp(nloc, nd)
   REAL ftd(nloc, nd), fqd(nloc, nd)
   REAL fu(nloc, nd), fv(nloc, nd)
@@ -528,4 +532,5 @@
   REAL qtc(nloc, nd)         ! cld
   REAL sigt(nloc, nd)        ! cld
+  REAL detrain(nloc, nd)     ! cld
  
 ! RomP >>>
@@ -593,4 +598,5 @@
   ft1(:, :) = 0.0
   fq1(:, :) = 0.0
+  fqcomp1(:, :) = 0.0
   fu1(:, :) = 0.0
   fv1(:, :) = 0.0
@@ -646,4 +652,5 @@
 ! RomP >>>
   sigt1(:, :) = 0.
+  detrain1(:, :) = 0.
   qtc1(:, :) = 0.
   wdtrainA1(:, :) = 0.
@@ -1139,10 +1146,10 @@
                      nent, elij, traent, sig, &
                      tv, tvp, wghti, &
-                     iflag, precip, Vprecip, Vprecipi, ft, fq, fu, fv, ftra, &      ! jyg
+                     iflag, precip, Vprecip, Vprecipi, ft, fq, fqcomp, fu, fv, ftra, &      ! jyg
                      cbmf, upwd, dnwd, dnwd0, ma, mip, &
 !!                     tls, tps, &                            ! useless . jyg
                      qcondc, wd, &
 !!                     ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
-                     ftd, fqd, qta, qtc, sigt, tau_cld_cv, coefw_cld_cv)         !!jygprl
+                     ftd, fqd, qta, qtc, sigt, detrain, tau_cld_cv, coefw_cld_cv)         !!jygprl
 !
 !         Test conseravtion de l'eau
@@ -1204,5 +1211,5 @@
                            iflag, icb, inb, &
                            precip, cbmf, plcl, plfc, wbeff, sig, w0, ptop2, &
-                           ft, fq, fu, fv, ftra, &
+                           ft, fq, fqcomp, fu, fv, ftra, &
                            sigd, ma, mip, vprecip, vprecipi, upwd, dnwd, dnwd0, &
                            qcondc, wd, cape, cin, &
@@ -1214,8 +1221,8 @@
                            qta, clw, elij, evap, ep, epmlmMm, eplaMm, &  ! RomP
                            wdtrainA, wdtrainS, wdtrainM, &                         ! RomP
-                           qtc, sigt, epmax_diag, & ! epmax_cape
+                           qtc, sigt, detrain, epmax_diag, & ! epmax_cape
                            iflag1, kbas1, ktop1, &
                            precip1, cbmf1, plcl1, plfc1, wbeff1, sig1, w01, ptop21, &
-                           ft1, fq1, fu1, fv1, ftra1, &
+                           ft1, fq1, fqcomp1, fu1, fv1, ftra1, &
                            sigd1, ma1, mip1, vprecip1, vprecipi1, upwd1, dnwd1, dnwd01, &
                            qcondc1, wd1, cape1, cin1, &
@@ -1227,5 +1234,5 @@
                            qta1, clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, & ! RomP
                            wdtrainA1, wdtrainS1, wdtrainM1,                  & ! RomP
-                           qtc1, sigt1, epmax_diag1) ! epmax_cape
+                           qtc1, sigt1, detrain1, epmax_diag1) ! epmax_cape
 !   
       IF (prt_level >= 10) THEN
Index: LMDZ6/trunk/libf/phylmd/lmdz_ratqs_ini.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_ratqs_ini.F90	(revision 4613)
+++ LMDZ6/trunk/libf/phylmd/lmdz_ratqs_ini.F90	(revision 4613)
@@ -0,0 +1,80 @@
+MODULE lmdz_ratqs_ini
+
+IMPLICIT NONE
+
+save
+
+integer :: lunout
+
+INTEGER, SAVE, PROTECTED :: nbsrf,is_lic,is_ter
+REAL, SAVE, PROTECTED :: RG,RV,RD,RCPD,RLSTT,RLVTT,RTT
+REAL, SAVE, PROTECTED :: a_ratqs_cv
+REAL, SAVE, PROTECTED :: tau_var
+REAL, SAVE, PROTECTED :: fac_tau
+REAL, SAVE, PROTECTED :: tau_cumul
+REAL, SAVE, PROTECTED :: a_ratqs_wake
+INTEGER, SAVE, PROTECTED :: dqimpl
+
+real, allocatable, SAVE :: povariance(:,:)
+!$OMP THREADPRIVATE(povariance)
+real, allocatable, SAVE :: var_conv(:,:)
+!$OMP THREADPRIVATE(var_conv)
+
+!$OMP THREADPRIVATE(a_ratqs_cv,tau_var,fac_tau,tau_cumul,a_ratqs_wake)
+
+
+CONTAINS
+SUBROUTINE ratqs_ini(klon,klev,iflag_thermals,lunout_in,nbsrf_in,is_lic_in,is_ter_in,RG_in,RV_in,RD_in,RCPD_in,RLSTT_in,RLVTT_in,RTT_in)
+
+   USE ioipsl_getin_p_mod, ONLY : getin_p
+
+integer, intent(in) :: klon,klev,iflag_thermals,lunout_in
+integer, intent(in) :: nbsrf_in,is_lic_in,is_ter_in
+real, intent(in) :: RG_in,RV_in,RD_in,RCPD_in,RLSTT_in,RLVTT_in,RTT_in
+
+integer :: i,k
+
+nbsrf   =nbsrf_in
+is_lic  =is_lic_in
+is_ter  =is_ter_in
+RG      =RG_in    
+RV      =RV_in
+RD      =RD_in
+RCPD    =RCPD_in
+RLSTT   =RLSTT_in
+RLVTT   =RLVTT_in
+RTT     =RTT_in
+
+print*,'ratqs_ini'
+if (iflag_thermals==15.or.iflag_thermals==16) then
+   dqimpl=-1
+else
+   dqimpl=1
+endif
+lunout=lunout_in
+
+allocate(povariance(klon,klev))
+allocate(var_conv(klon,klev))
+
+
+tau_var = 700. ! temps de relaxation de la variabilité pour les thermiques
+fac_tau = 2. 
+tau_cumul = 1000 
+a_ratqs_wake = 3.    ! paramètre pilotant l'importance du terme dépendant des poches froides
+a_ratqs_cv = 1.
+CALL getin_p('tau_var', tau_var)
+CALL getin_p('fac_tau', fac_tau)
+CALL getin_p('tau_cumul', tau_cumul)
+CALL getin_p('a_ratqs_wake', a_ratqs_wake) 
+CALL getin_p('a_ratqs_cv', a_ratqs_cv)         
+do k=1,klev
+  do i=1,klon
+   povariance(i,k) = 0.00001
+   var_conv(i,k) = 0.
+  enddo
+enddo
+
+ RETURN
+
+END SUBROUTINE ratqs_ini
+END MODULE lmdz_ratqs_ini
Index: LMDZ6/trunk/libf/phylmd/lmdz_ratqs_main.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_ratqs_main.F90	(revision 4613)
+++ LMDZ6/trunk/libf/phylmd/lmdz_ratqs_main.F90	(revision 4613)
@@ -0,0 +1,285 @@
+MODULE lmdz_ratqs_main
+
+CONTAINS
+
+SUBROUTINE ratqs_main(klon,klev,nbsrf,prt_level,lunout,       &
+           iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
+           ratqsbas,ratqshaut,ratqsp0,ratqsdp, &
+           pctsrf,s_pblh,zstd, &
+           tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
+           ptconv,ptconvth,clwcon0th, rnebcon0th,       &
+           paprs,pplay,t_seri,q_seri,                   &
+           qtc_cv, sigt_cv,detrain_cv,fm_cv,fqd,fqcomp,sigd,zqsat,             &
+           omega,tke,tke_dissip,lmix,wprime, &
+           t2m,q2m,fm_therm,entr_therm,detr_therm,cell_area,&
+           ratqs,ratqsc,ratqs_inter_)
+
+
+USE lmdz_ratqs_multi,   ONLY: ratqs_inter, ratqs_oro, ratqs_hetero, ratqs_tke
+
+implicit none
+
+!========================================================================
+! Computation of ratqs, the width of the subrid scale water distribution
+! (normalized by the mean value)
+! Various options controled by flags iflag_con and iflag_ratqs
+! F Hourdin 2012/12/06
+!========================================================================
+
+! Declarations
+
+! Input
+integer,intent(in) :: klon,klev,nbsrf,prt_level,lunout
+integer,intent(in) :: iflag_con,iflag_cld_th,iflag_ratqs
+real,intent(in) :: pdtphys,ratqsbas,ratqshaut,fact_cldcon,tau_ratqs
+real,intent(in) :: ratqsp0, ratqsdp
+real, dimension(klon,klev),intent(in) :: omega
+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
+real, dimension(klon,klev),intent(in) :: entr_therm,detr_therm,qtc_cv, sigt_cv
+real, dimension(klon,klev) :: detrain_cv,fm_cv,fqd,fqcomp
+real, dimension(klon) :: sigd
+
+real, dimension(klon,klev+1),intent(in) :: fm_therm
+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
+real, dimension(klon), intent(in) :: cell_area
+real, dimension(klon,nbsrf),intent(in) :: pctsrf
+real, dimension(klon),intent(in) :: s_pblh
+real, dimension(klon),intent(in) :: zstd
+
+! Output
+real, dimension(klon,klev),intent(inout) :: ratqs,ratqsc,ratqs_inter_
+
+logical, dimension(klon,klev),intent(inout) :: ptconvth
+
+! local
+integer i,k
+real, dimension(klon,klev) :: ratqss
+real facteur,zfratqs1,zfratqs2
+real, dimension(klon,klev) :: ratqs_hetero_,ratqs_oro_,ratqs_tke_
+real resol,resolmax,fact
+
+!-------------------------------------------------------------------------
+!  Caclul des ratqs
+!-------------------------------------------------------------------------
+
+!      print*,'calcul des ratqs'
+!   ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
+!   ----------------
+!   on ecrase le tableau ratqsc calcule par clouds_gno
+      if (iflag_cld_th.eq.1) then
+         do k=1,klev
+         do i=1,klon
+            if(ptconv(i,k)) then
+              ratqsc(i,k)=ratqsbas &
+              +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
+            else
+               ratqsc(i,k)=0.
+            endif
+         enddo
+         enddo
+
+!-----------------------------------------------------------------------
+!  par nversion de la fonction log normale
+!-----------------------------------------------------------------------
+      else if (iflag_cld_th.eq.4) then
+         ptconvth(:,:)=.false.
+         ratqsc(:,:)=0.
+         if(prt_level.ge.9) print*,'avant clouds_gno thermique'
+         call clouds_gno &
+         (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
+         if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
+       
+       endif
+
+!   ratqs stables
+!   -------------
+
+      if (iflag_ratqs.eq.0) then
+
+! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
+         do k=1,klev
+            do i=1, klon
+               ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)* &
+               min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.) 
+            enddo 
+         enddo
+
+! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de 
+! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
+! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
+! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
+! Il s'agit de differents tests dans la phase de reglage du modele
+! avec thermiques.
+
+      else if (iflag_ratqs.eq.1) then
+
+         do k=1,klev
+            do i=1, klon
+               if (pplay(i,k).ge.60000.) then
+                  ratqss(i,k)=ratqsbas
+               else if ((pplay(i,k).ge.30000.).and.(pplay(i,k).lt.60000.)) then
+                  ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*(60000.-pplay(i,k))/(60000.-30000.)
+               else
+                  ratqss(i,k)=ratqshaut
+               endif
+            enddo
+         enddo
+
+      else if (iflag_ratqs.eq.2) then
+
+         do k=1,klev
+            do i=1, klon
+               if (pplay(i,k).ge.60000.) then
+                  ratqss(i,k)=ratqsbas*(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
+               else if ((pplay(i,k).ge.30000.).and.(pplay(i,k).lt.60000.)) then
+                    ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*(60000.-pplay(i,k))/(60000.-30000.)
+               else
+                    ratqss(i,k)=ratqshaut
+               endif
+            enddo
+         enddo
+
+      else if (iflag_ratqs==3) then
+         do k=1,klev
+           ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas) &
+           *min( ((paprs(:,1)-pplay(:,k))/70000.)**2 , 1. )
+         enddo
+
+      else if (iflag_ratqs==4) then 
+         do k=1,klev
+           ratqss(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) &
+!          *( tanh( (50000.-pplay(:,k))/20000.) + 1.)
+           *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.)
+         enddo
+
+
+      else if (iflag_ratqs==5) then
+! Dependency of ratqs on model resolution
+! Audran, Meryl, Lea, Gwendal and Etienne
+! April 2023 
+        resolmax=sqrt(maxval(cell_area))
+         do k=1,klev
+            do i=1,klon
+              resol=sqrt(cell_area(i))
+              fact=sqrt(resol/resolmax)
+              ratqss(i,k)=ratqsbas*fact+0.5*(ratqshaut-ratqsbas)*fact &
+              *( 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 ratqs_inter(klon,klev,iflag_ratqs,pdtphys,paprs, &
+                       ratqsbas,wake_deltaq,wake_s,q_seri,qtc_cv, sigt_cv, &
+                       fm_therm,entr_therm,detr_therm,detrain_cv,fm_cv,fqd,fqcomp,sigd, &
+                       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
+            print*,'avant ratqs_inter'
+            ! interactive ratqs with several sources
+             call ratqs_inter(klon,klev,iflag_ratqs,pdtphys,paprs, &
+                       ratqsbas,wake_deltaq,wake_s,q_seri,qtc_cv, sigt_cv, &
+                       fm_therm,entr_therm,detr_therm,detrain_cv,fm_cv,fqd,fqcomp,sigd, &
+                       ratqs_inter_)
+             ratqss=ratqss+ratqs_inter_
+          else if (iflag_ratqs .EQ. 12) then
+             ! contribution of surface heterogeneities to ratqs
+             call ratqs_hetero(klon,klev,pctsrf,s_pblh,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 ratqs_oro(klon,klev,pctsrf,zstd,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 ratqs_tke(klon,klev,pdtphys,t_seri,q_seri,zqsat,pplay,paprs,omega,tke,tke_dissip,lmix,wprime,ratqs_tke_)     
+             ratqss=ratqss+ratqs_tke_
+          endif
+          
+      
+      endif
+
+
+!  ratqs final
+!  -----------
+
+      if (iflag_cld_th.eq.1 .or.iflag_cld_th.eq.2.or.iflag_cld_th.eq.4) then
+
+! On ajoute une constante au ratqsc*2 pour tenir compte de 
+! fluctuations turbulentes de petite echelle
+
+         do k=1,klev
+            do i=1,klon
+               if ((fm_therm(i,k)>1.e-10)) then
+                  ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
+               endif
+            enddo
+         enddo
+
+!   les ratqs sont une combinaison de ratqss et ratqsc
+       if(prt_level.ge.9) write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
+
+         if (tau_ratqs>1.e-10) then
+            facteur=exp(-pdtphys/tau_ratqs)
+         else
+            facteur=0.
+         endif
+         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH 22/09/2009
+! La ligne ci-dessous faisait osciller le modele et donnait une solution
+! assymptotique bidon et dépendant fortement du pas de temps.
+!        ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+         ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
+      else if (iflag_cld_th<=6) then
+!   on ne prend que le ratqs stable pour fisrtilp
+         ratqs(:,:)=ratqss(:,:)
+      else
+          zfratqs1=exp(-pdtphys/10800.)
+          zfratqs2=exp(-pdtphys/10800.)
+          do k=1,klev
+             do i=1,klon
+                if (ratqsc(i,k).gt.1.e-10) then
+                   ratqs(i,k)=ratqs(i,k)*zfratqs2+(iflag_cld_th/100.)*ratqsc(i,k)*(1.-zfratqs2)
+                endif
+                ratqs(i,k)=min(ratqs(i,k)*zfratqs1+ratqss(i,k)*(1.-zfratqs1),0.5)
+             enddo
+          enddo
+      endif
+
+
+return
+END SUBROUTINE ratqs_main
+
+END MODULE lmdz_ratqs_main
Index: LMDZ6/trunk/libf/phylmd/lmdz_ratqs_multi.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/lmdz_ratqs_multi.F90	(revision 4613)
+++ LMDZ6/trunk/libf/phylmd/lmdz_ratqs_multi.F90	(revision 4613)
@@ -0,0 +1,483 @@
+MODULE lmdz_ratqs_multi
+
+!=============================================
+! A FAIRE :
+! Traiter le probleme de USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF
+!=============================================
+
+!=============================================
+! 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"
+
+
+      CONTAINS
+
+
+!========================================================================    
+SUBROUTINE ratqs_inter(klon,klev,iflag_ratqs,pdtphys,paprs, &
+           ratqsbas, wake_deltaq, wake_s, q_seri,qtc_cv, sigt_cv,     &
+           fm_therm,entr_therm,detr_therm,detrain_cv,fm_cv,fqd,fqcomp,sigd, &
+           ratqs_inter_)
+
+USE lmdz_ratqs_ini, ONLY : a_ratqs_cv,tau_var,fac_tau,tau_cumul,a_ratqs_wake, dqimpl
+USE lmdz_ratqs_ini, ONLY : RG
+USE lmdz_ratqs_ini, ONLY : povariance, var_conv
+USE lmdz_thermcell_dq,  ONLY : thermcell_dq
+
+implicit none
+
+!========================================================================
+! L. d'Alençon, 25/02/2021
+! Cette subroutine calcule une valeur de ratqsbas interactive  
+! Elle est appelée par la subroutine ratqs lorsque iflag_ratqs = 11. 
+!========================================================================
+
+! Declarations
+! Input
+integer,intent(in) :: klon,klev,iflag_ratqs
+real,intent(in) :: pdtphys,ratqsbas
+real, dimension(klon,klev+1),intent(in) :: paprs
+real, dimension(klon,klev),intent(in) :: wake_deltaq, q_seri,qtc_cv, sigt_cv
+real, dimension(klon),intent(in) :: wake_s
+real, dimension(klon,klev+1),intent(in) :: fm_therm
+real, dimension(klon,klev),intent(in) :: entr_therm,detr_therm,detrain_cv,fm_cv,fqd,fqcomp
+real, dimension(klon),intent(in) :: sigd
+
+! Output
+real, dimension(klon,klev),intent(inout) :: ratqs_inter_
+
+! local
+LOGICAL :: klein = .false.
+LOGICAL :: klein_conv = .true.
+REAL :: taup0 = 70000
+REAL :: taudp = 500
+integer :: lev_out=10
+REAL, DIMENSION (klon,klev) :: zmasse,entr0,detr0,detraincv,dqp,detrain_p,q0,qd0,tau_diss
+REAL, DIMENSION (klon,klev+1) :: fm0
+integer i,k
+real, dimension(klon,klev) :: wake_dq
+
+
+real, dimension(klon) :: max_sigd, max_dqconv,max_sigt
+real, dimension(klon,klev) :: zoa,zocarrea,pdocarreadj,pocarre,po,pdoadj,varq_therm
+
+
+lev_out=0.
+
+print*,'ratqs_inter'
+
+!-----------------------------------------------------------------------
+!   Calcul des masses
+!-----------------------------------------------------------------------
+
+      do k=1,klev
+         zmasse(:,k)=(paprs(:,k)-paprs(:,k+1))/RG
+      enddo
+!-------------------------------------------------------------------------
+!  Caclul du terme de détrainement de la variance pour les thermiques
+!-------------------------------------------------------------------------
+
+! initialisations 
+
+
+      do k=1,klev
+         do i=1,klon
+            tau_diss(i,k)=tau_var +0.5*fac_tau*tau_var*(tanh((taup0-paprs(i,k))/taudp) + 1.)
+         enddo
+      enddo
+       
+                    
+      
+      entr0(:,:) = entr_therm(:,:) 
+      fm0(:,:) = fm_therm(:,:)  
+      detr0(:,:) = detr_therm(:,:) 
+
+! calcul du carré de l'humidité spécifique
+      po(:,:) = q_seri(:,:)
+      call thermcell_dq(klon,klev,dqimpl,pdtphys,fm0,entr0,zmasse,  &
+     &                   po,pdoadj,zoa,lev_out)
+      do k=1,klev
+         do i=1,klon
+            pocarre(i,k)=po(i,k)*po(i,k) + povariance(i,k)
+         enddo
+      enddo
+      call thermcell_dq(klon,klev,dqimpl,pdtphys,fm0,entr0,zmasse,  &
+      &                   pocarre,pdocarreadj,zocarrea,lev_out) 
+
+      
+      do k=1,klev
+         do i=1,klon      
+            varq_therm(i,k)=zocarrea(i,k)-zoa(i,k)*zoa(i,k)
+         enddo
+      enddo
+
+      
+      if (klein) then
+         do k=1,klev-1
+            do i=1,klon
+              qd0(:,:) = 0.0
+              if (sigd(i).ne.0) then
+                qd0(i,k) = fqd(i,k)*tau_cumul/sigd(i) 
+              endif
+            enddo
+         enddo
+         do k=1,klev-1
+            do i=1,klon      
+               povariance(i,k)= (detr0(i,k)*((zoa(i,k)-po(i,k))**2 + (varq_therm(i,k)-povariance(i,k)))/zmasse(i,k) + fm0(i,k+1)*povariance(i,k+1)/zmasse(i,k) -  &
+               fm0(i,k)*povariance(i,k)/zmasse(i,k) +  &
+               a_ratqs_cv*(detrain_cv(i,k)/zmasse(i,k)) +  sigd(i)*(1-sigd(i))*qd0(i,k)**2/tau_cumul  &   
+                + ((povariance(i,k+1)-povariance(i,k))*(fm_cv(i,k)/zmasse(i,k))))*pdtphys + povariance(i,k)
+               povariance(i,k)= povariance(i,k)*exp(-pdtphys/tau_diss(i,k))
+            enddo
+         enddo
+         povariance(:,klev) = povariance(:,klev-1)
+         
+      else ! calcul direct     
+         qd0(:,:) = 0.0
+         q0(:,:) = 0.0
+         do k=1,klev-1
+            do i=1,klon
+               if (sigd(i).ne.0) then
+                 qd0(i,k) = fqd(i,k)*tau_cumul/sigd(i) 
+               endif
+               if (sigt_cv(i,k).ne.0) then
+                 q0(i,k) = fqcomp(i,k)*tau_cumul/sigt_cv(i,k)
+               endif
+            enddo
+         enddo
+         do k=1,klev-1
+            do i=1,klon      
+               povariance(i,k)= (pdocarreadj(i,k)-2.*po(i,k)*pdoadj(i,k) +  &
+               a_ratqs_cv*(sigt_cv(i,k)*(1-sigt_cv(i,k))*q0(i,k)**2/tau_cumul + sigd(i)*(1-sigd(i))*qd0(i,k)**2/tau_cumul) +  &
+               ((povariance(i,k+1)-povariance(i,k))*(fm_cv(i,k)/zmasse(i,k))))*pdtphys + povariance(i,k)
+               povariance(i,k)=povariance(i,k)*exp(-pdtphys/tau_diss(i,k))
+            enddo
+         enddo
+         povariance(:,klev) = povariance(:,klev-1)
+!         fqd(:,:)=sigt_cv(:,:)*(1-sigt_cv(:,:))*q0(:,:)**2/tau_cumul 
+      endif
+      
+!-------------------------------------------------------------------------
+!  Caclul du terme de détrainement de la variance pour la convection (version fausse avec deux calculs de variance indépendants)
+!-------------------------------------------------------------------------
+
+!      if (klein_conv) then
+!        detrain_p(:,:) = 0.
+!        detraincv(:,:) = 0.
+!        dqp(:,:)    = 0. 
+      
+!        do k=1,klev-1
+!          do i=1,klon
+!           dqp(i,k) = q_seri(i,k) - qp(i,k)
+!           detraincv(i,k) = abs(detrain_cv(i,k))
+         
+!         if ((mp(i,k)-mp(i,k+1)).le.0) then
+!            detrain_p(i,k) = (mp(i,k+1)-mp(i,k))
+!         endif
+!          enddo
+!        enddo
+
+!        do k=1,klev-1
+!          do i=1,klon
+!            qd0(:,:) = 0.0
+!            if (sigd(i).ne.0) then
+!              qd0(i,k) = fqd(i,k)*tau_cumul/sigd(i) 
+!            endif
+!          enddo
+!        enddo
+
+!        do k=1,klev-1
+!          do i=1,klon            
+!            var_conv(i,k)= var_conv(i,k)*exp(-pdtphys/tau_conv) +  &
+!                (a_ratqs_cv*(detraincv(i,k)/zmasse(i,k)) +  sigd(i)*(1-sigd(i))*qd0(i,k)**2/tau_cumul  &   
+!                + ((var_conv(i,k+1)-var_conv(i,k))*(fm_cv(i,k)/zmasse(i,k))))* &
+!              ((dqp(i,k)**2)*(detrain_p(i,k)/zmasse(i,k))))* & ! les termes de descentes précipitantes qui seront traités autrement
+!              (((mp(i,k)-mp(i,k+1))/zmasse(i,k))*dqp(i,k)**2) + (2*mp(i,k)*dqp(i,k)*(qp(i,k+1)-qp(i,k))/zmasse(i,k)))* &
+!                (1.-exp(-pdtphys/tau_conv))/(1/tau_conv)      
+!          enddo
+!        enddo
+!        var_conv(:,klev) = var_conv(:,klev-1)
+!        fqd(:,:) = var_conv(:,:)
+!      else
+
+
+!        do k=1,klev-1
+!          do i=1,klon            
+!            var_conv(i,k)= var_conv(i,k)*exp(-pdtphys/tau_conv) +  &
+!                (a_ratqs_cv*(sigt_cv(i,k)*(1-sigt_cv(i,k))*q0(i,k)**2/pdtphys + sigd(i)*(1-sigd(i))*qd0(i,k)**2/tau_cumul) +  &   ! somme des contributions des descentes précipitantes et des flux mélangés
+!                + ((var_conv(i,k+1)-var_conv(i,k))*(fm_cv(i,k)/zmasse(i,k))))* & ! flux compensatoires dans l'environnement
+!                (1.-exp(-pdtphys/tau_conv))/(1/tau_conv) 
+!             ((var_conv(i,k+1)-var_conv(i,k))*(fm_cv(i,k)/zmasse(i,k)))  
+!          enddo
+!        enddo
+!        var_conv(:,klev) = var_conv(:,klev-1) 
+!      endif
+
+!-------------------------------------------------------------------------
+!  Caclul du ratqs_inter_
+!-------------------------------------------------------------------------
+
+      do k=1,klev
+        do i=1,klon
+           if(q_seri(i,k).ge.1E-7) then
+               ratqs_inter_(i,k) = abs(povariance(i,k))**0.5/q_seri(i,k) 
+           else 
+               ratqs_inter_(i,k) = 0.  
+           endif
+        enddo
+      enddo
+      
+return
+end             
+
+!------------------------------------------------------------------
+SUBROUTINE ratqs_oro(klon,klev,pctsrf,zstd,qsat,temp,pplay,paprs,ratqs_oro_)
+
+! Etienne Vignon, November 2021: effect of subgrid orography on ratqs
+
+USE lmdz_ratqs_ini, ONLY : RG,RV,RD,RLSTT,RLVTT,RTT,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,nbsrf) :: pctsrf
+REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat    ! saturation specific humidity [kg/kg]
+REAL, DIMENSION(klon), INTENT(IN) :: zstd    ! sub grid orography standard deviation
+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 ratqs_oro
+
+!=============================================
+
+SUBROUTINE ratqs_hetero(klon,klev,pctsrf,s_pblh,t2m,q2m,temp,q,pplay,paprs,ratqs_hetero_)
+
+! Etienne Vignon, November 2021
+! Effect of subgrid surface heterogeneities on ratqs
+
+USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF
+
+USE lmdz_ratqs_ini, ONLY : RG,RD,RTT,nbsrf
+
+IMPLICIT NONE
+
+! INPUTS
+
+
+INTEGER, INTENT(IN) :: klon                       ! number of horizontal grid points
+INTEGER, INTENT(IN) :: klev                       ! number of vertical layers
+REAL, DIMENSION(klon)                   :: s_pblh ! height of the planetary boundary layer(HPBL)
+REAL, DIMENSION(klon,nbsrf)             :: pctsrf ! Fractional cover of subsurfaces
+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, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT
+REAL, DIMENSION (klon,klev) :: zlay
+
+
+
+! Calculation of near-surface surface ratqs
+!-------------------------------------------
+
+    
+    ratiom(:)=0.
+    xsi0(:)=0.
+    
+    DO nsrf=1,nbsrf
+    CALL CALC_QSAT_ECMWF(klon,t2m(:,nsrf),q2m(:,nsrf),paprs(:,1),RTT,0,.false.,qsat2m,dqsatdT)
+    ratiom(:)=ratiom(:)+pctsrf(:,nsrf)*(q2m(:,nsrf)/qsat2m(:))
+    xsi0(:)=xsi0(:)+pctsrf(:,nsrf)*((q2m(:,nsrf)/qsat2m(:)-ratiom(:))**2)
+    END DO
+    
+    xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6)
+
+
+
+! 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 ratqs_hetero
+
+!=============================================
+
+SUBROUTINE ratqs_tke(klon,klev,pdtphys,temp,q,qsat,pplay,paprs,omega,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 lmdz_ratqs_ini, ONLY : RG,RV,RD,RCPD,RLSTT,RLVTT,RTT
+
+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), INTENT(IN) :: omega       ! 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 ratqs_tke
+
+END MODULE lmdz_ratqs_multi
Index: LMDZ6/trunk/libf/phylmd/nuage.h
===================================================================
--- LMDZ6/trunk/libf/phylmd/nuage.h	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/nuage.h	(revision 4613)
@@ -6,5 +6,5 @@
       REAL tau_cld_cv,coefw_cld_cv
       REAL tmax_fonte_cv
-      INTEGER iflag_cld_cv
+      INTEGER iflag_cld_cv, iflag_ratqs
       INTEGER iflag_rei
       LOGICAL ok_icefra_lscp
@@ -13,5 +13,5 @@
      &                  rei_min,rei_max,                                &
      &                  tau_cld_cv,coefw_cld_cv,                        &
-     &                  tmax_fonte_cv,                                  &
+     &                  tmax_fonte_cv,iflag_ratqs,                      &
      &                  iflag_cld_cv,                                   &
      &                  ok_icefra_lscp,                                 &
Index: LMDZ6/trunk/libf/phylmd/phyetat0_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyetat0_mod.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/phyetat0_mod.F90	(revision 4613)
@@ -28,5 +28,5 @@
        zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, u10m, v10m, treedrg, &
        ale_wake, ale_bl_stat, ds_ns, dt_ns, delta_sst, delta_sal, dter, dser, &
-       dt_ds, ratqs_inter
+       dt_ds, ratqs_inter_
 !FC
   USE geometry_mod,     ONLY: longitude_deg, latitude_deg
@@ -491,5 +491,5 @@
 
 ! fisrtilp/Clouds 0.002 could be ratqsbas. But can stay like this as well
-  found=phyetat0_get(ratqs_inter,"RATQS_INTER","Relative width of the lsc sugrid scale water",0.002)
+  found=phyetat0_get(ratqs_inter_,"RATQS_INTER","Relative width of the lsc sugrid scale water",0.002)
 
 !===========================================
Index: LMDZ6/trunk/libf/phylmd/phyredem.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/phyredem.F90	(revision 4613)
@@ -30,5 +30,5 @@
                                 du_gwd_rando, du_gwd_front, u10m, v10m, &
                                 treedrg, solswfdiff, delta_sal, ds_ns, dt_ns, &
-                                delta_sst, ratqs_inter, dter, dser, dt_ds
+                                delta_sst, ratqs_inter_, dter, dser, dt_ds
 
   USE geometry_mod, ONLY : longitude_deg, latitude_deg
@@ -334,5 +334,5 @@
 
     ! fisrtilp/clouds
-    CALL put_field(pass,"RATQS_INTER","Relative width of the lsc sugrid scale water",ratqs_inter)
+    CALL put_field(pass,"RATQS_INTER","Relative width of the lsc sugrid scale water",ratqs_inter_)
 
 
Index: LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/phys_state_var_mod.F90	(revision 4613)
@@ -102,6 +102,6 @@
       REAL, ALLOCATABLE, SAVE :: rneb_ancien(:,:)
 !$OMP THREADPRIVATE(rneb_ancien)
-      REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:)
-!$OMP THREADPRIVATE(qtc_cv,sigt_cv)
+      REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:),detrain_cv(:,:),fm_cv(:,:)
+!$OMP THREADPRIVATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
       REAL, ALLOCATABLE, SAVE :: ratqs(:,:)
 !$OMP THREADPRIVATE(ratqs)
@@ -229,6 +229,6 @@
 !$OMP THREADPRIVATE(ftd)
 ! fqd : convective moistening due to unsaturated downdraughts
-      REAL,ALLOCATABLE,SAVE :: fqd(:,:)     
-!$OMP THREADPRIVATE(fqd)
+      REAL,ALLOCATABLE,SAVE :: fqd(:,:),fqcomp(:,:)     
+!$OMP THREADPRIVATE(fqd,fqcomp)
 #ifdef ISO
       REAL, ALLOCATABLE, SAVE :: fxtd(:,:,:)
@@ -454,6 +454,6 @@
 !$OMP THREADPRIVATE(ale_bl_trig)
 
-      REAL,SAVE,ALLOCATABLE :: ratqs_inter(:,:)
-!$OMP THREADPRIVATE(ratqs_inter)
+      REAL,SAVE,ALLOCATABLE :: ratqs_inter_(:,:)
+!$OMP THREADPRIVATE(ratqs_inter_)
 
 #ifdef ISO
@@ -567,5 +567,5 @@
       ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
       ALLOCATE(rneb_ancien(klon,klev))
-      ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev))
+      ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev),detrain_cv(klon,klev),fm_cv(klon,klev))
       ALLOCATE(ratqs(klon,klev))
       ALLOCATE(pbl_tke(klon,klev+1,nbsrf+1))
@@ -638,5 +638,5 @@
       ALLOCATE(sigd(klon))
       ALLOCATE(cin(klon), ALE(klon), ALP(klon))
-      ALLOCATE(ftd(klon,klev), fqd(klon,klev))
+      ALLOCATE(ftd(klon,klev), fqd(klon,klev),fqcomp(klon,klev))
       ALLOCATE(Ale_bl(klon))
       ALLOCATE(ale_wake(klon))
@@ -733,5 +733,5 @@
 
       ALLOCATE(ale_bl_trig(klon))
-      ALLOCATE(ratqs_inter(klon,klev))
+      ALLOCATE(ratqs_inter_(klon,klev))
       IF (ok_gwd_rando) THEN
         ALLOCATE(du_gwd_rando(klon, klev))
@@ -773,5 +773,5 @@
       DEALLOCATE(qs_ancien, ql_ancien, qbs_ancien, rneb_ancien)
       DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien)
-      DEALLOCATE(qtc_cv,sigt_cv)
+      DEALLOCATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
       DEALLOCATE(u_ancien, v_ancien)
       DEALLOCATE(tr_ancien)                           !RomP
@@ -821,5 +821,5 @@
       DEALLOCATE(wd, sigd)
       DEALLOCATE(cin, ALE, ALP)
-      DEALLOCATE(ftd, fqd)
+      DEALLOCATE(ftd, fqd, fqcomp)
       DEALLOCATE(Ale_bl, Alp_bl)
       DEALLOCATE(ale_wake)
@@ -886,5 +886,5 @@
       if (.not. ok_hines .and. ok_gwd_rando) DEALLOCATE(du_gwd_front)
       DEALLOCATE(ale_bl_trig)
-      DEALLOCATE(ratqs_inter)
+      DEALLOCATE(ratqs_inter_)
 
       if (activate_ocean_skin >= 1) then
Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4608)
+++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90	(revision 4613)
@@ -87,4 +87,6 @@
     USE blowing_snow_ini_mod, ONLY : blowing_snow_ini , qbst_bs 
     USE lscp_ini_mod, ONLY : lscp_ini
+    USE lmdz_ratqs_main, ONLY : ratqs_main
+    USE lmdz_ratqs_ini, ONLY : ratqs_ini
 
     !USE cmp_seri_mod
@@ -458,7 +460,7 @@
     SAVE seuil_inversion
     !$OMP THREADPRIVATE(seuil_inversion)
-    INTEGER iflag_ratqs
-    SAVE iflag_ratqs
-    !$OMP THREADPRIVATE(iflag_ratqs)
+    
+    
+    
     real facteur
 
@@ -1334,5 +1336,5 @@
             solarlong0,seuil_inversion, &
             fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
-            iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
+            iflag_cld_th,ratqsbas,ratqshaut,tau_ratqs, &
             ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat, aerosol_couple, &
             chemistry_couple, flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
@@ -1801,4 +1803,5 @@
        CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
    &    RG,RD,RCPD,RKAPPA,RLVTT,RETV)
+       CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT)
        CALL lscp_ini(pdtphys,ok_ice_sursat, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
        CALL blowing_snow_ini(prt_level,lunout, &
@@ -3126,5 +3129,5 @@
                ALE,ALP, &
                sig1,w01, &
-               d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
+               d_t_con,d_q_con,fqcomp,d_u_con,d_v_con,d_tr, &
                rain_con, snow_con, ibas_con, itop_con, sigd, &
                ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
@@ -3137,5 +3140,5 @@
                ftd,fqd,lalim_conv,wght_th, &
                ev, ep,epmlmMm,eplaMm, &
-               wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
+               wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv,detrain_cv, &
                tau_cld_cv,coefw_cld_cv,epmax_diag)
 
@@ -3149,4 +3152,5 @@
           clwcon0=qcondc
           pmfu(:,:)=upwd(:,:)+dnwd(:,:)
+          fm_cv(:,:)=upwd(:,:)+dnwd(:,:)+dnwd0(:,:)
           !
           !jyg<
@@ -3699,12 +3703,15 @@
     ENDDO
 
-    CALL  calcratqs(klon,klev,prt_level,lunout,        &
+    CALL ratqs_main(klon,klev,nbsrf,prt_level,lunout,        &
          iflag_ratqs,iflag_con,iflag_cld_th,pdtphys,  &
          ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
+         pctsrf,s_pblh,zstd, &
          tau_ratqs,fact_cldcon,wake_s, wake_deltaq,   &
          ptconv,ptconvth,clwcon0th, rnebcon0th,     &
-         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,cell_area, &
-         ratqs,ratqsc,ratqs_inter)
+         paprs,pplay,t_seri,q_seri, &
+         qtc_cv, sigt_cv,detrain_cv,fm_cv,fqd,fqcomp,sigd,zqsat, &
+         omega,pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave, &
+         t2m,q2m,fm_therm,entr_therm,detr_therm,cell_area, &
+         ratqs,ratqsc,ratqs_inter_)
 
     !
@@ -5213,6 +5220,4 @@
     !ENDIF    ! (iflag_phytrac=1)
 
-
-#ifndef INCA 
     IF (offline) THEN
 
@@ -5230,5 +5235,5 @@
 
     ENDIF
-#endif
+
     !
     ! Calculer le transport de l'eau et de l'energie (diagnostique)
@@ -5549,5 +5554,5 @@
          ok_sync, ptconv, read_climoz, clevSTD,          &
          ptconvth, d_u, d_t, qx, d_qx, zmasse,           &
-         flag_aerosol, flag_aerosol_strat, ok_cdnc, t, u1, v1 )
+         flag_aerosol, flag_aerosol_strat, ok_cdnc,t, u1, v1)
 #endif
 
