!=====================================================================
! $BB.EYJ,I[7W;;%b%8%e!<%k(B
! 2001/12/08 $B1tD>J,I[%k!<%A%s40N;(B
! 2003/11/02 $B=$O@$+$i0\?"(B
! 2003/11/26 $BI=LLB.EY$b=PNO(B KMAX+1 $BHVL\(B
! 2003/11/27 $B1tD>B.EY=PNO(B
! 2003/12/06 $BDlLL?eJ?B.EY$N%5%V%k!<%A%s:n@.(B($B$?$@$7$9$Y$F(B 0 $BCM(B)
! 2003/01/01 $B?eJ?<ANL%U%i%C%/%9$r3J;RKh$K=PNO(B($B%(%M%k%.!<J]B8$r9M$($k$?$a(B)
! 2004/02/03 $BITMW$J%W%m%0%i%`$N:o=|(B, $B$*$h$S(B Hthk $B$NF3F~(B
!=====================================================================
module u_pack
  
  use prm, only : IMAX, JMAX, KMAX,   &
       &          dx, dt
  use prmphys, only :     rho, grav, Npw, A_con, s_a
  use com_var, only : x, zk, zk_half, wgt, Apara


  implicit none

  private

  public get_U                  ! $B?eJ?B.EY$NF3=P(B
  public get_Ubnd                  ! $B?eJ?B.EY$NF3=P(B
  public get_U_ana                  ! $B?eJ?B.EY$NF3=P(B
  public get_w                  ! $B1tD>B.EY$NF3=P(B
  public get_F_ana              ! $B?eJ?B.EY%U%i%C%/%9$NF3=P(B
  public get_Fgrd               ! $B3F3J;R$K$*$1$k?eJ?<ANL%U%i%C%/%9$NF3=P(B


contains
!=====================================================================
! $B?eJ?B.EYJ,I[(B($BBf7AB'$rMQ$$$?J}K!(B)
!=====================================================================
subroutine get_U(Hsfc1, Hbnd, Hthk1, Usfc1, Uvec1)
  implicit none

  real(8), parameter   :: one = 1.0D0

  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(OUT) :: Usfc1(0:IMAX, 1:JMAX)
  real(8), INTENT(OUT) :: Uvec1(0:IMAX, 1:JMAX, 1:KMAX)

  real(8)              :: Ubnd( 0:IMAX, 1:JMAX)
  real(8)              :: Uscle(0:IMAX, 1:JMAX)

  real(8)              :: Uvec1R(1:IMAX/2+1, 1:JMAX, 1:KMAX)
  real(8)              :: Uvec1L(1:IMAX/2+1, 1:JMAX, 1:KMAX)

  integer              :: i, j, k 

  ! x $BJ}8~$N6-3&$NCM$O$9$Y$F(B 0 
  Uvec1(0,    :, :) = 0.0D0
  Uvec1(IMAX, :, :) = 0.0D0
  Usfc1(:, :) = 0.D0
  
  ! $BDlLLB.EY(B
  call get_Ubnd(Hsfc1, Hbnd, Ubnd)
  call monit_add_2('Hbnd      ', Hbnd,   one)
 
  do j = 1, JMAX
     do i = 1,IMAX-1
        Uscle(i, j) = - 2.0D0 * (rho *grav)** Npw     &
        &             * ((Hsfc1(i+1, j) - Hsfc1(i, j))/dx)**Npw    &
        &             * ((Hthk1(i+1, j) +Hthk1(i, j))/2.0D0)**(Npw+1.0D0)

     ! $B:G2<AX(B
        Uvec1(i, j, 1) = Ubnd(i, j) + Uscle(i, j) &
             &    * (Apara(i+1, j, 1)+Apara(i, j, 1))/2.0D0 &
          &       * (zk_half(1) - zk_half(0))/2.0D0 &
          &       * ((1 - zk(1))** Npw  + 1.0D0) /2.0D0

     ! $B:G2<AX$+$i:G>eAX(B
        do k = 2, KMAX
           Uvec1(i, j, k) = Uvec1(i, j, k-1)                    &
                &  + Uscle(i, j)  * (zk(k) - zk(k-1))           &
                &   * (Apara(i+1, j, k) +Apara(i, j, k))/2.0D0  &
                &  * ((1 - zk(k-1))**Npw + (1 -zk(k))**Npw)/2.0D0
        end do

     ! $BI=LLB.EY(B( 2003/11/26 $BDI5-(B)
        Usfc1(i, j) = Uvec1(i, j, KMAX)                        &
             & + Uscle(i, j) * (Apara(i+1, j, KMAX)+Apara(i, j, KMAX))/2.0D0              &
             & * (zk_half(KMAX) - zk_half(KMAX-1))/2.0D0       &
             & * ((1 - zk(KMAX))**Npw)/2.0D0
     end do
  end do
  call monit_add_2('Usfc1     ', Usfc1(1:IMAX, :)*s_a, one )
  call monit_add_3('Uvec1     ', Uvec1(1:IMAX, :, :)*s_a, one )


  Uvec1L(1:IMAX/2+1, :, :) = Uvec1(0:IMAX/2, :, :)
    do  i = 1, IMAX/2+1
       Uvec1R(i, :, :) = Uvec1(IMAX+1-i, :, :)
    end do


    do j = 1, JMAX
       do i = 1, IMAX/2
          do k = 1, KMAX
             if((abs(Uvec1L(i, j, k)/Uvec1R(i, j, k))  < 0.9999).or.  &
                  (abs(Uvec1L(i, j, k)/Uvec1R(i, j, k))  > 1.0001))then
                write(6, *)i, j, k, 'Uvec1', Uvec1L(i, j, k), Uvec1R(i, j, k)
                
             end if
          end do
       end do
    end do



end subroutine get_u
!=====================================================================
! $B3F3J;R$K$*$1$k?eJ?<ANL%U%i%C%/%9J,I[(B($BBf7AB'$rMQ$$$?J}K!(B)
!=====================================================================


subroutine get_Fgrd(Hsfc1, Hbnd, Hthk1, Usfc1, Uvec1, Ufl, Uflgrd)

  implicit none

  real(8), INTENT(IN)  :: Hsfc1(1:IMAX, 1:JMAX)            !m
  real(8), INTENT(IN)  :: Hbnd( 1:IMAX, 1:JMAX)            !m
  real(8), INTENT(IN)  :: Hthk1( 1:IMAX, 1:JMAX)            !m
  real(8), INTENT(IN)  :: Usfc1(0:IMAX, 1:JMAX)            !m/s
  real(8), INTENT(IN)  :: Uvec1(0:IMAX, 1:JMAX, 1:KMAX)    !m/s
  real(8), INTENT(OUT) :: Uflgrd(0:IMAX, 1:JMAX, 1:KMAX)   !m^2/s
  real(8), INTENT(OUT) :: Ufl(0:IMAX, 1:JMAX)              !m^2/s

  real(8)              :: Ubnd(0:IMAX, 1:JMAX)
  integer              :: i, j, k

  !$B%U%i%C%/%9=i4|2=(B
  Ufl = 0.0D0

  ! x $BJ}8~$N6-3&$NCM$O(B 0 $B$K8GDj(B
  Uflgrd(0, :, :)    = 0.0D0
  uflgrd(IMAX, :, :) = 0.0D0
  ! $BDlLLB.EY(B
  call get_Ubnd(Hsfc1, Hbnd, Ubnd)

  do j = 1, JMAX
     do i = 1, IMAX-1
        ! $B:G2<AX(B (U_h(0) + U_h(1))*wgt(1)*(h-b)_i+1/2
        Uflgrd(i, j, 1)    = (Ubnd(i, j)   &
             &                +(Uvec1(i, j, 1)*wgt(2)+Uvec1(i, j, 2)*wgt(1))&
             &                  /(wgt(2)+wgt(1)))                   & 
             &         * (Hthk1(i+1, j)+Hthk1(i, j))/2.0D0*wgt(1)/2.0D0
        ! 2 $BAXL\$+$i:G>eAX2<$^$G(B
        do k = 2, KMAX-1
           Uflgrd(i, j, k) =  &
                &      + ((Uvec1(i, j, k+1)*wgt(k)                      &
                &          + Uvec1(i, j, k)*wgt(k+1))/(wgt(k+1)+wgt(k)) &
                &        + (Uvec1(i, j, k)*wgt(k-1)                     &
                &          + Uvec1(i, j, k-1)*wgt(k))/(wgt(k)+wgt(k-1)))&
                &      * (Hthk1(i+1, j) + Hthk1(i, j))/2.0D0 * wgt(k)/2.0D0
        end do
        ! $B:G>eAX(B
        
        Uflgrd(i, j, KMAX) = & 
             &       ((Uvec1(i, j, KMAX)*wgt(KMAX-1)                   &
             &        + Uvec1(i, j, KMAX-1)*wgt(KMAX))/(wgt(KMAX)+wgt(KMAX-1))&
             &          + Usfc1(i, j))                                  &
             &      * (Hthk1(i+1, j) + Hthk1(i, j))/2.0D0 * wgt(k)/2.0D0
     end do
  end do


     

  do k =1, KMAX
!     write(6, *)k,  Uflgrd(1, 2, k), Uflgrd(2, 2, k)
     Ufl(:, :) = Ufl(:, :) + Uflgrd(:, :, k)
  end do

end subroutine get_Fgrd

!=====================================================================
! $B?eJ?B.EYJ,I[(B($B2r@OE*$K5a$a$kJ}K!(B)
! $BB.EY$,29EY$K0MB8$9$k>l9g$O;H$($J$$(B
!=====================================================================
subroutine get_U_ana(Hsfc1, Hbnd, Hthk1, Usfc1, Uvec1)
 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(OUT) :: Usfc1(0:IMAX, 1:JMAX)
  real(8), INTENT(OUT) :: Uvec1(0:IMAX, 1:JMAX, 1:KMAX)
  real(8)              :: Ubnd (0:IMAX, 1:JMAX)
 

  integer              :: i, j,k
  ! x $BJ}8~$N6-3&$NCM$O(B 0 $B$K8GDj(B
  Uvec1(0, :, :)    = 0.0D0
  Uvec1(IMAX, :, :) = 0.0D0
 

  call get_Ubnd(Hsfc1, Hbnd, Ubnd)
  
  do j = 1, JMAX
     do i =1, IMAX-1

        ! $BI=LLB.EY(B
        Usfc1(i, j) = - 2.0D0 * (rho *grav)** Npw                           &
                &     * ((Hsfc1(i+1, j) - Hsfc1(i, j))/dx)**Npw             &
                &     * ((Hthk1(i+1, j) + Hthk1(i, j))/2.0D0)**(Npw+1.0D0)  &
                &     *  A_con /(Npw + 1.0D0)                

        ! $B3FAX$N7W;;(B
        do k = 1, KMAX
           Uvec1(i, j, k) = Usfc1(i, j)                         &
        &             * (1 - (1 - zk(k))**(Npw+1.0D0)) + Ubnd(i, j)
        end do

        
     end do
  end do

end subroutine get_U_ana

!=====================================================================
! $B?eJ?<ANL%U%i%C%/%9J,I[(B($B2r@OE*$K5a$a$kJ}K!(B)
!    $B$?$@$7(B, $B29EY0MB8%Q%i%a!<%?$ODj?t$H$7$F7W;;(B
!=====================================================================
subroutine get_F_ana(Hsfc1, Hbnd, Hthk1, Ufl)

  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(OUT) :: Ufl(  0:IMAX, 1:JMAX)

  real(8)              :: Ubnd( 0:IMAX, 1:JMAX)

  integer              :: i, j


  ! x $BJ}8~$N6-3&$NCM$O(B 0 $B$K8GDj(B
  Ufl(0, :)    = 0.0D0
  Ufl(IMAX, :) = 0.0D0

  call get_Ubnd(Hsfc1, Hbnd, Ubnd)
  
  do j = 1, JMAX
     do i = 1, IMAX-1

        ! $BDlLL$9$Y$j$K$h$k4sM?(B($B$9$Y$j$,$J$1$l$P(B Ufl = 0 )
        Ufl(i, j) =   (Hthk1(i+1, j) + Hthk1(i, j))/2.0D0 * Ubnd(i, j)

        Ufl(i, j) =   Ufl(i, j) - 2.0D0 * (rho *grav)** Npw     &
        &             * ((Hsfc1(i+1, j) - Hsfc1(i, j))/dx)**Npw    &
        &             * ((Hthk1(i+1, j) + Hthk1(i, j))/2.0D0)**(Npw+2.0D0) &
        &             * A_con /(Npw + 2.0D0)

        
     end do
  end do


end subroutine get_F_ana


!=====================================================================
! $B1tD>B.EYJ,I[(B
! $B%3%a%s%H(B: 
!         $B:BI8JQ49$K$h$C$F(B
!            w(x, z)  = - \int^{z}_{b} du/dx dz 
!            w(x, $B&F(B) = - h \int^{$B&F(B}_{0} du/dx d$B&F(B
!                        + dh/dx (- \int^{$B&F(B}_{0} u d$B&F(B + $B&F(B u )
!                        + db/dx (u($B&F(B) - u(0))
!         $B$H$J$k(B.
! half level $B$NCM$OMQ$$$J$$(B
!=====================================================================
subroutine get_W(Hsfc1, Hbnd, Hthk1, Uvec1, Usfc1,  Wvec1)
  
  implicit none
  
  real(8), parameter   :: one=1.0D0

  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), INTENT(OUT) :: Wvec1(1:IMAX, 1:JMAX, 0:KMAX)

  real(8)              ::  Ubnd(0:IMAX, 1:JMAX)
  real(8)              ::  Wbnd(1:IMAX, 1:JMAX)


! W $B$rF3=P$9$k$?$a$N9`(B(w $B$N<0$,D9$$$N$G(B)
  real(8)              :: dU(        1:IMAX, 1:JMAX, 0:KMAX)
  real(8)              :: Uvec1_w(   1:IMAX, 1:JMAX, 0:KMAX)
  real(8)              :: Wvec1_dudx(1:IMAX, 1:JMAX, 0:KMAX) !
  real(8)              :: Wvec1_dbdx(1:IMAX, 1:JMAX, 0:KMAX) !
  real(8)              :: Wvec1_intU(1:IMAX, 1:JMAX, 0:KMAX) !
  real(8)              :: Wvec1_dHdx(1:IMAX, 1:JMAX, 0:KMAX) !



  integer              :: i, j, k

  ! $BDlLLB.EY(B
  call get_Ubnd(Hsfc1, Hbnd, Ubnd)
  call get_Wbnd(Hbnd,  Ubnd, Wbnd)
  Wvec1(:, :, 0) = Wbnd(:, :)


  ! \int^0_0 f($B&F(B)d$B&F(B
  Wvec1_dudx(:, :, 0) = 0.0D0
  Wvec1_intU(:, :, 0) = 0.0D0


  ! $B=i4|2=$7$F$*$/(B
  Uvec1_w = 0.0D0
  dU      = 0.0D0

  do k = 1, KMAX
     Wvec1(:, :, k) = Wvec1(:, :, 0)
  end do

  do j = 1, JMAX
     do i = 2, IMAX-1

        dU(i, j, 0     ) = (Ubnd( i, j        ) - Ubnd( i-1, j        ))/dx
        dU(i, j, 1:KMAX-1) = (Uvec1(i, j, 1:KMAX-1) - Uvec1(i-1, j, 1:KMAX-1))/dx
        ! $BI=LL$NCM$O(B Usfc1 $B$+$i(B
        dU(i, j, KMAX:KMAX)   = (Usfc1(i, j        ) - Usfc1(i-1, j))/dx
        Uvec1_w(i, j, KMAX:KMAX) = (Usfc1(i, j) + Usfc1(i-1, j))/2.0D0

        do k =1, KMAX-1
           ! W $B$N3J;R>e$N(B dU/dx 
           dU(i, j, k) = (dU(i, j, k)*wgt(k+1) + dU(i, j, k+1)*wgt(k))&
                &            /(wgt(k) +wgt(k+1))
           ! W $B$N3J;R>e$N(B U $B$NCM(B
           Uvec1_w(i, j, k) = &
                &    ((Uvec1(i-1, j, k) + Uvec1(i, j, k))/2.0D0 * wgt(k+1)  &
                &   + (Uvec1(i-1, j, k+1) + Uvec1(i, j, k+1))/2.0D0 *wgt(k))&
                &        /(wgt(k)+wgt(k+1))
        end do

        do k = 1, KMAX
           
           ! - H \int^{$B&F(B}_{0} du/dx d$B&F(B
           Wvec1_dudx(i, j, k) = Wvec1_dudx(i, j, k-1)  &
                & - Hthk1(i, j)* (dU(i, j, k-1) + dU(i, j, k)) * wgt(k)/2.0D0
           !  db/dx (u($B&F(B) - u(0))
           Wvec1_dbdx(i, j, k) = (Hbnd(i+1, j) - Hbnd(i-1,j))/(2.0D0*dx) &
                & * (Uvec1_w(i, j, k)  - (Ubnd(i-1, j)+Ubnd(i, j))/2.0D0)
           ! \int^$B&F(B_0 u d$B&F(B
           Wvec1_intU(i, j, k) = Wvec1_intU(i, j, k-1) &
                & + (Uvec1_w(i, j, k-1) + Uvec1_w(i, j, k))*wgt(k)/2.0D0
           ! dH/dx($B&F(B u($B&F(B) - \int^$B&F(B_0 u d$B&F(B)
           Wvec1_dHdx(i, j, k) = (Hthk1(i+1,j)  - Hthk1(i-1, j))/(2.0D0*dx) &
                &  * (zk_half(k) * Uvec1_w(i, j, k) - Wvec1_intU(i, j, k))
        end do
     end do
  end do
!  write(6, *)'Wvec1_dudx', Wvec1_dudx(16, 2, :)*s_a
!  write(6, *)'Wvec1_dbdx', Wvec1_dbdx(16, 2, :)*s_a
!  write(6, *)'Wvec1_dHdx', Wvec1_dHdx(16, 2, :)*s_a




  Wvec1 =  Wvec1      &
  &      + Wvec1_dudx &
  &      + Wvec1_dbdx &
  &      + Wvec1_dHdx

  call monit_add_3('Wvec1     ', Wvec1(:, :, 1:KMAX), one)

end subroutine get_W

!$BDlLLB.EYJ,I[(B

subroutine get_Ubnd(Hsfc1, Hbnd, Ubnd)

  implicit none

  real(8), INTENT(IN)  :: Hsfc1(1:IMAX, 1:JMAX)
  real(8), INTENT(IN)  :: Hbnd (1:IMAX, 1:JMAX)
  real(8), INTENT(OUT) :: Ubnd (0:IMAX, 1:JMAX)
  
  Ubnd(1:IMAX, :) = Hsfc1-Hbnd  ! $B%@%_!<(B
  Ubnd = 0.0D0/s_a
  
  
end subroutine get_Ubnd


subroutine get_Wbnd(Hsfc1, Hbnd, Wbnd)

  implicit none

  real(8), INTENT(IN)  :: Hsfc1(1:IMAX, 1:JMAX)
  real(8), INTENT(IN)  :: Hbnd (1:IMAX, 1:JMAX)
  real(8), INTENT(OUT) :: Wbnd (1:IMAX, 1:JMAX)
  
  Wbnd = Hsfc1-Hbnd  ! $B%@%_!<(B
  Wbnd = 0.0D0/s_a
  
  
end subroutine get_Wbnd





end module u_pack

