subroutine get_kinetic(Hsfc1, Hbnd, Hthk1, Uvec1, Usfc1, Wvec1, energy_k, sigmaXX, sigmaXZ)
  use prm    , only : IMAX, JMAX, KMAX, dx, dt
  use prmphys, only : rho, grav
  use com_var, only : x, zk, zk_half, wgt
  use u_pack,  only : get_Ubnd
  implicit none
  
  real(8), INTENT(in) :: Hsfc1(1:IMAX, 1:JMAX)
  real(8), INTENT(in) :: Hbnd(1:IMAX, 1:JMAX)
  real(8), INTENT(in) :: Hthk1(1:IMAX, 1:JMAX)
  real(8), INTENT(in) :: Uvec1(0:IMAX, 1:JMAX, 1:KMAX)
  real(8), INTENT(in) :: Usfc1( 0:IMAX, 1:JMAX)
  real(8)             :: Ubnd( 0:IMAX, 1:JMAX)
  real(8)            :: WORK_U(0:IMAX, 1:JMAX, 0:KMAX+1)

  real(8), INTENT(in) :: Wvec1(1:IMAX, 1:JMAX, 0:KMAX)
  real(8), INTENT(OUT):: energy_k(1:IMAX, 1:JMAX, 1:KMAX)
  real(8),INTENT(OUT) :: sigmaXX(0:IMAX, 1:JMAX, 0:KMAX+1)
  real(8),INTENT(OUT) :: sigmaXZ(0:IMAX, 1:JMAX, 0:KMAX+1)
  real(8)             :: WORK_wgt(0:KMAX+1)
  integer           :: i, j, k 


  call get_Ubnd(Hsfc1, Hbnd, Ubnd)

  ! $B:n6HJQ?t(B
  WORK_U(:, :, 0) = Ubnd(:, :)
  WORK_U(:, :, 1:KMAX) = Uvec1(:, :, 1:KMAX)
  WORK_U(:, :, KMAX+1) = Usfc1(:, :)
  
  WORK_wgt(0) = wgt(1)/2.0D0
  WORK_wgt(1:KMAX) = wgt(1:KMAX)
  WORK_wgt(KMAX+1) = wgt(KMAX)/2.0D0

  ! $B&R(Bxx $B!A(B $B&R(Bzz
  sigmaXX(0,:, :) =  0.0D0
  sigmaXX(IMAX,:, :) =  0.0D0
  do i = 1, IMAX-1
     sigmaXX(i, :, 0) = rho*grav*(Hthk1(i+1, :)+Hthk1(i, :))/2.0D0 
  end do
  sigmaXX(1:IMAX-1, :, KMAX+1) = 0.0D0
  do k =1, KMAX
     do i = 1, IMAX-1
        sigmaXX(i, :, k) =  &
             & rho * grav*(Hsfc1(i+1, :)+Hsfc1(i, :))/2.0D0 * (1.0D0 - zk(k))
     end do
  end do
  ! $B&R(Bxz
  sigmaXZ(0, :, :) = 0.0D0
  sigmaXZ(IMAX, :, :) = 0.0D0
  do i = 1, IMAX-1
     sigmaXZ(i, :, 0) = &
          & sigmaXX(i,:,0)*(Hsfc1(i+1, :) - Hsfc1(i, :))/dx
  end do
  sigmaXZ(1:IMAX-1, :, KMAX+1) = 0.0D0
  
  do k = 1, KMAX
     do i = 1, IMAX-1
        do j = 1, JMAX
           sigmaXZ(i, j, k) = &
                & - sigmaXX(i, j, k)*(Hsfc1(i+1, j) - Hsfc1(i, j))/dx
        end do
     end do
  end do
  
  energy_k(1, :, :)=0.0D0
  energy_k(IMAX, :, :)=0.0D0

  do i = 2,IMAX-1
     do j =1, JMAX
!        energy_k(i, j, 1) =   dt*&  ! x $BLLJ}8~(B
!                & (Uvec1(i, j, 1)*sigmaXX(i, j, 1)      &
!                &    Uvec1(i-1, j, 1)*sigmaXX(i-1, j, 1))&
!                &   * Hthk1(i, j)*wgt(1)  

        do k=1, KMAX
           energy_k(i, j, k) =   dt*&
                & (Uvec1(i, j, k)*sigmaXX(i, j, k)      &
                &    -Uvec1(i-1, j, k)*sigmaXX(i-1, j, k)) & 
                &* Hthk1(i, j)*wgt(k)  !& 
!                & +dt *(Wvec1(i, j, k)*(sigmaXX(i, j, k)+sigmaXX(i-1, j, k))/2.0D0      &
           !                &    -Wvec1(i, j, k-1)*(sigmaXX(i, j, k-1)+sigmaXX(i-1, j,k-1)/2.0D0)) * dx  
           



!           energy_k(i, j, k)=  dt* & ! $B2<It(B
!                &(((WORK_U(i, j, k-1)*WORK_wgt(k) &
!                &      +WORK_U(i, j, k)*WORK_wgt(k-1))&
!                &           /(WORK_wgt(k)+WORK_wgt(k-1))&
!                & * (sigmaXX(i, j, k-1)*WORK_wgt(k)                         &
!                &         +sigmaXX(i, j, k)*WORK_wgt(k-1))                  &
!                &         /(WORK_wgt(k)+WORK_wgt(k-1))                        &
!                &   ! $B>eIt(B
!                & + (WORK_U(i, j, k)*WORK_wgt(k+1)&
!                & +WORK_U(i, j, k+1)*WORK_wgt(k))&
!                & /(WORK_wgt(k+1)+WORK_wgt(k))&
!                & * (sigmaXX(i, j, k)*WORK_wgt(k+1)                         &
!                &              +sigmaXX(i, j, k+1)*WORK_wgt(k))             &
!                &         /(WORK_wgt(k+1)+WORK_wgt(k)))/2.0D0               &
!                &  * (Hthk1(i, j)+Hthk1(i-1, j))/2.0D0*WORK_wgt(k)          &
!                & ! $B3J;R$N1&B&(B($B2<It(B)
!                & - ((WORK_U(i+1, j, k-1)*WORK_wgt(k)&
!                & +WORK_U(i+1, j, k)*WORK_wgt(k-1))    &
!                &         /(WORK_wgt(k)+WORK_wgt(k-1))                       &
!                & * (sigmaXX(i+1, j, k-1)*WORK_wgt(k)&
!                & +sigmaXX(i+1, j, k)*WORK_wgt(k-1))&
!                &         /(WORK_wgt(k)+WORK_wgt(k-1))                       &
!                &    ! $B>eIt(B
!                & + (WORK_U(i+1, j, k)*WORK_wgt(k+1)&
!                & +WORK_U(i+1, j, k+1)*WORK_wgt(k))    &
!                &         /(WORK_wgt(k+1)+WORK_wgt(k))                       &
!                & * (sigmaXX(i+1, j, k)*WORK_wgt(k+1)&
!                &  +sigmaXX(i+1, j, k+1)*WORK_wgt(k))&
!                &         /(WORK_wgt(k+1)+WORK_wgt(k)))/2.0D0                & 
!                &  * (Hthk1(i+1, j)+Hthk1(i, j))/2.0D0*WORK_wgt(k))
!           energy_k(i, j, k) = energy_k(i, j, k)+    dt*                &
!                &  (((Wvec1(i, j, k-1)+Wvec1(i-1, k, k-1))/2.0D0          &
!                & * (sigmaXZ(i+1, j, k-1)*WORK_wgt(k)                         &
!                &         +sigmaXZ(i+1, j, k)*WORK_wgt(k-1))                  &
!                &         /(WORK_wgt(k)+WORK_wgt(k-1))                      &
!                &   ! $B>eIt(B
!                & + (Wvec1(i, j, k)+Wvec1(i-1, k, k))/2.0D0                 &
!                & * (sigmaXZ(i+1, j, k)*WORK_wgt(k+1)                       &
!                &              +sigmaXZ(i+1, j, k+1)*WORK_wgt(k))           &
!                &         /(WORK_wgt(k+1)+WORK_wgt(k)))/2.0D0               &
!                &  * (Hthk1(i, j)+Hthk1(i-1, j))/2.0D0*WORK_wgt(k)          &
!                & ! $B3J;R$N1&B&(B($B2<It(B)
!                & + ((Wvec1(i+1, j, k-1)+Wvec1(i, k, k-1))/2.0D0            &
!                & * (sigmaXZ(i, j, k-1)*WORK_wgt(k)&
!                &          +sigmaXZ(i, j, k)*WORK_wgt(k-1))&
!                &         /(WORK_wgt(k)+WORK_wgt(k-1))                            &
!                &    ! $B>eIt(B
!                & + (Wvec1(i+1, j, k)+Wvec1(i, k, k))/2.0D0          &
!                & * (sigmaXZ(i, j, k)*WORK_wgt(k+1)                  &
!                &    +sigmaXZ(i, j, k+1)*WORK_wgt(k))                &
!                &         /(WORK_wgt(k+1)+WORK_wgt(k)))/2.0D0        &  
!                &  * (Hthk1(i+1, j)+Hthk1(i, j))/2.0D0*WORK_wgt(k))
        end do
        
     end do
  end do
!  do k = 1, KMAX
!     write(6, *)'energy_kinetic',k,  energy_k(1:IMAX, 2, k)
!  end do
end subroutine get_kinetic
