Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_cine.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_cine.F	(revision 1126)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_cine.F	(revision 1127)
@@ -33,8 +33,10 @@
       integer ifst(nloc),isublcl(nloc)
       logical lswitch(nloc),lswitch1(nloc),lswitch2(nloc)
+      logical exist_lfc(nloc)
+      real plfc(nloc)
       real dpmax
       real deltap,dcin
       real buoylcl(nloc),tvplcl(nloc),tvlcl(nloc)
-      real plfc(nloc),p0(nloc)
+      real p0(nloc)
       real buoyz(nloc), buoy(nloc,nd)
 c
@@ -50,32 +52,10 @@
 c      Recompute buoyancies
 c--------------------------------------------------------------
-      DO k = 1,nl
+      DO k = 1,nd
         DO il = 1,ncum
+!      print*,'tvp tv=',tvp(il,k),tv(il,k)
           buoy(il,k) = tvp(il,k) - tv(il,k)
         ENDDO
       ENDDO
-c
-c---------------------------------------------------------------
-c premiere couche contenant un  niveau de flotabilite positive
-c et premiere couche contenant un  niveau de flotabilite negative
-c  au dessus du niveau de condensation
-c---------------------------------------------------------------
-      do il = 1,ncum
-        itop(il) =nl-1
-        ineg(il) = nl-1
-      enddo
-      do 100 k=nl,1,-1
-       do 110 il=1,ncum
-        if (k .ge. icb(il)) then
-         if (buoy(il,k) .gt. 0.) then
-          itop(il)=k
-         else
-          ineg(il)=k
-         endif
-        endif
-110    continue
-100   continue
-c      print *,' itop, ineg, icb ',itop(1),ineg(1), icb(1)
-c
 c---------------------------------------------------------------
 c
@@ -109,4 +89,40 @@
 c
 c---------------------------------------------------------------
+c premiere couche contenant un  niveau de flotabilite positive
+c et premiere couche contenant un  niveau de flotabilite negative
+c  au dessus du niveau de condensation
+c---------------------------------------------------------------
+      do il = 1,ncum
+        itop(il) =nl-1
+        ineg(il) = nl-1
+        exist_lfc(il) = .FALSE.
+      enddo
+      do 100 k=nl-1,1,-1
+       do 110 il=1,ncum
+        if (k .ge. ifst(il)) then
+         if (buoy(il,k) .gt. 0.) then
+          itop(il)=k
+          exist_lfc(il) = .TRUE.
+         else
+          ineg(il)=k
+         endif
+        endif
+110    continue
+100   continue
+c
+c---------------------------------------------------------------
+c When there is no positive buoyancy level, set Plfc, Cina and Cinb
+c to arbitrary extreme values.
+c---------------------------------------------------------------
+      DO il = 1,ncum
+       IF (.NOT.exist_lfc(il)) THEN
+         Plfc(il) = 1.111
+         Cinb(il) = -1111.
+         Cina(il) = -1112.
+       ENDIF
+      ENDDO
+c
+c
+c---------------------------------------------------------------
 c -- Two cases : BUOYlcl >= 0 and BUOYlcl < 0.
 c---------------------------------------------------------------
@@ -118,5 +134,5 @@
       DPMAX = 50.
       DO il = 1,ncum
-        lswitch1(il)=BUOYlcl(il) .GE. 0.
+        lswitch1(il)=BUOYlcl(il) .GE. 0. .AND. exist_lfc(il)
         lswitch(il) = lswitch1(il)
       ENDDO
@@ -233,5 +249,5 @@
 C
       DO il = 1,ncum
-        lswitch1(il)=BUOYlcl(il) .LT. 0.
+        lswitch1(il)=BUOYlcl(il) .LT. 0. .AND. exist_lfc(il)
         lswitch(il) = lswitch1(il)
       ENDDO
@@ -239,5 +255,5 @@
 c 2.0.1 Premiere  couche ou la flotabilite est negative au dessus du sol
 c ----------------------------------------------------
-c    au cas ou il existe  sinon ilow=1 (nk apres)
+c    au cas ou elle existe  sinon ilow=1 (nk apres)
 c      on suppose que la parcelle part de la premiere couche
 c
@@ -248,5 +264,5 @@
       ENDDO
 c
-      do 200 i=nl,1,-1
+      do 200 k=nl,1,-1
         DO il = 1,ncum
         IF (lswitch(il) .AND. k .LE.icb(il)-1) THEN
@@ -292,9 +308,6 @@
         dcin = RD*(BUOYz(il)+BUOYlcl(il))*deltap/(P0(il)+Plcl(il))
         CINB(il) = min(0.,dcin)
-cc        print *,'buoyz(il),buoylcl(il),deltap,p0(il),plcl(il),dcin ',
-cc     $           buoyz(il),buoylcl(il),deltap,p0(il),plcl(il),dcin
-      ENDIF
-      ENDDO
-c        print*, 'CINB ',CINB(1),'DCIN ',DCIN,I,BUOYz(1),BUOYlcl(1)
+      ENDIF
+      ENDDO
 c
       DO il = 1,ncum
@@ -316,8 +329,4 @@
       ENDDO
 c
-      IF (lswitch(1)) THEN
-c        print*,'ilow= ',ilow(1),'DCIN0 ',DCIN,P0(1),P(1,ilow(1))
-c        print*,'buoy',(BUOY(1,k),k=1,itop(1))
-      ENDIF
 c
 C  Middle part of CINB : integral from P(ilow) to P(isublcl)
@@ -332,5 +341,4 @@
         ENDIF
         ENDDO
