!=======================================
! 2D cumulus model - kaminari
!  - subroutine vel_u_tb
! 
! Author    : TAKAHASHI Koko
! Date      : 2003/11/21 $B:G=*99?7(B
!             2003/11/11 $B?75,:n@.(B
! Note      : x $BJ}8~1?F0J}Dx<0$ND9$$;~4V(B
!             $B4V3V$G7W;;$9$kItJ,(B
! 
!=======================================
subroutine vel_u_tb(i,k,im,km,nuh,nuv,cm,ml,dx,dz,          &
     &              g_sqrt_x,g_sqrt_xz,g13,g13z,g13xz,u,omg,&
     &              u_old,w_old,km_sub,u_adv,u_cori,u_trb)

  implicit none

  integer(8), intent(in) :: i, k
  integer(8), intent(in) :: im, km

  real(8), intent(in) :: nuh,nuv,cm,ml
  real(8), intent(in) :: dx, dz

  real(8), intent(in) :: g_sqrt_x(-2:im+1,-2:km+1)
  real(8), intent(in) :: g_sqrt_xz(-2:im+1,-2:km+1)
  real(8), intent(in) :: g13(-2:im+1,-2:km+1)
  real(8), intent(in) :: g13z(-2:im+1,-2:km+1)
  real(8), intent(in) :: g13xz(-2:im+1,-2:km+1)

  real(8), intent(in) :: u(-2:im+1,-2:km+1)
  real(8), intent(in) :: omg(-2:im+1,-2:km+1)
  real(8), intent(in) :: u_old(-2:im+1,-2:km+1)
  real(8), intent(in) :: w_old(-2:im+1,-2:km+1)
  real(8), intent(in) :: km_sub(-2:im+1,-2:km+1)

  real(8), intent(out) :: u_adv(-2:im+1,-2:km+1)
  real(8), intent(out) :: u_cori(-2:im+1,-2:km+1)
  real(8), intent(out) :: u_trb(-2:im+1,-2:km+1)

  real(8)              :: u_dif(-2:im+1,-2:km+1)

  !--- $B2;GH8:?j9`(B
  u_dif(i,k) = nuh*(u(i+1,k) - 2*u(i,k) + u(i-1,k))         &
       &       /dx**2.0D0                                   &
       &       + nuv*(u(i,k+1) - 2*u(i,k) + u(i,k-1))       &
       &         /dz**2.0D0

  !--- $B0\N.9`(B
  u_adv(i,k) = u(i,k)*(u(i+1,k) - u(i-1,k))/(4.0D0*dx)      &
       &      + (                                           &
       &         omg(i+1,k+1) + omg(i+1,k-1)                &
       &         + omg(i-1,k+1) + omg(i-1,k-1)              &
       &        )/4.0D0                                     &
       &        *(u(i,k+1) - u(i,k-1))/(2.0D0*dz)           &
       &      - u_dif(i,k)

  !--- $B%3%j%*%j9`(B
  u_cori(i,k) = 0.0D0

  !--- $B3H;69`(B
  !    $BM?$($kCM$O(B ts - dts $B$N$H$-$NCM(B  ***_old(i,k)
  u_trb(i,k) = 2.0D0                                        &
       &       *(                                           &
       &         km_sub(i,k)                                &
       &         *(                                         &
       &           (u_old(i+1,k) - u_old(i,k))/dx           &
       &           + g13xz(i,k)                             &
       &             *(                                     &
       &               u_old(i+1,k+1) + u_old(i-1,k+1)      &
       &               + u_old(i+1,k-1) + u_old(i-1,k-1)    &
       &              )/(4.0d0*dz)                          &
       &          )                                         &
       &         - (                                        &
       &            km_sub(i-1,k)                           &
       &            *(                                      &
       &              (u_old(i,k) - u_old(i-1,k))/dx        &
       &               + g13xz(i-1,k)                       &
       &             *(                                     &
       &               u_old(i,k+1) + u_old(i-2,k+1)        &
       &               + u_old(i,k-1) + u_old(i-2,k-1)      &
       &              )/(4.0d0*dz)                          &
       &             )                                      &
       &            )                                       &
       &         )/dx                                       &
       &       + 2.0d0*g13z(i,k)                            &
       &         *(                                         &
       &           (                                        &
       &             km_sub(i+1,k+1) + km_sub(i,k+1)        &
       &             + km_sub(i+1,k-1) + km_sub(i,k-1)      &
       &           )/4.0d0                                  &
       &           *(                                       &
       &             (                                      &
       &              u_old(i+1,k+1) + u_old(i+1,k-1)       &
       &              - u_old(i-1,k+1) - u_old(i-1,k-1)     &
       &             )/(4.0d0*dx)                           &
       &             + g13(i,k+1)                           &
       &               *(u_old(i,k+1) - u_old(i,k))/dz      &
       &            )                                       &
       &           - (                                      &
       &               km_sub(i+1,k) + km_sub(i,k)          &
       &               + km_sub(i+1,k-2) + km_sub(i,k-2)    &
       &             )/4.0d0                                &
       &             *(                                     &
       &               (                                    &
       &                u_old(i+1,k) + u_old(i+1,k-2)       &
       &                - u_old(i-1,k) - u_old(i-1,k-2)     &
       &               )/(4.0d0*dx)                         &
       &               + g13(i,k)                           &
       &                 *(u_old(i,k) - u_old(i,k-1))/dz    &
       &              )                                     &
       &          )/dz                                      &
       &       + 1.0d0/g_sqrt_x(i,k)                        &
       &         *(                                         &
       &           (                                        &
       &             km_sub(i+1,k+1) + km_sub(i-1,k+1)      &
       &             + km_sub(i+1,k) + km_sub(i-1,k)        &
       &           )/4.0d0                                  &
       &           *(                                       &
       &             (w_old(i+1,k+1) -w_old(i,k+1))/dx      &
       &             + 1.0d0/g_sqrt_xz(i,k+1)               &
       &               *(u_old(i,k+1) - u_old(i,k))/dz      &
       &             + g13(i,k+1)                           &
       &               *(                                   &
       &                 w_old(i+1,k+1) + w_old(i-1,k+1)    &
       &                 - w_old(i+1,k-1) - w_old(i-1,k-1)  &
       &                )/(4.0d0*dz)                        &
       &            )                                       &
       &           - (                                      &
       &               km_sub(i+1,k) + km_sub(i-1,k)        &
       &               + km_sub(i+1,k-1) + km_sub(i-1,k-1)  &
       &             )/4.0d0                                &
       &             *(                                     &
       &               (w_old(i+1,k) -w_old(i,k))/dx        &
       &               + 1.0d0/g_sqrt_xz(i,k)               &
       &                 *(u_old(i,k) - u_old(i,k-1))/dz    &
       &               + g13(i,k)                           &
       &                 *(                                 &
       &                   w_old(i+1,k) + w_old(i-1,k)      &
       &                   - w_old(i+1,k-2) + w_old(i-1,k-2)&
       &                  )/(4.0d0*dz)                      &
       &              )                                     &
       &          )/dz                                      &
       &       - 2.0d0/(3.0d0*(cm*ml)**2.0d0)               &
       &         *(                                         &
       &           (km_sub(i,k)**2.0d0 - km_sub(i,k-1)**2.0d0)&
       &           /dx                                      &
       &           + g13z(i,k)                              &
       &             *(                                     &
       &                km_sub(i+1,k+1) + km_sub(i-1,k+1)   &
       &                - km_sub(i+1,k-1) - km_sub(i-1,k-1) &
       &              )/(4.0d0*dz)                          &
       &          )

end subroutine vel_u_tb
