!=======================================
! 2D cumulus model - kaminari
!  - subroutine kener_sub_1stp
! 
! Author    : TAKAHASHI Koko
! Date      : 2003/11/05
! Note      : 1$B%9%F%C%WL\$N%5%V%0%j%C%I(B
!             $B1?F0%(%M%k%.!<(B
! 
!=======================================
subroutine kener_sub_1stp(i,k,im,km,nuh,nuv,eps,lv,rd,cp,cm,&
     &                    ml,grv,dtb,dx,dz,ptemp_bs,u,w,omg,&
     &                    temp,ptemp,qv,qc,g12,g12z,        &
     &                    g13x,g13z,g13xz,e_sub,km_sub)

  implicit none

  integer, intent(in) :: i, k
  integer, 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) :: u(0:im+1,0:km+1)
  real(8), intent(in) :: w(0:im+1,0:km+1)
  real(8), intent(in) :: omg(0:im+1,0:km+1)
  real(8), intent(in) :: ptemp_bs(0:im+1,0:km+1)
  real(8), intent(in) :: temp(0:im+1,0:km+1)
  real(8), intent(in) :: ptemp(0:im+1,0:km+1)
  real(8), intent(in) :: qv(0:im+1,0:km+1)
  real(8), intent(in) :: qc(0:im+1,0:km+1)
  real(8), intent(in) :: g12(0:im+1,0:km+1)
  real(8), intent(in) :: g12z(0:im+1,0:km+1)
  real(8), intent(in) :: g13x(0:im+1,0:km+1)
  real(8), intent(in) :: g13z(0:im+1,0:km+1)
  real(8), intent(in) :: g13xz(0:im+1,0:km+1)
!  real(8), intent(in) :: (0:im+1,0:km+1)
!  real(8), intent(in) :: 

  real(8), intent(inout) :: e_sub(0:im+1,0:km+1)
  real(8), intent(inout) :: km_sub(0:im+1,0:km+1)

  real(8)             :: e_sub_dif(0:im+1,0:km+1)
  real(8)             :: e_sub_adv(0:im+1,0:km+1)
  real(8)             :: a(0:im+1,0:km+1)
  real(8)             :: ptemp_e(0:im+1,0:km+1)
  real(8)             :: ql(0:im+1,0:km+1)
!  real(8)             :: (0:im+1,0:km+1)

  !--- $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*e_sub(i,k)            &
       &             + e_sub(i-1,k)                         &
       &             )                                      &
       &             /(4*dx**2)                             &
       &           + nuv                                    &
       &           *(                                       &
       &             e_sub(i,k+1) - 2*e_sub(i,k)            &
       &              + e_sub(i,k-1)                        &
       &           )/(4*dz**2)
  e_sub_dif(i,k) = 0.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                    &
       &           *(e_sub(i+1,k) - e_sub(i,k))/(4*dx)      &
       &           + (omg(i,k+1) + omg(i,k))/2              &
       &             *(e_sub(i,k+1) - e_sub(i,k))/(4*dz)    &
       &           + e_sub_dif(i,k)

  e_sub_adv(i,k) = 0.0D0

  !--- $B78?t(B A
  a(i,k) = 1/ptemp_bs(i,k)                                  &
       &   *(1 + 1.61*eps*lv*qv(i,k)/(rd*temp(i,k)))        &
       &   /(1 + eps*(lv**2)*qv(i,k)/(cp*rd*temp(i,k)**2))

  !--- $BAjEv290L(B
  ptemp_e(i,k) = ptemp(i,k)*(1 + 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*grv*cm*ml*e_sub(i,k)**(1/2)/g12(i,k)  &
       &          *(- a(i,k)*ptemp_e(i,k+1)                 &
       &            - ptemp_e(i,k-1)/(4*dz)                 &
       &            + (ql(i,k+1) - ql(i,k-1))/(4*dz)        &
       &           )                                        &
       &          + 2*cm*ml*e_sub(i,k)**(1/2)               &
       &            *(                                      &
       &              ((u(i+1,k) - u(i,k))/(2*dx)           &
       &               + g13xz(i,k)                         &
       &                 *(u(i+1,k+1) + u(i,k+1)            &
       &                   - u(i+1,k-1)- u(i,k-1)           &
       &                  )/(8*dz)                          &
       &               )**2                                 &
       &               + (                                  &
       &                  1/g12(i,k)                        &
       &                  *(w(i,k+1) - w(i,k))              &
       &                  /(2*dz)                           &
       &                 )**2                               &
       &               + (w(i+1,k+1) + w(i+1,k-1)           &
       &                  - w(i-1,k+1) - w(i-1,k-1)         &
       &                 )/(8*dx)                           &
       &               + g13xz(i,k)                         &
       &                 *(w(i,k+1) - w(i,k))               &
       &                 /(2*dz)                            &
       &               + 1/g12(i,k)                         &
       &                 *(                                 &
       &                   u(i+1,k+1) + u(i-1,k+1)          &
       &                   - u(i+1,k-1) - u(i-1,k-1)        &
       &                  )/(8*dz)                          &
       &             )                                      &
       &          + cm*ml                                   &
       &            /(2*e_sub(i,k)**(1/2))                  &
       &            *(                                      &
       &              (                                     &
       &               e_sub(i+1,k)**2                      &
       &               - e_sub(i,k)**2                      &
       &              )/dx**2                               &
       &              + (                                   &
       &                 g13z(i+1,k)                        &            
       &                 *(                                 &
       &                   e_sub(i+2,k+1)**2                &
       &                   + e_sub(i,k+1)**2                &
       &                   + e_sub(i+2,k-1)**2              &
       &                   + e_sub(i,k-1)**2                &
       &                  )/(8*dz)                          &
       &                 - g13z(i,k)                        &
       &                   *(                               &
       &                     e_sub(i+1,k+1)**2              &
       &                     + e_sub(i-1,k+1)**2            &
       &                     + e_sub(i+1,k-1)**2            &
       &                     + e_sub(i-1,k-1)**2            &
       &                    )/(8*dz)                        &
       &                )/(2*dx)                            &
       &              + (                                   &
       &                 g13x(i,k+1)                        &
       &                 *(                                 &
       &                   e_sub(i+1,k+2)**2                &
       &                   + e_sub(i+1,k)**2                &
       &                   + e_sub(i-1,k+2)**2              &
       &                   + e_sub(i-1,k)**2                &
       &                  )/(8*dx)                          &
       &                 - g13x(i,k)                        &
       &                   *(                               &
       &                     e_sub(i+1,k+1)**2              &
       &                     + e_sub(i+1,k-1)**2            &
       &                     + e_sub(i-1,k+1)**2            &
       &                     + e_sub(i-1,k-1)**2            &
       &                    )/(8*dx)                        &
       &                )/(2*dz)                            &
       &              + g13xz(i,k)                          &
       &                *(                                  &
       &                  g13x(i+1,k)                       &
       &                  *(                                &
       &                    e_sub(i,k+2)**2                 &
       &                    - e_sub(i,k+1)**2               &
       &                   )/(2*dz)                         &
       &                  - g13x(i,k)                       &
       &                    *(                              &
       &                      e_sub(i,k+1)**2               &
       &                      - e_sub(i,k)**2               &
       &                     )/(2*dz)                       &
       &                 )/(2*dz)                           &
       &              + 1/g12(i,k)                          &
       &                *(                                  &
       &                  g12z(i,k+1)                       &
       &                    *(                              &
       &                      e_sub(i,k+2)**2               &
       &                      - e_sub(i,k+1)**2                &
       &                     )/(2*dz)                         &
       &                  - g12z(i,k+1)         &            
       &                      *(                              &
       &                        e_sub(i,k+1)**2                &
       &                        - e_sub(i,k)**2                &
       &                       )/(2*dz)                       &
       &                 )                                    &
       &              - (                                     &
       &                 (                                    &
       &                  e_sub(i+1,k)**2                      &
       &                  - e_sub(i-1,k)                       &
       &                 )/(4*dx)                             &
       &                 + g13xz(i,k)           &            
       &                   *(                                 &
       &                     e_sub(i,k+1)**2                   &
       &                     - e_sub(i,k-1)                    &
       &                 )/(4*dz)                             &
       &                )**2                                  &
       &              - (                                     &
       &                 1/g12z(i,k)            &            
       &                 *(                                   &
       &                   e_sub(i,k+1)                        &
       &                   - e_sub(i,k-1)                      &
       &                  )/(2*dz)                            &
       &                 )**2                                 &
       &             )                                        &
       &          + cm/ml*e_sub(i,k)**(3/2)                               &
       &         )

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

end subroutine kener_sub_1stp

