!=======================================
! 2D cumulus model - kaminari
!  - subroutine div_diff
! 
! Author    : TAKAHASHI Koko
! Date      : 2003/11/13
! Note      : $B2;GH8:?j9`(B div
!             x$BJ}8~1?F0J}Dx<0Cf$N(B
!              - $B2;GH8:?j9`$N(B x $BHyJ,(B:     div_x_u
!              - $B2;GH8:?j9`$N(B z^{*} $BHyJ,(B: div_z_u
!             z^{*}$BJ}8~1?F0J}Dx<0Cf$N(B
!              - $B2;GH8:?j9`$N(B z^{*} $BHyJ,(B: div_z_w
! 
!=======================================
subroutine div_diff(i,k,im,km,dx,dz,u,w,u_old,w_old,g12,    &
     &              g13xz,div_x_u,div_z_u,div_z_w,div_z_w_old)

  implicit none

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

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

  real(8), intent(in) :: u(0:im+1,0:km+1)
  real(8), intent(in) :: w(0:im+1,0:km+1)
  real(8), intent(in) :: u_old(0:im+1,0:km+1)
  real(8), intent(in) :: w_old(0:im+1,0:km+1)
  real(8), intent(in) :: g12(0:im+1,0:km+1)
  real(8), intent(in) :: g13xz(0:im+1,0:km+1)

  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)
  real(8), intent(out) :: div_z_w_old(0:im+1,0:km+1)

  real(8)              :: div(0:im+1,0:km+1)
  real(8)              :: div_old(0:im+1,0:km+1)


  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/g12(i,k)*(w(i,k+1) - w(i,k))/dz

  div_x_u(i,k) = (div(i+1,k) - div(i,k))/dx

  div_z_u(i,k) = (div(i,k+1) - div(i,k-1))/(2.0D0*dz)

  div_z_w(i,k) = (div(i,k+1) - div(i,k))/dz

  div_old(i,k) = (u_old(i+1,k) - u_old(i,k))/dx             &
       &     + g13xz(i,k)                                   &
       &       *(                                           &
       &         u_old(i+1,k+1) + u_old(i,k+1)              &
       &         - u_old(i+1,k-1) - u_old(i,k+1)            &
       &         )/(4.0D0*dz)                               &
       &     + 1/g12(i,k)*(w_old(i,k+1) - w_old(i,k))/dz

  div_z_w_old(i,k) = (div_old(i,k+1) - div_old(i,k))/dz
  
end subroutine div_diff
