!=====================================================================
! $BB.EYJ,I[7W;;%5%V%k!<%A%s(B
! 2001/12/08 $B1tD>J,I[%k!<%A%s40N;(B
! 
!=====================================================================

module u_pack
  use prm
  implicit none

  private                       ! defult $B$r(B private $BB0@-(B

  public get_u                  ! $B?eJ?B.EY$NF3=P(B
  public get_w                  ! $B1tD>B.EY$NF3=P(B
  public get_F                  ! $B?eJ?B.EY%U%i%C%/%9$NF3=P(B



contains
!=====================================================================
! $B?eJ?B.EYJ,I[(B 
!=====================================================================
subroutine get_u(h, zk, u)

  real(8), intent(in)       :: h(0:im-1)          ! $BI98|JQ?t(B
 !h(im) = 0 $B$H$7$F9M$($?J}$,$h$$$+$b(B.
  real(8), intent(in)       :: zk(0:km)           ! $B&F:BI8(B
  real(8), intent(out)      :: u(0:im-1, 0:km)    ! $B?eJ?B.EY(B
  integer                   :: i, k               ! do $BJQ?t(B

  real(8)                   :: u_0(0:im-1)        ! $B?eJ?B.EY%9%1!<%k(B
  ! $BNN0h$N6-3&>r7o(B
  ! z = 0 $B$G$O>o$K(B u $B$O(B 0
  u(:, 0) = 0.
  
  ! x $B:BI8$N6-3&$G$N>r7o(B
  ! $B>o$K1&6-3&$G$O(B h = 0 $B$H$9$k(B. 

  u_0(im-1) = - (2 * (rho * g)**n)                             &
                & * (abs((- h(im-1))/(2 *dx)))**(n-1)          &
                & * (- h(im-1)) / (2 * dx)                     &
                & * (h(im-1)/2) ** (n+1) 
  do k =1, km
     u(im-1, k)  = u(im-1, k-1)                                 &
          &  + u_0(im-1) * ((zk(k) - zk(k-1))/2)                &
          &  * A * ((1 - zk(k))**n + (1 - zk(k-1))**n)
  end do

  !$BNN0hFb$G$N7W;;(B
  do i = 0, im-2
     do k = 1, km
        ! $B?eJ?B.EY%9%1!<%k$N7W;;(B
        u_0(i)= - (2 * (rho * g)**n)                     &
             & * (abs((h(i+1) - h(i))/dx))**(n-1)        &
             & * ((h(i+1) - h(i)) / dx)                  &
             & * ((h(i+1) + h(i))/2) ** (n+1) 
        ! $B?eJ?B.EYJ,I[$N7W;;(B
        u(i, k)  = u(i, k-1)                                     &
             &   + u_0(i) * ((zk(k) - zk(k-1))/2)                &
             &   * A * ((1 - zk(k))**n + (1 - zk(k-1))**n)
     enddo
  enddo

end subroutine get_u

!=====================================================================
! $B?eJ?%U%i%C%/%9J,I[(B 
! F = \int^{h}_{0} u dz 
!=====================================================================
subroutine get_F(h, zk, u, df)
  
  real(8) , intent(in)      :: h(0:im-1)                ! $BI98|JQ?t(B
  real(8) , intent(in)      :: zk(0:km)                 ! $B&F:BI8(B
  real(8) , intent(in)      :: u(0:im-1, 0:km)          ! $B?eJ?B.EY(B
  real(8) , intent(out)     :: df(0:im-1)               ! $B?eJ?B.EY:9J,(B

  real(8)                   :: f(-1:im-1)               ! $B%U%i%C%/%9(B
  real(8)                   :: f_0(0:im-1, 0:km)        ! u $B$N&F@QJ,G[Ns(B
  integer                   :: i                        ! do $BJQ?t(B

  !$B6-3&>r7o(B( i = 1/2 $B$G$N(B df $B$r5a$a$k$?$a(B)
  f(-1) = - f(0)

  !$BG[NsMWAG$r07$&$?$a$K(B f_0 $B$H$$G[Ns$X$NCV$-49$((B($B%@%5%$(B)
  !$B3F!9G[NsMWAG$r07$o$J$$$G(B($BG[Ns(B)$B%3!<%I$r=q$/J}K!$,$$$$$N$@$,(B...
  f_0 = int_f(u, zk)   

  !$BNN0hFb$N7W;;(B
  do i = 0, im-1
     f(i) = (h(i+1) + h(i))/2 * f_0(i, km)              ! 
     df(i) = f(i) - f(i-1)                              ! $B%U%i%C%/%9$N:9J,(B
  end do
  
end subroutine get_F



!=====================================================================
! $B1tD>B.EYJ,I[(B 
! 
! $B:BI8JQ49$K$h$C$F(B
! w(x, z)  = - \int^{z}_{0} 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
!            + dh/dx $B&F(B u
! $B$H$J$k(B.
!=====================================================================
subroutine get_w(h, zk, u, w)

  real(8), intent(in)       :: h(0:im-1)               ! $BI98|JQ?t(B
  real(8), intent(in)       :: zk(0:km)                ! $B&F:BI8(B
  real(8), intent(in)       :: u(0:im-1, 0:km)         ! $B?eJ?B.EY(B
  real(8), intent(out)      :: w(0:im-1, 0:km)         ! $B1tD>B.EY(B
  integer                   :: i, k                    ! do $BJQ?t(B

  real(8)                   :: F(0:im-1, 0:km)         ! $B?eJ?B.EY$N@QJ,JQ?t(B
  real(8)                   :: F_du_dx(0:im-1, 0:km)   ! du_dx $B$N@QJ,JQ?t(B



  F       = int_f(u, zk)             ! $BG[NsMWAG$r07$$$?$$$N$G(B
  F_du_dx = int_f(du_dx(u), zk)      ! $BJLG[Ns(B(F, F_du_dx)$B$KCV$-49$($?(B.


  do i = 0, im-1
     w(i, 0) = 0.
     do k = 1, km
        w(i, k) = - (h(i+1) + h(i))/2 * f_du_dx(i, k)               &
             &   + (h(i+1) - h(i))/dx                               &
             &   * (zk(k) * u(i, k)  - f(i, k))
     end do
  end do

end subroutine get_w

!=====================================================================
! $B$"$kJQ?t$r&FJ}8~$N@QJ,4X?t(B
! int_f(i, k) = \int^{km}_{0} f d$B&F(B' $B4X?t(B
!=====================================================================
function int_f(f, zk)
  real(8)    , intent(in)  :: f(0:im-1, 0:km)     ! $BHo@QJ,4X?t(B
  real(8)    , intent(in)  :: zk(0:km)            ! $B&F:BI8(B
  real(8)                  :: int_f(0:im-1, 0:km) ! $B@QJ,4X?t(B

  integer                  :: i, k                ! do $BJQ?t(B


  ! -- $B6-3&>r7o(B --- 
  int_f(:, 0) = 0.

  do i = 0, im-1
     do k = 1, km
        int_f(i, k) = int_f(i, k-1)                               & 
             &    + (f(i, k) + f(i, k-1)) * (zk(k) - zk(k-1)) /2
     end do
  end do
end function int_f

!=====================================================================
!u $B$N(B x $BHyJ,7W;;(B 
!=====================================================================
function du_dx(u)

  real(8)    , intent(in)   :: u(0:im-1, 0:km)      ! $B?eJ?B.EY(B
  real(8)                   :: du_dx(0:im-1, 0:km)  ! $BHyJ,JQ?t(B
  integer                   :: i, k                 ! do $BJQ?t(B

  du_dx(:, 0) = 0.
  do k = 1, km
     do i = 0, im-1
        du_dx(i, k) = (u(i+1, k) - u(i, k))/dx
     end do
  end do

end function du_dx

end module u_pack
