!=======================================
! 2D cumulus model - kaminari
!  - subroutine div_diff
!
! Author    : TAKAHASHI Koko
! Date      : 2003/12/04  最終更新
!             2003/11/13  新規作成
! Note      : 音波減衰項 div
!             x方向運動方程式中の
!              - 音波減衰項の x 微分:     div_x_u
!              - 音波減衰項の z^{*} 微分: div_z_u
!             z^{*}方向運動方程式中の
!              - 音波減衰項の z^{*} 微分: div_z_w
!
!=======================================
!--- 1行目入力変数, 2行目出力変数
subroutine div_diff(im,km,dx,dz,u,w,g_sqrt,g13xz,           &
     &              div_x_u,div_z_u,div_z_w)

  implicit none

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

  real(8), intent(in) :: dx, dz

  real(8), intent(in) :: u(-2:im+2,-2:km+2)
  real(8), intent(in) :: w(-2:im+2,-2:km+3)

  real(8), intent(in) :: g_sqrt(-2:im+2,-2:km+2)
  real(8), intent(in) :: g13xz(-2:im+2,-2:km+2)

  real(8), intent(out) :: div_x_u(0:im+1,0:km+1)
  real(8), intent(out) :: div_z_u(0:im+1,0:km+1)
  real(8), intent(out) :: div_z_w(0:im+1,0:km+1)

  integer(8)          :: i, k

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

  do k = -1,km+2
     do i = -1,im+2
        !--- div: スカラー量と同じ位置
        div(i,k) = (u(i+1,k) - u(i,k))/dx                   &
             &     + g13xz(i,k)                             &
             &       *(                                     &
             &         u(i+1,k+1) + u(i,k+1)                &
             &         - u(i+1,k-1) - u(i,k-1)              &
             &         )/(4.0d0*dz)                         &
             &     + 1.0d0/g_sqrt(i,k)*(w(i,k+1) - w(i,k))/dz
     end do
  end do

  do k = 0,km+1
     do i = 0,im+1
        !--- x方向運動方程式のdivのx微分: u(i,k)と同じ位置
        div_x_u(i,k) = (div(i,k) - div(i-1,k))/dx

        !--- x方向運動方程式のdivのz微分: u(i,k)と同じ位置
        div_z_u(i,k) = (                                    &
             &          div(i,k+1) + div(i-1,k+1)           &
             &          - div(i,k-1) - div(i-1,k-1)         &
             &         )/(4.0d0*dz)

        !--- z方向運動方程式のdivのz微分: w(i,k)と同じ位置
        div_z_w(i,k) = (div(i,k) - div(i,k-1))/dz
     end do
  end do

end subroutine div_diff