MODULE lmdz_wake_pkupper
PUBLIC wake_pkupper
CONTAINS

SUBROUTINE wake_pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
                    dth, hw_, rho, delta_t_min_in, &
                    ktop, wk_adv, h_zzz, ptop1, ktop1)

USE lmdz_wake_ini , ONLY : wk_pupper
USE lmdz_wake_ini , ONLY : RG
USE lmdz_wake_ini , ONLY : hwmin
USE lmdz_wake_ini , ONLY : iflag_wk_new_ptop, wk_delta_t_min, wk_frac_int_delta_t
USE lmdz_wake_ini , ONLY : wk_int_delta_t_min

IMPLICIT NONE

INTEGER,                              INTENT(IN) :: klon,klev
REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: ph, p
REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: rho
LOGICAL,    DIMENSION (klon)        , INTENT(IN) :: wk_adv
REAL,       DIMENSION (klon,klev+1) , INTENT(IN) :: dth
REAL,                                 INTENT(IN) :: delta_t_min_in


REAL,       DIMENSION (klon)  , INTENT(OUT)        :: hw_
REAL,       DIMENSION (klon)  , INTENT(OUT)        :: ptop
INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: Ktop
REAL,       DIMENSION (klon)  , INTENT(OUT)        :: pupper
INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: kupper
REAL,       DIMENSION (klon)  , INTENT(OUT)        :: h_zzz       !!
REAL,       DIMENSION (klon)  , INTENT(OUT)        :: Ptop1      !!
INTEGER,    DIMENSION (klon)  , INTENT(OUT)        :: ktop1      !!

INTEGER :: i,k

LOGICAL,    DIMENSION (klon)       :: wk_active
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

INTEGER,     DIMENSION (klon)                     :: k_ptop_provis
REAL,     DIMENSION (klon)                     :: zk_ptop_provis
REAL,     DIMENSION (klon)                        :: omega        !!
REAL,     DIMENSION (klon,klev+1)                 :: int_dth      !! 
REAL,     DIMENSION (klon,klev+1)                 :: dzz          !!
REAL,     DIMENSION (klon,klev+1)                 :: zzz          !!
REAL,     DIMENSION (klon)                 :: frac_int_dth          !!
REAL                                              :: ddd!!

REAL :: www

INTEGER, SAVE :: ipas=0
!$OMP THREADPRIVATE(ipas) 



!INTEGER, SAVE :: compte=0

! LJYF : a priori z, dz sum_dth sont aussi des variables internes
! Les eliminer apres verification convergence numerique

!compte=compte+1
!print*,'compte=',compte

    ! Determine Ptop from buoyancy integral
    ! ---------------------------------------

    ! -     1/ Pressure of the level where dth changes sign.
    !print*,'WAKE LJYF'


if (iflag_wk_new_ptop==0) then
    delta_t_min=delta_t_min_in
else
    delta_t_min=wk_delta_t_min
endif

    DO i = 1, klon
        ptop_provis(i) = ph(i, 1)
        k_ptop_provis(i) = 1
    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))
          k_ptop_provis(i) = k
        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
        ptop_new(i) = ptop(i)
    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
        ptop(i) = ptop_new(i)
    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

    ! -----------------------------------------------------------------------
    ! nouveau calcul de hw et ptop
    ! -----------------------------------------------------------------------
!if (iflag_wk_new_ptop>0) then
do i=1,klon
   ptop1(i)=ph(i,1)
   ktop1(i)=1
   h_zzz(i)=0.
enddo
    