-c      print*, 'CINB ', CINB(1), 'DCIN',DCIN,k,BUOY(1,k),BUOY(1,k+1)
       ENDDO
 c
@@ -345,5 +353,4 @@
       ENDDO
 C
-c        print*, ' CINB ', CINB(1), 'Dcin ',dcin
 c
 cc      ENDIF
@@ -439,8 +446,4 @@
       ENDDO
 cc      ENDIF
-c       Print *,' Plcl,P(itop-1),P(itop),PLFC,BUOYlcl'
-c     $      ,Plcl(1),P(1,itop(1)-1),P(1,itop(1)),PLFC(1),BUOYlcl(1)
-C
-c        print*, 'CIN above', CINA(1), 'CIN below',CINB(1)
 c
 
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_routines.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_routines.F	(revision 1126)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_routines.F	(revision 1127)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Header: /home/cvsroot/LMDZ4/libf/phylmd/cv3_routines.F,v 1.16 2008-11-06 16:29:35 lmdzadmin Exp $
 !
 c
@@ -120,5 +120,7 @@
 
 c ori      do 110 k=1,nlp
-      do 110 k=1,nl ! convect3
+! abderr     do 110 k=1,nl ! convect3
+       do 110 k=1,nlp
+      
         do 100 i=1,len
 cdebug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
@@ -2256,5 +2258,5 @@
       SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
      :                    ,icb,inb,delt
-     :                    ,t,rr,t_wake,rr_wake,u,v,tra
+     :                    ,t,rr,t_wake,rr_wake,s_wake,u,v,tra
      :                    ,gz,p,ph,h,hp,lv,cpn,th,th_wake
      :                    ,ep,clw,m,tp,mp,rp,up,vp,trap
@@ -2283,4 +2285,5 @@
       real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
       real t_wake(nloc,nd), rr_wake(nloc,nd)
+      real s_wake(nloc)
       real tra(nloc,nd,ntra), sig(nloc,nd)
       real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
@@ -2327,4 +2330,5 @@
       real esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
       real th_wake(nloc,nd)
+      real alpha_qpos(nloc)
       real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld
       real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld
@@ -2961,8 +2965,9 @@
 c   ***          integrated enthalpy and water tendencies         ***
 c
+c Correction bug le 18-03-09
       do 503 il=1,ncum
       IF (iflag(il) .le. 1) THEN
-      ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il))
-     : +t(il,inb(il))*(cpv-cpd)
+       ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))
+     : -h(il,inb(il))+t(il,inb(il))*(cpv-cpd)
      : *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
      :  /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
@@ -3065,4 +3070,67 @@
        enddo
       enddo
