Index: /LMDZ6/trunk/libf/phylmd/lmdz_wake.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lmdz_wake.F90	(revision 4815)
+++ /LMDZ6/trunk/libf/phylmd/lmdz_wake.F90	(revision 4816)
@@ -227,5 +227,5 @@
 
   ! Variables liees au calcul de hw
-  REAL, DIMENSION (klon)                                :: ptop_provis, ptop, ptop_new
+  REAL, DIMENSION (klon)                                :: ptop
   REAL, DIMENSION (klon)                                :: sum_dth
   REAL, DIMENSION (klon)                                :: dthmin
@@ -311,4 +311,5 @@
 
   LOGICAL                                               :: first_call=.true.
+
 
   ! -------------------------------------------------------------------------
@@ -697,113 +698,13 @@
   ! Determine Wake top pressure (Ptop) from buoyancy integral
   ! --------------------------------------------------------
-
-  ! -1/ Pressure of the level where dth becomes less than delta_t_min.
-
-  DO i = 1, klon
-    ptop_provis(i) = ph(i, 1)
-  END DO
-  DO k = 2, klev
-    DO i = 1, klon
-
-      ! IM v3JYG; ptop_provis(i).LT. ph(i,1)
-
-! LJYF : mettre cette ligne à la place de la suivante      IF (dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
-      IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. &
-          ptop_provis(i)==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))
-      END IF
-    END DO
-  END DO
-
-  ! -2/ dth integral
-
-  DO i = 1, klon
-    sum_dth(i) = 0.
-    dthmin(i) = -delta_t_min
-    z(i) = 0.
-  END DO
-
-  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)>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))
-      END IF
-    END DO
-  END DO
-
-  ! -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))
-  END DO
-
-  ! -4/ now, get Ptop
-
-  DO i = 1, klon
-    z(i) = 0.
-    ptop(i) = ph(i, 1)
-  END DO
-
-  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)>0) THEN
-        z(i) = z(i) + dz(i)
-        ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i)
-      END IF
-    END DO
-  END DO
-
-  IF (prt_level>=10) THEN
-    PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout)
-  ENDIF
-
-
-  ! -5/ Determination de ktop et kupper
-
-  DO k = klev, 1, -1
-    DO i = 1, klon
-      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
-    END DO
-  END DO
-  !print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper
- 
-
-
-  ! -6/ Correct ktop and ptop
-
-  DO i = 1, klon
-    ptop_new(i) = ptop(i)
-  END DO
-  DO k = klev, 2, -1
-    DO i = 1, klon
-      IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
-! LJYF changer par cette ligne          dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
-          dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-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))
-      END IF
-    END DO
-  END DO
-
-  DO i = 1, klon
-    ptop(i) = ptop_new(i)
-  END DO
-
-  DO k = klev, 1, -1
-    DO i = 1, klon
-      IF (ph(i,k+1)<ptop(i)) ktop(i) = k
-    END DO
-  END DO
-
-  CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
-  
-  IF (prt_level>=10) THEN
-    PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
+   Do i=1, klon
+       wk_adv(i) = .True.
+   Enddo
+   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
+                    dth, hw0, rho, &
+                    ktop, wk_adv)
+   
+   IF (prt_level>=10) THEN
+     PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
   ENDIF
 
@@ -848,4 +749,5 @@
       dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
       IF (dz(i)>0) THEN
+              ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau
         z(i) = z(i) + dz(i)
         sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
@@ -988,4 +890,5 @@
   ! ------------------------------------------------------------------------
   !
+  wk_adv = .True.
   DO isubstep = 1, nsub
   !
@@ -993,4 +896,5 @@
     ! 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) >= 1.
@@ -1826,116 +1730,8 @@
 
 
-    ! Determine Ptop from buoyancy integral
-    ! ---------------------------------------
-
-    ! -     1/ Pressure of the level where dth changes sign.
-
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN
-        ptop_provis(i) = ph(i, 1)
-      END IF
-    END DO
-
-    DO k = 2, klev
-      DO i = 1, klon
-        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
-! LJYF changer :           dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
-            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-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))
-        END IF
-      END DO
-    END DO
-
-    ! -     2/ dth integral
-
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN !!! nrlmd
-        sum_dth(i) = 0.
-        dthmin(i) = -delta_t_min
-        z(i) = 0.
-      END IF
-    END DO
-
-    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)>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))
-          END IF
-        END IF
-      END DO
-    END DO
-
-    ! -     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))
-      END IF
-    END DO
-
-    ! -     4/ now, get Ptop
-
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN !!! nrlmd
-        ktop(i) = 0
-        z(i) = 0.
-      END IF
-    END DO
-
-    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)>0) THEN
-            z(i) = z(i) + dz(i)
-            ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i)
-            ktop(i) = k
-          END IF
-        END IF
-      END DO
-    END DO
-
-    ! 4.5/Correct ktop and ptop
-
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN
-        ptop_new(i) = ptop(i)
-      END IF
-    END DO
-
-    DO k = klev, 2, -1
-      DO i = 1, klon
-        ! IM v3JYG; IF (k .GE. ktop(i)
-        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
-! LJYF changer :           dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
-            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-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))
-        END IF
-      END DO
-    END DO
-
-
-    DO i = 1, klon
-      IF (wk_adv(i)) THEN
-        ptop(i) = ptop_new(i)
-      END IF
-    END DO
-
-    DO k = klev, 1, -1
-      DO i = 1, klon
-        IF (wk_adv(i)) THEN !!! nrlmd
-          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
-        END IF
-      END DO
-    END DO
-
-    CALL pkupper (klon, klev, ptop, ph, pupper, kupper)
+
+   Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
+                    dth, hw, rho, &
+                    ktop, wk_adv)
  
 
@@ -2651,16 +2447,153 @@
 
 
-SUBROUTINE pkupper (klon, klev, ptop, ph, pupper, kupper)
+SUBROUTINE pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
+                    dth, hw_, rho,  &
+                    ktop, wk_adv)
 
 USE lmdz_wake_ini , ONLY : wk_pupper
