Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/etat0_netcdf.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/etat0_netcdf.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/etat0_netcdf.F	(revision 1322)
@@ -144,4 +144,5 @@
       REAL      :: solarlong0
       real :: seuil_inversion
+      real :: alp_offset
 
       integer  read_climoz ! read ozone climatology
@@ -184,5 +185,6 @@
      &                 iflag_thermals,nsplit_thermals,tau_thermals,     &
      &                 iflag_thermals_ed,iflag_thermals_optflux,        &
-     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz )
+     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz,  &
+     &                 alp_offset)
 
 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F	(revision 1322)
@@ -144,4 +144,5 @@
       REAL      :: solarlong0
       real :: seuil_inversion
+      REAL :: alp_offset
 
       integer  read_climoz ! read ozone climatology
@@ -184,5 +185,6 @@
      &                 iflag_thermals,nsplit_thermals,tau_thermals,     &
      &                 iflag_thermals_ed,iflag_thermals_optflux,        &
-     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz )
+     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz,  &
+     &                 alp_offset)
 
 ! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F	(revision 1322)
@@ -677,6 +677,8 @@
          call suspend_timer(timer_caldyn)
 
+        if (prt_level >= 10) then
          write(lunout,*)
      &   'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
+        endif
 c$OMP END MASTER
 
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calwake.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calwake.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calwake.F	(revision 1322)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
       SUBROUTINE CALWAKE(paprs,pplay,dtime
      :             ,t,q,omgb
@@ -189,8 +192,8 @@
      $                ,Cstar,d_deltat_gw
      $                ,d_deltatw,d_deltaqw)
-
-      DO i=1,klon 
-       IF (ktopw(i) .GT. 0) THEN
-         DO l=1,klev
+c
+      DO l=1,klev
+       DO i=1,klon 
+        IF (ktopw(i) .GT. 0) THEN
            wake_deltat(i,l)= dtw(i,l)
            wake_deltaq(i,l)= dqw(i,l)
@@ -212,7 +215,5 @@
            wake_ddeltat(i,l) = d_deltatw(i,l)
            wake_ddeltaq(i,l) = d_deltaqw(i,l)
-         ENDDO
-       ELSE
-         DO l = 1,klev
+        ELSE
            wake_deltat(i,l)= 0.
            wake_deltaq(i,l)= 0.
@@ -222,4 +223,6 @@
            wake_dtKE(i,l) = 0.
            wake_dqKE(i,l) = 0.
+ 	   wake_dtPBL(i,l) = 0.
+	   wake_dqPBL(i,l) = 0.
            wake_omg(i,l) = 0.
            wake_dp_deltomg(i,l) = 0.
@@ -230,7 +233,11 @@
            undi_t(i,l)=te(i,l)
            undi_q(i,l)=qe(i,l)
-         ENDDO
-       ENDIF
-
+           wake_ddeltat(i,l) = 0.
+           wake_ddeltaq(i,l) = 0.
+        ENDIF
+       ENDDO
+      ENDDO
+c
+      DO i=1,klon 
        wake_h(i)= hw(i)
        wake_s(i)= sigmaw(i)
@@ -241,10 +248,9 @@
        wake_Cstar(i) = Cstar(i)
        wake_dens(i) = wdens(i)
-c
-cIM 290108 999  CONTINUE
-c
-      ENDDO
+      ENDDO
+c
       RETURN
       END
+
       SUBROUTINE CALWAKE_scal(paprs,pplay,dtime
      :             ,t,q,omgb
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/conf_phys.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/conf_phys.F90	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/conf_phys.F90	(revision 1322)
@@ -22,5 +22,6 @@
                        iflag_thermals,nsplit_thermals,tau_thermals, &
                        iflag_thermals_ed,iflag_thermals_optflux, &
-                       iflag_coupl,iflag_clos,iflag_wake, read_climoz)
+                       iflag_coupl,iflag_clos,iflag_wake, read_climoz, &
+                       alp_offset)
 
    use IOIPSL
@@ -110,4 +111,6 @@
   integer :: iflag_clos
   integer :: iflag_wake
+  real :: alp_offset
+  REAL, SAVE :: alp_offset_omp
   integer,SAVE :: iflag_coupl_omp,iflag_clos_omp,iflag_wake_omp
   integer,SAVE :: iflag_cvl_sigd_omp
@@ -1032,4 +1035,12 @@
   iflag_wake_omp = 0
   call getin('iflag_wake',iflag_wake_omp)
+
+!Config Key  = alp_offset
+!Config Desc =  
+!Config Def  = 0
+!Config Help = 
+!
+  alp_offset_omp = 0.01
+  call getin('alp_offset',alp_offset_omp)
 
 !
@@ -1449,4 +1460,5 @@
     iflag_clos = iflag_clos_omp
     iflag_wake = iflag_wake_omp
+    alp_offset = alp_offset_omp
     iflag_cvl_sigd = iflag_cvl_sigd_omp
     type_run = type_run_omp
@@ -1616,4 +1628,6 @@
   write(numout,*)' Fmax = ', Fmax
   write(numout,*)' alphas = ', alphas
+  write(numout,*)' iflag_wake = ', iflag_wake
+  write(numout,*)' alp_offset = ', alp_offset
 
   write(numout,*)' lonmin lonmax latmin latmax bilKP_ins =',&
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3_cine.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3_cine.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3_cine.F	(revision 1322)
@@ -1,5 +1,8 @@
+!
+! $Id$
+!
         SUBROUTINE cv3_cine(nloc,ncum,nd,icb,inb
      :                      ,pbase,plcl,p,ph,tv,tvp
-     :                      ,cina,cinb)
+     :                      ,cina,cinb,plfc)
 
 ***************************************************************
@@ -26,5 +29,5 @@
 c
 c output
-      real cina(nloc),cinb(nloc)
+      real cina(nloc),cinb(nloc),plfc(nloc)
 c
 c local variables
@@ -34,5 +37,4 @@
       logical lswitch(nloc),lswitch1(nloc),lswitch2(nloc)
       logical exist_lfc(nloc)
-      real plfc(nloc)
       real dpmax
       real deltap,dcin
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3_routines.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3_routines.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3_routines.F	(revision 1322)
@@ -74,7 +74,6 @@
 c      dtcrit = -5.0
 c      tau    = 3000.
-cc      tau = 1800.
-c     tau= 2800.
-      tau=8000.
+      tau = 1800.
+cc      tau=8000.
       beta   = 1.0 - delt/tau
       alpha1 = 1.5e-3
@@ -2094,5 +2093,5 @@
 cc---end jyg---
 c
-c--------retour à la formulation originale d''Emanuel.
+c--------retour à la formulation originale d'Emanuel.
       b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
       c6=water(il,i+1)+bfac*wdtrain(il)
@@ -2100,6 +2099,12 @@
       if(c6.gt.0.0)then
        revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
-       water(il,i)=revap*revap      !equation de conservation
+cjyg    Dans sa formulation originale, Emanuel calcule l'evaporation par:
+cc             evap(il,i)=sigt*afac*revap
+c     ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
+c     Ici,l'evaporation evap est simplement calculee par l'equation de
+c     conservation.
+       water(il,i)=revap*revap
       else
+cjyg----   Correction : si c6 <= 0, water(il,i)=0.
        water(il,i) = 0.
       endif
@@ -2345,5 +2350,5 @@
       real esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
       real th_wake(nloc,nd)
-      real alpha_qpos(nloc)
+      real alpha_qpos(nloc),alpha_qpos1(nloc)
       real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld
       real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld
@@ -3050,4 +3055,13 @@
        do il=1,ncum
         IF (iflag(il) .le. 1) THEN
+	IF (cvflag_grav) then
+        ex=0.01*grav*ment(il,inb(il),inb(il))
+     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
+     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
+        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
+        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
+     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
+     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
+	else
         ex=0.1*ment(il,inb(il),inb(il))
      :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
@@ -3057,4 +3071,5 @@
      :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
      :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
+        ENDIF   !cvflag grav
         ENDIF    !iflag
        enddo
@@ -3122,7 +3137,11 @@
 c        in order to ensure moisture positivity
       DO il = 1,ncum
+      alpha_qpos(il)=1.
        IF (iflag(il) .le. 1) THEN