+
+c
+c   ***   Check that moisture stays positive. If not, scale tendencies
+c        in order to ensure moisture positivity
+      DO il = 1,ncum
+       IF (iflag(il) .le. 1) THEN
+        alpha_qpos(il) = max(1. , -delt*fr(il,1)/
+     :     (s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
+       ENDIF
+      ENDDO
+      DO i = 2,nl
+       DO il = 1,ncum
+        IF (iflag(il) .le. 1) THEN
+        alpha_qpos(il) = max(alpha_qpos(il) , -delt*fr(il,i)/
+     :     (s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
+        ENDIF
+       ENDDO
+      ENDDO
+      DO il = 1,ncum
+       IF (iflag(il) .le. 1 .and. alpha_qpos(il) .gt. 1.001) THEN
+        alpha_qpos(il) = alpha_qpos(il)*1.1
+       ENDIF
+      ENDDO
+      DO il = 1,ncum
+       IF (iflag(il) .le. 1) THEN
+        sigd(il) = sigd(il)/alpha_qpos(il)
+        precip(il) = precip(il)/alpha_qpos(il)
+       ENDIF
+      ENDDO
+      DO i = 1,nl
+       DO il = 1,ncum
+        IF (iflag(il) .le. 1) THEN
+         fr(il,i) = fr(il,i)/alpha_qpos(il)
+         ft(il,i) = ft(il,i)/alpha_qpos(il)
+         fqd(il,i) = fqd(il,i)/alpha_qpos(il)
+         ftd(il,i) = ftd(il,i)/alpha_qpos(il)
+         fu(il,i) = fu(il,i)/alpha_qpos(il)
+         fv(il,i) = fv(il,i)/alpha_qpos(il)
+         m(il,i) = m(il,i)/alpha_qpos(il)
+         mp(il,i) = mp(il,i)/alpha_qpos(il)
+         Vprecip(il,i) = Vprecip(il,i)/alpha_qpos(il)
+        ENDIF
+       ENDDO
+      ENDDO
+      DO i = 1,nl
+      DO j = 1,nl
+       DO il = 1,ncum
+        IF (iflag(il) .le. 1) THEN
+         ment(il,i,j) = ment(il,i,j)/alpha_qpos(il)
+        ENDIF
+       ENDDO
+      ENDDO
+      ENDDO
+      DO j = 1,ntra
+      DO i = 1,nl
+       DO il = 1,ncum
+        IF (iflag(il) .le. 1) THEN
+         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
+        ENDIF
+       ENDDO
+      ENDDO
+      ENDDO
+
 c
 c   ***           reset counter and return           ***
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3a_compress.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3a_compress.F	(revision 1126)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3a_compress.F	(revision 1127)
@@ -3,5 +3,6 @@
      :    ,plcl1,tnk1,qnk1,gznk1,hnk1,unk1,vnk1
      :    ,wghti1,pbase1,buoybase1
-     :    ,t1,q1,qs1,t1_wake,q1_wake,qs1_wake,u1,v1,gz1,th1,th1_wake
+     :    ,t1,q1,qs1,t1_wake,q1_wake,qs1_wake,s1_wake
+     :    ,u1,v1,gz1,th1,th1_wake
      :    ,tra1
      :    ,h1     ,lv1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
@@ -12,5 +13,6 @@
      o    ,plcl,tnk,qnk,gznk,hnk,unk,vnk
      o    ,wghti,pbase,buoybase
-     o    ,t,q,qs,t_wake,q_wake,qs_wake,u,v,gz,th,th_wake
+     o    ,t,q,qs,t_wake,q_wake,qs_wake,s_wake
+     o    ,u,v,gz,th,th_wake
      o    ,tra
      o    ,h     ,lv     ,cpn    ,p,ph,tv    ,tp,tvp,clw
@@ -39,4 +41,5 @@
       real t1(len,nd),q1(len,nd),qs1(len,nd)
       real t1_wake(len,nd),q1_wake(len,nd),qs1_wake(len,nd)
+      real s1_wake(len)
       real u1(len,nd),v1(len,nd)
       real gz1(len,nd),th1(len,nd),th1_wake(len,nd)
@@ -58,4 +61,5 @@
       real t(len,nd),q(len,nd),qs(len,nd)
       real t_wake(len,nd),q_wake(len,nd),qs_wake(len,nd)
+      real s_wake(len)
       real u(len,nd),v(len,nd)
       real gz(len,nd),th(len,nd),th_wake(len,nd)
@@ -131,4 +135,5 @@
       if(iflag1(i).eq.0)then
       nn=nn+1
+      s_wake(nn)=s1_wake(i)
       iflag(nn)=iflag1(i)
       nk(nn)=nk1(i)
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cva_driver.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cva_driver.F	(revision 1126)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/cva_driver.F	(revision 1127)
@@ -2,5 +2,5 @@
      &                   iflag_con,iflag_mix,
      &                   iflag_clos,delt,
-     &                   t1,q1,qs1,t1_wake,q1_wake,qs1_wake,
+     &                   t1,q1,qs1,t1_wake,q1_wake,qs1_wake,s1_wake,
      &                   u1,v1,tra1,
      &                   p1,ph1,
@@ -50,4 +50,5 @@
 C      q1_wake       Real           Input        specific hum(unsat draught envt)
 C      qs1_wake      Real           Input        sat specific hum(unsat draughts envt)
+C      s1_wake       Real           Input        fractionnal area covered by wakes
 C      u1            Real           Input        u-wind
 C      v1            Real           Input        v-wind
@@ -121,4 +122,5 @@
       real q1_wake(len,nd)
       real qs1_wake(len,nd)
+      real s1_wake(len)
       real u1(len,nd)
       real v1(len,nd)
@@ -198,4 +200,6 @@
 !       Must be defined at same grid levels as T.
 !
+!s_wake: Array of fractionnal area occupied by the wakes.
+!
 !  u:   Array of zonal wind velocity (m/s) of dimension ND, witth first
 !       index corresponding with the lowest model level. Defined at
@@ -358,4 +362,5 @@
       real t(nloc,klev),q(nloc,klev),qs(nloc,klev)
       real t_wake(nloc,klev),q_wake(nloc,klev),qs_wake(nloc,klev)
+      real s_wake(nloc)
       real u(nloc,klev),v(nloc,klev)
       real gz(nloc,klev),h(nloc,klev)     ,hm(nloc,klev)
@@ -531,5 +536,5 @@
         print*,'Emanuel version 3 nouvelle'
        endif 
-
+!       print*,'t1, q1 ',t1,q1
        CALL cv3_prelim(len,nd,ndp1,t1,q1,p1,ph1      ! nd->na
      o               ,lv1,cpn1,tv1,gz1,h1,hm1,th1)
@@ -668,10 +673,12 @@
 
       if (iflag_con.eq.3) then
-     
+!       print*,'ncum tv1 ',ncum,tv1 
+!       print*,'tvp1 ',tvp1
        CALL cv3a_compress( len,nloc,ncum,nd,ntra
      :    ,iflag1,nk1,icb1,icbs1
      :    ,plcl1,tnk1,qnk1,gznk1,hnk1,unk1,vnk1
      :    ,wghti1,pbase1,buoybase1
-     :    ,t1,q1,qs1,t1_wake,q1_wake,qs1_wake,u1,v1,gz1,th1,th1_wake
+     :    ,t1,q1,qs1,t1_wake,q1_wake,qs1_wake,s1_wake
+     :    ,u1,v1,gz1,th1,th1_wake
      :    ,tra1
      :    ,h1     ,lv1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
@@ -682,5 +689,6 @@
      o    ,plcl,tnk,qnk,gznk,hnk,unk,vnk
      o    ,wghti,pbase,buoybase
-     o    ,t,q,qs,t_wake,q_wake,qs_wake,u,v,gz,th,th_wake
+     o    ,t,q,qs,t_wake,q_wake,qs_wake,s_wake
+     o    ,u,v,gz,th,th_wake
      o    ,tra
      o    ,h     ,lv     ,cpn    ,p,ph,tv    ,tp,tvp,clw
@@ -688,4 +696,7 @@
      o    ,sig,w0,ptop2
      o    ,Ale,Alp  ) 
+
+!       print*,'tv ',tv
+!       print*,'tvp ',tvp
 
       endif
@@ -856,5 +867,5 @@
        CALL cv3_yield(nloc,ncum,nd,nd,ntra            ! na->nd
      :                     ,icb,inb,delt
-     :                     ,t,q,t_wake,q_wake,u,v,tra
+     :                     ,t,q,t_wake,q_wake,s_wake,u,v,tra
      :                     ,gz,p,ph,h,hp,lv,cpn,th,th_wake
      :                     ,ep,clw,m,tp,mp,qp,up,vp,trap
Index: /LMDZ4/branches/LMDZ4-dev/libf/phylmd/wake.F
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/phylmd/wake.F	(revision 1126)
+++ /LMDZ4/branches/LMDZ4-dev/libf/phylmd/wake.F	(revision 1127)
@@ -12,4 +12,5 @@
      o                ,Cstar,d_deltat_gw
      o                ,d_deltatw2,d_deltaqw2)
+
 
 ***************************************************************
@@ -157,8 +158,8 @@
       REAL alpk
       REAL delta_t_min
-      REAL Pupper
       INTEGER nsub
       REAL dtimesub
       REAL sigmad, hwmin
+      REAL :: sigmaw_max
 cIM 080208
       LOGICAL, dimension(klon) :: gwake
@@ -183,4 +184,13 @@
       INTEGER, DIMENSION(klon) :: ktop, kupper
 
+c Sub-timestep tendencies and related variables
+       REAL d_deltatw(klon,klev),d_deltaqw(klon,klev)
+       REAL d_te(klon,klev),d_qe(klon,klev)
+       REAL d_sigmaw(klon),alpha(klon)
+       REAL q0_min(klon),q1_min(klon)
+       LOGICAL wk_adv(klon), OK_qx_qw(klon)
+       REAL epsilon
+       DATA epsilon/1.e-9/
+
 c Autres variables internes
       INTEGER isubstep, k, i
@@ -202,7 +212,8 @@
       REAL, DIMENSION(klon,klev) :: the, thu
 
-      REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
+!      REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
 
       REAL, DIMENSION(klon,klev+1) :: omgbw 
+      REAL, DIMENSION(klon) :: pupper
       REAL, DIMENSION(klon) :: omgtop
       REAL, DIMENSION(klon,klev) :: dp_omgbw
@@ -279,4 +290,8 @@
 	dqls(i,k) = 0.
         d_deltat_gw(i,k)=0.
+        d_te(i,k) = 0.
+        d_qe(i,k) = 0.
+        d_deltatw(i,k) = 0.
+        d_deltaqw(i,k) = 0.
 !IM 060508 beg
         d_deltatw2(i,k)=0.
@@ -294,4 +309,8 @@
       sigmaw(i) = amin1(sigmaw(i),0.99)
       sigmaw0(i) = sigmaw(i)
+      wape(i) = 0.
+      wape2(i) = 0.
+      d_sigmaw(i) = 0.
+      ktopw(i) = 0
       ENDDO
 C
@@ -406,6 +425,11 @@
 c
 C       Pupper = 50000.  ! melting level
-       Pupper = 60000.
+c       Pupper = 60000.
 c       Pupper = 80000.  ! essais pour case_e
+       DO i = 1,klon
+ccc       Pupper(i) = 0.6*ph(i,1)
+        Pupper(i) = 60000.
+       ENDDO
+
 C
 C    Determine Wake top pressure (Ptop) from buoyancy integral
@@ -481,5 +505,5 @@
       DO i=1,klon
         IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
-        IF (ph(i,k+1) .lt. pupper) kupper(i)=k
+        IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
       ENDDO
       ENDDO
@@ -622,6 +646,27 @@
       ENDIF
       ENDDO
-c
-C
+
+c
+c Check qx and qw positivity
+c --------------------------
+      DO i = 1,klon
+       q0_min(i)=min(  (qe(i,1)-sigmaw(i)*deltaqw(i,1)),
+     $              (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1))  )
+      ENDDO
+      DO k = 2,klev
+      DO i = 1,klon
+        q1_min(i)=min(  (qe(i,k)-sigmaw(i)*deltaqw(i,k)),
+     $              (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k))  )
+        IF (q1_min(i).le.q0_min(i)) THEN
+          q0_min(i)=q1_min(i)
+        ENDIF
+      ENDDO
+      ENDDO
+c
+      DO i = 1,klon
+       OK_qx_qw(i) = q0_min(i) .GE. 0.
+       alpha(i) = 1.
+      ENDDO
+c
 CC -----------------------------------------------------------------
 C    Sub-time-stepping
@@ -634,13 +679,24 @@
       DO isubstep = 1,nsub
 c------------------------------------------------------------
-      DO i=1,klon
+c
+c wk_adv is the logical flag enabling wake evolution in the time advance loop
+      DO i = 1,klon
+       wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1.
+      ENDDO
+c
+      DO i=1,klon
+        IF (wk_adv(i)) THEN
         gfl(i) = 2.*sqrt(3.14*wdens*sigmaw(i))
-      ENDDO
-      DO i=1,klon
-        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
-        sigmaw(i) =amin1(sigmaw(i),0.99)     !!!!!!!!
+        ENDIF
+      ENDDO
+      DO i=1,klon
+        IF (wk_adv(i)) THEN
+         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
+c        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
+c        sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
 c        wdens = wdens0/(10.*sigmaw)
 c        sigmaw =max(sigmaw,sigd_con)
 c        sigmaw =max(sigmaw,sigmad)
+        ENDIF
       ENDDO
 C
@@ -650,5 +706,9 @@
 cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
 cIM 060208 au niveau k=1..?
-      dp_deltomg(1:klon,1:klev)=0.
+      DO k= 1,klev
+      DO i = 1,klon
+        dp_deltomg(i,k)=0.
+      ENDDO
+      ENDDO
       DO k= 1,klev+1
       DO i = 1,klon
@@ -658,12 +718,14 @@
 c
       DO i=1,klon
+        IF (wk_adv(i)) THEN
         z(i)= 0.
         omg(i,1) = 0.
         dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i)))
