!=======================================
! 2D cumulus model - kaminari
!  - subroutine eddymx_sub
! 
! Author    : TAKAHASHI Koko
! Date      : 2004/01/07  $B:G=*99?7(B
!             2003/11/05  $B?75,:n@.(B
! Note      : 1$B%9%F%C%WL\$N%5%V%0%j%C%I(B
!             $B1?F0%(%M%k%.!<$r5a$a$F$+$i(B
!             $B12G4@-78?t$r5a$a$k(B
!             $B=PNO$O12G4@-78?t(B
! 
!=======================================
!--- 1,2,3,4$B9TL\F~NOJQ?t(B, 5$B9TL\F~=PNOJQ?t(B
subroutine eddymx_sub(im,km,dx,dz,nuh,nuv,cm,ml,grv,dtb,    &
     &                g_sqrt,g_sqrt_z,g13x,g13z,g13xz,      &
     &                ptemp_bs,u_old,w_old,omg_old,         &
     &                ptemp,km_sub_old,km_sub_old2,         &
     &                km_sub)

  implicit none

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

  real(8), intent(in) :: dx, dz
  real(8), intent(in) :: nuh,nuv,cm,ml !, rd,cp,lv,eps
  real(8), intent(in) :: grv, dtb

  real(8), intent(in) :: g_sqrt(-2:im+2,-2:km+2)
  real(8), intent(in) :: g_sqrt_z(-2:im+2,-2:km+2)
  real(8), intent(in) :: g13x(-2:im+2,-2:km+2)
  real(8), intent(in) :: g13z(-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) :: 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) :: pi(-2:im+2,-2:km+2)
!  real(8), intent(in) :: temp(-2:im+2,-2:km+2)
  real(8), intent(in) :: ptemp(-2:im+2,-2:km+2)
  real(8), intent(in) :: km_sub_old(-2:im+2,-2:km+2)
  real(8), intent(in) :: km_sub_old2(-2:im+2,-2:km+2)

  real(8), intent(inout) :: km_sub(-2:im+2,-2:km+2)

  integer(8)          :: i, k

  real(8)             :: km_sub_dif(-2:im+2,-2:km+2)
  real(8)             :: km_sub_adv(-2:im+2,-2:km+2)
  real(8)             :: a(-2:im+2,-2:km+2)
  real(8)             :: ptemp_e(-2:im+2,-2:km+2)

  do k = 0,km
     do i = 0,im
        !--- $B%5%V%0%j%C%I1?F0%(%M%k%.!<$N?tCMG4@-9`(B
        km_sub_dif(i,k) = nuh                               &
             &           *(                                 &
             &             km_sub_old(i,k)                  &
             &             - 2.0d0*km_sub_old(i-1,k)        &
             &             + km_sub_old(i-2,k)              &
             &             )/(dx**2.0d0)                    &
             &           + nuv                              &
             &           *(                                 &
             &             km_sub_old(i,k)                  &
             &             - 2.0d0*km_sub_old(i,k-1)        &
             &             + km_sub_old(i,k-2)              &
             &           )/(dz**2.0d0)
      
        !--- $B%5%V%0%j%C%I1?F0%(%M%k%.!<$N0\N.9`(B
        km_sub_adv(i,k) = (u_old(i+1,k) + u_old(i,k))/2.0d0 &
             &           *(                                 &
             &             km_sub_old(i+1,k)                &
             &             - km_sub_old(i-1,k)              &
             &            )/(2.0d0*dx)                      &
             &           + (omg_old(i,k+1) + omg_old(i,k))  &
             &             /2.0d0                           &
             &             *(                               &
             &               km_sub_old(i,k+1)              &
             &               - km_sub_old(i,k-1)            &
             &              )/(2.0d0*dz)                    &
             &           + km_sub_dif(i,k)

        !--- $B78?t(B A
        a(i,k) = 1.0d0/ptemp_bs(i,k)  !                      &
!             &   *(1.0d0 + 1.61d0*eps*lv*qv(i,k)            &
!             &   /(rd*temp(i,k)))                           &
!             &   /(1.0d0                                    &
!             &     + (                                      &
!             &        eps*(lv**2.0d0)*qv(i,k)               &
!             &        /(cp*rd*temp(i,k)**2.0d0)             &
!             &       )                                      &
!             &     )
      
        !--- $BAjEv290L(B
        ptemp_e(i,k) = ptemp(i,k)   !                        &
!             &         *(1.0d0 + lv*qv(i,k)/(cp*temp(i,k))) &

        !-- $B%5%V%0%j%C%I1?F0%(%M%k%.!<(B
        km_sub(i,k) = km_sub_old2(i,k)                      &
             &      + 2.0d0*dtb                             &
             &        *(                                    &
             &          - km_sub_adv(i,k)                   &
             &          + 3.0d0*grv*(cm*ml)**2.0d0          &
             &            /(2.0d0*g_sqrt(i,k))              &
             &            *(- a(i,k)                        &
             &                *(                            &
             &                  ptemp_e(i,k+1)              &
             &                  - ptemp_e(i,k-1)            &
             &                 )/(2.0d0*dz)                 &
!             &            + (                               &
!             &               ql(i+1,k+1) + ql(i-1,k+1)      &
!             &               - ql(i+1,k-1) - ql(i-1,k-1)    &
!             &              )/(4.0d0*dz)                    &
             &             )                                &
             &          + (cm*ml)**2.0d0                    &
             &            *(                                &
             &              ((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)                &
             &              )**2.0d0                        &
             &              + (                             &
             &                 1.0d0/g_sqrt(i,k)            &
             &                 *(w_old(i,k+1) - w_old(i,k)) &
             &                 /dz                          &
             &                 )**2.0d0                     &
             &              )                               &
             &          + (cm*ml)**2.0d0/2.0d0              &
             &            *(                                &
             &              (                               &
             &               w_old(i+1,k+1) + w_old(i+1,k)  &
             &               - w_old(i-1,k+1) - w_old(i-1,k)&
             &              )/(4.0d0*dx)                    &
             &              + g13xz(i,k)                    &
             &                *(w_old(i,k+1) - w_old(i,k))  &
             &                /dz                           &
             &               + 1.0d0/g_sqrt(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)                &
             &             )**2.0d0                         &
             &          + 1.0d0/2.0d0                       &
             &            *(                                &
             &              (                               &
             &               km_sub_old(i+1,k)**2.0d0       &
             &               - 2.0d0*km_sub_old(i,k)**2.0d0 &
             &               + km_sub_old(i-1,k)**2.0d0     &
             &              )/(dx**2.0d0)                   &
             &              + (                             &
             &                 g13z(i+1,k)                  &            
             &                 *(                           &
             &                   km_sub_old(i+1,k+1)**2.0d0  &
             &                   + km_sub_old(i,k+1)**2.0d0  &
             &                   - km_sub_old(i+1,k-1)**2.0d0&
             &                   - km_sub_old(i,k-1)**2.0d0  &
             &                  )/(4.0d0*dz)                &
             &                 - g13z(i,k)                  &
             &                   *(                         &
             &                     km_sub_old(i,k+1)**2.0d0    &
             &                     + km_sub_old(i-1,k+1)**2.0d0&
             &                     - km_sub_old(i,k-1)**2.0d0  &
             &                     - km_sub_old(i-1,k-1)**2.0d0&
             &                    )/(4.0d0*dz)              &
             &                )/dx                          &
             &              + (                             &
             &                 g13x(i,k+1)                  &
             &                 *(                           &
             &                   km_sub_old(i+1,k+1)**2.0d0         &
             &                   + km_sub_old(i+1,k)**2.0d0         &
             &                   - km_sub_old(i-1,k+1)**2.0d0       &
             &                   - km_sub_old(i-1,k)**2.0d0         &
             &                  )/(4.0d0*dx)                &
             &                 - g13x(i,k)                  &
             &                   *(                         &
             &                     km_sub_old(i+1,k)**2.0d0         &
             &                     + km_sub_old(i+1,k-1)**2.0d0     &
             &                     - km_sub_old(i-1,k)**2.0d0       &
             &                     - km_sub_old(i-1,k-1)**2.0d0     &
             &                    )/(4.0d0*dx)              &
             &                )/dz                          &
             &              + g13xz(i,k)                    &
             &                *(                            &
             &                  g13x(i+1,k)                 &
             &                  *(                          &
             &                    km_sub_old(i,k+1)**2.0d0  &
             &                    - km_sub_old(i,k)**2.0d0  &
             &                   )/dz                       &
             &                  - g13x(i,k)                 &
             &                    *(                        &
             &                      km_sub_old(i,k)**2.0d0  &
             &                      - km_sub_old(i,k-1)**2.0d0&
             &                     )/dz                     &
             &                 )/dz                         &
             &              + 1.0d0/g_sqrt(i,k)             &
             &                *(                            &
             &                  1.0d0/g_sqrt_z(i,k+1)       &
             &                  *(                          &
             &                    km_sub_old(i,k+1)**2.0d0  &
             &                    - km_sub_old(i,k)**2.0d0  &
             &                   )/dz                       &
             &                  - 1.0d0/g_sqrt_z(i,k)       &
             &                    *(                        &
             &                      km_sub_old(i,k)**2.0d0  &
             &                      - km_sub_old(i,k-1)**2.0d0&
             &                     )/dz                     &
             &                 )/dz                         &
             &             )                                &
             &          + (                                 &
             &             (                                &
             &              km_sub_old(i+1,k)               &
             &              - km_sub_old(i-1,k)             &
             &             )/(2.0d0*dx)                     &
             &             + g13xz(i,k)                     &
             &               *(                             &
             &                 km_sub_old(i,k+1)            &
             &                 - km_sub_old(i,k-1)          &
             &                 )/(2.0d0*dz)                 &
             &            )**2.0d0                          &
             &          + (                                 &
             &             1.0d0/g_sqrt_z(i,k)              & 
             &             *(                               &
             &               km_sub_old(i,k+1)              &
             &               - km_sub_old(i,k-1)            &
             &              )/(2.0d0*dz)                    &
             &            )**2.0d0                          &
             &          - km_sub_old(i,k)**2.0d0            &
             &            /(2.0d0*ml**2.0d0)                &
             &         )
     end do
  end do

end subroutine eddymx_sub