-        alpha_qpos(il) = max(1. , -delt*fr(il,1)/
+        if (fr(il,1) .le. 0.) then
+            alpha_qpos(il) = max(alpha_qpos(il) ,
+     :	   (-delt*fr(il,1))/
      :     (s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
+        end if
        ENDIF
       ENDDO
@@ -3130,6 +3149,10 @@
        DO il = 1,ncum
         IF (iflag(il) .le. 1) THEN
-        alpha_qpos(il) = max(alpha_qpos(il) , -delt*fr(il,i)/
+          IF (fr(il,i) .le. 0.) THEN
+           alpha_qpos1(il)=max(1. , (-delt*fr(il,i))/
      :     (s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
+             IF (alpha_qpos1(il) .ge. alpha_qpos(il))
+     :           alpha_qpos(il)=alpha_qpos1(il)
+          ENDIF
         ENDIF
        ENDDO
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3p1_closure.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3p1_closure.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3p1_closure.F	(revision 1322)
@@ -63,5 +63,5 @@
       integer nsupmax(nloc)
       real supcrit,temp(nloc,nd)
-      real P1(nloc),Pmin(nloc)
+      real P1(nloc),Pmin(nloc),plfc(nloc)
       real asupmax0(nloc)
       logical ok(nloc)
@@ -385,5 +385,5 @@
       CALL cv3_cine (nloc,ncum,nd,icb,inb
      :                      ,pbase,plcl,p,ph,tv,tvp
-     :                      ,cina,cinb)
+     :                      ,cina,cinb,plfc)
 c
       DO il = 1,ncum
@@ -495,6 +495,8 @@
       do k= 1,nl
        do il = 1,ncum
-!IM       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
-        IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
+!old       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
+!IM        IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
+       IF (k .ge. icb(il) .and. k .le. inb(il)         !cor jyg
+     $     .and. icb(il)+1 .le. inb(il)) THEN          !cor jyg
          cbmflim(il) = cbmflim(il)+MLIM(il,k)
         ENDIF
@@ -548,8 +550,10 @@
         do il = 1,ncum
          IF ( k .ge. icb(il)+1 .AND. k .le. inb(il)) THEN
-         sig(il,k) = beta*sig(il,k)+(1.-beta)*coef(il)*siglim(il,k)
-cc         sig(il,k) = beta*sig(il,k)+siglim(il,k)
-         w0(il,k) = beta*w0(il,k)  +(1.-beta)*wlim(il,k)
-         AMU=SIG(il,k)*W0(il,k)
+         amu=beta*sig(il,k)*w0(il,k)+
+     :   (1.-beta)*coef(il)*siglim(il,k)*wlim(il,k)
+         w0(il,k) = wlim(il,k)
+         w0(il,k) =max(w0(il,k),1.e-10)
+         sig(il,k)=amu/w0(il,k)
+         sig(il,k)=min(sig(il,k),1.)
 cc         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
          M(il,k)=AMU*0.007*P(il,k)*(PH(il,k)-PH(il,k+1))/TV(il,k)
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/ini_wake.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/ini_wake.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/ini_wake.F	(revision 1322)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
       SUBROUTINE INI_WAKE(wape,fip,it_wape_prescr,
      :     wape_prescr, fip_prescr, alp_bl_prescr, ale_bl_prescr)
@@ -23,4 +26,6 @@
 c   wape_prescr    : valeur prescrite de la WAPE.
 c   fip_prescr     : valeur prescrite de la FIP.
+c   ale_bl_prescr  : valeur prescrite de la Ale de PBL.
+c   alp_bl_prescr  : valeur prescrite de la Alp de PBL.
 c
 c Variables internes
@@ -29,17 +34,12 @@
 c   w  = WAPE lue
 c   f  = FIP lue
+c   alebl  = Ale de PBL lue
+c   alpbl  = Alp de PBL lue
 c
+      include 'iniprint.h'
 cdeclarations
       real ale_bl_prescr
       real alp_bl_prescr
       real it
-cCR: on rajoute ale et alp de la PBL precrits
-c     open (99,file='wake.data',form='formatted')
-c     read (99,*) it
-c     read (99,*) w
-c     read (99,*) f
-c     read (99,*) u
-c     read (99,*) p
-c     close (99)
 
 ! FH A mettre si besoin dans physiq.def
@@ -48,8 +48,18 @@
       w=4.
       f=0.1
-      u=0.1
-      p=4.
+      alebl=4.
+      alpbl=0.1
 c
-      print *,' it,w ',it,w
+cCR: on rajoute ale et alp de la PBL precrits
+      open (99,file='wake.data',form='formatted',status='old',err=902)
+      read (99,*) it
+      read (99,*) w
+      read (99,*) f
+      read (99,*,end=901) alebl
+      read (99,*,end=901) alpbl
+901   close (99)
+902   continue
+c
+      write(lunout,*)' it,wape ',it,wape
       it_wape_prescr = it
       if (w .lt. 0) then
@@ -61,7 +71,7 @@
       endif
 c
-      print *,' u,p ',u,p
-      alp_bl_prescr=u
-      ale_bl_prescr=p
+      write(lunout,*)' alebl, alpbl ',alebl,alpbl
+      ale_bl_prescr=alebl
+      alp_bl_prescr=alpbl
       print *,'Initialisation de la poche : WAPE, FIP imposees ='
      $               ,wape_prescr, fip_prescr
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyetat0.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyetat0.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyetat0.F	(revision 1322)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id$
 !
 c
@@ -228,4 +228,5 @@
      $          'coherente ', i, zmasq(i), pctsrf(i, is_ter)
      $          ,pctsrf(i, is_lic)
+            WRITE(*,*) 'Je force la coherence zmasq=fractint'
             zmasq(i) = fractint(i)
         ENDIF 
@@ -238,4 +239,5 @@
      $          'coherente ', i, zmasq(i) , pctsrf(i, is_oce)
      $          ,pctsrf(i, is_sic)
+            WRITE(*,*) 'Je force la coherence zmasq=fractint'
             zmasq(i) = fractint(i)
         ENDIF 
@@ -987,4 +989,18 @@
       PRINT*,'(ecart-type) wake_cstar:', xmin, xmax
 c
+c wake_pe
+c
+      CALL get_field("WAKE_PE",wake_pe,found)
+      IF (.NOT. found) THEN
+         PRINT*, "phyetat0: Le champ <WAKE_PE> est absent"
+         PRINT*, "Depart legerement fausse. Mais je continue"
+         wake_pe=0.
+      ENDIF
+      xmin = 1.0E+20
+      xmax = -1.0E+20
+      xmin = MINval(wake_pe)
+      xmax = MAXval(wake_pe)
+      PRINT*,'(ecart-type) wake_pe:', xmin, xmax
+c
 c wake_fip
 c
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyredem.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyredem.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyredem.F	(revision 1322)
@@ -299,5 +299,5 @@
 !!!!!!!!!!!!!!!!!!!! FIN TKE PBL !!!!!!!!!!!!!!!!!!!!!!!!!
 cIM ajout zmax0, f0, ema_work1, ema_work2
-cIM wake_deltat, wake_deltaq, wake_s, wake_cstar, wake_fip
+cIM wake_deltat, wake_deltaq, wake_s, wake_cstar, wake_pe, wake_fip
       
       CALL put_field("ZMAX0","ZMAX0",zmax0)
@@ -318,4 +318,6 @@
       CALL put_field("WAKE_CSTAR","WAKE_CSTAR",wake_cstar)
       
+      CALL put_field("WAKE_PE","WAKE_PE",wake_pe)
+
       CALL put_field("WAKE_FIP","WAKE_FIP",wake_fip)
 
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_state_var_mod.F90	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_state_var_mod.F90	(revision 1322)
@@ -1,2 +1,5 @@
+!
+! $Id$
+!
       MODULE phys_state_var_mod
 ! Variables sauvegardees pour le startphy.nc
@@ -169,4 +172,5 @@
 ! wake_Cstar  : vitesse d'etalement de la poche
 ! wake_s      : fraction surfacique occupee par la poche froide
+! wake_pe     : wake potential energy - WAPE
 ! wake_fip    : Gust Front Impinging power - ALP
 ! dt_wake, dq_wake: LS tendencies due to wake
@@ -179,4 +183,6 @@
       REAL,ALLOCATABLE,SAVE :: wake_s(:)
 !$OMP THREADPRIVATE(wake_s)
+      REAL,ALLOCATABLE,SAVE :: wake_pe(:)
+!$OMP THREADPRIVATE(wake_pe)
       REAL,ALLOCATABLE,SAVE :: wake_fip(:)
 !$OMP THREADPRIVATE(wake_fip)
@@ -370,5 +376,6 @@
       ALLOCATE(wght_th(klon,klev))
       ALLOCATE(wake_deltat(klon,klev), wake_deltaq(klon,klev))
-      ALLOCATE(wake_Cstar(klon), wake_s(klon), wake_fip(klon))
+      ALLOCATE(wake_Cstar(klon), wake_s(klon))
+      ALLOCATE(wake_pe(klon), wake_fip(klon))
       ALLOCATE(dt_wake(klon,klev), dq_wake(klon,klev))
       ALLOCATE(pfrac_impa(klon,klev), pfrac_nucl(klon,klev))
@@ -464,5 +471,5 @@
       deallocate(lalim_conv, wght_th)
       deallocate(wake_deltat, wake_deltaq)
-      deallocate(wake_Cstar, wake_s, wake_fip)
+      deallocate(wake_Cstar, wake_s, wake_pe, wake_fip)
       deallocate(dt_wake, dq_wake)
       deallocate(pfrac_impa, pfrac_nucl)
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F	(revision 1322)
@@ -625,5 +625,6 @@
       REAL q_undi(klon,klev)               ! humidite moyenne dans la zone non perturbee
 c
-      REAL wake_pe(klon)              ! Wake potential energy - WAPE 
+cjyg
+ccc      REAL wake_pe(klon)              ! Wake potential energy - WAPE 
 
       REAL wake_gfl(klon)             ! Gust Front Length
@@ -641,4 +642,7 @@
       REAL dt_a(klon,klev)
       REAL dq_a(klon,klev)
+      REAL, SAVE :: alp_offset
+c$OMP THREADPRIVATE(alp_offset)
+
 c
 cRR:fin declarations poches froides
@@ -1188,5 +1192,6 @@
      .     iflag_thermals_ed,iflag_thermals_optflux,
 c     nv flags pour la convection et les poches froides
-     .     iflag_coupl,iflag_clos,iflag_wake, read_climoz)
+     .     iflag_coupl,iflag_clos,iflag_wake, read_climoz,
+     &     alp_offset)
       call phys_state_var_init(read_climoz)
       print*, '================================================='
@@ -2021,17 +2026,29 @@
 c
 ccalcul de ale_wake et alp_wake
-       do i = 1,klon
-          if (iflag_wake.eq.1) then
-          ale_wake(i) = 0.5*wake_cstar(i)**2
-          alp_wake(i) = wake_fip(i)
-          else
-          ale_wake(i) = 0.
-          alp_wake(i) = 0.
-          endif
-       enddo
+       if (iflag_wake.eq.1) then
+         if (itap .le. it_wape_prescr) then
+          do i = 1,klon
+           ale_wake(i) = wape_prescr
+           alp_wake(i) = fip_prescr
+          enddo
+         else
+          do i = 1,klon
+cjyg  ALE=WAPE au lieu de ALE = 1/2 Cstar**2
+ccc           ale_wake(i) = 0.5*wake_cstar(i)**2
+           ale_wake(i) = wake_pe(i)
+           alp_wake(i) = wake_fip(i)
+          enddo
+         endif
+       else
+         do i = 1,klon
+           ale_wake(i) = 0.
+           alp_wake(i) = 0.
+         enddo
+       endif
 ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
 cdans le thermique sinon
        if (iflag_coupl.eq.0) then
-          if (debut) print*,'ALE et ALP imposes'
+          if (debut.and.prt_level.gt.9)
+     $                     WRITE(lunout,*)'ALE et ALP imposes'
           do i = 1,klon
 con ne couple que ale
@@ -2046,5 +2063,6 @@
           do i = 1,klon
               ALE(i) = max(ale_wake(i),Ale_bl(i))
-              ALP(i) = alp_wake(i) + Alp_bl(i)
+c avant        ALP(i) = alp_wake(i) + Alp_bl(i)
+              ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
 c         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
 c         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
@@ -2289,4 +2307,19 @@
 
       endif
+c
+c===================================================================
+cJYG
+      IF (ip_ebil_phy.ge.2) THEN 
+        ztit='after wake'
+        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
+     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
+     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
+      END IF 
+
 c      print*,'apres callwake iflag_cldcon=', iflag_cldcon
 c
@@ -2402,4 +2435,9 @@
      e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
      s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
       END IF 
 
@@ -2450,5 +2488,7 @@
          enddo
          tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
-         print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
+         if(prt_level.ge.9)
+     &     write(lunout,*)'TAU TH OK ',
+     &     tau_overturning_th(1),detr_therm(1,3)
 
 c On impose que l'air autour de la fraction couverte par le thermique
@@ -2561,5 +2601,6 @@
 
 !   les ratqs sont une combinaison de ratqss et ratqsc
-!       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
+       if(prt_level.ge.9)
+     $       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
 
          if (tau_ratqs>1.e-10) then
@@ -2785,4 +2826,9 @@
      e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
      s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+        call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
       END IF 
 c
@@ -2992,5 +3038,5 @@
       itaprad = itaprad + 1
 
-      if (iflag_radia.eq.0) then
+      if (iflag_radia.eq.0 .and. prt_level.ge.9) then
       print *,'--------------------------------------------------'
       print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
@@ -3195,4 +3241,9 @@
      e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
      s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
+         call diagphy(airephy,ztit,ip_ebil_phy
+     e      , zero_v, zero_v, zero_v, zero_v, zero_v
+     e      , zero_v, zero_v, zero_v, ztsol
+     e      , d_h_vcol, d_qt, d_ec
+     s      , fs_bound, fq_bound )
       END IF 
 c
@@ -3262,6 +3313,7 @@
       IF (offline) THEN
 
-         print*,'Attention on met a 0 les thermiques pour phystoke'
-         call phystokenc (
+       IF (prt_level.ge.9)
+     $    print*,'Attention on met a 0 les thermiques pour phystoke'
+	 call phystokenc (
      I                   nlon,klev,pdtphys,rlon,rlat,
      I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/wake.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/wake.F	(revision 1321)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/wake.F	(revision 1322)
@@ -1,3 +1,6 @@
-      Subroutine WAKE (p,ph,ppi,dtime,sigd_con
+!
+! $Id$
+!
+      Subroutine WAKE (p,ph,pi,dtime,sigd_con
      :                ,te0,qe0,omgb
      :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
@@ -24,4 +27,1830 @@
 c
       use dimphy
+      IMPLICIT none
+c============================================================================
+C
+C
+C   But : Decrire le comportement des poches froides apparaissant dans les
+C        grands systemes convectifs, et fournir l'energie disponible pour
+C        le declenchement de nouvelles colonnes convectives.
+C
+C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
+C                      deltaqw    : ecart d'humidite wake-undisturbed area
+C                      sigmaw     : fraction d'aire occupee par la poche.
+C
+C   Variable de sortie : 
+c
+c			 wape : WAke Potential Energy
+c                        fip  : Front Incident Power (W/m2) - ALP
+c                        gfl  : Gust Front Length per unit area (m-1)
+C                        dtls : large scale temperature tendency due to wake
+C                        dqls : large scale humidity tendency due to wake
+C                        hw   : hauteur de la poche
+C                     dp_omgb : vertical gradient of large scale omega
+C                     wdens   : densite de poches
+C                      omgbdth: flux of Delta_Theta transported by LS omega
+C                      dtKE   : differential heating (wake - unpertubed)
+C                      dqKE   : differential moistening (wake - unpertubed)
+C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
+C                 dp_deltomg  : vertical gradient of omg (s-1)
+C                     spread  : spreading term in dt_wake and dq_wake
+C                 deltatw     : updated temperature difference (T_w-T_u).
+C                 deltaqw     : updated humidity difference (q_w-q_u).
+C                 sigmaw      : updated wake fractional area.
+C                 d_deltat_gw : delta T tendency due to GW
+c
+C   Variables d'entree : 
+c
+c		         aire : aire de la maille
+c			 te0  : temperature dans l'environnement  (K)
+C                        qe0  : humidite dans l'environnement     (kg/kg)
+C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
+C                        dtdwn: source de chaleur due aux descentes (K/s)
+C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
+C			 dta  : source de chaleur due courants satures et detrain  (K/s)
+C			 dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
+C                        amdwn: flux de masse total des descentes, par unite de
+C                                surface de la maille (kg/m2/s)
+C                        amup : flux de masse total des ascendances, par unite de
+C                                surface de la maille (kg/m2/s)
+C                        p    : pressions aux milieux des couches (Pa)
+C                        ph   : pressions aux interfaces (Pa)
+C                        pi  : (p/p_0)**kapa (adim)
+C                        dtime: increment temporel (s)
+c
+C   Variables internes :
+c
+c			 rhow : masse volumique de la poche froide
+C                        rho  : environment density at P levels
+C                        rhoh : environment density at Ph levels
+C                        te   : environment temperature | may change within
+C                        qe   : environment humidity    | sub-time-stepping
+C                        the  : environment potential temperature
+C                        thu  : potential temperature in undisturbed area
+C                        tu   :  temperature  in undisturbed area
+C                        qu   : humidity in undisturbed area
+C                      dp_omgb: vertical gradient og LS omega
+C                      omgbw  : wake average vertical omega
+C                     dp_omgbw: vertical gradient of omgbw
+C                      omgbdq : flux of Delta_q transported by LS omega
+C                        dth  : potential temperature diff. wake-undist.
+C                        th1  : first pot. temp. for vertical advection (=thu)
+C                        th2  : second pot. temp. for vertical advection (=thw)
+C                        q1   : first humidity for vertical advection
+C                        q2   : second humidity for vertical advection
+C                     d_deltatw   : terme de redistribution pour deltatw
+C                     d_deltaqw   : terme de redistribution pour deltaqw
+C                      deltatw0   : deltatw initial
+C                      deltaqw0   : deltaqw initial
+C                      hw0    : hw initial
+C                      sigmaw0: sigmaw initial
+C                      amflux : horizontal mass flux through wake boundary
+C                      wdens_ref: initial number of wakes per unit area (3D) or per
+C                               unit length (2D), at the beginning of each time step
+C                      Tgw    : 1 sur la période de onde de gravité
+c                      Cgw    : vitesse de propagation de onde de gravité
+c                      LL     : distance entre 2 poches
+
+c-------------------------------------------------------------------------
+c          Déclaration de variables
+c-------------------------------------------------------------------------
+
+#include "dimensions.h"
+#include "YOMCST.h"
+#include "cvthermo.h"
+#include "iniprint.h"
+
+c Arguments en entree
+c--------------------
+
+      REAL, dimension(klon,klev) :: p, pi
+      REAL, dimension(klon,klev+1) :: ph, omgb
+      REAL dtime
+      REAL, dimension(klon,klev) :: te0,qe0
+      REAL, dimension(klon,klev) :: dtdwn, dqdwn
+      REAL, dimension(klon,klev) :: wdtPBL,wdqPBL
+      REAL, dimension(klon,klev) :: udtPBL,udqPBL
+      REAL, dimension(klon,klev) :: amdwn, amup
+      REAL, dimension(klon,klev) :: dta, dqa
+      REAL, dimension(klon) :: sigd_con
+
+c Sorties
+c--------
+
+      REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
+      REAL, dimension(klon,klev) :: tu, qu
+      REAL, dimension(klon,klev) :: dtls, dqls
+      REAL, dimension(klon,klev) :: dtKE, dqKE
+      REAL, dimension(klon,klev) :: dtPBL, dqPBL
+      REAL, dimension(klon,klev) :: spread
+      REAL, dimension(klon,klev) :: d_deltatgw
+      REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
+      REAL, dimension(klon,klev+1) :: omgbdth, omg
+      REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
+      REAL, dimension(klon,klev) :: d_deltat_gw
+      REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar
+      REAL, dimension(klon) :: wdens
+      INTEGER, dimension(klon) :: ktopw
+
+c Variables internes
+c-------------------
+
+c Variables à fixer
+      REAL ALON
+      REAL coefgw
+      REAL :: wdens_ref
+      REAL stark
+      REAL alpk
+      REAL delta_t_min
+      INTEGER nsub
+      REAL dtimesub
+      REAL sigmad, hwmin,wapecut
+      REAL :: sigmaw_max
+      REAL :: dens_rate
+      REAL wdens0
+cIM 080208
+      LOGICAL, dimension(klon) :: gwake
+
+c Variables de sauvegarde
+      REAL, dimension(klon,klev) :: deltatw0
+      REAL, dimension(klon,klev) :: deltaqw0
+      REAL, dimension(klon,klev) :: te, qe
+      REAL, dimension(klon) :: sigmaw0, sigmaw1
+
+c Variables pour les GW
+      REAL, DIMENSION(klon) :: LL
+      REAL, dimension(klon,klev) :: N2
+      REAL, dimension(klon,klev) :: Cgw
+      REAL, dimension(klon,klev) :: Tgw
+
+c Variables liées au calcul de hw
+      REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new
+      REAL, DIMENSION(klon) :: sum_dth
+      REAL, DIMENSION(klon) :: dthmin
+      REAL, DIMENSION(klon) :: z, dz, hw0
+      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-15/
+
+c Autres variables internes
+      INTEGER isubstep, k, i
+
+      REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu
+      REAL, DIMENSION(klon) :: sum_dq, sum_rho
+      REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn
+      REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu
+      REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho
+      REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn
+
+      REAL, DIMENSION(klon,klev) :: rho, rhow
+      REAL, DIMENSION(klon,klev+1) :: rhoh
+      REAL, DIMENSION(klon,klev) :: rhow_moyen
+      REAL, DIMENSION(klon,klev) :: zh
+      REAL, DIMENSION(klon,klev+1) :: zhh
+      REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2
+
+      REAL, DIMENSION(klon,klev) :: the, thu
+
+!      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
+      REAL, DIMENSION(klon) :: ztop, dztop
+      REAL, DIMENSION(klon,klev) :: alpha_up
+      
+      REAL, dimension(klon) :: RRe1, RRe2
+      REAL :: RRd1, RRd2
+      REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2
+      REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth
+      REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq
+      REAL, DIMENSION(klon,klev) :: omgbdq
+
+      REAL, dimension(klon) :: ff, gg
+      REAL, dimension(klon) :: wape2, Cstar2, heff
+
+      REAL, DIMENSION(klon,klev) :: Crep
+      REAL Crep_upper, Crep_sol
+
+      REAL, DIMENSION(klon,klev) :: ppi
+
+ccc nrlmd
+      real, dimension(klon) :: death_rate,nat_rate
+      real, dimension(klon,klev) :: entr
+      real, dimension(klon,klev) :: detr
+
+C-------------------------------------------------------------------------
+c         Initialisations
+c-------------------------------------------------------------------------
+
+c      print*, 'wake initialisations'
+
+c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
+c-------------------------------------------------------------------------
+
+      DATA wapecut,sigmad, hwmin /5.,.02,10./
+ccc nrlmd
+      DATA sigmaw_max /0.4/
+      DATA dens_rate /0.1/
+ccc
+C Longueur de maille (en m)
+c-------------------------------------------------------------------------
+
+c      ALON = 3.e5
+      ALON = 1.e6
+
+
+C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
+c
+c      coefgw : Coefficient pour les ondes de gravité
+c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
+c       wdens : Densité de poche froide par maille
+c-------------------------------------------------------------------------
+
+ccc nrlmd      coefgw=10
+c      coefgw=1
+c      wdens0 = 1.0/(alon**2)
+ccc nrlmd      wdens = 1.0/(alon**2)
+ccc nrlmd      stark = 0.50
+cCRtest
+ccc nrlmd      alpk=0.1
+c      alpk = 1.0
+c      alpk = 0.5
+c      alpk = 0.05
+c
+       stark  = 0.33
+       Alpk   = 0.25
+       wdens_ref  = 8.e-12
+       coefgw = 4.
+      Crep_upper=0.9
+      Crep_sol=1.0
+
+ccc nrlmd Lecture du fichier wake_param.data
+c$$$      OPEN(99,file='wake_param.data',status='old',
+c$$$     $          form='formatted',err=9999)
+c$$$      READ(99,*,end=9998) stark
+c$$$      READ(99,*,end=9998) Alpk
+c$$$      READ(99,*,end=9998) wdens_ref
+c$$$      READ(99,*,end=9998) coefgw
+c$$$9998  Continue
+c$$$      CLOSE(99)
+c$$$9999  Continue
+c
+c   Initialisation de toutes des densites a wdens_ref.
+c   Les densites peuvent evoluer si les poches debordent
+c   (voir au tout debut de la boucle sur les substeps)
+      wdens = wdens_ref
+c
+c      print*,'stark',stark
+c      print*,'alpk',alpk
+c      print*,'wdens',wdens
+c      print*,'coefgw',coefgw
+ccc
+C Minimum value for |T_wake - T_undist|. Used for wake top definition
+c-------------------------------------------------------------------------
+
+      delta_t_min = 0.2
+
+C 1. - Save initial values and initialize tendencies
+C --------------------------------------------------
+
+      DO k=1,klev
+      DO i=1, klon
+        ppi(i,k)=pi(i,k)
+	deltatw0(i,k) = deltatw(i,k)
+	deltaqw0(i,k)= deltaqw(i,k)
+	te(i,k) = te0(i,k)
+	qe(i,k) = qe0(i,k)
+	dtls(i,k) = 0.
+	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.
+        d_deltaqw2(i,k)=0.
+!IM 060508 end
+      ENDDO
+      ENDDO
+c      sigmaw1=sigmaw
+c      IF (sigd_con.GT.sigmaw1) THEN
+c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
+c      ENDIF
+      DO i=1, klon
+cc      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
+      sigmaw(i) = amax1(sigmaw(i),sigmad)
+      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
+C
+C 2. - Prognostic part
+C --------------------
+C
+C
+C 2.1 - Undisturbed area and Wake integrals
+C ---------------------------------------------------------
+
+      DO i=1, klon
+      z(i) = 0.
+      ktop(i)=0
+      kupper(i) = 0
+      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 Distance between wakes
+       DO i = 1,klon
+        LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
+       ENDDO
+C Potential temperatures and humidity
+c----------------------------------------------------------
+      DO k =1,klev
+       DO i=1, klon 
+!        write(*,*)'wake 1',i,k,rd,te(i,k)
+        rho(i,k) = p(i,k)/(rd*te(i,k))
+!        write(*,*)'wake 2',rho(i,k)
+        IF(k .eq. 1) THEN
+!        write(*,*)'wake 3',i,k,rd,te(i,k)
+          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
+!        write(*,*)'wake 4',i,k,rd,te(i,k)
+          zhh(i,k)=0
+        ELSE
+!          write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
+          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
+!          write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
+          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
+        ENDIF
+!          write(*,*)'wake 7',ppi(i,k)
+        the(i,k) = te(i,k)/ppi(i,k)
+        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
+        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
+        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
+!          write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
+        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
+        dth(i,k) = deltatw(i,k)/ppi(i,k)
+       ENDDO
+      ENDDO
+        
+      DO k = 1, klev-1
+      DO i=1, klon 
+        IF(k.eq.1) THEN
+          N2(i,k)=0
+        ELSE
+          N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-
+     $            the(i,k-1))/(p(i,k+1)-p(i,k-1)))
+        ENDIF
+        ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2
+
+        Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k)
+        Tgw(i,k)=coefgw*Cgw(i,k)/LL(i)
+      ENDDO
+      ENDDO
+
+      DO i=1, klon
+      N2(i,klev)=0
+      ZH(i,klev)=0
+      Cgw(i,klev)=0
+      Tgw(i,klev)=0
+      ENDDO
+
+c  Calcul de la masse volumique moyenne de la colonne   (bdlmd)
+c-----------------------------------------------------------------
+
+      DO k=1,klev
+       DO i=1, klon
+        epaisseur1(i,k)=0.
+        epaisseur2(i,k)=0.
+       ENDDO
+      ENDDO
+
+      DO i=1, klon
+      epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
+      epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
+      rhow_moyen(i,1) = rhow(i,1)
+      ENDDO
+
+      DO k = 2, klev
+      DO i=1, klon
+        epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.
+        epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)
+        rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+
+     $                 rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)
+      ENDDO
+      ENDDO
+
+C
+C Choose an integration bound well above wake top
+c-----------------------------------------------------------------
+c
+C       Pupper = 50000.  ! melting level
+c       Pupper = 60000.
+c       Pupper = 80000.  ! essais pour case_e
+       DO i = 1,klon
+        Pupper(i) = 0.6*ph(i,1)
+        Pupper(i) = max(Pupper(i), 45000.)
+ccc        Pupper(i) = 60000.
+       ENDDO
+
+C
+C    Determine Wake top pressure (Ptop) from buoyancy integral
+C    --------------------------------------------------------
+c
+c-1/ Pressure of the level where dth becomes less than delta_t_min.
+
+      DO i=1,klon
+      ptop_provis(i)=ph(i,1)
+      ENDDO
+      DO k= 2,klev
+      DO i=1,klon
+c
+cIM v3JYG; ptop_provis(i).LT. ph(i,1)
+c
+        IF (dth(i,k) .GT. -delta_t_min .and.
+     $      dth(i,k-1).LT. -delta_t_min .and.
+     $      ptop_provis(i).EQ. ph(i,1)) THEN
+          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
+     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
+     $          (dth(i,k) - dth(i,k-1))
+        ENDIF
+      ENDDO
+      ENDDO
+
+c-2/ dth integral
+
+      DO i=1,klon
+      sum_dth(i) = 0.
+      dthmin(i) = -delta_t_min
+      z(i) = 0.
+      ENDDO
+
+      DO k = 1,klev
+      DO i=1,klon
+        dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
+        IF (dz(i) .gt. 0) THEN
+          z(i) = z(i)+dz(i)
+          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
+          dthmin(i) = amin1(dthmin(i),dth(i,k))
+        ENDIF
+      ENDDO
+      ENDDO
+
+c-3/ height of triangle with area= sum_dth and base = dthmin
+
+      DO i=1,klon
+      hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
+      hw0(i) = amax1(hwmin,hw0(i))
+      ENDDO
+
+c-4/ now, get Ptop
+
+      DO i=1,klon
+      z(i) = 0.
+      ptop(i) = ph(i,1)
+      ENDDO
+
+      DO k = 1,klev
+      DO i=1,klon
+        dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))
+        IF (dz(i) .gt. 0) THEN
+         z(i) = z(i)+dz(i)
+         ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
+        ENDIF
+      ENDDO
+      ENDDO
+
+
+C-5/ Determination de ktop et kupper
+
+      DO k=klev,1,-1
+      DO i=1,klon
+        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
+        IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
+      ENDDO
+      ENDDO
+
+c      On evite kupper = 1
+      DO i=1,klon
+        kupper(i) = max(kupper(i),2)
+      ENDDO
+
+
+c-6/ Correct ktop and ptop
+
+      DO i = 1,klon
+        ptop_new(i)=ptop(i)
+      ENDDO
+      DO k= klev,2,-1
+      DO i=1,klon
+        IF (k .LE. ktop(i) .and.
+     $      ptop_new(i) .EQ. ptop(i) .and.
+     $      dth(i,k) .GT. -delta_t_min .and.
+     $      dth(i,k-1).LT. -delta_t_min) THEN
+          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
+     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
+     $          (dth(i,k) - dth(i,k-1))
+        ENDIF
+      ENDDO
+      ENDDO
+
+      DO i=1,klon
+        ptop(i) = ptop_new(i)
+      ENDDO
+
+      DO k=klev,1,-1
+      DO i=1,klon
+        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
+      ENDDO
+      ENDDO
+c
+c-5/ Set deltatw & deltaqw to 0 above kupper
+c
+      DO k = 1,klev
+      DO i=1,klon
+       IF (k.GE. kupper(i)) THEN
+        deltatw(i,k) = 0.
+        deltaqw(i,k) = 0.
+       ENDIF
+      ENDDO
+      ENDDO
+c
+C
+C Vertical gradient of LS omega
+C
+      DO k = 1,klev
+      DO i=1,klon
+       IF (k.LE. kupper(i)) THEN
+        dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
+       ENDIF
+      ENDDO
+      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) = -(amax1(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 2.1 - 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)
+c      av_thve = sum_thve/hw0
+      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)
+
+      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 2.2 Prognostic variable update
+C ------------------------------
+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) = amax1(sigmad,sigd_con(i))
+        fip(i) = 0.
+        gwake(i) = .FALSE.
+      ELSE
+        Cstar(i) = stark*sqrt(2.*wape(i))
+        gwake(i) = .TRUE.
+      ENDIF
+      ENDDO
+
+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
+C    -----------------
+C
+      nsub=10
+      dtimesub=dtime/nsub
+c
+c------------------------------------------------------------
+      DO isubstep = 1,nsub
+c------------------------------------------------------------
+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
+ccc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement négatif de ktop à kupper --------
+ccc           On calcule pour cela une densité wdens0 pour laquelle on aurait un entrainement nul ---
+      DO i=1,klon
+cc       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
+cc     $           isubstep,wk_adv(i),cstar(i),wape(i)
+        IF (wk_adv(i) .AND. cstar(i).GT.0.01) THEN
+           omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
+     $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
+           wdens0 = ( sigmaw(i) / (4.*3.14) ) *
+     $     ( (1.-sigmaw(i)) * omg(i,kupper(i)+1) /
+     $     ( (ph(i,1)-pupper(i)) * cstar(i) )  ) **(2)
+         IF ( wdens(i) .LE. wdens0*1.1 ) THEN
+            wdens(i) = wdens0
+         ENDIF
+cc	   print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
+cc     $     ,ph(i,1)-pupper(i)',
+cc     $             omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
+cc     $     ,ph(i,1)-pupper(i)
+        ENDIF
+      ENDDO
+
+ccc nrlmd
+
+      DO i=1,klon
+       IF (wk_adv(i)) THEN
+        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
+        sigmaw(i)=amin1(sigmaw(i),sigmaw_max)
+       ENDIF
+      ENDDO
+      DO i=1,klon
+        IF (wk_adv(i)) THEN
+ccc nrlmd          Introduction du taux de mortalité des poches et test sur sigmaw_max=0.4
+ccc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
+           IF (sigmaw(i).ge.sigmaw_max) THEN
+           death_rate(i)=gfl(i)*Cstar(i)/sigmaw(i)
+           ELSE
+             death_rate(i)=0.
+           END IF
+        d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
+     $               - death_rate(i)*sigmaw(i)*dtimesub
+c     $              - nat_rate(i)*sigmaw(i)*dtimesub
+cc        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
+cc     $  death_rate(i),ktop(i),kupper(i)',
+cc     $	         d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
+cc     $  death_rate(i),ktop(i),kupper(i)
+
+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
+C
+c calcul de la difference de vitesse verticale poche - zone non perturbee
+cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
+cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
+cIM 060208 au niveau k=1..?
+      DO k= 1,klev
+      DO i = 1,klon
+      if (wk_adv(i)) THEN !!! nrlmd
+        dp_deltomg(i,k)=0.
+      end if
+      ENDDO
+      ENDDO
+      DO k= 1,klev+1
+      DO i = 1,klon
+      if (wk_adv(i)) THEN !!! nrlmd
+        omg(i,k)=0.
+      end if
+      ENDDO
+      ENDDO
+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( 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)
+          dp_deltomg(i,k)= dp_deltomg(i,1)
+          omg(i,k)= dp_deltomg(i,1)*z(i)
+       ENDIF
+      ENDDO
+      ENDDO
+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
+c        -----------------
+c        From m/s to Pa/s
+c        -----------------
+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( 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)
+       ENDIF
+      ENDDO
+      ENDDO
+c
+c   raccordement lineaire de omg de ptop a pupper
+
+      DO i=1,klon
+      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(i))
+      ENDIF
+      ENDDO
+c
+cc      DO i=1,klon
+cc        print*,'Pente entre 0 et kupper (référence)'
+cc     $   	,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
+cc        print*,'Pente entre ktop et kupper'
+cc     $  	,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
+cc      ENDDO
+cc
+      DO k= 1,klev
+      DO i = 1,klon
+       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))
+       ENDIF
+      ENDDO
+      ENDDO
+ccc nrlmd
+cc      DO i=1,klon
+cc      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
+cc      END DO
+ccc
+c
+c
+c--    Compute wake average vertical velocity omgbw
+c
+c
+      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
+c--    and its vertical gradient dp_omgbw
+c
+      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
+C
+c--    Upstream coefficients for omgb velocity
+c--    (alpha_up(k) is the coefficient of the value at level k)
+c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
+      DO k = 1,klev
+      DO i=1,klon
+        IF ( wk_adv(i)) THEN
+         alpha_up(i,k) = 0.
+         IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
+        ENDIF
+      ENDDO
+      ENDDO
+
+c  Matrix expressing [The,deltatw] from  [Th1,Th2]
+
+      DO i=1,klon
+        IF ( wk_adv(i)) THEN
+         RRe1(i) = 1.-sigmaw(i)
+         RRe2(i) = sigmaw(i)
+        ENDIF
+      ENDDO
+      RRd1 = -1.
+      RRd2 = 1.
+c
+c--    Get [Th1,Th2], dth and [q1,q2]
+c
+      DO k= 1,klev
+      DO i = 1,klon
+       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
+        Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k)   ! wake
+        q1(i,k) = qe(i,k) - sigmaw(i)     *deltaqw(i,k) ! undisturbed area
+        q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake
+       ENDIF
+      ENDDO
+      ENDDO
+
+      DO i=1,klon
+       if (wk_adv(i)) then !!! nrlmd
+       D_Th1(i,1) = 0.
+       D_Th2(i,1) = 0.
+       D_dth(i,1) = 0.
+       D_q1(i,1) = 0.
+       D_q2(i,1) = 0.
+       D_dq(i,1) = 0.
+       end if
+      ENDDO
+
+      DO k= 2,klev
+      DO i = 1,klon
+       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)
+        D_dth(i,k) = dth(i,k-1)-dth(i,k)
+        D_q1(i,k) = q1(i,k-1)-q1(i,k)
+        D_q2(i,k) = q2(i,k-1)-q2(i,k)
+        D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)
+       ENDIF
+      ENDDO
+      ENDDO
+
+      DO i=1,klon
+        IF( wk_adv(i)) THEN
+         omgbdth(i,1) = 0.
+         omgbdq(i,1) = 0.
+        ENDIF
+      ENDDO
+
+      DO k= 2,klev
+      DO i = 1,klon
+       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))
+       ENDIF
+      ENDDO
+      ENDDO
+c
+c-----------------------------------------------------------------
+      DO k= 1,klev
+      DO i = 1,klon
+       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
+c-----------------------------------------------------------------
+c
+c   Compute redistribution (advective) term
+c
+         d_deltatw(i,k) =
+     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
+     $       RRd1*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
+     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1)
+     $      -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)*
+     $      omgbdth(i,k+1))*ppi(i,k)
+c         print*,'d_deltatw=',d_deltatw(i,k)
+c
+         d_deltaqw(i,k) =
+     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
+     $       RRd1*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
+     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1)
+     $      -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)*
+     $      omgbdq(i,k+1))
+c         print*,'d_deltaqw=',d_deltaqw(i,k)
+c
+c   and increment large scale tendencies
+c
+
+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) )
+     $               /(Ph(i,k)-Ph(i,k+1))
+ccc nrlmd     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
+     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)
+     $         *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1))
+ccc
+     $                      )*ppi(i,k)
+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) )
+     $               /(Ph(i,k)-Ph(i,k+1))
+ccc nrlmd     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
+     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)
+     $         *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1))
+ccc
+     $                      )
+ccc nrlmd
+       ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN
+         d_te(i,k) =  dtimesub*(
+     $       ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
+     $        /(Ph(i,k)-Ph(i,k+1)))
+     $                       )*ppi(i,k)
+
+         d_qe(i,k) =  dtimesub*(
+     $       ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)       
+     $        /(Ph(i,k)-Ph(i,k+1)))
+     $                       )
+
+       ENDIF
+ccc
+      ENDDO
+      ENDDO
+c------------------------------------------------------------------
+C
+C   Increment state variables
+
+      DO k= 1,klev
+      DO i = 1,klon
+ccc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
+        IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
+ccc
+
+
+c
+c Coefficient de répartition
+
+        Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i))
+     $          -ph(i,1))
+        Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-
+     $          ph(i,kupper(i)))
+        
+
+c Reintroduce compensating subsidence term.
+
+c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
+c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
+c     .                   /(1-sigmaw)
+c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
+c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
+c     .                   /(1-sigmaw)
+c
+c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
+c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
+c     .                   /(1-sigmaw)
+c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
+c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
+c     .                   /(1-sigmaw)
+
+        dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))
+        dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))
+c        print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
+c
+        dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i)))
+        dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i)))
+c        print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
+c
+ccc nrlmd          Prise en compte du taux de mortalité
+ccc               Définitions de entr, detr
+        detr(i,k)=0.
+
+        entr(i,k)=detr(i,k)+gfl(i)*cstar(i)+
+     $          sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i,k)
+
+        spread(i,k) = (entr(i,k)-detr(i,k))/sigmaw(i)
+ccc        spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
+ccc     $  sigmaw(i)
+
+
+c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
+
+!      write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
+!     &  Tgw(i,k),deltatw(i,k)
+        d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)*
+     $  dtimesub
+!      write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
+        ff(i)=d_deltatw(i,k)/dtimesub
+
+c Sans GW
+c
+c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
+c
+c GW formule 1
+c
+c        deltatw(k) = deltatw(k)+dtimesub*
+c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
+c
+c GW formule 2
+
+        IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
+          d_deltatw(i,k) = dtimesub*
+     $       (ff(i)+dtKE(i,k)+dtPBL(i,k)
+ccc     $       -spread(i,k)*deltatw(i,k)
+     $       - entr(i,k)*deltatw(i,k)/sigmaw(i)
+     $       - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)
+     $       / (1.-sigmaw(i))
+ccc
+     $       -Tgw(i,k)*deltatw(i,k))
+        ELSE
+           d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub*
+     $       Tgw(i,k)))*
+     $       (ff(i)+dtKE(i,k)+dtPBL(i,k)
+ccc     $       -spread(i,k)*deltatw(i,k)
+     $       - entr(i,k)*deltatw(i,k)/sigmaw(i)
+     $       - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)
+     $       / (1.-sigmaw(i))
+ccc
+     $       -Tgw(i,k)*deltatw(i,k))
+        ENDIF
+
+        dth(i,k) = deltatw(i,k)/ppi(i,k)
+
+        gg(i)=d_deltaqw(i,k)/dtimesub
+
+       d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k)
+ccc     $     -spread(i,k)*deltaqw(i,k))
+     $        - entr(i,k)*deltaqw(i,k)/sigmaw(i)
+     $        - (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)
+     $        /(1.-sigmaw(i)))
+ccc
+
+ccc nrlmd
+ccc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
+ccc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
+ccc
+       ENDIF
+      ENDDO
+      ENDDO
+
+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
+ccc nrlmd
+cc      print*,'alpha'
+cc      do i=1,klon
+cc         print*,alpha(i)
+cc      end do
+ccc
+      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( 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)
+ccc nrlmd
+        d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
+        d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
+ccc
+       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)
+        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)
+cc      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
+cc     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
+       ENDIF
+      ENDDO
+      ENDDO
+      DO i = 1,klon
+       IF( wk_adv(i)) THEN
+        sigmaw(i) = sigmaw(i)+d_sigmaw(i)
+       ENDIF
+      ENDDO
+c
+C
+c     Determine Ptop from buoyancy integral
+c     ---------------------------------------
+c
+c-     1/ Pressure of the level where dth changes sign.
+c
+      DO i=1,klon
+       IF ( wk_adv(i)) THEN
+        Ptop_provis(i)=ph(i,1)
+       ENDIF
+      ENDDO
+c
+      DO k= 2,klev
+      DO i=1,klon
+        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
+          Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
+     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
+     $          - dth(i,k-1))
+        ENDIF
+      ENDDO
+      ENDDO
+c
+c-     2/ dth integral
+c
+      DO i=1,klon
+       if (wk_adv(i)) then !!! nrlmd
+       sum_dth(i) = 0.
+       dthmin(i) = -delta_t_min
+       z(i) = 0.
+       end if
+      ENDDO
+
+      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
+         z(i) = z(i)+dz(i)
+         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
+         dthmin(i) = amin1(dthmin(i),dth(i,k))
+        ENDIF
+       ENDIF
+      ENDDO
+      ENDDO
+c
+c-     3/ height of triangle with area= sum_dth and base = dthmin
+
+      DO i=1,klon
+       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
+c-     4/ now, get Ptop
+c
+      DO i=1,klon
+       if (wk_adv(i)) then !!! nrlmd
+       ktop(i) = 0
+       z(i)=0.
+       end if
+      ENDDO
+c
+      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
+         z(i) = z(i)+dz(i)
+         Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i)
+         ktop(i) = k
+        ENDIF
+       ENDIF
+      ENDDO
+      ENDDO
+c
+c      4.5/Correct ktop and ptop
+c
+      DO i=1,klon
+       IF ( wk_adv(i)) THEN
+        Ptop_new(i)=ptop(i)
+       ENDIF
+      ENDDO
+c
+      DO k= klev,2,-1
+      DO i=1,klon
+cIM v3JYG; IF (k .GE. ktop(i)
+       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.
+     $      dth(i,k-1).LT. -delta_t_min) THEN
+          Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
+     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
+     $          - dth(i,k-1))
+        ENDIF
+      ENDDO
+      ENDDO
+c
+c
+      DO i=1,klon
+       IF ( wk_adv(i)) THEN
+        ptop(i) = ptop_new(i)
+       ENDIF
+      ENDDO
+
+      DO k=klev,1,-1
+      DO i=1,klon
+      if (wk_adv(i)) then !!! nrlmd
+        IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
+      end if
+      ENDDO
+      ENDDO
+c
+c      5/ Set deltatw & deltaqw to 0 above kupper
+c
+      DO k = 1,klev
+      DO i=1,klon
+        IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
+         deltatw(i,k) = 0.
+         deltaqw(i,k) = 0.
+        ENDIF
+      ENDDO
+      ENDDO
+c
+C
+c-------------Cstar computation---------------------------------
+      DO i=1, klon
+       if (wk_adv(i)) then !!! nrlmd
+      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.
+       end if
+      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
+       if (wk_adv(i)) then !!! nrlmd
+      z(i) = 1.
+      dz(i) = 1.
+      sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
+      sum_dth(i) = 0.
+       end if
+      ENDDO
+
+      DO k = 1,klev
+      DO i=1,klon
+       if (wk_adv(i)) then !!! nrlmd
+        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
+       end if
+      ENDDO
+      ENDDO
+c
+      DO i=1,klon
+       if (wk_adv(i)) then !!! nrlmd
+        hw0(i) = z(i)
+       end if
+      ENDDO
+c
+C
+C - WAPE and mean forcing computation
+C ---------------------------------------
+C
+C ---------------------------------------
+C
+C Means
+
+      DO i=1,klon
+       if (wk_adv(i)) then !!! nrlmd
+       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)
+       end if
+      ENDDO
+C
+C Filter out bad wakes
+
+      DO k = 1,klev
+       DO i=1,klon
+        if (wk_adv(i)) then !!! nrlmd
+        IF ( wape(i) .LT. 0.) THEN
+          deltatw(i,k) = 0.
+          deltaqw(i,k) = 0.
+          dth(i,k) = 0.
+        ENDIF
+	end if
+       ENDDO
+      ENDDO
+c
+      DO i=1,klon
+      if (wk_adv(i)) then !!! nrlmd
+      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
+      end if
+      ENDDO
+
+       ENDDO      ! end sub-timestep loop
+C
+C -----------------------------------------------------------------
+c   Get back to tendencies per second
+c
+      DO k = 1,klev
+      DO i=1,klon
+
+ccc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
+        IF ( OK_qx_qw(i) .AND. k .LE. kupper(i)) THEN
+ccc
+         dtls(i,k) = dtls(i,k)/dtime
+         dqls(i,k) = dqls(i,k)/dtime
+         d_deltatw2(i,k)=d_deltatw2(i,k)/dtime
+         d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime
+         d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime
+cc      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
+cc     $         ,death_rate(i)*sigmaw(i)
+        ENDIF
+      ENDDO
+      ENDDO
+
+c
+c----------------------------------------------------------
+c   Determine wake final state; recompute wape, cstar, ktop;
+c   filter out bad wakes.
+c----------------------------------------------------------
+c
+C 2.1 - Undisturbed area and Wake integrals
+C ---------------------------------------------------------
+
+      DO i=1,klon
+ccc nrlmd       if (wk_adv(i)) then !!! nrlmd
+      if (OK_qx_qw(i)) then
+ccc
+        z(i) = 0.
+        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.
+       end if	
+      ENDDO
+C Potential temperatures and humidity
+c----------------------------------------------------------
+
+      DO k =1,klev
+      DO i=1,klon
+ccc nrlmd       IF ( wk_adv(i)) THEN
+       if (OK_qx_qw(i)) then
+ccc
+        rho(i,k) = p(i,k)/(rd*te(i,k))
+        IF(k .eq. 1) THEN
+          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
+          zhh(i,k)=0
+        ELSE
+          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
+          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
+        ENDIF
+        the(i,k) = te(i,k)/ppi(i,k)
+        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
+        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
+        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
+        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
+
+C Integrals (and wake top level number)
+C -----------------------------------------------------------
+
+C Initialize sum_thvu to 1st level virt. pot. temp.
+
+      DO i=1,klon
+ccc nrlmd       IF ( wk_adv(i)) THEN
+      if (OK_qx_qw(i)) then
+ccc
+        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
+ccc nrlmd       IF ( wk_adv(i)) THEN
+       if (OK_qx_qw(i)) then
+ccc
+        dz(i) = -(amax1(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
+       ENDIF
+      ENDDO
+      ENDDO
+c
+      DO i=1,klon
+ccc nrlmd       IF ( wk_adv(i)) THEN
+       if (OK_qx_qw(i)) then
+ccc
+        hw0(i) = z(i)
+       ENDIF
+      ENDDO
+c
+C - WAPE and mean forcing computation
+C-------------------------------------------------------------
+
+C Means
+
+      DO i=1, klon
+ccc nrlmd       IF ( wk_adv(i)) THEN
+      if (OK_qx_qw(i)) then
+ccc
+        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)
+
+        wape2(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)
+       ENDIF
+      ENDDO
+
+C Prognostic variable update
+C ------------------------------------------------------------
+
+C Filter out bad wakes
+c
+      DO k = 1,klev
+      DO i=1,klon
+ccc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
+      if (OK_qx_qw(i) .AND. wape2(i) .LT. 0.) then
+ccc
+          deltatw(i,k) = 0.
+          deltaqw(i,k) = 0.
+          dth(i,k) = 0.
+        ENDIF
+      ENDDO
+      ENDDO
+c
+
+      DO i=1, klon
+ccc nrlmd       IF ( wk_adv(i)) THEN
+      if (OK_qx_qw(i)) then
+ccc
+       IF ( wape2(i) .LT. 0.) THEN
+        wape2(i) = 0.
+        Cstar2(i) = 0.
+        hw(i) = hwmin
+        sigmaw(i) = amax1(sigmad,sigd_con(i))
+        fip(i) = 0.
+        gwake(i) = .FALSE.
+      ELSE
+        if(prt_level.ge.10) print*,'wape2>0'
+        Cstar2(i) = stark*sqrt(2.*wape2(i))
+        gwake(i) = .TRUE.
+      ENDIF
+      ENDIF
+      ENDDO
+c
+      DO i=1, klon
+ccc nrlmd       IF ( wk_adv(i)) THEN
+       if (OK_qx_qw(i)) then
+ccc
+        ktopw(i) = ktop(i)
+       ENDIF
+      ENDDO
+c
+      DO i=1, klon
+ccc nrlmd       IF ( wk_adv(i)) THEN
+       if (OK_qx_qw(i)) then
+ccc
+       IF (ktopw(i) .gt. 0 .and. gwake(i)) then
+
+Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
+ccc       heff = 600.
+C      Utilisation de la hauteur hw
+cc       heff = 0.7*hw
+       heff(i) = hw(i)
+
+       FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2*
+     $      sqrt(sigmaw(i)*wdens(i)*3.14)
+       FIP(i) = alpk * FIP(i)
+Cjyg2
+       ELSE
+         FIP(i) = 0.
+       ENDIF
+       ENDIF
+      ENDDO
+c
+C   Limitation de sigmaw
+
+ccc nrlmd
+c       DO i=1,klon
+c         IF (OK_qx_qw(i)) THEN
+c	   IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
+c	 ENDIF
+c       ENDDO
+ccc
+      DO k = 1,klev
+       DO i=1, klon
+
+ccc nrlmd      On maintient désormais constant sigmaw en régime permanent
+ccc      IF ((sigmaw(i).GT.sigmaw_max).or.
+        IF     ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
+     $         (ktopw(i).le.2) .OR.
+     $         .not. OK_qx_qw(i) ) THEN
+ccc
+          dtls(i,k) = 0.
+          dqls(i,k) = 0.
+          deltatw(i,k) = 0.
+          deltaqw(i,k) = 0.
+        ENDIF
+       ENDDO
+      ENDDO
+c
+ccc nrlmd      On maintient désormais constant sigmaw en régime permanent
+      DO i=1, klon
+        IF  ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
+     $        (ktopw(i).le.2) .OR.
+     $        .not. OK_qx_qw(i)   ) THEN
+         wape(i) = 0.
+         cstar(i)=0.
+         hw(i) = hwmin
+         sigmaw(i) = sigmad
+         fip(i) = 0.
+        ELSE
+         wape(i) = wape2(i)
+         cstar(i)=cstar2(i)
+        ENDIF
+cc        print*,'wape wape2 ktopw OK_qx_qw =',
+cc     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
+      ENDDO
+c
+c
+      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 zeta(nlon,nl)
+      REAL alpha1(nlon)
+      REAL x,a,b,c,discrim
+      REAL epsilon
+!      DATA epsilon/1.e-15/
+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,k)=0.
+        ELSE
+         zeta(i,k)=1.
+        END IF
+       ENDIF
+      ENDDO
+      DO i = 1,nlon
+       IF (wk_adv(i)) THEN
+        x = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k)
+     $    + d_qe(i,k)+(zeta(i,k)-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,k)-sigmaw(i))*d_deltaqw(i,k)
+     $    - deltaqw(i,k)*d_sigmaw(i)
+        c = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k)+epsilon
+        discrim = b*b-4.*a*c
+c      print*, 'x, a, b, c, discrim', x, a, b, c, discrim
+        IF (a+b .GE. 0.) THEN !! Condition suffisante pour la positivité de ovap
+         alpha1(i)=1.
+        ELSE
+         IF (x .GE. 0.) THEN
+            alpha1(i)=1.
+         ELSE
+              IF (a .GT. 0.) THEN
+                 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
+c         print*,'a,b,c discrim',a,b,c discrim
+                 alpha1(i)=0.9*max(   (2.*c)/(-b+sqrt(discrim)),
+     $                        (-b+sqrt(discrim))/(2.*a)   )
+              ENDIF
+         ENDIF
+        ENDIF
+       alpha(i) = min(alpha(i),alpha1(i))
+       ENDIF
+      ENDDO
+      ENDDO
+!
+      return
+      end
+
+      Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
+     :                ,te0,qe0,omgb
+     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
+     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
+     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
+     o                ,dtls,dqls
+     o                ,ktopw,omgbdth,dp_omgb,wdens
+     o                ,tu,qu
+     o                ,dtKE,dqKE
+     o                ,dtPBL,dqPBL
+     o                ,omg,dp_deltomg,spread
+     o                ,Cstar,d_deltat_gw
+     o                ,d_deltatw2,d_deltaqw2)
+
+***************************************************************
+*                                                             *
+* WAKE                                                        *
+*      retour a un Pupper fixe                                *
+*                                                             *
+* written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
+* modified by :   ROEHRIG Romain        01/29/2007            *
+***************************************************************
+c
+      USE dimphy
       IMPLICIT none
 c============================================================================