+        ENDIF
       ENDDO
 c
       DO k= 2,klev
       DO i = 1,klon
-       IF( k .LE. ktop(i)) THEN
+       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
           dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
           z(i) = z(i)+dz(i)
@@ -675,7 +737,9 @@
 c
       DO i = 1,klon
+        IF (wk_adv(i)) THEN
         dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
         ztop(i) = z(i)+dztop(i)
         omgtop(i)=dp_deltomg(i,1)*ztop(i)
+        ENDIF
       ENDDO
 c
@@ -685,11 +749,13 @@
 c
        DO i=1,klon
+        IF (wk_adv(i)) THEN
         omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)
         dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
+        ENDIF
        ENDDO
 c
       DO k= 1,klev
       DO i = 1,klon
-       IF( k .LE. ktop(i)) THEN
+       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
           omg(i,k) = - rho(i,k)*rg*omg(i,k)
           dp_deltomg(i,k) = dp_deltomg(i,1)
@@ -701,9 +767,9 @@
 
       DO i=1,klon
-      IF (kupper(i) .GT. ktop(i)) THEN
+      IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN
         omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
      $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
         dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/
-     $                     (ptop(i)-pupper)
+     $                     (ptop(i)-pupper(i))
       ENDIF
       ENDDO
@@ -711,5 +777,5 @@
       DO k= 1,klev
       DO i = 1,klon
