!=======================================
! 2D cumulus model - kaminari
!  - subroutine eddymx_sub
!
! Author    : TAKAHASHI Koko
! Date      : 2003/11/25  最終更新
!             2003/11/05  新規作成
! Note      : 1ステップ目のサブグリッド
!             運動エネルギーを求めてから
!             渦粘性係数を求める
!             出力は渦粘性係数
!
!=======================================
!--- 1,2,3,4,5,6行目入力変数, 7行目入出力変数
subroutine eddymx_sub(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_old,  &
     &                ptemp_old,qv_old,qc_old,              &
     &                km_sub_old,km_sub_old2,               &
     &                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_old(-2:im+1,-2:km+1)
  real(8), intent(in) :: ptemp_old(-2:im+1,-2:km+1)
  real(8), intent(in) :: qv_old(-2:im+1,-2:km+1)
  real(8), intent(in) :: qc_old(-2:im+1,-2:km+1)
  real(8), intent(in) :: km_sub_old(-2:im+1,-2:km+1)
  real(8), intent(in) :: km_sub_old2(-2:im+1,-2:km+1)

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

  real(8)             :: temp_old(-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)

  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)             :: e_sub_old(-2:im+1,-2:km+1)
  real(8)             :: e_sub_old2(-2:im+1,-2:km+1)

  !--- サブグリッド運動エネルギー
  e_sub(i,k) = (km_sub(i,k)/(cm*ml))**2.0d0
  e_sub_old(i,k) = (km_sub_old(i,k)/(cm*ml))**2.0d0
  e_sub_old2(i,k) = (km_sub_old2(i,k)/(cm*ml))**2.0d0

  !--- サブグリッド運動エネルギーの数値粘性項
  e_sub_dif(i,k) = nuh                                      &
       &           *(                                       &
       &             e_sub_old(i+1,k) - 2.0d0*e_sub_old(i,k)&
       &             + e_sub_old(i-1,k)                     &
       &             )/dx**2.0d0                            &
       &           + nuv                                    &
       &           *(                                       &
       &             e_sub_old(i,k+1) - 2.0d0*e_sub_old(i,k)&
       &              + e_sub_old(i,k-1)                    &
       &           )/dz**2.0d0

  !--- サブグリッド運動エネルギーの移流項
  e_sub_adv(i,k) = (u_old(i+1,k) + u_old(i,k))/2.0d0        &
       &           *(                                       &
       &             e_sub_old(i+1,k+1) + e_sub_old(i+1,k-1)&
       &             - e_sub_old(i-1,k+1)                   &
       &             - e_sub_old(i-1,k-1)                   &
       &            )/(4.0d0*dx)                            &
       &           + (omg_old(i,k+1) + omg_old(i,k))/2.0d0  &
       &             *(                                     &
       &               e_sub_old(i+1,k+1)                   &
       &               + e_sub_old(i-1,k+1)                 &
       &               - e_sub_old(i+1,k-1)                 &
       &               - e_sub_old(i-1,k-1)                 &
       &             )/(4.0d0*dz)                           &
       &           + e_sub_dif(i,k)

  !--- 温位 => 温度
  temp_old(i,k) = ptemp_old(i,k)*pi_old(i,k)

  !--- 係数 A
  a(i,k) = 1.0d0/ptemp_bs(i,k)                              &
       &   *(                                               &
       &     1.0d0                                          &
       &     + 1.61d0*eps*lv*qv_old(i,k)/(rd*temp_old(i,k)) &
       &    )/(1.0d0                                        &
       &     + (                                            &
       &        eps*(lv**2.0d0)*qv_old(i,k)                 &
       &        /(cp*rd*temp_old(i,k)**2.0d0)               &
       &       )                                            &
       &     )

  !--- 相当温位
  ptemp_e(i,k) = ptemp_old(i,k)                             &
       &         *(1.0d0 + lv*qv_old(i,k)/(cp*temp_old(i,k)))

  !--- 水蒸気と雲水混合比の和
  ql(i,k) = qv_old(i,k) + qc_old(i,k)

  !-- サブグリッド運動エネルギー
  e_sub(i,k) = e_sub_old2(i,k)                              &
       &      + 2.0d0*dtb                                   &
       &        *(- e_sub_adv(i,k)                          &
       &          + 3.0d0*grv*cm*ml*sqrt(e_sub_old(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_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-1,k+1)  &
       &                   - u_old(i+1,k-1) - u_old(i-1,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*sqrt(e_sub_old(i,k))              &
       &            *(                                      &
       &              (                                     &
       &               w_old(i+1,k+1) + w_old(i+1,k-1)      &
       &               - w_old(i-1,k+1) - w_old(i-1,k-1)    &
       &              )/(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-1,k+1)  &
       &                   - u_old(i+1,k-1) - u_old(i-1,k-1)&
       &                  )/(4.0d0*dz)                      &
       &             )                                      &
       &          + cm*ml*sqrt(e_sub_old(i,k))              &
       &            *(                                      &
       &              (                                     &
       &               e_sub_old(i+1,k)                     &
       &               - 2.0d0*e_sub_old(i,k)               &
       &               + e_sub_old(i-1,k)                   &
       &              )/dx**2.0d0                           &
       &              + (                                   &
       &                 g13z(i+1,k)                        &
       &                 *(                                 &
       &                   e_sub_old(i+1,k+1)               &
       &                   + e_sub_old(i,k+1)               &
       &                   - e_sub_old(i+1,k-1)             &
       &                   - e_sub_old(i,k-1)               &
       &                  )/(4.0d0*dz)                      &
       &                 - g13z(i,k)                        &
       &                   *(                               &
       &                     e_sub_old(i,k+1)               &
       &                     + e_sub_old(i-1,k+1)           &
       &                     - e_sub_old(i,k-1)             &
       &                     - e_sub_old(i-1,k-1)           &
       &                    )/(4.0d0*dz)                    &
       &                )/dx                                &
       &              + (                                   &
       &                 g13x(i,k+1)                        &
       &                 *(                                 &
       &                   e_sub_old(i+1,k+1)               &
       &                   + e_sub_old(i+1,k)               &
       &                   - e_sub_old(i-1,k+1)             &
       &                   - e_sub_old(i-1,k)               &
       &                  )/(4.0d0*dx)                      &
       &                 - g13x(i,k)                        &
       &                   *(                               &
       &                     e_sub_old(i+1,k)               &
       &                     + e_sub_old(i+1,k-1)           &
       &                     - e_sub_old(i-1,k)             &
       &                     - e_sub_old(i-1,k-1)           &
       &                    )/(4.0d0*dx)                    &
       &                )/dz                                &
       &              + g13xz(i,k)                          &
       &                *(                                  &
       &                  g13x(i+1,k)                       &
       &                  *(                                &
       &                    e_sub_old(i,k+1)                &
       &                    - e_sub_old(i,k)                &
       &                   )/dz                             &
       &                  - g13x(i,k)                       &
       &                    *(                              &
       &                      e_sub_old(i,k)                &
       &                      - e_sub_old(i,k-1)            &
       &                     )/dz                           &
       &                 )/dz                               &
       &              + 1.0d0/g_sqrt(i,k)                   &
       &                *(                                  &
       &                  1.0d0/g_sqrt_z(i,k+1)             &
       &                  *(                                &
       &                    e_sub_old(i,k+1)                &
       &                    - e_sub_old(i,k)                &
       &                   )/dz                             &
       &                  - 1.0d0/g_sqrt_z(i,k)             &
       &                    *(                              &
       &                      e_sub_old(i,k)                &
       &                      - e_sub_old(i,k-1)            &
       &                     )/dz                           &
       &                 )/dz                               &
       &             )                                      &
       &          + 2.0d0*cm*ml*sqrt(e_sub_old(i,k))        &
       &            *(                                      &
       &              (                                     &
       &               (                                    &
       &                sqrt(e_sub_old(i+1,k))              &
       &                - sqrt(e_sub_old(i-1,k))            &
       &               )/(2.0d0*dx)                         &
       &               + g13xz(i,k)                         &
       &                 *(                                 &
       &                   sqrt(e_sub_old(i,k+1))           &
       &                   - sqrt(e_sub_old(i,k-1))         &
       &                   )/(2.0d0*dz)                     &
       &               )**2.0d0                             &
       &               + (                                  &
       &                  1.0d0/g_sqrt_z(i,k)               &
       &                  *(                                &
       &                    sqrt(e_sub_old(i,k+1))          &
       &                    - sqrt(e_sub_old(i,k-1))        &
       &                   )/(2.0d0*dz)                     &
       &                 )**2.0d0                           &
       &              )                                     &
       &          - cm/ml*e_sub_old(i,k)**(3.0d0/2.0d0)     &
       &         )

  !--- 渦混合係数
  km_sub(i,k) = cm*ml*sqrt(e_sub(i,k))

  write(*,*) i,k,e_sub(i,k),e_sub_old(i,k),e_sub_old2(i,k)
  write(*,*) e_sub_dif(i,k),e_sub_adv(i,k)

end subroutine eddymx_sub