@@ -113,1643 +1942,4 @@
 
 #include "dimensions.h"
-#include "YOMCST.h"
-#include "cvthermo.h"
-#include "iniprint.h"
-
-c Arguments en entree
-c--------------------
-
-      REAL, dimension(klon,klev) :: p, ppi
-      REAL, dimension(klon,klev+1) :: ph, omgb
-      REAL dtime
-      REAL, dimension(klon,klev) :: te0,qe0
-      REAL, dimension(klon,klev) :: dtdwn, dqdwn
-      REAL, dimension(klon,klev) :: wdtPBL,wdqPBL
-      REAL, dimension(klon,klev) :: udtPBL,udqPBL
-      REAL, dimension(klon,klev) :: amdwn, amup
-      REAL, dimension(klon,klev) :: dta, dqa
-      REAL, dimension(klon) :: sigd_con
-
-c Sorties
-c--------
-
-      REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
-      REAL, dimension(klon,klev) :: tu, qu
-      REAL, dimension(klon,klev) :: dtls, dqls
-      REAL, dimension(klon,klev) :: dtKE, dqKE
-      REAL, dimension(klon,klev) :: dtPBL, dqPBL
-      REAL, dimension(klon,klev) :: spread
-      REAL, dimension(klon,klev) :: d_deltatgw
-      REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
-      REAL, dimension(klon,klev+1) :: omgbdth, omg
-      REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
-      REAL, dimension(klon,klev) :: d_deltat_gw
-      REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar
-      INTEGER, dimension(klon) :: ktopw
-
-c Variables internes
-c-------------------
-
-c Variables à fixer
-      REAL ALON
-      REAL coefgw
-      REAL :: wdens0, wdens
-      REAL stark
-      REAL alpk
-      REAL delta_t_min
-      INTEGER nsub
-      REAL dtimesub
-      REAL sigmad, hwmin
-      REAL :: sigmaw_max
-cIM 080208
-      LOGICAL, dimension(klon) :: gwake
-
-c Variables de sauvegarde
-      REAL, dimension(klon,klev) :: deltatw0
-      REAL, dimension(klon,klev) :: deltaqw0
-      REAL, dimension(klon,klev) :: te, qe
-      REAL, dimension(klon) :: sigmaw0, sigmaw1
-
-c Variables pour les GW
-      REAL, DIMENSION(klon) :: LL
-      REAL, dimension(klon,klev) :: N2
-      REAL, dimension(klon,klev) :: Cgw
-      REAL, dimension(klon,klev) :: Tgw
-
-c Variables liées au calcul de hw
-      REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new
-      REAL, DIMENSION(klon) :: sum_dth
-      REAL, DIMENSION(klon) :: dthmin
-      REAL, DIMENSION(klon) :: z, dz, hw0
-      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)
-
-c Autres variables internes
-      INTEGER isubstep, k, i
-
-      REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu
-      REAL, DIMENSION(klon) :: sum_dq, sum_rho
-      REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn
-      REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu
-      REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho
-      REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn
-
-      REAL, DIMENSION(klon,klev) :: rho, rhow
-      REAL, DIMENSION(klon,klev+1) :: rhoh
-      REAL, DIMENSION(klon,klev) :: rhow_moyen
-      REAL, DIMENSION(klon,klev) :: zh
-      REAL, DIMENSION(klon,klev+1) :: zhh
-      REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2
-
-      REAL, DIMENSION(klon,klev) :: the, thu
-
-!      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
-      REAL, DIMENSION(klon) :: ztop, dztop
-      REAL, DIMENSION(klon,klev) :: alpha_up
-      
-      REAL, dimension(klon) :: RRe1, RRe2
-      REAL :: RRd1, RRd2
-      REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2, T1 
-      REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth
-      REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq
-      REAL, DIMENSION(klon,klev) :: omgbdq
-
-      REAL, dimension(klon) :: ff, gg
-      REAL, dimension(klon) :: wape2, Cstar2, heff
-
-      REAL, DIMENSION(klon,klev) :: Crep
-      REAL Crep_upper, Crep_sol
-
-C-------------------------------------------------------------------------
-c         Initialisations
-c-------------------------------------------------------------------------
-
-c      print*, 'wake initialisations'
-
-c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
-c-------------------------------------------------------------------------
-
-      DATA sigmad, hwmin /.02,10./
-
-C Longueur de maille (en m)
-c-------------------------------------------------------------------------
-
-c      ALON = 3.e5
-      ALON = 1.e6
-
-
-C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
-c
-c      coefgw : Coefficient pour les ondes de gravité
-c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
-c       wdens : Densité de poche froide par maille
-c-------------------------------------------------------------------------
-
-      coefgw=10
-c      coefgw=1
-c      wdens0 = 1.0/(alon**2)   
-      wdens = 1.0/(alon**2)       
-      stark = 0.50
-cCRtest
-      alpk=0.1
-c      alpk = 1.0 
-c      alpk = 0.5
-c      alpk = 0.05
-      Crep_upper=0.9
-      Crep_sol=1.0
-
-C Minimum value for |T_wake - T_undist|. Used for wake top definition
-c-------------------------------------------------------------------------
-
-      delta_t_min = 0.2
-
-C 1. - Save initial values and initialize tendencies
-C --------------------------------------------------
-
-      DO k=1,klev
-      DO i=1, klon
-	deltatw0(i,k) = deltatw(i,k)
-	deltaqw0(i,k)= deltaqw(i,k)
-	te(i,k) = te0(i,k)
-	qe(i,k) = qe0(i,k)
-	dtls(i,k) = 0.
-	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.
-        d_deltaqw2(i,k)=0.
-!IM 060508 end
-      ENDDO
-      ENDDO
-c      sigmaw1=sigmaw
-c      IF (sigd_con.GT.sigmaw1) THEN
-c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
-c      ENDIF
-      DO i=1, klon
-cc      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
-      sigmaw(i) = amax1(sigmaw(i),sigmad)
-      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
-C
-C 2. - Prognostic part
-C --------------------
-C
-C
-C 2.1 - Undisturbed area and Wake integrals
-C ---------------------------------------------------------
-
-      DO i=1, klon
-      z(i) = 0.
-      ktop(i)=0
-      kupper(i) = 0
-      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 Distance between wakes
-       DO i = 1,klon
-        LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens)
-       ENDDO
-C Potential temperatures and humidity
-c----------------------------------------------------------
-      DO k =1,klev
-       DO i=1, klon 
-        rho(i,k) = p(i,k)/(rd*te(i,k))
-        IF(k .eq. 1) THEN
-          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
-          zhh(i,k)=0
-        ELSE
-          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
-          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
-        ENDIF
-        the(i,k) = te(i,k)/ppi(i,k)
-        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
-        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
-        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
-        rhow(i,k) = p(i,k)/(rd*(tu(i,k)+deltatw(i,k)))
-        dth(i,k) = deltatw(i,k)/ppi(i,k)
-       ENDDO
-      ENDDO
-        
-      DO k = 1, klev-1
-      DO i=1, klon 
-        IF(k.eq.1) THEN
-          N2(i,k)=0
-        ELSE
-          N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-
-     $            the(i,k-1))/(p(i,k+1)-p(i,k-1)))
-        ENDIF
-        ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2
-
-        Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k)
-        Tgw(i,k)=coefgw*Cgw(i,k)/LL(i)
-      ENDDO
-      ENDDO
-
-      DO i=1, klon
-      N2(i,klev)=0
-      ZH(i,klev)=0
-      Cgw(i,klev)=0
-      Tgw(i,klev)=0
-      ENDDO
-
-c  Calcul de la masse volumique moyenne de la colonne   (bdlmd)
-c-----------------------------------------------------------------
-
-      DO k=1,klev
-       DO i=1, klon
-        epaisseur1(i,k)=0.
-        epaisseur2(i,k)=0.
-       ENDDO
-      ENDDO
-
-      DO i=1, klon
-      epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
-      epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
-      rhow_moyen(i,1) = rhow(i,1)
-      ENDDO
-
-      DO k = 2, klev
-      DO i=1, klon
-        epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.
-        epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)
-        rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+
-     $                 rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)
-      ENDDO
-      ENDDO
-
-C
-C Choose an integration bound well above wake top
-c-----------------------------------------------------------------
-c
-C       Pupper = 50000.  ! melting level
-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
-C    --------------------------------------------------------
-c
-c-1/ Pressure of the level where dth becomes less than delta_t_min.
-
-      DO i=1,klon
-      ptop_provis(i)=ph(i,1)
-      ENDDO
-      DO k= 2,klev
-      DO i=1,klon
-c
-cIM v3JYG; ptop_provis(i).LT. ph(i,1)
-c
-        IF (dth(i,k) .GT. -delta_t_min .and.
-     $      dth(i,k-1).LT. -delta_t_min .and.
-     $      ptop_provis(i).EQ. ph(i,1)) THEN
-          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
-     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
-     $          (dth(i,k) - dth(i,k-1))
-        ENDIF
-      ENDDO
-      ENDDO
-
-c-2/ dth integral
-
-      DO i=1,klon
-      sum_dth(i) = 0.
-      dthmin(i) = -delta_t_min
-      z(i) = 0.
-      ENDDO
-
-      DO k = 1,klev
-      DO i=1,klon
-        dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
-        IF (dz(i) .gt. 0) THEN
-          z(i) = z(i)+dz(i)
-          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
-          dthmin(i) = amin1(dthmin(i),dth(i,k))
-        ENDIF
-      ENDDO
-      ENDDO
-
-c-3/ height of triangle with area= sum_dth and base = dthmin
-
-      DO i=1,klon
-      hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
-      hw0(i) = amax1(hwmin,hw0(i))
-      ENDDO
-
-c-4/ now, get Ptop
-
-      DO i=1,klon
-      z(i) = 0.
-      ptop(i) = ph(i,1)
-      ENDDO
-
-      DO k = 1,klev
-      DO i=1,klon
-        dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))
-        IF (dz(i) .gt. 0) THEN
-         z(i) = z(i)+dz(i)
-         ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
-        ENDIF
-      ENDDO
-      ENDDO
-
-
-C-5/ Determination de ktop et kupper
-
-      DO k=klev,1,-1
-      DO i=1,klon
-        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
-        IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
-      ENDDO
-      ENDDO
-
-c-6/ Correct ktop and ptop
-
-      DO i = 1,klon
-        ptop_new(i)=ptop(i)
-      ENDDO
-      DO k= klev,2,-1
-      DO i=1,klon
-        IF (k .LE. ktop(i) .and.
-     $      ptop_new(i) .EQ. ptop(i) .and.
-     $      dth(i,k) .GT. -delta_t_min .and.
-     $      dth(i,k-1).LT. -delta_t_min) THEN
-          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
-     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
-     $          (dth(i,k) - dth(i,k-1))
-        ENDIF
-      ENDDO
-      ENDDO
-
-      DO i=1,klon
-        ptop(i) = ptop_new(i)
-      ENDDO
-
-      DO k=klev,1,-1
-      DO i=1,klon
-        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
-      ENDDO
-      ENDDO
-c
-c-5/ Set deltatw & deltaqw to 0 above kupper
-c
-      DO k = 1,klev
-      DO i=1,klon
-       IF (k.GE. kupper(i)) THEN
-        deltatw(i,k) = 0.
-        deltaqw(i,k) = 0.
-       ENDIF
-      ENDDO
-      ENDDO
-c
-C
-C Vertical gradient of LS omega
-C
-      DO k = 1,klev
-      DO i=1,klon
-       IF (k.LE. kupper(i)) THEN
-        dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
-       ENDIF
-      ENDDO
-      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) = -(amax1(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 2.1 - 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)
-c      av_thve = sum_thve/hw0
-      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)
-
-      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 2.2 Prognostic variable update
-C ------------------------------
-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) = amax1(sigmad,sigd_con(i))
-        fip(i) = 0.
-        gwake(i) = .FALSE.
-      ELSE
-        Cstar(i) = stark*sqrt(2.*wape(i))
-        gwake(i) = .TRUE.
-      ENDIF
-      ENDDO
-
-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
-C    -----------------
-C
-      nsub=10
-      dtimesub=dtime/nsub
-c
-c------------------------------------------------------------
-      DO isubstep = 1,nsub
-c------------------------------------------------------------
-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))
-        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
-C
-c calcul de la difference de vitesse verticale poche - zone non perturbee
-cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
-cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
-cIM 060208 au niveau k=1..?
-      DO k= 1,klev
-      DO i = 1,klon
-        dp_deltomg(i,k)=0.
-      ENDDO
-      ENDDO
-      DO k= 1,klev+1
-      DO i = 1,klon
-        omg(i,k)=0.
-      ENDDO
-      ENDDO
-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( 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)
-          dp_deltomg(i,k)= dp_deltomg(i,1)
-          omg(i,k)= dp_deltomg(i,1)*z(i)
-       ENDIF
-      ENDDO
-      ENDDO
-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
-c        -----------------
-c        From m/s to Pa/s
-c        -----------------
-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( 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)
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c   raccordement lineaire de omg de ptop a pupper
-
-      DO i=1,klon
-      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(i))
-      ENDIF
-      ENDDO
-c
-      DO k= 1,klev
-      DO i = 1,klon
-       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))
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c
-c--    Compute wake average vertical velocity omgbw
-c
-c
-      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
-c--    and its vertical gradient dp_omgbw
-c
-      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
-C
-c--    Upstream coefficients for omgb velocity
-c--    (alpha_up(k) is the coefficient of the value at level k)
-c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
-      DO k = 1,klev
-      DO i=1,klon
-        IF ( wk_adv(i)) THEN
-         alpha_up(i,k) = 0.
-         IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
-        ENDIF
-      ENDDO
-      ENDDO
-
-c  Matrix expressing [The,deltatw] from  [Th1,Th2]
-
-      DO i=1,klon
-        IF ( wk_adv(i)) THEN
-         RRe1(i) = 1.-sigmaw(i)
-         RRe2(i) = sigmaw(i)
-        ENDIF
-      ENDDO
-      RRd1 = -1.
-      RRd2 = 1.
-c
-c--    Get [Th1,Th2], dth and [q1,q2]
-c
-      DO k= 1,klev
-      DO i = 1,klon
-       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
-        Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k)   ! wake
-        q1(i,k) = qe(i,k) - sigmaw(i)     *deltaqw(i,k) ! undisturbed area
-        q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake
-        T1(i,k) = te(i,k) - sigmaw(i)*deltatw(i,k)! undisturb itlmd
-       ENDIF
-      ENDDO
-      ENDDO
-
-      DO i=1,klon
-       D_Th1(i,1) = 0.   !!!itlmd : ne pas mettre if wk_adv cf nrlmd?
-       D_Th2(i,1) = 0.
-       D_dth(i,1) = 0.
-       D_q1(i,1) = 0.
-       D_q2(i,1) = 0.
-       D_dq(i,1) = 0.
-      ENDDO
-
-      DO k= 2,klev
-      DO i = 1,klon
-       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)
-        D_dth(i,k) = dth(i,k-1)-dth(i,k)
-        D_q1(i,k) = q1(i,k-1)-q1(i,k)
-        D_q2(i,k) = q2(i,k-1)-q2(i,k)
-        D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)
-       ENDIF
-      ENDDO
-      ENDDO
-
-      DO i=1,klon
-        IF( wk_adv(i)) THEN
-         omgbdth(i,1) = 0.
-         omgbdq(i,1) = 0.
-        ENDIF
-      ENDDO
-
-      DO k= 2,klev
-      DO i = 1,klon
-       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))
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c-----------------------------------------------------------------
-      DO k= 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
-c-----------------------------------------------------------------
-c
-c   Compute redistribution (advective) term
-c
-         d_deltatw(i,k) =
-     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
-     $       RRd1*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
-     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1)
-     $      -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)*
-     $      omgbdth(i,k+1))*ppi(i,k)
-c         print*,'d_deltatw=',d_deltatw(i,k)
-c
-         d_deltaqw(i,k) =
-     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
-     $       RRd1*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
-     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1)
-     $      -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)*
-     $      omgbdq(i,k+1))
-c         print*,'d_deltaqw=',d_deltaqw(i,k)
-c
-c   and increment large scale tendencies
-c
-
-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) )
-     $               /(Ph(i,k)-Ph(i,k+1))
-     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*
-     $            (omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) !instead of dp_deltomg(i,k) 
-     $                      )*ppi(i,k)
-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) )
-     $               /(Ph(i,k)-Ph(i,k+1))
-     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*
-     $           (omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1))!instead of dp_deltomg(i,k)
-     $                      )
-        ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN ! corr pour conserver l'eau
-
-         d_te(i,k) =  dtimesub*(
-     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k))
-     $               /(Ph(i,k)-Ph(i,k+1))
-     $                      )*ppi(i,k)
-
-         d_qe(i,k) =  dtimesub*(
-     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k))
-     $               /(Ph(i,k)-Ph(i,k+1))
-     $                      )
-       ENDIF
-
-c-------------------------------------------------------------------
-      ENDDO
-      ENDDO
-c------------------------------------------------------------------
-C
-C   Increment state variables
-
-      DO k= 1,klev
-      DO i = 1,klon
-       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
-c
-c Coefficient de répartition
-
-        Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i))
-     $          -ph(i,1))
-        Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-
-     $          ph(i,kupper(i)))
-        
-
-c Reintroduce compensating subsidence term.
-
-c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
-c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
-c     .                   /(1-sigmaw)
-c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
-c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
-c     .                   /(1-sigmaw)
-c
-c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
-c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
-c     .                   /(1-sigmaw)
-c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
-c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
-c     .                   /(1-sigmaw)
-
-        dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))
-        dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))
-c        print*,'dtKE=',dtKE(k)
-c        print*,'dqKE=',dqKE(k)
-c
-        dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i)))
-        dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i)))
-c
-        spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
-     $  sigmaw(i)
-
-
-c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
-
-        d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)* 
-     $  dtimesub
-        ff(i)=d_deltatw(i,k)/dtimesub
-
-c Sans GW
-c
-c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) 
-c
-c GW formule 1
-c
-c        deltatw(k) = deltatw(k)+dtimesub*
-c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
-c
-c GW formule 2
-
-        IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
-          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
-           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
-
-       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)
-       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
-       ENDIF
-      ENDDO
-      ENDDO
-
-C
-C   Scale tendencies so that water vapour remains positive in w and x.
-C
-      call wake_vec_modulation(klon,klev,wk_adv,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( 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)
-        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
-C
-c     Determine Ptop from buoyancy integral
-c     ---------------------------------------
-c
-c-     1/ Pressure of the level where dth changes sign.
-c
-      DO i=1,klon
-       IF ( wk_adv(i)) THEN
-        Ptop_provis(i)=ph(i,1)
-       ENDIF
-      ENDDO
-c
-      DO k= 2,klev
-      DO i=1,klon
-        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
-          Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
-     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
-     $          - dth(i,k-1))
-        ENDIF
-      ENDDO
-      ENDDO
-c
-c-     2/ dth integral
-c
-      DO i=1,klon
-       sum_dth(i) = 0.
-       dthmin(i) = -delta_t_min
-       z(i) = 0.
-      ENDDO
-
-      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
-         z(i) = z(i)+dz(i)
-         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
-         dthmin(i) = amin1(dthmin(i),dth(i,k))
-        ENDIF
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c-     3/ height of triangle with area= sum_dth and base = dthmin
-
-      DO i=1,klon
-       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
-c-     4/ now, get Ptop
-c
-      DO i=1,klon
-       ktop(i) = 0
-       z(i)=0.
-      ENDDO
-c
-      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
-         z(i) = z(i)+dz(i)
-         Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i)
-         ktop(i) = k
-        ENDIF
-       ENDIF
-      ENDDO
-      ENDDO
-c
-c      4.5/Correct ktop and ptop
-c
-      DO i=1,klon
-       IF ( wk_adv(i)) THEN
-        Ptop_new(i)=ptop(i)
-       ENDIF
-      ENDDO
-c
-      DO k= klev,2,-1
-      DO i=1,klon
-cIM v3JYG; IF (k .GE. ktop(i)
-       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.
-     $      dth(i,k-1).LT. -delta_t_min) THEN
-          Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
-     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
-     $          - dth(i,k-1))
-        ENDIF
-      ENDDO
-      ENDDO
-c
-c
-      DO i=1,klon
-       IF ( wk_adv(i)) THEN
-        ptop(i) = ptop_new(i)
-       ENDIF
-      ENDDO
-
-      DO k=klev,1,-1
-      DO i=1,klon
-        IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
-      ENDDO
-      ENDDO
-c
-c      5/ Set deltatw & deltaqw to 0 above kupper
-c
-      DO k = 1,klev
-      DO i=1,klon
-        IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
-         deltatw(i,k) = 0.
-         deltaqw(i,k) = 0.
-        ENDIF
-      ENDDO
-      ENDDO
-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) + th1(i,k)*dz(i)
-         sum_tu(i) = sum_tu(i) + t1(i,k)*dz(i)
-         sum_qu(i) = sum_qu(i) + q1(i,k)*dz(i)
-         sum_thvu(i) = sum_thvu(i) + th1(i,k)*(1.+eps*q1(i,k))*dz(i)!itlmd
-
-         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
-C -----------------------------------------------------------------
-c   Get back to tendencies per second
-c
-      DO k = 1,klev
-      DO i=1,klon
-       IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN  !! corr conservation eau
-         dtls(i,k) = dtls(i,k)/dtime
-         dqls(i,k) = dqls(i,k)/dtime
-         d_deltatw2(i,k)=d_deltatw2(i,k)/dtime
-         d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime
-         d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime
-        ENDIF
-      ENDDO
-      ENDDO
-c
-c
-c----------------------------------------------------------
-c   Determine wake final state; recompute wape, cstar, ktop;
-c   filter out bad wakes.
-c----------------------------------------------------------
-c
-C 2.1 - Undisturbed area and Wake integrals
-C ---------------------------------------------------------
-
-      DO i=1,klon
-        z(i) = 0.
-        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 Potential temperatures and humidity
-c----------------------------------------------------------
-
-      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
-          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
-          zhh(i,k)=0
-        ELSE
-          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
-          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
-        ENDIF
-        the(i,k) = te(i,k)/ppi(i,k)
-        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
-        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
-        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
-        rhow(i,k) = p(i,k)/(rd*(tu(i,k)+deltatw(i,k)))
-        dth(i,k) = deltatw(i,k)/ppi(i,k)
-       ENDIF
-      ENDDO
-      ENDDO
-
-C Integrals (and wake top level number)
-C -----------------------------------------------------------
-
-C Initialize sum_thvu to 1st level virt. pot. temp.
-
-      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
-         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
-       ENDIF
-      ENDDO
-      ENDDO
-c
-      DO i=1,klon
-       IF ( wk_adv(i)) THEN
-        hw0(i) = z(i)
-       ENDIF
-      ENDDO
-c
-C - WAPE and mean forcing computation
-C-------------------------------------------------------------
-
-C Means
-
-      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)
-        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)
-
-        wape2(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)
-       ENDIF
-      ENDDO
-
-C Prognostic variable update
-C ------------------------------------------------------------
-
-C Filter out bad wakes
-c
-      DO k = 1,klev
-      DO i=1,klon
-        IF ( wk_adv(i) .AND. wape2(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 ( wk_adv(i)) THEN
-       IF ( wape2(i) .LT. 0.) THEN
-        wape2(i) = 0.
-        Cstar2(i) = 0.
-        hw(i) = hwmin
-        sigmaw(i) = amax1(sigmad,sigd_con(i))
-        fip(i) = 0.
-        gwake(i) = .FALSE.
-      ELSE
-        if(prt_level.ge.10) print*,'wape2>0'
-        Cstar2(i) = stark*sqrt(2.*wape2(i))
-        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 ( wk_adv(i)) THEN
-       IF (ktopw(i) .gt. 0 .and. gwake(i)) then
-
-Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
-ccc       heff = 600.
-C      Utilisation de la hauteur hw
-cc       heff = 0.7*hw
-       heff(i) = hw(i)
-
-       FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2*
-     $      sqrt(sigmaw(i)*wdens*3.14)
-       FIP(i) = alpk * FIP(i)
-Cjyg2
-       ELSE
-         FIP(i) = 0.
-       ENDIF
-       ENDIF
-      ENDDO
-c
-C   Limitation de sigmaw
-c
-C   sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
-C              alors il disparait en se mélangeant à la partie undisturbed
-c
-      sigmaw_max = 0.9
-      DO k = 1,klev
-       DO i=1, klon
-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) .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
-          dtls(i,k) = 0.
-          dqls(i,k) = 0.
-          deltatw(i,k) = 0.
-          deltaqw(i,k) = 0.
-        ENDIF
-       ENDDO
-      ENDDO
-c
-      DO i=1, klon
-         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.
-         cstar(i)= 0.  !!corr itlmd
-         hw(i) = hwmin
-         sigmaw(i) = sigmad
-         fip(i) = 0.
-        ELSE
-         wape(i) = wape2(i)
-         cstar(i)= cstar2(i) !!corr itlmd
-        ENDIF
-      ENDDO
-c
-c
-      RETURN
-      END
-
-      SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,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-15/
-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
-     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
-     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
-     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
-     o                ,dtls,dqls
-     o                ,ktopw,omgbdth,dp_omgb,wdens
-     o                ,tu,qu
-     o                ,dtKE,dqKE
-     o                ,dtPBL,dqPBL
-     o                ,omg,dp_deltomg,spread
-     o                ,Cstar,d_deltat_gw
-     o                ,d_deltatw2,d_deltaqw2)
-
-***************************************************************
-*                                                             *
-* WAKE                                                        *
-*      retour a un Pupper fixe                                *
-*                                                             *
-* written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
-* modified by :   ROEHRIG Romain        01/29/2007            *
-***************************************************************
-c
-      USE dimphy
-      IMPLICIT none
-c============================================================================
-C
-C
-C   But : Decrire le comportement des poches froides apparaissant dans les
-C        grands systemes convectifs, et fournir l'energie disponible pour
-C        le declenchement de nouvelles colonnes convectives.
-C
-C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
-C                      deltaqw    : ecart d'humidite wake-undisturbed area
-C                      sigmaw     : fraction d'aire occupee par la poche.
-C
-C   Variable de sortie : 
-c
-c			 wape : WAke Potential Energy
-c                        fip  : Front Incident Power (W/m2) - ALP
-c                        gfl  : Gust Front Length per unit area (m-1)
-C                        dtls : large scale temperature tendency due to wake
-C                        dqls : large scale humidity tendency due to wake
-C                        hw   : hauteur de la poche
-C                     dp_omgb : vertical gradient of large scale omega
-C                      omgbdth: flux of Delta_Theta transported by LS omega
-C                      dtKE   : differential heating (wake - unpertubed)
-C                      dqKE   : differential moistening (wake - unpertubed)
-C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
-C                 dp_deltomg  : vertical gradient of omg (s-1)
-C                     spread  : spreading term in dt_wake and dq_wake
-C                 deltatw     : updated temperature difference (T_w-T_u).
-C                 deltaqw     : updated humidity difference (q_w-q_u).
-C                 sigmaw      : updated wake fractional area.
-C                 d_deltat_gw : delta T tendency due to GW
-c
-C   Variables d'entree : 
-c
-c		         aire : aire de la maille
-c			 te0  : temperature dans l'environnement  (K)
-C                        qe0  : humidite dans l'environnement     (kg/kg)
-C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
-C                        dtdwn: source de chaleur due aux descentes (K/s)
-C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
-C			 dta  : source de chaleur due courants satures et detrain  (K/s)
-C			 dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
-C                        amdwn: flux de masse total des descentes, par unite de
-C                                surface de la maille (kg/m2/s)
-C                        amup : flux de masse total des ascendances, par unite de
-C                                surface de la maille (kg/m2/s)
-C                        p    : pressions aux milieux des couches (Pa)
-C                        ph   : pressions aux interfaces (Pa)
-C                        ppi  : (p/p_0)**kapa (adim)
-C                        dtime: increment temporel (s)
-c
-C   Variables internes :
-c
-c			 rhow : masse volumique de la poche froide
-C                        rho  : environment density at P levels
-C                        rhoh : environment density at Ph levels
-C                        te   : environment temperature | may change within
-C                        qe   : environment humidity    | sub-time-stepping
-C                        the  : environment potential temperature
-C                        thu  : potential temperature in undisturbed area
-C                        tu   :  temperature  in undisturbed area
-C                        qu   : humidity in undisturbed area
-C                      dp_omgb: vertical gradient og LS omega
-C                      omgbw  : wake average vertical omega
-C                     dp_omgbw: vertical gradient of omgbw
-C                      omgbdq : flux of Delta_q transported by LS omega
-C                        dth  : potential temperature diff. wake-undist.
-C                        th1  : first pot. temp. for vertical advection (=thu)
-C                        th2  : second pot. temp. for vertical advection (=thw)
-C                        q1   : first humidity for vertical advection
-C                        q2   : second humidity for vertical advection
-C                     d_deltatw   : terme de redistribution pour deltatw
-C                     d_deltaqw   : terme de redistribution pour deltaqw
-C                      deltatw0   : deltatw initial
-C                      deltaqw0   : deltaqw initial
-C                      hw0    : hw initial
-C                      sigmaw0: sigmaw initial
-C                      amflux : horizontal mass flux through wake boundary
-C                      wdens  : number of wakes per unit area (3D) or per
-C                               unit length (2D)
-C                      Tgw    : 1 sur la période de onde de gravité
-c                      Cgw    : vitesse de propagation de onde de gravité
-c                      LL     : distance entre 2 poches
-
-c-------------------------------------------------------------------------
-c          Déclaration de variables
-c-------------------------------------------------------------------------
-
-#include "dimensions.h"
 cccc#include "dimphy.h"
 #include "YOMCST.h"