IF (iflag_wk_new_ptop/=0) THEN
    
    int_dth(1:klon,1:klev+1)=0.
    DO i = 1, klon
       IF (wk_adv(i)) THEN 
          int_dth(i,1) = 0.
      END IF
    END DO
    
    if (abs(iflag_wk_new_ptop) == 1 ) then
        DO k = 2, klev+1
           Do i = 1, klon
              IF (wk_adv(i)) THEN
                 if (k<=k_ptop_provis(i)) then
                      ddd=dth(i,k-1)*(ph(i,k-1) - max(ptop_provis(i),ph(i,k)))
                      !ddd=dth(i,k-1)*(ph(i,k-1) - ph(i,k))
                 else
                      ddd=0.
                 endif             
                 int_dth(i,k) = int_dth(i,k-1) + ddd
              !ELSE
              !   int_dth(i,k) = 0.
              END IF
           END DO
        END DO
    else
        k_ptop_provis(:)=klev+1
        dthmin(:)=dth(:,1)
        ! calcul de l'int??grale de dT * dP jusqu'au dernier
        ! niveau avec dT<0. (en s'assurant qu'on a bien un 
        ! dT negatif plus bas)
        DO k = 1, klev
           DO i = 1, klon
              dthmin(i)=min(dthmin(i),dth(i,k))
              ddd=dth(i,k)*(ph(i,k)-ph(i,k+1))
              if (dthmin(i)<0.) then
                  if (k>=k_ptop_provis(i)) then
                      ddd=0.
                  else if (dth(i,k)>=0.) then
                      ddd=0.
                      k_ptop_provis(i)=k+1
                  endif
              endif
              int_dth(i,k+1) = int_dth(i,k)+ ddd
           ENDDO
        ENDDO

        DO i = 1, klon
           if ( k_ptop_provis(i)==klev+1 .or. .not. wk_adv(i)) then
                k_ptop_provis(i)=1
           endif
        ENDDO
    endif ! (abs(iflag_wk_new_ptop) == 1 )
   ! print*, 'xxx, int_dth', (k,int_dth(1,k),k=1,klev)
   ! print*, 'xxx, k_ptop_provis', k_ptop_provis(1)
   

 
    ! On se limite ?? des poches avec integrale dT * dp < -wk_int_delta_t_min
    do i=1,klon
          if (int_dth(i,k_ptop_provis(i)) > -wk_int_delta_t_min .or. k_ptop_provis(i)==1) then
          !if (1==0) then
             wk_active(i)=.false.
             ptop(i)=ph(i,1)
             ktop(i)=1
             hw_(i)=0.
          else
             wk_active(i)=wk_adv(i)
          endif
    enddo

    DO i=1,klon
       IF (wk_active(i)) THEN
        frac_int_dth(i)=wk_frac_int_delta_t*int_dth(i,k_ptop_provis(i))
       ENDIF
    ENDDO
    DO k = 1,klev
       DO i =1, klon
          IF (wk_active(i)) THEN
            IF (int_dth(i,k)>=frac_int_dth(i)) THEN
              ktop1(i) = min(k, k_ptop_provis(i))
              !ktop1(i) = k
              !print*,ipas,'yyy ktop1= ',ktop1
            ENDIF
          ENDIF
       END DO
    END DO
    !print*, 'LAMINE'
    
    DO i = 1, klon
       IF (wk_active(i)) THEN
           !print*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop1
           ddd=int_dth(i,ktop1(i)+1)-int_dth(i,ktop1(i))
           if (ddd==0.) then
              omega(i)=0.
           else
              omega(i) = (frac_int_dth(i) - int_dth(i,ktop1(i)))/ddd
           endif
           !! print*,'OMEGA ',omega(i)
       END IF
    END DO
    
    !! print*, 'xxx'
    DO i = 1, klon
       IF (wk_active(i)) THEN
      ! print*, 'xxx, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ', &
      !               int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1)
      ! print*, 'xxx, omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ',  &
      !e               omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1)
          ptop1(i) = min((1 - omega(i))*ph(i,ktop1(i)) + omega(i)*ph(i,ktop1(i)+1), ph(i,1))
      END IF
    END DO
    
    DO i=1, klon
       IF (wk_active(i)) THEN 
           zzz(i, 1) = 0
       END IF
     END DO
    DO k = 1, klev
       DO i = 1, klon 
           IF (wk_active(i)) THEN         
              dzz(i,k) = (ph(i,k) - ph(i,k+1))/(rho(i,k)*RG)
              zzz(i,k+1) = zzz(i,k) + dzz(i,k)
           END IF
       END DO
    END DO
    
    DO i =1, klon
       IF (wk_active(i)) THEN
           h_zzz(i) = max((1- omega(i))*zzz(i,ktop1(i)) + omega(i)*zzz(i,ktop1(i)+1), hwmin)
       END IF
    END DO


ENDIF ! (iflag_wk_new_ptop/=0)

!if (iflag_wk_new_ptop==2) then
IF (iflag_wk_new_ptop>0) THEN 
   do i=1,klon
      ptop(i)=ptop1(i)
      ktop(i)=ktop1(i)
      hw_(i)=h_zzz(i)
   enddo

!endif
ENDIF

 kupper = 0
 
IF (0.<wk_pupper .and. wk_pupper<1.) THEN
 ! Choose an integration bound well above wake top
  ! -----------------------------------------------------------------

  ! Pupper = 50000.  ! melting level
  ! Pupper = 60000.
  ! Pupper = 80000.  ! essais pour case_e
  DO i = 1, klon
  !  pupper(i) = 0.6*ph(i, 1)
    pupper(i) = wk_pupper*ph(i, 1)
    pupper(i) = max(pupper(i), 45000.)
    ! cc        Pupper(i) = 60000.
  END DO

ELSE IF (1.<=wk_pupper) THEN
  DO i=1, klon
     ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) 
     !  pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.) 
      pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-5000.) 
  END DO
ELSE
  www=abs(wk_pupper)
  DO i=1, klon
     ! pupper(i) = wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) 
     !  pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-50.) 
      pupper(i) = min(   exp(www*log(ptop(i))+(1.-www)*log(ph(i, 1)) ) , ptop(i)-5000.) 
  END DO
END IF
 
  ! -5/ Determination de kupper

  DO k = klev, 1, -1
    DO i = 1, klon
      IF (ph(i,k+1)<pupper(i)) kupper(i) = k
    END DO
  END DO

  ! On evite kupper = 1 et kupper = klev
  DO i = 1, klon
    kupper(i) = max(kupper(i), 2)
    kupper(i) = min(kupper(i), klev-1)
  END DO 
  !---------- FIN nouveau calcul hw et ptop -------------------------------------

IF (iflag_wk_new_ptop==999) THEN
    DO i = 1, klon
    hw_(i)=0.
    ptop(i)=ph(i,1)
    Ktop(i)=1
    pupper(i)=ph(i,2)
    kupper(i)=2
    h_zzz(i)=0.
    Ptop1(i)=ph(i,1)
    ENDDO
ENDIF

zk_ptop_provis=k_ptop_provis

    RETURN
END SUBROUTINE wake_pkupper
END MODULE lmdz_wake_pkupper