!=======================================
! 2D cumulus model - kaminari
!  - subroutine eddymx_sub_1stp
! 
! Author    : TAKAHASHI Koko
! Date      : 2003/11/25  $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_1stp(i,k,im,km,nuh,nuv,eps,lv,rd,cp,   &
     &                     cm,ml,grv,dtb,dx,dz,              &
     &                     g_sqrt,g_sqrt_z,g13x,g13z,g13xz,  &
     &                     ptemp_bs,u_old,w_old,omg_old,pi,ptemp,qv,qc,  &
     &                     km_sub)

  implicit none

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

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

  real(8), intent(in) :: g_sqrt(-2:im+1,-2:km+1)
  real(8), intent(in) :: g_sqrt_z(-2:im+1,-2:km+1)
  real(8), intent(in) :: g13x(-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) :: ptemp_bs(-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) :: omg_old(-2:im+1,-2:km+1)
  real(8), intent(in) :: pi(-2:im+1,-2:km+1)
  real(8), intent(in) :: ptemp(-2:im+1,-2:km+1)
  real(8), intent(in) :: qv(-2:im+1,-2:km+1)
  real(8), intent(in) :: qc(-2:im+1,-2:km+1)

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

  real(8)             :: e_sub(-2:im+1,-2:km+1)
  real(8)             :: e_sub_dif(-2:im+1,-2:km+1)
  real(8)             :: e_sub_adv(-2:im+1,-2:km+1)
  real(8)             :: temp(-2:im+1,-2:km+1)
  real(8)             :: a(-2:im+1,-2:km+1)
  real(8)             :: ptemp_e(-2:im+1,-2:km+1)
  real(8)             :: ql(-2:im+1,-2:km+1)

  !--- $B%5%V%0%j%C%I1?F0%(%M%k%.!<(B
  e_sub(i,k) = (km_sub(i,k)/(cm*ml))**2.0d0

  !--- $B%5%V%0%j%C%I1?F0%(%M%k%.!<$N?tCMG4@-9`(B
  e_sub_dif(i,k) = nuh                                      &
       &           *(                                       &
       &             e_sub(i+1,k) - 2.0d0*e_sub(i,k)        &
       &             + e_sub(i-1,k)                         &
       &             )/dx**2.0d0                            &
       &           + nuv                                    &
       &           *(                                       &
       &             e_sub(i,k+1) - 2.0d0*e_sub(i,k)        &
       &              + e_sub(i,k-1)                        &
       &           )/dz**2.0d0

  !--- $B%5%V%0%j%C%I1?F0%(%M%k%.!<$N0\N.9`(B
  e_sub_adv(i,k) = (u(i+1,k) + u(i,k))/2.0d0                &
       &           *(                                       &
       &             e_sub(i+1,k+1) + e_sub(i+1,k-1)        &
       &             - e_sub(i-1,k+1) - e_sub(i-1,k-1)      &
       &            )/(4.0d0*dx)                            &
       &           + (omg(i,k+1) + omg(i,k))/2.0d0          &
       &             *(                                     &
       &               e_sub(i+1,k+1) + e_sub(i-1,k+1)      &
       &               - e_sub(i+1,k-1) - e_sub(i-1,k-1)    &
       &             )/(4.0d0*dz)                           &
       &           + e_sub_dif(i,k)

  !--- $B290L(B => $B29EY(B
  temp(i,k) = ptemp(i,k)*pi(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?e>x5$$H1@?e:.9gHf$NOB(B
  ql(i,k) = qv(i,k) + qc(i,k)

  !-- $B%5%V%0%j%C%I1?F0%(%M%k%.!<(B
  e_sub(i,k) = e_sub(i,k)                                   &
       &      + dtb                                         &
       &        *(- e_sub_adv(i,k)                          &
       &          + 3.0d0*grv*cm*ml*sqrt(e_sub(i,k))        &
       &            /g_sqrt(i,k)                            &
       &          *(- a(i,k)                                &
       &              *(                                    &
       &                ptemp_e(i+1,k+1) + ptemp_e(i-1,k+1) &
       &                - ptemp_e(i+1,k-1)                  &
       &                - ptemp_e(i-1,k-1)                  &
       &                )/(4.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)                          &
       &           )                                        &
       &          + 2.0d0*cm*ml*sqrt(e_sub(i,k))            &
       &            *(                                      &
       &              ((u(i+1,k) - u(i,k))/dx               &
       &               + g13xz(i,k)                         &
       &                 *(u(i+1,k+1) + u(i-1,k+1)          &
       &                   - u(i+1,k-1) - u(i-1,k-1)        &
       &                  )/(4.0d0*dz)                      &
       &              )**2.0d0                              &
       &              + (                                   &
       &                 1.0d0/g_sqrt(i,k)                  &
       &                 *(w(i,k+1) - w(i,k))/dz            &
       &                 )**2.0d0                           &
       &              )                                     &
       &          + cm*ml*sqrt(e_sub(i,k))                  &
       &            *(                                      &
       &              (                                     &
       &               w(i+1,k+1) + w(i+1,k-1)              &
       &               - w(i-1,k+1) - w(i-1,k-1)            &
       &              )/(4.0d0*dx)                          &
       &              + g13xz(i,k)                          &
       &                *(w(i,k+1) - w(i,k))/dz             &
       &               + 1.0d0/g_sqrt(i,k)                  &
       &                 *(                                 &
       &                   u(i+1,k+1) + u(i-1,k+1)          &
       &                   - u(i+1,k-1) - u(i-1,k-1)        &
       &                  )/(4.0d0*dz)                      &
       &             )                                      &
       &          + cm*ml*sqrt(e_sub(i,k))                  &
       &            *(                                      &
       &              (                                     &
       &               e_sub(i+1,k) - 2.0d0*e_sub(i,k)      &
       &               + e_sub(i-1,k)                       &
       &              )/dx**2.0d0                           &
       &              + (                                   &
       &                 g13z(i+1,k)                        &            
       &                 *(                                 &
       &                   e_sub(i+1,k+1) + e_sub(i,k+1)    &
       &                   - e_sub(i+1,k-1) - e_sub(i,k-1)  &
       &                  )/(4.0d0*dz)                      &
       &                 - g13z(i,k)                        &
       &                   *(                               &
       &                     e_sub(i,k+1) + e_sub(i-1,k+1)  &
       &                     - e_sub(i,k-1) - e_sub(i-1,k-1)&
       &                    )/(4.0d0*dz)                    &
       &                )/dx                                &
       &              + (                                   &
       &                 g13x(i,k+1)                        &
       &                 *(                                 &
       &                   e_sub(i+1,k+1) + e_sub(i+1,k)    &
       &                   - e_sub(i-1,k+1) - e_sub(i-1,k)  &
       &                  )/(4.0d0*dx)                      &
       &                 - g13x(i,k)                        &
       &                   *(                               &
       &                     e_sub(i+1,k) + e_sub(i+1,k-1)  &
       &                     - e_sub(i-1,k) - e_sub(i-1,k-1)&
       &                    )/(4.0d0*dx)                    &
       &                )/dz                                &
       &              + g13xz(i,k)                          &
       &                *(                                  &
       &                  g13x(i+1,k)                       &
       &                  *(e_sub(i,k+1) - e_sub(i,k))/dz   &
       &                  - g13x(i,k)                       &
       &                    *(e_sub(i,k) - e_sub(i,k-1))/dz &
       &                 )/dz                               &
       &              + 1.0d0/g_sqrt(i,k)                   &
       &                *(                                  &
       &                  1.0d0/g_sqrt_z(i,k+1)             &
       &                  *(e_sub(i,k+1) - e_sub(i,k))/dz   &
       &                  - 1.0d0/g_sqrt_z(i,k)             &
       &                    *(e_sub(i,k) - e_sub(i,k-1))/dz &
       &                 )/dz                               &
       &             )                                      &
       &          + 2.0d0*cm*ml*sqrt(e_sub(i,k))            &
       &            *(                                      &
       &              (                                     &
       &               (                                    &
       &                sqrt(e_sub(i+1,k))                  &
       &                - sqrt(e_sub(i-1,k))                &
       &               )/(2.0d0*dx)                         &
       &               + g13xz(i,k)                         &
       &                 *(                                 &
       &                   sqrt(e_sub(i,k+1))               &
       &                   - sqrt(e_sub(i,k-1))             &
       &                   )/(2.0d0*dz)                     &
       &               )**2.0d0                             &
       &               + (                                  &
       &                  1.0d0/g_sqrt_z(i,k)               & 
       &                  *(                                &
       &                    sqrt(e_sub(i,k+1))              &
       &                    - sqrt(e_sub(i,k-1))            &
       &                   )/(2.0d0*dz)                     &
       &                 )**2.0d0                           &
       &              )                                     &
       &          - cm/ml*e_sub(i,k)**(3.0d0/2.0d0)         &
       &         )

  !--- $B12:.9g78?t(B
  km_sub(i,k) = cm*ml*sqrt(e_sub(i,k))

end subroutine eddymx_sub_1stp