+USE lmdz_wake_ini , ONLY : RG
+USE lmdz_wake_ini , ONLY : hwmin
 IMPLICIT NONE
 
 INTEGER,  INTENT(IN)                              :: klon,klev
-REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: ph
-REAL,     INTENT(IN),   DIMENSION (klon)          :: ptop
-REAL,     INTENT(OUT),  DIMENSION (klon)          :: pupper
-INTEGER,  INTENT(OUT),  DIMENSION (klon)          :: kupper
+REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: ph, p
+REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: rho
+LOGICAL,  INTENT(IN),   DIMENSION (klon)          :: wk_adv
+REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: dth
+
+REAL,     INTENT(OUT),   DIMENSION (klon)         :: hw_
+REAL,     INTENT(OUT),   DIMENSION (klon)         :: ptop
+INTEGER,  INTENT(OUT),   DIMENSION (klon)         :: Ktop 
+REAL,     INTENT(OUT),   DIMENSION (klon)         :: pupper
+INTEGER,  INTENT(OUT),   DIMENSION (klon)         :: kupper
 INTEGER                                           :: i,k
-
+REAL                                              :: delta_t_min
+
+
+REAL,     DIMENSION (klon)         :: dthmin
+REAL,     DIMENSION (klon)         :: ptop_provis,ptop_new
+REAL,     DIMENSION (klon)         :: z, dz
+REAL,     DIMENSION (klon)         :: sum_dth
+
+! LJYF : a priori z, dz sum_dth sont aussi des variables internes
+! Les eliminer apres verification convergence numerique
+
+  delta_t_min = 0.2
+ 
+    ! Determine Ptop from buoyancy integral
+    ! ---------------------------------------
+
+    ! -     1/ Pressure of the level where dth changes sign.
+    print*,'WAKE LJYF'
+
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN
+        ptop_provis(i) = ph(i, 1)
+      END IF
+    END DO
+
+    DO k = 2, klev
+      DO i = 1, klon
+        IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. &
+! LJYF changer :           dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
+            dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-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))
+        END IF
+      END DO
+    END DO
+
+    ! -     2/ dth integral
+
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN !!! nrlmd
+        sum_dth(i) = 0.
+        dthmin(i) = -delta_t_min
+        z(i) = 0.
+      END IF
+    END DO
+
+    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)>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))
+          END IF
+        END IF
+      END DO
+    END DO
+
+    ! -     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))
+      END IF
+    END DO
+
+    ! -     4/ now, get Ptop
+
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN !!! nrlmd
+        ktop(i) = 0
+        z(i) = 0.
+      END IF
+    END DO
+
+    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)>0) THEN
+            z(i) = z(i) + dz(i)
+            ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i)
+            ktop(i) = k
+          END IF
+        END IF
+      END DO
+    END DO
+
+    ! 4.5/Correct ktop and ptop
+
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN
+        ptop_new(i) = ptop(i)
+      END IF
+    END DO
+
+    DO k = klev, 2, -1
+      DO i = 1, klon
+        ! IM v3JYG; IF (k .GE. ktop(i)
+        IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. &
+! LJYF changer :           dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN
+            dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-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))
+        END IF
+      END DO
+    END DO
+
+
+    DO i = 1, klon
+      IF (wk_adv(i)) THEN
+        ptop(i) = ptop_new(i)
+      END IF
+    END DO
+
+    DO k = klev, 1, -1
+      DO i = 1, klon
+        IF (wk_adv(i)) THEN !!! nrlmd
+          IF (ph(i,k+1)<ptop(i)) ktop(i) = k
+        END IF
+      END DO
+    END DO
+  
+!  IF (prt_level>=10) THEN
+!    PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout)
+!  ENDIF
 
  kupper = 0