-       IF( k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
+       IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
           dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
           omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
@@ -718,4 +784,5 @@
       ENDDO
 c
+c
 c--    Compute wake average vertical velocity omgbw
 c
@@ -723,5 +790,7 @@
       DO k = 1,klev+1
       DO i=1,klon
+        IF ( wk_adv(i)) THEN
         omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
+        ENDIF
       ENDDO
       ENDDO
@@ -730,5 +799,7 @@
       DO k = 1,klev
       DO i=1,klon
+        IF ( wk_adv(i)) THEN
         dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
+        ENDIF
       ENDDO
       ENDDO
@@ -739,6 +810,8 @@
       DO k = 1,klev
       DO i=1,klon
-       alpha_up(i,k) = 0.
-       IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
+        IF ( wk_adv(i)) THEN
+         alpha_up(i,k) = 0.
+         IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
+        ENDIF
       ENDDO
       ENDDO
@@ -747,6 +820,8 @@
 
       DO i=1,klon
-        RRe1(i) = 1.-sigmaw(i)
-        RRe2(i) = sigmaw(i)
+        IF ( wk_adv(i)) THEN
+         RRe1(i) = 1.-sigmaw(i)
+         RRe2(i) = sigmaw(i)
+        ENDIF
       ENDDO
       RRd1 = -1.
@@ -757,5 +832,5 @@
       DO k= 1,klev
       DO i = 1,klon
-       IF(k .LE. kupper(i)+1) THEN
+       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
         dth(i,k) = deltatw(i,k)/ppi(i,k)
         Th1(i,k) = the(i,k) - sigmaw(i)     *dth(i,k)   ! undisturbed area
@@ -778,5 +853,5 @@
       DO k= 2,klev
       DO i = 1,klon
-       IF(k .LE. kupper(i)+1) THEN
+       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
         D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)
         D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)
@@ -790,11 +865,13 @@
 
       DO i=1,klon
-        omgbdth(i,1) = 0.
-        omgbdq(i,1) = 0.
+        IF( wk_adv(i)) THEN
+         omgbdth(i,1) = 0.
+         omgbdq(i,1) = 0.
+        ENDIF
       ENDDO
 
       DO k= 2,klev
       DO i = 1,klon
-       IF(k .LE. kupper(i)+1) THEN  !   loop on interfaces
+       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN  !   loop on interfaces
         omgbdth(i,k) = omgb(i,k)*(    dth(i,k-1) -     dth(i,k))
         omgbdq(i,k)  = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
@@ -806,5 +883,5 @@
       DO k= 1,klev
       DO i = 1,klon
-       IF(k .LE. kupper(i)-1) THEN
+       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
 c-----------------------------------------------------------------
 c
@@ -829,6 +906,9 @@
 c   and increment large scale tendencies
 c
-         dtls(i,k) = dtls(i,k) +
-     $               dtimesub*(
+
+c
+C
+CC -----------------------------------------------------------------
+         d_te(i,k) =  dtimesub*(
      $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
      $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) )
@@ -836,8 +916,6 @@
      $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
      $                      )*ppi(i,k)
-c         print*,'dtls=',dtls(i,k)
-c
-         dqls(i,k) = dqls(i,k) +
-     $               dtimesub*(
+c
+         d_qe(i,k) =  dtimesub*(
      $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
      $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) )
@@ -845,6 +923,6 @@
      $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
      $                      )
-c         print*,'dqls=',dqls(k)
-       ENDIF
+       ENDIF
+
 c-------------------------------------------------------------------
       ENDDO
@@ -856,5 +934,5 @@
       DO k= 1,klev
       DO i = 1,klon
-       IF(k .LE. kupper(i)-1) THEN
+       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
 c
 c Coefficient de répartition
@@ -912,21 +990,20 @@
 
         IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
