!=======================================
! 2D cumulus model - kaminari
!  - subroutine vel_w_tb
! 
! Author    : TAKAHASHI Koko
! Date      : 2004/01/06  $B:G=*99?7(B
!             2003/11/11  $B?75,:n@.(B
! Note      : z $BJ}8~1?F0J}Dx<0$ND9$$;~4V(B
!             $B4V3V$G7W;;$9$kItJ,(B
! 
!=======================================
!--- 1,2,3,4$B9TL\F~NOJQ?t(B, 5$B9TL\=PNOJQ?t(B
subroutine vel_w_tb(im,km,dx,dz,nuh,nuv,grv,cm,ml,ptemp_bs, &
     &              g_sqrt,g_sqrt_x,g_sqrt_z,g13,g13x,g13xz,&
     &              u_old,w_old,omg_old,u_old2,w_old2,      &
     &              ptemp,km_sub_old2,                      &
     &              w_adv,w_buoy,w_trb)

  implicit none

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

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

  real(8), intent(in) :: g_sqrt(-2:im+2,-2:km+2)
  real(8), intent(in) :: g_sqrt_x(-2:im+2,-2:km+2)
  real(8), intent(in) :: g_sqrt_z(-2:im+2,-2:km+2)
  real(8), intent(in) :: g13(-2:im+2,-2:km+2)
  real(8), intent(in) :: g13x(-2:im+2,-2:km+2)
  real(8), intent(in) :: g13xz(-2:im+2,-2:km+2)

  real(8), intent(in) :: ptemp_bs(-2:im+2,-2:km+2)
!  real(8), intent(in) :: qv_bs(-2:im+2,-2:km+2)

  real(8), intent(in) :: u_old(-2:im+2,-2:km+2)
  real(8), intent(in) :: w_old(-2:im+2,-2:km+3)
  real(8), intent(in) :: omg_old(-2:im+2,-2:km+3)
  real(8), intent(in) :: u_old2(-2:im+2,-2:km+2)
  real(8), intent(in) :: w_old2(-2:im+2,-2:km+3)
  real(8), intent(in) :: ptemp(-2:im+2,-2:km+2)
!  real(8), intent(in) :: qv(-2:im+2,-2:km+2)
!  real(8), intent(in) :: qc(-2:im+2,-2:km+2)
!  real(8), intent(in) :: qr(-2:im+2,-2:km+2)
  real(8), intent(in) :: km_sub_old2(-2:im+2,-2:km+2)

  real(8), intent(out) :: w_adv(-2:im+2,-2:km+3)
  real(8), intent(out) :: w_buoy(-2:im+2,-2:km+3)
  real(8), intent(out) :: w_trb(-2:im+2,-2:km+3)

  integer(8)          :: i, k

  real(8)              :: w_dif(-2:im+2,-2:km+3)

  do k = 1,km
     do i = 0,im
        !--- $B?tCMG4@-9`(B
        w_dif(i,k) = nuh*(                                  &
             &            w_old(i,k) - 2.0d0*w_old(i-1,k)   &
             &            + w_old(i-2,k)                    &
             &           )/(dx**2.0d0)                      &
             &       + nuv*(                                &
             &              w_old(i,k) - 2.0d0*w_old(i,k-1) &
             &              + w_old(i,k-2)                  &
             &             )/(dz**2.0d0)

        !--- $B0\N.9`(B
        w_adv(i,k) = (                                      &
             &        u_old(i,k) + u_old(i,k-1)             &
             &        + u_old(i+1,k) + u_old(i+1,k-1)       &
             &       )/4.0d0                                &
             &       *(w_old(i+1,k) - w_old(i-1,k))         &
             &       /(2.0d0*dx)                            &
             &       + omg_old(i,k)                         &
             &         *(w_old(i,k+1) - w_old(i,k-1))       &
             &         /(2.0d0*dz)                          &
             &       + w_dif(i,k)

        !--- $BIbNO9`(B
        w_buoy(i,k) = grv                                   &
             &        *(                                    &
             &          (ptemp(i,k) + ptemp(i,k-1))/2.0d0   &
             &          /(                                  &
             &            (ptemp_bs(i,k) + ptemp_bs(i,k-1)) &
             &            /2.0d0                            &
             &           )                                  &
             &          - 1.0d0                             &
!             &          + 0.61d0                            &
!             &            *(                                &
!             &              (qv(i,k+1) + qv(i,k))/2.0d0     &
!             &              - (qv_bs(i,k+1) - qv_bs(i,k))   &
!             &                /2.0d0                        &
!             &             )                                &
!             &          - (qc(i,k+1) + qc(i,k))/2.0d0       &
!             &          - (qr(i,k+1) + qr(i,k))/2.0d0       &
             &            )

        !--- $B3H;69`(B
        w_trb(i,k) = (                                      &
             &        (                                     &
             &         (                                    &
             &          km_sub_old2(i+1,k+1)                &
             &          + km_sub_old2(i,k+1)                &
             &          + km_sub_old2(i+1,k)                &
             &          + km_sub_old2(i,k)                  &
             &         )/4.0d0                              &
             &         *(                                   &
             &           (w_old2(i+1,k) - w_old2(i,k))/dx   &
             &           + 1.0d0/g_sqrt_x(i+1,k)            &
             &             *(                               &
             &               u_old2(i+1,k) + u_old2(i,k)    &
             &               - u_old2(i+1,k-1)              &
             &               - u_old2(i,k-1)                &
             &              )/(4.0d0*dz)                    &
             &           + g13(i+1,k)                       &
             &             *(                               &
             &               w_old2(i+1,k+1) + w_old2(i,k+1)&
             &               - w_old2(i+1,k-1)              &
             &               - w_old2(i,k-1)                &
             &              )/(4.0d0*dz)                    &
             &           )                                  &
             &        )                                     &
             &        - (                                   &
             &           (                                  &
             &            km_sub_old2(i,k+1)                &
             &            + km_sub_old2(i-1,k+1)            &
             &            + km_sub_old2(i,k)                &
             &            + km_sub_old2(i-1,k)              &
             &           )/4.0d0                            &
             &           *(                                 &
             &             (w_old2(i,k) - w_old2(i-1,k))/dx &
             &             + 1.0d0/g_sqrt_x(i,k)            &
             &               *(                             &
             &                 u_old2(i,k) + u_old2(i-1,k)  &
             &                 - u_old2(i,k-1)              &
             &                 - u_old2(i-1,k-1)            &
             &                )/(4.0d0*dz)                  &
             &             + g13(i,k)                       &
             &               *(                             &
             &                 w_old2(i,k) + w_old2(i-1,k)  &
             &                 - w_old2(i,k-1)              &
             &                 - w_old2(i-1,k-1)            &
             &                )/(4.0d0*dz)                  &
             &            )                                 &
             &           )                                  &
             &       )/dx                                   &
             &       + g13x(i,k)                            &
             &         *(                                   &
             &           km_sub_old2(i,k)                   &
             &           *(                                 &
             &             (                                &
             &              w_old2(i,k+1)                   &
             &              + w_old2(i,k)                   &
             &              - w_old2(i-1,k+1)               &
             &              - w_old2(i-1,k)                 &
             &             )/(4.0d0*dx)                     &
             &             + 1.0d0/g_sqrt(i,k)              &
             &               *(                             &
             &                 u_old2(i+1,k+1)              &
             &                 + u_old2(i,k+1)              &
             &                 - u_old2(i+1,k)              &
             &                 - u_old2(i,k)                &
             &                )/(4.0d0*dz)                  &
             &             + g13xz(i,k+1)                   &
             &               *(w_old2(i,k+1) - w_old2(i,k)) &
             &               /dz                            &
             &            )                                 &
             &           - km_sub_old2(i,k-1)               &
             &             *(                               &
             &               (                              &
             &                w_old2(i,k)                   &
             &                + w_old2(i,k-1)               &
             &                - w_old2(i-1,k)               &
             &                - w_old2(i-1,k-1)             &
             &               )/(4.0d0*dx)                   &
             &               + 1.0d0/g_sqrt(i,k)            &
             &                 *(                           &
             &                   u_old2(i+1,k)              &
             &                   + u_old2(i,k)              &
             &                   - u_old2(i+1,k-1)          &
             &                   - u_old2(i,k-1)            &
             &                  )/(4.0d0*dz)                &
             &               + g13xz(i,k)                   &
             &                 *(w_old2(i,k) - w_old2(i,k-1))&
             &                 /dz                          &
             &              )                               &
             &          )/dz                                &
             &       + 2.0d0/g_sqrt_z(i,k)                  &
             &         *(                                   &
             &           km_sub_old2(i,k)/g_sqrt(i,k)       &
             &           *(w_old2(i,k+1) -w_old2(i,k))/dz   &
             &           - km_sub_old2(i,k-1)/g_sqrt(i,k-1) &
             &             *(w_old2(i,k) -w_old2(i,k-1))/dz &
             &          )/dz                                &
             &       - 2.0d0/(3.0d0*(cm*ml)**2.0d0)         &
             &         *(                                   &
             &           (                                  &
             &            km_sub_old2(i,k)**2.0d0           &
             &            + km_sub_old2(i,k-1)**2.0d0       &
             &            - km_sub_old2(i-1,k)**2.0d0       &
             &            - km_sub_old2(i-1,k-1)**2.0d0     &
             &            )/(4.0d0*dx)                      &
             &            + g13x(i,k)                       &
             &              *(                              &
             &                km_sub_old2(i,k)**2.0d0       &
             &                - km_sub_old2(i,k-1)**20.d0   &
             &               )/dz                           &
             &            + 1.0d0/g_sqrt_z(i,k)             &
             &              *(                              &
             &                km_sub_old2(i,k)**2.0d0       &
             &                - km_sub_old2(i,k-1)**2.0d0   &
             &               )/dz                           &
             &          )
        w_trb(i,k) = 0.0d0
     end do
  end do

end subroutine vel_w_tb