-          deltatw(i,k) = deltatw(i,k)+dtimesub*
-     $          (ff(i)+dtKE(i,k)+dtPBL(i,k) 
+          d_deltatw(i,k) = dtimesub*
+     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
      $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
         ELSE
-           deltatw(i,k) = deltatw(i,k)+1/Tgw(i,k)*(1-exp(-dtimesub*
+           d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub*
      $          Tgw(i,k)))*
      $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
      $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
         ENDIF
-   
+
         dth(i,k) = deltatw(i,k)/ppi(i,k)
 
         gg(i)=d_deltaqw(i,k)/dtimesub
 
-       deltaqw(i,k) = deltaqw(i,k) +
-     $         dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k) - spread(i,k)*
-     $         deltaqw(i,k))
+       d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k)
+     $                            - spread(i,k)*deltaqw(i,k))
 
        d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
@@ -936,15 +1013,54 @@
       ENDDO
 
-C   And update large scale variables
+C
+C   Scale tendencies so that water vapour remains positive in w and x.
+C
+      call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw,
+     $                d_deltaqw,sigmaw,d_sigmaw,alpha)
+c
+      DO k = 1,klev
+      DO i = 1,klon
+       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
+        d_te(i,k)=alpha(i)*d_te(i,k)
+        d_qe(i,k)=alpha(i)*d_qe(i,k)
+        d_deltatw(i,k)=alpha(i)*d_deltatw(i,k)
+        d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k)
+        d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k)
+       ENDIF
+      ENDDO
+      ENDDO
+      DO i = 1,klon
+       IF( wk_adv(i)) THEN
+        d_sigmaw(i)=alpha(i)*d_sigmaw(i)
+       ENDIF
+      ENDDO
+
+C   Update large scale variables and wake variables
 cIM 060208 manque DO i + remplace DO k=1,kupper(i)
 cIM 060208     DO k = 1,kupper(i)
       DO k= 1,klev
       DO i = 1,klon
-       IF(k .LE. kupper(i)) THEN
+       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
+        dtls(i,k)=dtls(i,k)+d_te(i,k)
+        dqls(i,k)=dqls(i,k)+d_qe(i,k)
+       ENDIF
+      ENDDO
+      ENDDO
+      DO k= 1,klev
+      DO i = 1,klon
+       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
         te(i,k) = te0(i,k) + dtls(i,k)
         qe(i,k) = qe0(i,k) + dqls(i,k)
         the(i,k) = te(i,k)/ppi(i,k)
-       ENDIF
-      ENDDO
+        deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k)
+        deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k)
+        dth(i,k) = deltatw(i,k)/ppi(i,k)
+       ENDIF
+      ENDDO
+      ENDDO
+      DO i = 1,klon
+       IF( wk_adv(i)) THEN
+        sigmaw(i) = sigmaw(i)+d_sigmaw(i)
+       ENDIF
       ENDDO
 c
@@ -956,10 +1072,13 @@
 c
       DO i=1,klon
-      Ptop_provis(i)=ph(i,1)
+       IF ( wk_adv(i)) THEN
+        Ptop_provis(i)=ph(i,1)
+       ENDIF
       ENDDO
 c
       DO k= 2,klev
       DO i=1,klon
-        IF (Ptop_provis(i) .EQ. ph(i,1) .AND.
+        IF ( wk_adv(i) .AND.
+     $       Ptop_provis(i) .EQ. ph(i,1) .AND.
      $      dth(i,k) .GT. -delta_t_min .and.
      $      dth(i,k-1).LT. -delta_t_min) THEN
@@ -981,4 +1100,5 @@
       DO k = 1,klev
       DO i=1,klon
+       IF ( wk_adv(i)) THEN
         dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
         IF (dz(i) .gt. 0) THEN
@@ -987,4 +1107,5 @@
          dthmin(i) = amin1(dthmin(i),dth(i,k))
         ENDIF
+       ENDIF
       ENDDO
       ENDDO
@@ -993,6 +1114,8 @@
 
       DO i=1,klon
-       hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
-       hw(i) = amax1(hwmin,hw(i))
+       IF ( wk_adv(i)) THEN
+         hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
+         hw(i) = amax1(hwmin,hw(i))
+       ENDIF
       ENDDO
 c
@@ -1006,4 +1129,5 @@
       DO k = 1,klev
       DO i=1,klon
+       IF ( wk_adv(i)) THEN
         dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
         IF (dz(i) .gt. 0) THEN
@@ -1012,4 +1136,5 @@
          ktop(i) = k
         ENDIF
+       ENDIF
       ENDDO
       ENDDO
@@ -1018,5 +1143,7 @@
 c
       DO i=1,klon
+       IF ( wk_adv(i)) THEN
         Ptop_new(i)=ptop(i)
+       ENDIF
       ENDDO
 c
@@ -1024,5 +1151,6 @@
       DO i=1,klon
 cIM v3JYG; IF (k .GE. ktop(i)
-       IF (k .LE. ktop(i) .AND.
+       IF ( wk_adv(i) .AND.
+     $      k .LE. ktop(i) .AND.
      $      ptop_new(i) .EQ. ptop(i) .AND.
      $      dth(i,k) .GT. -delta_t_min .and.
@@ -1037,5 +1165,7 @@
 c
       DO i=1,klon
-      ptop(i) = ptop_new(i)
+       IF ( wk_adv(i)) THEN
+        ptop(i) = ptop_new(i)
+       ENDIF
       ENDDO
 
@@ -1050,5 +1180,5 @@
       DO k = 1,klev
       DO i=1,klon
-        IF (k .GE. kupper(i)) THEN
+        IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
          deltatw(i,k) = 0.
          deltaqw(i,k) = 0.
@@ -1058,4 +1188,111 @@
 c
 C
+c-------------Cstar computation---------------------------------
+      DO i=1, klon
+      sum_thu(i) = 0.
+      sum_tu(i) = 0.
+      sum_qu(i) = 0.
+      sum_thvu(i) = 0.
+      sum_dth(i) = 0.
+      sum_dq(i) = 0.
+      sum_rho(i) = 0.
+      sum_dtdwn(i) = 0.
+      sum_dqdwn(i) = 0.
+
+      av_thu(i) = 0.
+      av_tu(i) =0.
+      av_qu(i) =0.
+      av_thvu(i) = 0.
+      av_dth(i) = 0.
+      av_dq(i) = 0.
+      av_rho(i) =0.
+      av_dtdwn(i) =0.
+      av_dqdwn(i) = 0.
+      ENDDO
+C
+C Integrals (and wake top level number)
+C --------------------------------------
+C
+C Initialize sum_thvu to 1st level virt. pot. temp.
+
+      DO i=1,klon
+      z(i) = 1.
+      dz(i) = 1.
+      sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
+      sum_dth(i) = 0.
+      ENDDO
+
+      DO k = 1,klev
+      DO i=1,klon
+        dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
+        IF (dz(i) .GT. 0) THEN
+         z(i) = z(i)+dz(i)
+         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
+         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
+         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
+         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
+         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
+         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
+         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
+         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
+         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
+        ENDIF
+      ENDDO
+      ENDDO
+c
+      DO i=1,klon
+        hw0(i) = z(i)
+      ENDDO
+c
+C
+C - WAPE and mean forcing computation
+C ---------------------------------------
+C
+C ---------------------------------------
+C
+C Means
+
+      DO i=1,klon
+       av_thu(i) = sum_thu(i)/hw0(i)
+       av_tu(i) = sum_tu(i)/hw0(i)
+       av_qu(i) = sum_qu(i)/hw0(i)
+       av_thvu(i) = sum_thvu(i)/hw0(i)
+       av_dth(i) = sum_dth(i)/hw0(i)
+       av_dq(i) = sum_dq(i)/hw0(i)
+       av_rho(i) = sum_rho(i)/hw0(i)
+       av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
+       av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
+c
+       wape(i) = - rg*hw0(i)*(av_dth(i)
+     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
+     $     av_dq(i) ))/av_thvu(i)
+      ENDDO
+C
+C Filter out bad wakes
+
+      DO k = 1,klev
+       DO i=1,klon
+        IF ( wape(i) .LT. 0.) THEN
+          deltatw(i,k) = 0.
+          deltaqw(i,k) = 0.
+          dth(i,k) = 0.
+        ENDIF
+       ENDDO
+      ENDDO
+c
+      DO i=1,klon
+      IF ( wape(i) .LT. 0.) THEN
+        wape(i) = 0.
+        Cstar(i) = 0.
+        hw(i) = hwmin
+        sigmaw(i) = max(sigmad,sigd_con(i))
+        fip(i) = 0.
+        gwake(i) = .FALSE.
+      ELSE
+        Cstar(i) = stark*sqrt(2.*wape(i))
+        gwake(i) = .TRUE.
+      ENDIF
+      ENDDO
+
        ENDDO      ! end sub-timestep loop
 C
@@ -1065,5 +1302,5 @@
       DO k = 1,klev
       DO i=1,klon
-        IF (k .LE. kupper(i)-1) THEN
+        IF ( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
          dtls(i,k) = dtls(i,k)/dtime
          dqls(i,k) = dqls(i,k)/dtime
@@ -1111,4 +1348,5 @@
       DO k =1,klev
       DO i=1,klon
+       IF ( wk_adv(i)) THEN
         rho(i,k) = p(i,k)/(rd*te(i,k))
         IF(k .eq. 1) THEN
@@ -1125,4 +1363,5 @@
         rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
         dth(i,k) = deltatw(i,k)/ppi(i,k)
+       ENDIF
       ENDDO
       ENDDO
@@ -1134,12 +1373,15 @@
 
       DO i=1,klon
+       IF ( wk_adv(i)) THEN
         z(i) = 1.
         dz(i) = 1.
         sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
         sum_dth(i) = 0.
+      ENDIF
       ENDDO
 
       DO k = 1,klev
       DO i=1,klon
+       IF ( wk_adv(i)) THEN
         dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
         IF (dz(i) .GT. 0) THEN
@@ -1155,12 +1397,15 @@
          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
         ENDIF
-      ENDDO
-      ENDDO
-c
-      DO i=1,klon
+       ENDIF
+      ENDDO
+      ENDDO
+c
+      DO i=1,klon
+       IF ( wk_adv(i)) THEN
         hw0(i) = z(i)
-      ENDDO
-c
-C 2.1 - WAPE and mean forcing computation
+       ENDIF
+      ENDDO
+c
+C - WAPE and mean forcing computation
 C-------------------------------------------------------------
 
@@ -1168,4 +1413,5 @@
 
       DO i=1, klon
+       IF ( wk_adv(i)) THEN
         av_thu(i) = sum_thu(i)/hw0(i)
         av_tu(i) = sum_tu(i)/hw0(i)
@@ -1181,7 +1427,8 @@
      $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+
      $     av_dth(i)*av_dq(i) ))/av_thvu(i)
-      ENDDO
-
-C 2.2 Prognostic variable update
+       ENDIF
+      ENDDO
+
+C Prognostic variable update
 C ------------------------------------------------------------
 
@@ -1190,5 +1437,5 @@
       DO k = 1,klev
       DO i=1,klon
-        IF ( wape2(i) .LT. 0.) THEN
+        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
           deltatw(i,k) = 0.
           deltaqw(i,k) = 0.
@@ -1200,5 +1447,6 @@
 
       DO i=1, klon
-      IF ( wape2(i) .LT. 0.) THEN
+       IF ( wk_adv(i)) THEN
+       IF ( wape2(i) .LT. 0.) THEN
         wape2(i) = 0.
         Cstar2(i) = 0.
@@ -1212,12 +1460,16 @@
         gwake(i) = .TRUE.
       ENDIF
+      ENDIF
       ENDDO
 c
       DO i=1, klon
+       IF ( wk_adv(i)) THEN
         ktopw(i) = ktop(i)
+       ENDIF
       ENDDO
 c
       DO i=1, klon
-      IF (ktopw(i) .gt. 0 .and. gwake(i)) then
+       IF ( wk_adv(i)) THEN
+       IF (ktopw(i) .gt. 0 .and. gwake(i)) then
 
 Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
@@ -1234,4 +1486,5 @@
          FIP(i) = 0.
        ENDIF
+       ENDIF
       ENDDO
 c
@@ -1241,9 +1494,14 @@
 C              alors il disparait en se mélangeant à la partie undisturbed
 c
+      sigmaw_max = 0.9
       DO k = 1,klev
        DO i=1, klon
-         IF ((sigmaw(i).GT.0.9).or.
+c correction NICOLAS     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
+!         print*,'wape wape2 ktopw OK_qx_qw =',
+!     $           wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
+         IF ((sigmaw(i).GT.sigmaw_max).or.
      $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
-     $      (ktopw(i).le.2)) THEN
+     $      (ktopw(i).le.2) .OR.
+     $     .not. OK_qx_qw(i)) THEN
 cIM cf NR/JYG 251108  $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
 ccc      IF (sigmaw(i).GT.0.9) THEN
@@ -1257,6 +1515,9 @@
 c
       DO i=1, klon
-         IF ((sigmaw(i).GT.0.9).or.
-     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
+         IF ( (sigmaw(i).GT.sigmaw_max).or.
+     $      ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
+     $      (ktopw(i).le.2) .OR.
+     $     .not. OK_qx_qw(i)) THEN
+! correction NICOLAS     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
 ccc      IF (sigmaw(i).GT.0.9) THEN
          wape(i) = 0.
@@ -1272,4 +1533,96 @@
       RETURN
       END
+
+      SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe,
+     $           deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)
+c------------------------------------------------------
+cDtermination du coefficient alpha tel que les tendances
+c corriges alpha*d_G, pour toutes les grandeurs G, correspondent
+c a une humidite positive dans la zone (x) et dans la zone (w).
+c------------------------------------------------------
+c
+ 
+c  Input
+      REAL qe(nlon,nl),d_qe(nlon,nl)
+      REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl)
+      REAL sigmaw(nlon),d_sigmaw(nlon)
+      LOGICAL wk_adv(nlon)
+      INTEGER nl,nlon
+c  Output
+      REAL alpha(nlon)
+c  Internal variables
+      REAL alpha1(nlon)
+      REAL x,a,b,c,discrim,zeta(nlon)
+      REAL epsilon
+!      DATA epsilon/1.e-9/
+c
+      DO k=1,nl
+      DO i = 1,nlon
+       IF (wk_adv(i)) THEN
+        IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then
+         zeta(i)=0.
+        ELSE
+         zeta(i)=1.
+        END IF
+       ENDIF
+      ENDDO
+      DO i = 1,nlon
+       IF (wk_adv(i)) THEN
+        x = qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)
+     $   +d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)
+     $   -d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k))
+      a=-d_sigmaw(i)*d_deltaqw(i,k)
+      b=d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)
+     $           -deltaqw(i,k)*d_sigmaw(i)
+!      c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)-epsilon
+       c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)
+
+      discrim=b*b-4.*a*c
+!       print*,'ZETA *********************'  
+!       print*,'zeta sigmaw ',zeta(:)
+!       print*,'SIGMA *********************'
+!       print*,'sigmaw ',sigmaw(:)
+
+!       print*,' x ************************'
+!       print*,'x ',x
+!       print*,'  a+b ************************'
+!       print*,'a+b ',a+b
+
+!       print*,'a b c delta zeta ',a,b,c,discrim
+        IF (a+b .GE. 0.) THEN
+         alpha1(i)=1.
+        ELSE
+         IF (x .GE. 0.) THEN
+            alpha1(i)=1.
+         ELSE
+!              IF (a .GE. 0.) THEN
+              IF (a .GT. 0.) THEN
+!       print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)
+!       print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)
+                 alpha1(i)=0.9*min(   (2.*c)/(-b+sqrt(discrim)),
+     $                        (-b+sqrt(discrim))/(2.*a)   )
+              ELSE IF (a.eq.0.) THEN
+                 alpha1(i)=0.9*(-c/b)
+              ELSE
+!       print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)
+!       print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)
+                 alpha1(i)=0.9*max(   (2.*c)/(-b+sqrt(discrim)),
+     $                        (-b+sqrt(discrim))/(2.*a)   )
+              ENDIF
+         ENDIF
+        ENDIF
+       ENDIF
+      ENDDO
+      ENDDO
+c
+      DO i = 1,nlon
+       IF (wk_adv(i)) THEN
+        alpha(i) = min(alpha(i),alpha1(i))
+       ENDIF
+      ENDDO
+c
+      return
+      end
+
       Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
      :                ,te0,qe0,omgb
@@ -2315,4 +2668,5 @@
 C              alors il disparait en se mélangeant à la partie undisturbed
 
+! correction NICOLAS     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
       IF ((sigmaw.GT.0.9).or.
      .     ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN
