!= Dynamical core module
!
! Authors::   Yasuhiro MORIKAWA, Yukiko YAMADA
! Version::   $Id: dynamics.f90,v 1.9 2006/08/21 08:06:14 morikawa Exp $
! Tag Name::  $Name: dcpam3-20061118 $
! Copyright:: Copyright (C) GFD Dennou Club, 2004-2005. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module dynamics_mod
  !
  !== Overview
  !
  ! Calculate Dynamical Core.
  !
  ! 力学コア部分を演算するモジュール。演算している方程式系、
  ! および離散化の手法は以下の通り。
  !
  ! 支配方程式系 : Governing Equations
  ! * 球座標プリミティブ方程式系 :
  !   Primitive Equations in Spherical Coordinate
  !
  ! 水平離散化 : Horizontal Discretization
  ! * スペクトル法 : Spectral Method
  !   * 三角形切断   : Triangle Truncation
  !   * 変換法       : Transform Method
  !
  ! 鉛直離散化 : Vertical Discretization
  ! * σ座標系    : Sigma Coordinate (Arakawa and Suarez(1983))
  !   * Lorenz 格子 : Lorenz Grid
  !
  ! 時間積分 : Time Integral
  ! * leap frog : リープフロッグ
  !
  !== Reference
  !
  ! * Arakawa, A., Suarez, M. J., 1983:
  !   Vertical differencing of the primitive equations
  !   in sigma coordinates.
  !   Mon. Wea. Rev., 111, 34--35.
  !
  !== Error Handling
  !
  !== Known Bugs
  !
  ! * 粘性計算の wa_NumVis_wa と wa_NumVisScaler_wa をとりあえず SPMODEL
  !   から移植したため、名称などがそのままである。
  !   本来はモジュール内部のものとして名称その他の検討が必要であろう。
  !
  !== Note
  !
  !== Future Plans
  !
  !現在、力学過程の全てがこのモジュール内にある。本来ならば、
  !解く方程式や項の種類によってモジュール化がなされるべきである。
  !

  use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
  implicit none
  private
  public :: dynamics_init, dynamics_leapfrog  ! subroutines
  public :: dynamics_diagnostic, dynamics_end ! subroutines
  public :: dynamics_diffusion                ! subroutines

  !-------------------------------------------------------------------
  !   dynamics_init において値を設定するもの
  !-------------------------------------------------------------------
  real(DBKIND) :: Kappa               ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱)
  real(DBKIND), pointer, save:: xy_Coli(:,:) =>null()
                                      ! 格子点データ(コリオリ力)
  real(DBKIND), pointer, save:: xy_UVFact(:,:) =>null()
                                      ! u→U, v→V のファクター
  real(DBKIND), pointer, save:: z_DelSigma(:) =>null()
                                      ! Δσ(整数レベル)
  real(DBKIND), pointer, save:: wz_DiffVorDiv(:,:) =>null()
                                      ! 運動量 (渦度発散) 水平拡散係数
  real(DBKIND), pointer, save:: wz_DiffTherm(:,:) =>null()
                                      ! 熱、比湿 水平拡散係数

  real(DBKIND), pointer, save:: z_TempAvrXY(:) =>null()
                                      ! 平均温度 (整数レベル)
  real(DBKIND), pointer, save:: r_TempAvrXY(:) =>null()
                                      ! 平均温度 (半整数レベル)

  real(DBKIND), pointer, save:: z_HydroAlpha(:) =>null()
                                      ! 静水圧の式の係数 α
  real(DBKIND), pointer, save:: z_HydroBeta(:) =>null()
                                      ! 静水圧の式の係数 β
  real(DBKIND), pointer, save:: z_TempInpolKappa(:) =>null()
                                      ! 温度鉛直補間の係数κ
  real(DBKIND), pointer, save:: z_TempInpolA(:) =>null()
                                      ! 温度鉛直補間の係数ａ
  real(DBKIND), pointer, save:: z_TempInpolB(:) =>null()
                                      ! 温度鉛直補間の係数ｂ

  ! semi-implicit用行列
  real(DBKIND), pointer, save:: z_siMtxG(:)   =>null()
                                      ! Ｇ = ＣｐκＴ
  real(DBKIND), pointer, save:: zz_siMtxGCt(:,:)  =>null()
                                      ! ＧＣ^{Transposed}
                                      !   (Ｃ = Δσ)
  real(DBKIND), pointer, save:: zz_siMtxW(:,:)   =>null()
                                      ! Ｗ:発散の式での温度補正
                                      !    (線形重力波項の効果)
  real(DBKIND), pointer, save:: zz_siMtxWH(:,:)  =>null()
                                      ! Ｗｈ:
  real(DBKIND), pointer, save:: zz_siMtxH(:,:)  =>null()
                                      ! ｈ: ＱＳ (σ移流) - Ｒ
  real(DBKIND), pointer, save:: zz_siMtxQ(:,:) =>null()
                                      ! Ｑ: dＴ/dσ (線形重力波項の効果)
  real(DBKIND), pointer, save:: zz_siMtxS(:,:) =>null()
                                      ! Ｓ: dσ/dｔ (線形重力波項の効果)
  real(DBKIND), pointer, save:: zz_siMtxQS(:,:) =>null()
                                      ! ＱＳ:
  real(DBKIND), pointer, save:: zz_siMtxR(:,:) =>null()
                                      ! Ｒ: ＲＤ =
                                      !      κＴ（∂π/∂ｔ＋ｄσ/ｄｔ／σ）
                                      !     気圧変化の効果の係数
                                      !     (線形重力波項の効果)


  !-------------------------------------------------------------------
  !   dynamics_leapfrog において利用するもの
  !-------------------------------------------------------------------
  real(DBKIND), pointer, save:: xy_lnPs(:,:) =>null()
                                        ! π＝ ln Ｐs (t)
  real(DBKIND), pointer, save:: xy_DlnPsDtN(:,:) =>null()
                                        ! πの時間変化 dπ/dt (t)
  real(DBKIND), pointer, save:: w_DlnPsDtN(:) =>null()
                                        ! πの時間変化 (スペクトル) (t)
  real(DBKIND), pointer, save:: w_lnPsN(:) =>null()
                                        ! π (t)
  real(DBKIND), pointer, save:: xy_GradLonlnPsN(:,:) =>null()
                                        ! dπ/dx (t)
  real(DBKIND), pointer, save:: xy_GradLatlnPsN(:,:) =>null()
                                        ! dπ/dy (t)
  real(DBKIND), pointer, save:: wz_TempN(:,:) => null()
                                        ! 温度 (t)
  real(DBKIND), pointer, save:: w_lnPsA(:) =>null()
                                        ! π (t+Δt)

  real(DBKIND), pointer, save:: xyz_lnPsAdv(:,:,:) =>null()
                                        ! π の移流 ｖ・∇π
  real(DBKIND), pointer, save:: xyz_lnPsAdvSum(:,:,:) =>null()
                                        ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ
  real(DBKIND), pointer, save:: xyz_DivSum(:,:,:) =>null()
                                        ! 発散の積下げ Σ[k〜K] ＤΔσ
  real(DBKIND), pointer, save:: xyr_VSigma(:,:,:) =>null()
                                        ! 鉛直σ速度(半整数レベル) (t)
  real(DBKIND), pointer, save:: xyr_VSigmaNonG(:,:,:) => null()
                                        ! 鉛直σ速度 非重力波成分
                                        !         (半整数レベル) (t)

  real(DBKIND), pointer, save:: xyz_TempEdd(:,:,:) =>null()
                                        ! Ｔ' ：温度の擾乱 (整数レベル) (t)
  real(DBKIND), pointer, save:: xyr_TempEdd(:,:,:) =>null()
                                        ! Ｔ' ：温度の擾乱 (半整数レベル)(t)
  real(DBKIND), pointer, save:: xyz_TempVir(:,:,:) =>null()
                                        ! Ｔv ：仮温度       (t)
  real(DBKIND), pointer, save:: xyz_TempVirEdd(:,:,:) =>null()
                                        ! Ｔv'：仮温度の擾乱 (t)

  real(DBKIND), pointer, save:: xyz_UAN(:,:,:) =>null()
                                        ! 東西運動量変化項ＵＡ
  real(DBKIND), pointer, save:: xyz_VAN(:,:,:) =>null()
                                        ! 南北運動量変化項ＶＡ

  real(DBKIND), pointer, save:: wz_DVorDtN(:,:) =>null()
                                        ! 渦度変化のスペクトルデータ (Δt)
  real(DBKIND), pointer, save:: wz_VorA(:,:) =>null()
                                        ! 渦度のスペクトルデータ    (t+Δt)

  real(DBKIND), pointer, save:: xyz_KE(:,:,:) =>null()
                                        ! 運動エネルギー + 仮温度補正
                                        ! ｕ**2+ｖ**2    + ΣＷ(Ｔv-Ｔ)
  real(DBKIND), pointer, save:: xy_Phi(:,:) =>null()
                                        ! 地表ジオポテンシャルΦ

  real(DBKIND), pointer, save:: wz_DDivDtN(:,:) =>null()
                                        ! 発散変化のスペクトルデータ (Δt)
  real(DBKIND), pointer, save:: wz_DivA(:,:) =>null()
                                        ! 発散のスペクトルデータ (t+Δt)
  real(DBKIND), pointer, save:: wz_PresTendTemp(:,:) =>null()
                                        ! 温度による圧力傾度のスペクトルデータ
  real(DBKIND), pointer, save:: wz_PresTendPs(:,:) =>null()
                                        ! 地表圧力による圧力傾度のスペクトルデータ

  real(DBKIND), pointer, save:: xyz_DTempLocalDtN(:,:,:) => null()
                                        ! 温度の局所時間変化項  Ｈ
                                        ! (水平発散 Ｔ'Ｄ + σ移流
                                        !  + π変化の効果) (全て非重力波項)

  real(DBKIND), pointer, save:: wz_DTempDtN(:,:) =>null()
                                        ! 温度の変化スペクトルデータ (Δt)
  real(DBKIND), pointer, save:: wz_TempA(:,:) =>null()
                                        ! 温度のスペクトルデータ   (t+Δt)

  real(DBKIND), pointer, save:: xyz_QVapDivVSigmaAdv(:,:,:) => null()
                                        ! Ｒ＝ｑＤ＋σ移流
  real(DBKIND), pointer, save:: wz_DQVapDtN(:,:) =>null()
                                        ! 比湿の変化スペクトルデータ (Δt)
  real(DBKIND), pointer, save:: wz_QVapA(:,:) =>null()
                                        ! 比湿のスペクトルデータ   (t+Δt)

  !-------------------------------------------------------------------
  !   dynamics_diagnostic で計算する変数
  !-------------------------------------------------------------------
  real(DBKIND), pointer:: wz_PsiA(:,:) =>null() ! スペクトル(流線関数) (t+Δt)
  real(DBKIND), pointer:: wz_ChiA(:,:) =>null() ! スペクトル(ポテンシャル) (t+Δt)

  logical, save :: dynamics_initialized = .false.
  character(*), parameter:: version = &
    & '$Name: dcpam3-20061118 $' // &
    & '$Id: dynamics.f90,v 1.9 2006/08/21 08:06:14 morikawa Exp $'

contains

  subroutine dynamics_init( &
    & x_Lon, y_Lat, z_Sigma, r_Sigma & !(in)
    & )
    !
    ! Initialize module and Calculate Non-Predictional Value.
    !
    ! 以降のサブルーチンで用いる変数の allocate 、および
    ! 時間発展しない量の演算を行なう。
    ! また、変数データ出力のための初期設定も行なう。
    !
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use grid_3d_mod,         only: im, jm, km
    use grid_wavenumber_mod, only: nm
    use constants_mod, only: constants_init, R0, Omega, Cp, RAir, &
         &                   TempAve, VisOrder, EFoldTime
    use time_mod,    only: DelTime
    use spml_mod,    only: spml_init, xy_Lat, rn
    use io_gt4_out_mod,only: io_gt4_out_init, io_gt4_out_SetVars
    use dc_trace,    only: DbgMessage, BeginSub, EndSub, DataDump
    use dc_string,   only: toChar
    implicit none
    real(DBKIND), intent(in) :: x_Lon(:)   ! 経度座標
    real(DBKIND), intent(in) :: y_Lat(:)   ! 緯度座標
    real(DBKIND), intent(in) :: z_Sigma(:) ! σレベル(整数)座標
    real(DBKIND), intent(in) :: r_Sigma(:) ! σレベル(半整数)座標

    !----- 拡散係数計算用変数 -----
    real(DBKIND)          :: VisCoef             ! 超粘性係数
    real(DBKIND), pointer :: wz_rn(:,:) =>null() ! ラプラシアンの係数 -n*(n+1)

    !----- 作業用内部変数 -----
    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)   :: i, j, k, kk, l

    character(STRING),  parameter:: subname = "dynamics_init"
  continue
    !----------------------------------------------------------------
    !   Check Initialization
    !----------------------------------------------------------------
    call BeginSub(subname, version)
    if (dynamics_initialized) then
       call EndSub( subname, '%c is already called.', c1=trim(subname) )
       return
    else
       dynamics_initialized = .true.
    endif

    !-------------------------------------------------------------------
    !   物理定数の取得
    !-------------------------------------------------------------------
    call constants_init

    !-------------------------------------------------------------------
    !   SPMODEL の初期化
    !-------------------------------------------------------------------
    call spml_init

    !-------------------------------------------------------------------
    !   io_gt4_out の初期化
    !-------------------------------------------------------------------
    call io_gt4_out_init

    !-------------------------------------------------------------------
    !   コリオリ力の設定
    !-------------------------------------------------------------------
    allocate( xy_Coli(im, jm) )
    xy_Coli = 2.0 * Omega * sin(xy_Lat)

!--
!!$    call DataDump('xy_Coli', xy_Coli, strlen=80)
!!$    call DataDump('xy_Lat', xy_Lat, strlen=80)
!++

    !-------------------------------------------------------------------
    !   u→U, v→V のファクターの計算
    !-------------------------------------------------------------------
    allocate( xy_UVFact(im, jm) )
    xy_UVFact = cos( xy_Lat )

!--
!!$    call DataDump('xy_UVFact', xy_UVFact, strlen=80)
!++

    !-------------------------------------------------------------------
    !   Δσ (整数、半整数) の計算
    !-------------------------------------------------------------------
    allocate( z_DelSigma(km) )

    do k = 1, km
      z_DelSigma(k) = r_Sigma(k) - r_Sigma(k+1)
    enddo

!--
!!$    call DataDump('z_Delsigma', z_DelSigma, strlen=80)
!++

    !-------------------------------------------------------------------
    !   運動量 (渦度発散) 拡散係数 wz_DiffVorDiv, 
    !   熱拡散係数 wz_DiffTherm の計算
    !-------------------------------------------------------------------
    allocate( wz_DiffVorDiv((nm+1)*(nm+1), km) )
    allocate( wz_DiffTherm ((nm+1)*(nm+1), km) )
    allocate( wz_rn        ((nm+1)*(nm+1), km) )

    do k = 1, km
      wz_rn(:,k) = rn(:,1)
    enddo

    ! 粘性係数の計算 (最大波数の e-folding time が EFoldTime となるように)
    VisCoef = (real( (nm * (nm+1)), DBKIND ) / R0**2 )**(-VisOrder / 2) &
      &        / EFoldTime

    call DbgMessage('VisCoef=<%f>', d=(/VisCoef/))

    wz_DiffTherm  = - VisCoef * ( (-wz_rn / R0**2)**(VisOrder / 2) )

    wz_DiffVorDiv =   wz_DiffTherm  &
       &            - VisCoef * ( - (2.0d0 / R0**2)**(VisOrder / 2))

!--
!!$    call DataDump('wz_DiffTherm' , wz_DiffTherm,  strlen=80)
!!$    call DataDump('wz_DiffVorDiv', wz_DiffVorDiv, strlen=80)
!++

    deallocate(wz_rn)


    !-------------------------------------------------------------------
    !   静水圧の式の係数α、βの計算
    !-------------------------------------------------------------------
    allocate(z_HydroAlpha(km))
    allocate(z_HydroBeta(km))

    Kappa = RAir / Cp    ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱)
    do k = 1, km
      z_HydroAlpha(k) = ( r_Sigma(k) / z_Sigma(k) )**Kappa - 1.0
      z_HydroBeta(k)  = 1.0 - ( r_Sigma(k+1) / z_Sigma(k) )**Kappa
    enddo

!--
!!$    call DbgMessage('Kappa=<%f>',  d=(/Kappa/) )
!!$    call DataDump('z_HydroAlpha', z_HydroAlpha, strlen=80)
!!$    call DataDump('z_HydroBeta', z_HydroBeta, strlen=80)
!++

    !-------------------------------------------------------------------
    !   温度鉛直補間の係数ａ、ｂの計算
    !-------------------------------------------------------------------
    allocate(z_TempInpolA(km))
    allocate(z_TempInpolB(km))
    allocate(z_TempInpolKappa(km))

    do k = 1, km
      z_TempInpolA(k) = &
        & z_HydroAlpha(k) / ( 1.0 - (z_Sigma(k) / z_Sigma(k-1))**Kappa )

      z_TempInpolB(k) = &
        & z_HydroBeta(k) / ( (z_Sigma(k) / z_Sigma(k+1))**Kappa - 1.0 )

      z_TempInpolKappa(k) = &
        & ( r_Sigma(k) * z_HydroAlpha(k) &
        &   + r_Sigma(k+1) * z_HydroBeta(k) ) / z_DelSigma(k)
    enddo
    z_TempInpolA(1) = 0.0
    z_TempInpolB(km) = 0.0

!--
!!$    call DataDump('z_TempInpolA', z_TempInpolA, strlen=80)
!!$    call DataDump('z_TempInpolB', z_TempInpolB, strlen=80)
!!$    call DataDump('z_TempInpolKappa', z_TempInpolKappa, strlen=80)
!++

    !-------------------------------------------------------------------
    !   平均温度 (整数レベル、半整数レベル) の計算
    !-------------------------------------------------------------------
    allocate( z_TempAvrXY(km) )
    allocate( r_TempAvrXY(km+1) )

    z_TempAvrXY = TempAve

    !----- 下端の平均温度設定 -----
    ! ※ 正しいかは微妙だが一応変な値が入らないように
    r_TempAvrXY = 0.0d0

    do k = 2, km
      r_TempAvrXY(k) = &
        &   z_TempInpolA(k)   * z_TempAvrXY(k)   &
        & + z_TempInpolB(k-1) * z_TempAvrXY(k-1)
    enddo

!--
!!$    call DataDump('z_TempAvrXY', z_TempAvrXY, strlen=80)
!!$    call DataDump('r_TempAvrXY', r_TempAvrXY, strlen=80)
!++


    !-------------------------------------------------------------------
    !   semi-implicit 用行列の計算
    !-------------------------------------------------------------------
    allocate( z_siMtxG(km) )
    allocate( zz_siMtxGCt(km,km) )
    allocate( zz_siMtxW(km,km) )
    allocate( zz_siMtxWH(km,km) )
    allocate( zz_siMtxH(km,km) )
    allocate( zz_siMtxQ(km,km) )
    allocate( zz_siMtxS(km,km) )
    allocate( zz_siMtxQS(km,km) )
    allocate( zz_siMtxR(km,km) )

    ! Ｇ = ＣｐκＴ
    z_siMtxG = Cp * z_TempInpolKappa * z_TempAvrXY

    ! ＧＣ^{Transposed}     (Ｃ = Δσ)
    do k = 1, km
      do kk = 1, km
        zz_siMtxGCt(k, kk) = z_siMtxG(k) * z_DelSigma(kk)
      end do
    end do


    ! Ｗ:発散の式での温度補正 (線形重力波項の効果)

    ! 初期化
    zz_siMtxW = 0.0

    do k = 1, km
      do kk = 1, k
        zz_siMtxW(k, kk) = Cp * z_HydroAlpha(kk)
      enddo

      do kk = 1, k-1
        zz_siMtxW(k, kk) = zz_siMtxW(k, kk) + Cp * z_HydroBeta(kk)
      enddo
    enddo

    ! Ｓ: dσ/dｔ (線形重力波項の効果)

    ! 初期化
    zz_siMtxS = 0.0

    do k = 1, km
      do kk = 1, km
        zz_siMtxS(k,kk) = r_Sigma(k) * z_DelSigma(kk)
      enddo

      do kk = k, km
        zz_siMtxS(k,kk) = zz_siMtxS(k,kk) - z_DelSigma(kk)
      enddo

    enddo
!--
!!$    call DataDump('zz_siMtxS', zz_siMtxS, strlen=80)
!++

    ! Ｑ: dＴ/dσ (線形重力波項の効果)

    ! 初期化
    zz_siMtxQ = 0.0

    do k = 1, km
      zz_siMtxQ(k,k) = ( r_TempAvrXY(k) - z_TempAvrXY(k) ) / z_DelSigma(k)
    enddo
    do k = 1, km - 1
      zz_siMtxQ(k,k+1) = ( z_TempAvrXY(k) - r_TempAvrXY(k+1) ) / z_DelSigma(k)
    enddo
!--
!!$    call DataDump('zz_siMtxQ', zz_siMtxQ, strlen=80)
!++

    ! ＱＳ

    ! 初期化
    zz_siMtxQS = 0.0

    zz_siMtxQS = matmul(zz_siMtxQ, zz_siMtxS)


    ! Ｒ: ＲＤ = κＴ（∂π/∂ｔ＋ｄσ/ｄｔ／σ）
    !     気圧変化の効果の係数  (線形重力波項の効果)

    ! 初期化
    zz_siMtxR = 0.0

    do k = 1, km
      do kk = k, km
        zz_siMtxR(k,kk) = &
          & - z_HydroAlpha(k) / z_DelSigma(k) * z_DelSigma(kk) * z_TempAvrXY(k)
      enddo

      do kk = k+1, km
        zz_siMtxR(k,kk) = zz_siMtxR(k,kk)  &
          & - z_HydroBeta(k) / z_DelSigma(k) * z_DelSigma(kk) * z_TempAvrXY(k)
      enddo
    enddo
!--
!!$    call DataDump('zz_siMtxR', zz_siMtxR, strlen=80)
!++

    ! ｈ: ＱＳ (σ移流) - Ｒ

    ! 初期化
    zz_siMtxH = 0.0

    zz_siMtxH = zz_siMtxQS - zz_siMtxR

    ! Ｗｈ:

    ! 初期化
    zz_siMtxWH = 0.0

    zz_siMtxWH = matmul(zz_siMtxW, zz_siMtxH)

!--
!!$    call DataDump('zz_siMtxH', zz_siMtxH, strlen=80)
!++

    !-------------------------------------------------------------------
    !   データ出力用の初期設定
    !-------------------------------------------------------------------
    call io_gt4_out_SetVars('xyz_DivSum')
    call io_gt4_out_SetVars('xyz_lnPsAdv')
    call io_gt4_out_SetVars('xyz_lnPsAdvSum')
    call io_gt4_out_SetVars('xyz_UAN')
    call io_gt4_out_SetVars('xyz_VAN')
    call io_gt4_out_SetVars('xyz_KE')
    call io_gt4_out_SetVars('xyz_PresTendTemp')
    call io_gt4_out_SetVars('xyz_PresTendPs')
    call io_gt4_out_SetVars('xyr_VSigma')
    call io_gt4_out_SetVars('xyr_VSigmaNonG')
    call io_gt4_out_SetVars('xyz_DTempLocalDtN')
    call io_gt4_out_SetVars('xyz_TempAdv')
    call io_gt4_out_SetVars('xyr_TempEdd')
    call io_gt4_out_SetVars('xyz_Temp_T')
    call io_gt4_out_SetVars('xyz_DTempLocalDtN_lnPsAdvSum')

    call io_gt4_out_SetVars('TotalMass')

    call io_gt4_out_SetVars('xyz_Vor_Diffusion')
    call io_gt4_out_SetVars('xyz_Div_Diffusion')
    call io_gt4_out_SetVars('xyz_Temp_Diffusion')
    call io_gt4_out_SetVars('xyz_QVap_Diffusion')

    !-------------------------------------------------------------------
    !   dynamics_leapfrog で計算する変数の allocate 
    !-------------------------------------------------------------------
    allocate  & 
      & ( xy_lnPs (im, jm)          , & ! π＝ ln Ｐs           (t)
      &   xyz_lnPsAdv(im, jm, km)   , & ! π の移流 ｖ・∇π
      &   xyz_lnPsAdvSum(im, jm, km), & ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ
      &   xyz_DivSum(im, jm, km)    , & ! 発散の積下げ Σ[k〜K] ＤΔσ

      &   xy_DlnPsDtN( im, jm )       , & ! πの時間変化 dπ/dt   (Δt)
      &   w_DlnPsDtN( (nm+1)*(nm+1) ) , & ! πの時間変化 (スペクトル) (Δt)
      &   w_lnPsN( (nm+1)*(nm+1) ) , & ! π (t)
      &   xy_GradLonlnPsN( im, jm ), & ! dπ/dx (t)
      &   xy_GradLatlnPsN( im, jm ), & ! dπ/dy (t)
      &   wz_TempN( (nm+1)*(nm+1), km ) , & ! 温度 (t)
      &   w_lnPsA( (nm+1)*(nm+1) ) , & ! π (t+Δt)

      &   xyr_VSigma(im, jm, km+1) , & ! 鉛直σ速度(半整数レベル) (t)
      &   xyr_VSigmaNonG(im, jm, km+1) , & ! 鉛直σ速度 非重力波成分 (半整数レベル) (t)

      &   xyz_TempEdd(im, jm, km)   , & ! Ｔ' ：温度の擾乱 (整数レベル) (t)
      &   xyz_TempVir(im, jm, km)   , & ! Ｔv ：仮温度 (virtual temperature)(t)
      &   xyz_TempVirEdd(im, jm, km ), &! Ｔv'：仮温度の擾乱               (t)
      &   xyr_TempEdd(im, jm, km+1), &! Ｔ' ：温度の擾乱 (半整数レベル)   (t)

      &   xyz_UAN(im, jm, km)      , & ! 東西運動量変化項ＵＡ
      &   xyz_VAN(im, jm, km)      , & ! 南北運動量変化項ＶＡ

      &   wz_DVorDtN( (nm+1)*(nm+1), km ) , & ! 渦度変化のスペクトルデータ (Δt) 
      &   wz_VorA( (nm+1)*(nm+1), km ) , & ! 渦度のスペクトルデータ    (t+Δt)

      &   xyz_KE(im, jm, km)      , & ! 運動エネルギー + 仮温度補正
                                      ! ｕ**2+ｖ**2    + ΣＷ(Ｔv-Ｔ)

      &   wz_DDivDtN( (nm+1)*(nm+1), km) , & ! 発散のスペクトルデータ (Δt)
      &   wz_DivA( (nm+1)*(nm+1), km) , & ! 発散のスペクトルデータ (t+Δt)
      &   wz_PresTendTemp( (nm+1)*(nm+1), km) , & ! 温度による圧力傾度のスペクトルデータ
      &   wz_PresTendPs( (nm+1)*(nm+1), km) , & ! 地表圧力による圧力傾度のスペクトルデータ

      &   xyz_DTempLocalDtN(im, jm, km)  , & ! 温度の局所時間変化項  Ｈ
                                           ! (水平発散 Ｔ'Ｄ + σ移流
                                           !  + π変化の効果) (全て非重力波項)

      &   wz_DTempDtN( (nm+1)*(nm+1), km ) , & ! 温度の変化スペクトルデータ (Δt)
      &   wz_TempA( (nm+1)*(nm+1), km ) , & ! 温度のスペクトルデータ   (t+Δt)

      &   xyz_QVapDivVSigmaAdv( im, jm, km), & ! Ｒ＝ｑＤ＋σ移流
      &   wz_DQVapDtN( (nm+1)*(nm+1), km) , & ! 比湿の変化スペクトルデータ (Δt)
      &   wz_QVapA( (nm+1)*(nm+1), km)   & ! 比湿のスペクトルデータ   (t+Δt)
      & )


    !-------------------------------------------------------------------
    !   dynamics_diagnostic で計算する変数の allocate 
    !-------------------------------------------------------------------
    allocate  &
      & ( wz_PsiA( (nm+1)*(nm+1), km) , & ! スペクトル(流線関数)
      &   wz_ChiA( (nm+1)*(nm+1), km) )   ! スペクトル(ポテンシャル)


    call EndSub(subname)
  end subroutine dynamics_init


  subroutine dynamics_leapfrog &
    & ( x_Lon      , y_Lat      , z_Sigma , r_Sigma  , & !(in)

    &   xyz_VelLonB, xyz_VelLatB, xyz_VorB, xyz_DivB , & !(in)
    &   xyz_TempB  , xyz_QVapB  , xy_PsB  ,            & !(in)

    &   xyz_VelLonN, xyz_VelLatN, xyz_VorN, xyz_DivN , & !(in)
    &   xyz_TempN  , xyz_QVapN  , xy_PsN  ,            & !(in)

    &   xyz_VelLonA, xyz_VelLatA, xyz_VorA, xyz_DivA , & !(out)
    &   xyz_TempA  , xyz_QVapA  , xy_PsA             )   !(out)
    !
    ! Calculate Predictional Values.
    !
    ! 時間発展する量の演算を行なう。
    ! 演算するデータの出力も行なう。
    !
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use grid_3d_mod,         only: im, jm, km
    use grid_wavenumber_mod, only: nm
    use constants_mod, only: R0, Cp, EpsVT
    use time_mod,    only: DelTime, CurrentTime
    use spml_mod,    only: w_xy, xy_w , xy_GradLon_w, xy_GradLat_w, &
      &                    w_Div_xy_xy, w_LaplaInv_w,               &
      &                    wa_xya, xya_wa, wa_Div_xya_xya, wa_Lapla_wa
    use io_gt4_out_mod,only: io_gt4_out_Put
    use dc_trace,    only: DbgMessage, BeginSub, EndSub, DataDump
    use dc_string,   only: toChar
    implicit none

    real(DBKIND), intent(in) ::  x_Lon(:)           ! 経度座標
    real(DBKIND), intent(in) ::  y_Lat(:)           ! 緯度座標
    real(DBKIND), intent(in) ::  z_Sigma(:)         ! σレベル(整数)座標
    real(DBKIND), intent(in) ::  r_Sigma(:)         ! σレベル(半整数)座標
    
    real(DBKIND), intent(in) ::  xyz_VelLonB(:,:,:)  ! 速度経度成分 (t-Δt)
    real(DBKIND), intent(in) ::  xyz_VelLatB(:,:,:)  ! 速度緯度成分 (t-Δt)
    real(DBKIND), intent(in) ::  xyz_VorB(:,:,:)     ! 渦度         (t-Δt)
    real(DBKIND), intent(in) ::  xyz_DivB(:,:,:)     ! 発散         (t-Δt)
    real(DBKIND), intent(in) ::  xyz_TempB(:,:,:)    ! 温度         (t-Δt)
    real(DBKIND), intent(in) ::  xyz_QVapB(:,:,:)    ! 比湿         (t-Δt)
    real(DBKIND), intent(in) ::  xy_PsB(:,:)         ! 地表面気圧   (t-Δt)
    
    real(DBKIND), intent(in) ::  xyz_VelLonN(:,:,:)  ! 速度経度成分 (t)
    real(DBKIND), intent(in) ::  xyz_VelLatN(:,:,:)  ! 速度緯度成分 (t)
    real(DBKIND), intent(in) ::  xyz_VorN(:,:,:)     ! 渦度         (t)
    real(DBKIND), intent(in) ::  xyz_DivN(:,:,:)     ! 発散         (t)
    real(DBKIND), intent(in) ::  xyz_TempN(:,:,:)    ! 温度         (t)
    real(DBKIND), intent(in) ::  xyz_QVapN(:,:,:)    ! 比湿         (t)
    real(DBKIND), intent(in) ::  xy_PsN(:,:)         ! 地表面気圧   (t)

    real(DBKIND), intent(out) ::  xyz_VelLonA(:,:,:)! 速度経度成分 (t+Δt)
    real(DBKIND), intent(out) ::  xyz_VelLatA(:,:,:)! 速度緯度成分 (t+Δt)
    real(DBKIND), intent(out) ::  xyz_VorA(:,:,:)   ! 渦度         (t+Δt)
    real(DBKIND), intent(out) ::  xyz_DivA(:,:,:)   ! 発散         (t+Δt)
    real(DBKIND), intent(out) ::  xyz_TempA(:,:,:)  ! 温度         (t+Δt)
    real(DBKIND), intent(out) ::  xyz_QVapA(:,:,:)  ! 比湿         (t+Δt)
    real(DBKIND), intent(out) ::  xy_PsA(:,:)       ! 地表面気圧   (t+Δt)

    !----- 作業用内部変数 -----
    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)   :: i, j, k, kk, l


    character(STRING),  parameter:: subname = "dynamics_leapfrog"
  continue
    !----------------------------------------------------------------
    !   Check Initialization
    !----------------------------------------------------------------
    call BeginSub(subname, version)
    if (.not. dynamics_initialized) then
      call EndSub( subname, 'Call dynamics_init before call %c',  &
        &       c1=trim(subname) )
      return
    endif

    !-------------------------------------------------------------------
    !   地表気圧変化・鉛直σ速度の計算  (連続の式、熱の式で利用される)
    !-------------------------------------------------------------------

    xy_lnPs = log( xy_PsN )

    w_lnPsN = w_xy(xy_lnPs)
    xy_GradLonlnPsN = xy_GradLon_w( w_lnPsN )
    xy_GradLatlnPsN = xy_GradLat_w( w_lnPsN )

    ! 地表面気圧 π の移流 = ｖ・∇π
    do k = 1, km
      xyz_lnPsAdv(:,:,k) = &
        &  ( xyz_VelLonN(:,:,k) * xy_GradLonlnPsN &
        &    + xyz_VelLatN(:,:,k) * xy_GradLatlnPsN &
        &   ) / R0
    enddo
    call io_gt4_out_Put('xyz_lnPsAdv', xyz_lnPsAdv)
!--
!!$    call DataDump('xyz_lnPsAdv', xyz_lnPsAdv, strlen=80)
!++


    ! π移流の積み下げ Σ[k〜K] (ｖ・∇π) Δσ
    xyz_lnPsAdvSum(:,:,km) = xyz_lnPsAdv(:,:, km ) * z_DelSigma( km )

    do k = km-1 , 1, -1
      xyz_lnPsAdvSum(:,:,k) = &
        &   xyz_lnPsAdvSum(:,:,k+1) &
        & + xyz_lnPsAdv(:,:,k) * z_DelSigma(k)
    enddo
    call io_gt4_out_Put('xyz_lnPsAdvSum', xyz_lnPsAdvSum)
!--
!!$    call DataDump('xyz_lnPsAdvSum', xyz_lnPsAdvSum, strlen=80)
!++


    ! 発散の積下げ Σ[k〜K] ＤΔσ
    xyz_DivSum(:,:,km) = xyz_DivN(:,:, km ) * z_DelSigma( km )

    do k = km -1 , 1, -1
      xyz_DivSum(:,:,k) = &
        &   xyz_DivSum(:,:,k+1) + xyz_DivN(:,:,k) * z_DelSigma(k)
    enddo

    call io_gt4_out_Put('xyz_DivSum', xyz_DivSum)
!--
!!$    call DataDump('xyz_DivSum', xyz_DivSum, strlen=80)
!++

    !-------------------------------------------------------------------
    !   Ps の計算。連続の式を解く
    !-------------------------------------------------------------------

    !----- π = lnＰs の時間変化 dπ/dt の計算 -----
    !
    ! π移流の積み下げ Σ[1〜K] (ｖ・∇π) Δσ を格子点上で解く。
    !
    xy_DlnPsDtN = - xyz_lnPsAdvSum(:,:,1)   ! - Σ[1〜K](ｖ・∇π)Δσ

!--
!!$    call DataDump('xy_DlnPsDtN', xy_DlnPsDtN, strlen=80)
!++


    !----- 発散の積み下げ -Σ[1〜K] ＤΔσ をスペクトル空間で加える。-----
    w_DlnPsDtN = w_xy( xy_DlnPsDtN ) - w_xy( xyz_DivSum(:,:,1) )

    !----- A = B + 2Δt * T-----
    w_lnPsA = w_xy(  log( xy_PsB )  ) + 2.0*DelTime * ( w_DlnPsDtN )

!--
!!$    call DataDump('xy_PsB', xy_PsB, strlen=80)
!!$    call DataDump('w_lnPsB', w_xy(  log( xy_PsB )  ), strlen=80)
!!$    call DataDump('w_lnPsA', w_lnPsA, strlen=80)
!++

    !----- xy_PsA に代入-----
    xy_PsA = exp(  xy_w( w_lnPsA )  )

!--
!!$    call DataDump('xy_PsA', xy_PsA, strlen=80)
!++

    !-------------------------------------------------------------------
    !   鉛直σ速度 (整数、半整数レベル) の上昇流の計算。
    !-------------------------------------------------------------------

    ! σ速度 (3/2 〜 K-1/2)
    ! - σ[k-1/2] × ( Σ[1〜K](Ｄ + ｖ・∇π)Δσ )
    ! - Σ[k〜K] (Ｄ + ｖ・∇π) Δσ

    do k = 2, km+1 - 1
      xyr_VSigma(:,:,k) =  &
        & r_Sigma(k) * ( xyz_lnPsAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) &
        & - ( xyz_lnPsAdvSum(:,:,k) + xyz_DivSum(:,:,k) )

    ! - σ[k-1/2] × ( Σ[1〜K](ｖ・∇π)Δσ )
    ! - Σ[k〜K] (ｖ・∇π) Δσ

      xyr_VSigmaNonG(:,:,k) =  &
        & r_Sigma(k) * xyz_lnPsAdvSum(:,:,1) - xyz_lnPsAdvSum(:,:,k)
    enddo

    ! σ速度 (1/2, K+1/2)
    xyr_VSigma(:,:,1) = 0.0
    xyr_VSigma(:,:,km+1) = 0.0
    xyr_VSigmaNonG(:,:,1) = 0.0
    xyr_VSigmaNonG(:,:,km+1) = 0.0

    call io_gt4_out_Put('xyr_VSigma', xyr_VSigma)
    call io_gt4_out_Put('xyr_VSigmaNonG', xyr_VSigmaNonG)
!--
!!$    call DataDump('xyr_VSigma', xyr_VSigma, strlen=80)
!++

    !-------------------------------------------------------------------
    !   温度・仮温度の基本場からのずれ  (以降の ＵＡ、ＵＶ、Ｈなどの全てで利用)
    !-------------------------------------------------------------------
    do k = 1, km
      xyz_TempVir(:,:,k)   = &
        & xyz_TempN(:,:,k) * ( 1.0 + (EpsVT * xyz_QVapN(:,:,k)) )

      xyz_TempEdd(:,:,k)   = xyz_TempN(:,:,k) - z_TempAvrXY(k)
      xyz_TempVirEdd(:,:,k) = xyz_TempVir(:,:,k) - z_TempAvrXY(k)
    enddo

!--
!!$    call DataDump('xyz_Temp', xyz_Temp, strlen=80)
!!$    call DataDump('xyz_TempEdd', xyz_TempEdd, strlen=80)
!!$    call DataDump('xyz_TempVir', xyz_TempVir, strlen=80)
!!$    call DataDump('xyz_TempVirEdd', xyz_TempVirEdd, strlen=75)
!++

    !-------------------------------------------------------------------
    !   半整数レベルの温度の擾乱  (以降の ＵＡ、ＵＶ、Ｈなどの全てで利用)
    !-------------------------------------------------------------------

    do k = 2, km+1 - 1
      xyr_TempEdd(:,:,k) = &
        &   z_TempInpolA(k)   * xyz_TempN(:,:,k)   &
        & + z_TempInpolB(k-1) * xyz_TempN(:,:,k-1) &
        & - r_TempAvrXY(k)
    enddo

    call io_gt4_out_Put('xyr_TempEdd', xyr_TempEdd)
!--
!!$    call DataDump('xyr_TempEdd', xyr_TempEdd, strlen=80)
!++

    !-------------------------------------------------------------------
    !   ＵＡ＝Ｖζ＋σ移流、ＶＡ＝Ｕζ＋σ移流 (渦度と発散の式で利用)
    !   ※ AGCM5 のマニュアルの ＵＡ と異なり、UVFact = cosφは
    !      掛かっていないので注意。
    !-------------------------------------------------------------------

    do k = 1, km
       ! (ζ + ｆ) ｖ
       ! - (Ｃp κ (1/ａ) Ｔ'v (dπ/dλ) ) / cosφ

      xyz_UAN(:,:,k) = &
        & ( xyz_VorN(:,:,k) + xy_Coli ) * xyz_VelLatN(:,:,k) &
        & - Cp * z_TempInpolKappa(k) &
        &           * xyz_TempVirEdd(:,:,k) * xy_GradLonlnPsN / R0

       ! - (ζ + ｆ) ｕ
       ! - (Ｃp κ (1/ａ) Ｔ'v (dπ/dμ) ) / cosφ

      xyz_VAN(:,:,k) = &
        & - ( xyz_VorN(:,:,k) + xy_Coli ) * xyz_VelLonN(:,:,k) &
        & - Cp * z_TempInpolKappa(k) &
        &           * xyz_TempVirEdd(:,:,k) * xy_GradLatlnPsN / R0
    enddo

!--
!!$    do k = 1, km
!!$       call DataDump('xyz_Div(pi)', &
!!$            & xy_w( &
!!$            & w_Div_xy_xy( &
!!$            & - Cp * z_TempInpolKappa(k) &
!!$            &           * xyz_TempVirEdd(:,:,k)  &
!!$            &             * xy_GradLon_w(  w_xy( xy_lnPs )  ) &
!!$            &                / R0 &
!!$            & , &
!!$            &
!!$            & - Cp * z_TempInpolKappa(k) &
!!$            &           * xyz_TempVirEdd(:,:,k) &
!!$            &             * xy_GradLat_w(  w_xy( xy_lnPs )  ) &
!!$            &               / R0 &
!!$            &
!!$            &     )  &
!!$            &    ) / R0 &
!!$            &
!!$            & , strlen=80, multi=(/k/))
!!$    enddo
!++

!--
!!$    call DataDump('xyz_VelLon', xyz_VelLon, strlen=80)
!!$    call DataDump('xyz_VelLat', xyz_VelLat, strlen=80)
!!$    call DataDump('xyz_VAN', xyz_VAN(1,1,:), strlen=80, multi=(/-1/))
!!$    call DataDump('xyz_UAN', xyz_UAN(1,1,:), strlen=80, multi=(/-1/))
!!$    call DataDump('xyz_VAN', xyz_VAN(1,1,:), strlen=80, multi=(/-1/))
!++

    do k = 2, km

      ! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (ｕ[k-1] - ｕ[k])

      xyz_UAN(:,:,k) = xyz_UAN(:,:,k) &
        & - 1.0 / (2.0 * z_DelSigma(k)) &
        &    * xyr_VSigma(:,:,k) * ( xyz_VelLonN(:,:,k-1) - xyz_VelLonN(:,:,k) )

      ! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (ｖ[k-1] - ｖ[k])

      xyz_VAN(:,:,k) = xyz_VAN(:,:,k) &
        & - 1.0 / (2.0 * z_DelSigma(k)) &
        &    * xyr_VSigma(:,:,k) * ( xyz_VelLatN(:,:,k-1) - xyz_VelLatN(:,:,k) )
    enddo

!--
!!$    call DataDump('xyz_UAN', xyz_UAN(1,1,:), strlen=80, multi=(/-2/))
!!$    call DataDump('xyz_VAN', xyz_VAN(1,1,:), strlen=80, multi=(/-2/))
!++

    do k = 1, km - 1

      ! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (ｕ[k] - ｕ[k+1])

      xyz_UAN(:,:,k) = xyz_UAN(:,:,k) &
        & - 1.0 / (2.0 * z_DelSigma(k)) &
        &    * xyr_VSigma(:,:,k+1) * ( xyz_VelLonN(:,:,k) - xyz_VelLonN(:,:,k+1) )

      ! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (ｖ[k] - ｖ[k+1])

      xyz_VAN(:,:,k) = xyz_VAN(:,:,k) &
        & - 1.0 / (2.0 * z_DelSigma(k)) &
        &    * xyr_VSigma(:,:,k+1) * ( xyz_VelLatN(:,:,k) - xyz_VelLatN(:,:,k+1) )
    enddo

!--
!!$    call DataDump('xyz_UAN', xyz_UAN(1,1,:), strlen=80, multi=(/-3/))
!!$    call DataDump('xyz_VAN', xyz_VAN(1,1,:), strlen=80, multi=(/-3/))
!++

!--
!!$    call DataDump('xyz_UAN', xyz_UAN, strlen=80)
!!$    call DataDump('xyz_VAN', xyz_VAN, strlen=80)
!++

    call io_gt4_out_Put('xyz_UAN', xyz_UAN)
    call io_gt4_out_Put('xyz_VAN', xyz_VAN)

    !-------------------------------------------------------------------
    !   渦度 Vor の計算。渦度方程式を解く
    !-------------------------------------------------------------------

    ! dζ/dt の計算

    wz_DVorDtN = wa_Div_xya_xya( xyz_VAN, - xyz_UAN ) / R0
         !& + w_DiffV(k) * wa_xya(xyz_Vor(:,:,k))
         !  拡散係数の導出が厄介なので後回し

    !----- A = B + 2Δt * T-----
    wz_VorA = wa_xya( xyz_VorB ) + 2.0d0*DelTime * ( wz_DVorDtN )
!!!         & + 2.0d0*DelTime * wa_NumVis_wa(  wa_xya( xyz_VorB )  )


    !----- xyz_VorA に代入-----
    xyz_VorA = xya_wa( wz_VorA )

!--
!!$    call DataDump('wz_DVorDtN', wz_DVorDtN, strlen=70)
!!$    call DataDump('xyz_VorB', xyz_VorB, strlen=80)
!!$    call DataDump('xyz_Vor', xyz_Vor, strlen=80)
!!$    call DataDump('xyz_Vor_T', xya_wa(wz_DVorDtN), strlen=80)
!!$    call DataDump('xyz_VorA', xyz_VorA, strlen=80)
!++

    !-------------------------------------------------------------------
    !   Ｅ=ｕ**2+ｖ**2 + 仮温度補正 ( ΣＷ(Ｔv-Ｔ) )  (発散の式で利用)
    !-------------------------------------------------------------------

    ! 仮温度補正 ( ΣＷ(Ｔv-Ｔ) ) の計算。格子点上で積み上げ。

    ! 初期化
    xyz_KE = 0.0

    xyz_KE(:,:,1) = Cp * z_HydroAlpha(1) &
      & * ( xyz_TempVir(:,:,1) - xyz_TempN(:,:,1) )

    do k = 2, km
      xyz_KE(:,:,k) = xyz_KE(:,:,k-1)   &
        & + Cp * z_HydroAlpha(k) &
        &     * ( xyz_TempVir(:,:,k) - xyz_TempN(:,:,k) ) &
        & + Cp * z_HydroBeta(k-1) &
        &     * ( xyz_TempVir(:,:,k-1) - xyz_TempN(:,:,k-1) )
    enddo

!--
!!$    call DataDump('xyz_Temp', xyz_Temp, strlen=80)
!!$    call DataDump('xyz_TempVir', xyz_TempVir, strlen=80)
!!$    call DataDump('xyz_KE-w(Tv-T)', xyz_KE, strlen=80)
!++

    ! 運動エネルギー ｕ**2+ｖ**2 の加算

    xyz_KE = xyz_KE + ( xyz_VelLonN**2 + xyz_VelLatN**2 ) / 2.0

!--
!!$    call DataDump('xyz_KE-u**2+v**2', xyz_KE, strlen=80)
!++
    call io_gt4_out_Put('xyz_KE', xyz_KE)

    !-------------------------------------------------------------------
    !   発散 Div の計算。発散方程式を解く
    !-------------------------------------------------------------------

    ! ＷＴ = Ｃp α Ｔ[k] + Ｃp β Ｔ[k-1] をスペクトル空間で積み上げ

    ! 初期化
    wz_PresTendTemp = 0.0

    wz_TempN = wa_xya(xyz_TempN)

    wz_PresTendTemp(:,1) = Cp * z_HydroAlpha(1) * wz_TempN(:,1)

    do k = 2, km
      wz_PresTendTemp(:,k) = wz_PresTendTemp(:,k-1) &
        & + Cp * z_HydroAlpha(k)  * wz_TempN(:,k)   &
        & + Cp * z_HydroBeta(k-1) * wz_TempN(:,k-1)
    enddo

!--
!!$    call io_gt4_out_Put('xyz_PresTendTemp', xya_wa(wa_Lapla_wa(wz_PresTendTemp))/R0**2)
!++


!--
!!$    call DataDump('wz_PresTendTemp', wz_PresTendTemp, strlen=70)
!!$    call DataDump('xyz_Lapla(wz_PresTendTemp)/R0**2', &
!!$         & xya_wa(  wa_Lapla_wa( wz_PresTendTemp )  )/R0**2, strlen=60)
!++


    ! Ｇπ = (Ｃp κ Ｔ) π
    do k = 1, km
      wz_PresTendPs(:,k) = &
        & Cp * z_TempInpolKappa(k) * z_TempAvrXY(k) * w_xy(  log( xy_PsN )  )
    enddo

!--
!!$    call io_gt4_out_Put('xyz_PresTendPs', xya_wa(wa_Lapla_wa(wz_PresTendPs))/R0**2)
!++

!--
!!$    call DataDump('wz_PresTendPs', wz_PresTendPs, strlen=70)
!!$    call DataDump('xyz_Lapla(wz_PresTendPs)/R0**2', &
!!$         & xya_wa(  wa_Lapla_wa( wz_PresTendPs )  )/R0**2, strlen=60)
!++

!--
!!$    call DataDump('Div(wz_UAN - wz_VAN)', &
!!$         & wa_Div_xya_xya( xyz_UAN, - xyz_VAN ), &
!!$         & strlen=70)
!!$    call DataDump('Div(wz_UAN - wz_VAN)/R0', &
!!$         & wa_Div_xya_xya( xyz_UAN, - xyz_VAN ) /R0, &
!!$         & strlen=70)
!!$    call DataDump('xyz_Div(wz_UAN - wz_VAN)/R0', &
!!$         & xya_wa(  wa_Div_xya_xya( xyz_UAN, - xyz_VAN )  ) /R0, &
!!$         & strlen=70)
!!$    call DataDump('Lapla(wz_KE)', &
!!$         & wa_Lapla_wa( wa_xya(xyz_KE) ), strlen=70)
!++

    ! dＤ/dt の計算
    !  - ∇**2 (Φ + ＷＴ + Ｇπ)

    wz_DDivDtN = &
      & wa_Div_xya_xya( xyz_UAN, xyz_VAN ) / R0 &
      & - wa_Lapla_wa( wa_xya(xyz_KE) ) / R0**2 &
      !
      !& + w_DiffV(k) * wa_xya(xyz_Div(:,:,k)) !現在, 拡散の効果は
                                               !dynamics_diffusion で与えている.
      ! semi-implicitだとここには無い項達
      & - wa_Lapla_wa( &
      !&    w_xy( xy_Phi )  &  ! 地表ジオポテンシャルΦは後回し
      &    wz_PresTendTemp   & ! ＷＴ
      &    + wz_PresTendPs   & ! Ｇπ = (Ｃp κ Ｔ)π
      &   ) / R0**2

    !----- A = B + 2Δt * T-----
    wz_DivA = wa_xya( xyz_DivB ) + 2.0d0*DelTime * ( wz_DDivDtN )
!!!         & + 2.0d0*DelTime * wa_NumVis_wa(  wa_xya( xyz_DivB )  )

    !----- xyz_DivA に代入-----
    xyz_DivA = xya_wa( wz_DivA )

!--
!!$    call DataDump('xyz_DivB', xyz_DivB, strlen=80)
!!$    call DataDump('xyz_Div', xyz_Div, strlen=80)
!!$    call DataDump('xyz_Div_T', xya_wa(wz_DDivDtN), strlen=80)
!!$    call DataDump('xyz_DivA', xyz_DivA, strlen=80)
!++


    !-------------------------------------------------------------------
    !   Ｈ＝水平発散 Ｔ'Ｄ + σ移流 + π変化の効果  (熱力学の式)
    !-------------------------------------------------------------------

    ! 温度の局所時間変化項  Ｈ の計算

    ! Ｔ' Ｄ
    ! κＴv (ｖ・∇π)  : πの移流の効果
    ! - (α/Δσ) {Ｔv Σ[k〜K] (ｖ・∇π) + Ｔv' Σ[k〜K] (ＤΔσ)}

    do k = 1, km
      xyz_DTempLocalDtN(:,:,k) = &
        
        & xyz_TempEdd(:,:,k) * xyz_DivN(:,:,k) &
        
        & + z_TempInpolKappa(k) * xyz_TempVir(:,:,k)  &
        &    * xyz_lnPsAdv(:,:,k)                   &
        
        & - z_HydroAlpha(k) / z_DelSigma(k)                 &
        &    * (                                            &
        &       xyz_TempVir(:,:,k) * xyz_lnPsAdvSum(:,:,k)  &
        &       + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) &
        &      )
    enddo

    ! - (dσ/dt) * (１/Δσ) * (Ｔ'[k-1/2] − Ｔ'[k])
    ! - (dσ/dt)^NG * (１/Δσ) * (Ｔ￣[k-1/2] − Ｔ￣[k])

    do k = 2, km
      xyz_DTempLocalDtN(:,:,k) = xyz_DTempLocalDtN(:,:,k) &
        
        & - xyr_VSigma(:,:,k) / z_DelSigma(k)               &
        &     * ( xyr_TempEdd(:,:,k) - xyz_TempEdd(:,:,k) ) &
        
        & - xyr_VSigmaNonG(:,:,k) / z_DelSigma(k)           &
        &     * ( r_TempAvrXY(k) - z_TempAvrXY(k) )
    enddo

    ! - (dσ/dt) * (１/Δσ) * (Ｔ'[k] − Ｔ'[k+1/2])
    ! - (dσ/dt)^NG * (１/Δσ) * (Ｔ￣[k] − Ｔ￣[k+1/2])
    ! - (β/Δσ) {   Ｔv  Σ[k+1〜K] (ｖ・∇π)
    !               + Ｔv' Σ[k+1〜K] (ＤΔσ)}

    do k = 1, km - 1
      xyz_DTempLocalDtN(:,:,k) = xyz_DTempLocalDtN(:,:,k) &
        
        & - xyr_VSigma(:,:,k+1) / z_DelSigma(k)               &
        &     * ( xyz_TempEdd(:,:,k) - xyr_TempEdd(:,:,k+1) ) &
        
        & - xyr_VSigmaNonG(:,:,k+1) / z_DelSigma(k)           &
        &     * ( z_TempAvrXY(k) - r_TempAvrXY(k+1) )             &
        
        & - z_HydroBeta(k) / z_DelSigma(k)                        &
        &    * (                                                  &
        &       xyz_TempVir(:,:,k) * xyz_lnPsAdvSum(:,:,k+1)      &
        &       + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k+1)     &
        &      )
    enddo

!--
!!$    call io_gt4_out_Put('xyz_DTempLocalDtN_lnPsAdvSum',             &
!!$            & - z_HydroAlpha(2) / z_DelSigma(2)                   &
!!$            &    * (                                              &
!!$            &       xyz_TempVir(:,:,2) * xyz_lnPsAdvSum(:,:,2)    &
!!$            &      )                                              &
!!$            & - z_HydroBeta(2) / z_DelSigma(2)                    &
!!$            &    * (                                              &
!!$            &       xyz_TempVir(:,:,2) * xyz_lnPsAdvSum(:,:,2+1)  &
!!$            &      )                                              &
!!$            & )
!++


!--
!!$    call DataDump('xyz_DivSum*TempVirEdd*(Alpha+Beta)', &
!!$            & - z_HydroAlpha(k) / z_DelSigma(k)                    &
!!$            &    * ( xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) )   &
!!$            & - z_HydroBeta(k) / z_DelSigma(k)                     &
!!$            &    * ( xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k+1) ) &
!!$            & , strlen=70)
!++

!--
!!$    call DataDump('xyz_DTempLocalDtN', xyz_DTempLocalDtN, strlen=80)
!++
    call io_gt4_out_Put('xyz_DTempLocalDtN',  xyz_DTempLocalDtN)


!--
    !-------------------------------------------------------------------
    !   ＵＴ' 、ＶＴ'   (熱力学の式)
    !-------------------------------------------------------------------
!!$      DO 5100 K = 1, KMAX
!!$         DO 5100 IJ = 1, IDIM*JDIM
!!$            GTUT( IJ,K ) =  GAU ( IJ,K ) * GATED ( IJ,K )
!!$     &                                   * UVFACT ( IJ )
!!$            GTVT( IJ,K ) =  GAV ( IJ,K ) * GATED ( IJ,K )
!!$     &                                   * UVFACT ( IJ )
!!$ 5100 CONTINUE
!++

    !-------------------------------------------------------------------
    !   温度 Ｔ の計算。熱力学の式を解く
    !-------------------------------------------------------------------

    ! 非重力波項成分 (非線形項) の部分に関する温度変化を計算
    !
    ! dＵＴ'/dλ + dＶＴ'/dμ
    ! 温度の局所時間変化項 Ｈ

    do k = 1, km
      wz_DTempDtN(:,k) = &

        & - w_Div_xy_xy( &
        &               xyz_VelLonN(:,:,k) * xyz_TempEdd(:,:,k), &
        &               xyz_VelLatN(:,:,k) * xyz_TempEdd(:,:,k)  &
        &              ) / R0                                 &

        & + w_xy( xyz_DTempLocalDtN(:,:,k) )
      !&
      ! 非断熱加熱 Ｑ もまだ入ってないよ
      !&
      ! 摩擦熱および熱拡散は後回し
    enddo

!--
!!$    call DataDump('xyz_Temp_T_NonGrav', xya_wa(wz_DTempDtN), strlen=70)
!++


    ! 重力波項成分 (線形項) の部分に関する温度変化を計算
    !
    ! ｈＤ (重力波項成分によるπ変化の効果)

    do k = 1, km
      do kk = 1, km
        wz_DTempDtN(:,k) = wz_DTempDtN(:,k) &

          & - zz_siMtxH(k,kk) * w_xy( xyz_DivN(:,:,kk) )
        !&
        ! 拡散項は後回し
       enddo
    enddo


    !----- A = B + 2Δt * T-----
    wz_TempA = wa_xya( xyz_TempB ) + 2.0d0*DelTime * ( wz_DTempDtN )
!!!         & + 2.0d0*DelTime * wa_NumVisScaler_wa(  wa_xya( xyz_TempB )  )

    !----- xyz_TempA に代入-----
    xyz_TempA = xya_wa( wz_TempA )


    call io_gt4_out_Put(  'xyz_Temp_T', xya_wa( wz_DTempDtN )  )
!--
!!$    call DataDump('xyz_TempB', xyz_TempB, strlen=80)
!!$    call DataDump('xyz_Temp', xyz_Temp, strlen=80)
!!$    call DataDump('xyz_Temp_T', xya_wa(wz_DTempDtN), strlen=80)
!!$    call DataDump('xyz_TempA', xyz_TempA, strlen=80)
!++

    !-------------------------------------------------------------------
    !   比湿 QVap の計算。水蒸気の式を解く。
    !-------------------------------------------------------------------

    ! Ｒ＝ｑＤ＋σ移流 (非重力波項) を格子点上で計算
    xyz_QVapDivVSigmaAdv = xyz_QVapN * xyz_DivN

    do k = 2, km
      xyz_QVapDivVSigmaAdv(:,:,k) = xyz_QVapDivVSigmaAdv(:,:,k) &
        & - xyr_VSigma(:,:,k) / ( 2.0 * z_DelSigma(k) )  &
        &      * ( xyz_QVapN(:,:,k-1)  - xyz_QVapN(:,:,k) )
    enddo

    do k = 1, km - 1
      xyz_QVapDivVSigmaAdv(:,:,k) = xyz_QVapDivVSigmaAdv(:,:,k)  &
        & - xyr_VSigma(:,:,k+1) / ( 2.0 * z_DelSigma(k) ) &
        &      * ( xyz_QVapN(:,:,k) - xyz_QVapN(:,:,k+1) )
    enddo

    ! 水平移流 dＵｑ/dλ + dＶｑ/dμ と水平拡散項、水蒸気ソースを計算
    ! dＵｑ/dλ + dＶｑ/dμ
    ! Ｒ＝ｑＤ＋σ移流

    do k = 1, km
      wz_DQVapDtN(:,k) = &

        & - w_Div_xy_xy( &
        &               xyz_VelLonN(:,:,k) * xyz_QVapN(:,:,k), &
        &               xyz_VelLatN(:,:,k) * xyz_QVapN(:,:,k)  &
        &              ) / R0                              &

        & + w_xy( xyz_QVapDivVSigmaAdv(:,:,k) )
        !&
        ! 水平拡散項、水蒸気ソースは後回し
    enddo

    !----- A = B + 2Δt * T-----
    wz_QVapA = wa_xya( xyz_QVapB ) + 2.0d0*DelTime * ( wz_DQVapDtN )
!!!         & + 2.0d0*DelTime * wa_NumVisScaler_wa(  wa_xya( xyz_QVapB )  )

    !----- xyz_QVapA に代入-----
    xyz_QVapA = xya_wa( wz_QVapA )

!--
!!$    call DataDump('xyz_QVapB', xyz_QVapB, strlen=80)
!!$    call DataDump('xyz_QVap', xyz_QVap, strlen=80)
!!$    call DataDump('xyz_QVap_T', xya_wa(wz_DQVapDtN), strlen=80)
!!$    call DataDump('xyz_QVapA', xyz_QVapA, strlen=80)
!++

    call EndSub(subname)
  end subroutine dynamics_leapfrog

  subroutine dynamics_diffusion( &
    & xyz_VorB , xyz_DivB , xyz_TempB , xyz_QVapB  , &
    & xyz_VorA , xyz_DivA , xyz_TempA , xyz_QVapA        )

    ! Calculate Diffusion Term
    !
    ! t-Δt の値から水平拡散項を求め、それを t+Δt の値に加える。
    !
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use time_mod,    only: DelTime, CurrentTime
    use grid_3d_mod, only: km
    use grid_wavenumber_mod, only: nm
    use spml_mod,    only: wa_xya, xya_wa, l_nm
    use io_gt4_out_mod,only: io_gt4_out_Put
    use dc_trace,    only: DbgMessage, BeginSub, EndSub, DataDump
    implicit none

    real(DBKIND), intent(in) :: xyz_VorB(:,:,:)   ! 渦度 (t-Δt)
    real(DBKIND), intent(in) :: xyz_DivB(:,:,:)   ! 発散 (t-Δt)
    real(DBKIND), intent(in) :: xyz_TempB(:,:,:)  ! 温度 (t-Δt)
    real(DBKIND), intent(in) :: xyz_QVapB(:,:,:)  ! 比湿 (t-Δt)

    real(DBKIND), intent(inout) :: xyz_VorA(:,:,:)   ! 渦度 (t+Δt)
    real(DBKIND), intent(inout) :: xyz_DivA(:,:,:)   ! 発散 (t+Δt)
    real(DBKIND), intent(inout) :: xyz_TempA(:,:,:)  ! 温度 (t+Δt)
    real(DBKIND), intent(inout) :: xyz_QVapA(:,:,:)  ! 比湿 (t+Δt)

    real(DBKIND)    :: wz_VorA((nm+1)**2, km)   ! 渦度 (t+Δt)
    real(DBKIND)    :: wz_DivA((nm+1)**2, km)   ! 発散 (t+Δt)
    real(DBKIND)    :: wz_TempA((nm+1)**2, km)  ! 温度 (t+Δt)
    real(DBKIND)    :: wz_QVapA((nm+1)**2, km)  ! 比湿 (t+Δt)
    integer(INTKIND) :: i
    character(STRING),  parameter:: subname = "dynamics_diffusion"
  continue
    !----------------------------------------------------------------
    !   Check Initialization
    !----------------------------------------------------------------
    call BeginSub(subname, version)
    if (.not. dynamics_initialized) then
      call EndSub( subname, 'Call dynamics_init before call %c',  &
        &       c1=trim(subname) )
      return
    endif

    !-------------------------------------------------------------------
    !   それぞれの量の拡散項を計算
    !-------------------------------------------------------------------
    xyz_VorA = &
      & xya_wa( wa_xya(xyz_VorA) &
      &         + 2.0d0*DelTime * wz_DiffVorDiv * wa_xya( xyz_VorB ) )

!--
!!$    call io_gt4_out_Put('xyz_Vor_Diffusion', &
!!$         & xya_wa(   2.0d0*DelTime                     &
!!$         &             * wz_DiffVorDiv * wa_xya( xyz_VorB )  )    )
!++

    xyz_DivA = &
      & xya_wa( wa_xya(xyz_DivA) &
      &         + 2.0d0*DelTime * wz_DiffVorDiv * wa_xya( xyz_DivB ) )

!--
!!$    call io_gt4_out_Put('xyz_Div_Diffusion', &
!!$         & xya_wa(   2.0d0*DelTime                     &
!!$         &             * wz_DiffVorDiv * wa_xya( xyz_DivB )   )    )
!++

    xyz_TempA = &
      & xya_wa( wa_xya(xyz_TempA) &
      &         + 2.0d0*DelTime * wz_DiffTherm * wa_xya( xyz_TempB ) )

!--
!!$    call io_gt4_out_Put('xyz_Temp_Diffusion', &
!!$         & xya_wa(   2.0d0*DelTime                     &
!!$         &             * wz_DiffTherm * wa_xya( xyz_TempB )   )    )
!++

    xyz_QVapA  = &
      & xya_wa( wa_xya(xyz_QVapA) &
      &         + 2.0d0*DelTime * wz_DiffTherm * wa_xya( xyz_QVapB ) )

!--
!!$    call io_gt4_out_Put('xyz_QVap_Diffusion', &
!!$         & xya_wa(   2.0d0*DelTime                     &
!!$         &             * wz_DiffTherm * wa_xya( xyz_QVapB )   )   )
!++


    ! スペクトルデータとして見て拡散が効いているか確認するための
    ! データ出力を記述したソース。本計算にはまったく必要でない。
    !!wz_VorA  = wa_xya(xyz_VorA )
    !!wz_DivA  = wa_xya(xyz_DivA )
    !!wz_TempA = wa_xya(xyz_TempA)
    !!wz_QVapA = wa_xya(xyz_QVapA)
    !!
    !!write(80, *) CurrentTime, &
    !!     &       wz_VorA( lNm(21,0), 1) , wz_DivA( lNm(21,0), 1) , &
    !!     &       wz_TempA( lNm(21,0), 1) , wz_QVapA( lNm(21,0), 1)
    !!call flush(80)

    call EndSub(subname)
  end subroutine dynamics_diffusion


  subroutine dynamics_diagnostic &
    & ( x_Lon      , y_Lat      , z_Sigma  , r_Sigma  , &         !(in)

    &   xyz_VelLonA, xyz_VelLatA, &                               !(out)
    &   xyz_VorA   , xyz_DivA   , xyz_TempA, xyz_QVapA, xy_PsA  ) !(in)
    !
    ! Calculate Diagnostic Values.
    !
    ! 診断的に得られる量の演算を行なう。
    ! 現在は渦度発散から速度成分 (経度方向、緯度方向) の演算を行なうのみである。
    !

    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use grid_3d_mod,         only: im, jm, km
    use grid_wavenumber_mod, only: nm
    use constants_mod, only: R0, Grav
    use spml_mod,    only: wa_xya, xya_GradLon_wa, xya_GradLat_wa, &
      &                    wa_LaplaInv_wa, IntLonLat_xy
    use io_gt4_out_mod,only: io_gt4_out_Put
    use dc_trace,    only: DbgMessage, BeginSub, EndSub, DataDump
    use dc_string,   only: toChar
    implicit none

    
    real(DBKIND), intent(in):: x_Lon(:)        ! 経度座標
    real(DBKIND), intent(in):: y_Lat(:)        ! 緯度座標
    real(DBKIND), intent(in):: z_Sigma(:)      ! σレベル(整数)座標
    real(DBKIND), intent(in):: r_Sigma(:)      ! σレベル(半整数)座標

    real(DBKIND), intent(in):: xyz_VorA(:,:,:) ! 渦度         (t+Δt)
    real(DBKIND), intent(in):: xyz_DivA(:,:,:) ! 発散         (t+Δt)
    real(DBKIND), intent(in):: xyz_TempA(:,:,:)! 温度         (t+Δt)
    real(DBKIND), intent(in):: xyz_QVapA(:,:,:)! 比湿         (t+Δt)
    real(DBKIND), intent(in):: xy_PsA(:,:)     ! 地表面気圧   (t+Δt)

    real(DBKIND), intent(out) :: xyz_VelLonA(:,:,:)  ! 速度経度成分 (t+Δt)
    real(DBKIND), intent(out) :: xyz_VelLatA(:,:,:)  ! 速度緯度成分 (t+Δt)

    real(DBKIND) :: TotalMass      ! 全質量

    !----- 作業用内部変数 -----
    ! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
    integer(INTKIND)   :: i, j, k, l

    character(STRING),  parameter:: subname = "dynamics_diagnostic"
    !----------------------------------------------------------------
    !   Check Initialization
    !----------------------------------------------------------------
    call BeginSub(subname, version)
    if (.not. dynamics_initialized) then
      call EndSub( subname, 'Call dynamics_init before call %c',  &
        &          c1=trim(subname) )
      return
    endif

    !-------------------------------------------------------------------
    !   ｕ と ｖ を出力
    !-------------------------------------------------------------------
    wz_PsiA = wa_LaplaInv_wa(  wa_xya( xyz_VorA )  ) * R0**2
    wz_ChiA = wa_LaplaInv_wa(  wa_xya( xyz_DivA )  ) * R0**2

    xyz_VelLonA = &
      & ( xya_GradLon_wa( wz_ChiA ) - xya_GradLat_wa( wz_PsiA )  ) / R0

    xyz_VelLatA = &
      & ( xya_GradLon_wa( wz_PsiA ) + xya_GradLat_wa( wz_ChiA )  ) / R0

!--
!!$    call DataDump('xyz_VelLonA', xyz_VelLonA, strlen=80)
!!$    call DataDump('xyz_VelLatA', xyz_VelLatA, strlen=80)
!++

    !-------------------------------------------------------------------
    !   xy_Ps から全質量を求める。
    !-------------------------------------------------------------------
    TotalMass = IntLonLat_xy( xy_PsA / Grav )

    call io_gt4_out_Put('TotalMass', TotalMass)

    call EndSub(subname)
  end subroutine dynamics_diagnostic


  subroutine dynamics_end
    !
    ! Terminate module
    !
    ! dynamics_init で allocate した変数を deallocate し、
    ! 演算した値も全て破棄する。
    !

    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use dc_trace,    only: BeginSub, EndSub, DbgMessage
    implicit none

    !-----------------------------------------------------------------
    !   変数定義
    !-----------------------------------------------------------------
    !----- 作業用内部変数 -----
    character(STRING),  parameter:: subname = "dynamics_end"

  continue

    !-----------------------------------------------------------------
    !   Check Initialization
    !-----------------------------------------------------------------
    call BeginSub(subname, version)
    if ( .not. dynamics_initialized) then
      call EndSub( subname, 'dynamics_init was not called', &
        &       c1=trim(subname) )
      return
    else
      dynamics_initialized = .false.
    endif

    !-----------------------------------------------------------------
    !   allocate した変数の nullify とか
    !-----------------------------------------------------------------
    deallocate( xy_Coli       , &
      &      xy_UVFact        , &
      &      z_DelSigma       , &
      &      wz_DiffVorDiv    , &
      &      wz_DiffTherm     , &
      &      z_HydroAlpha     , &
      &      z_HydroBeta      , &
      &      z_TempInpolA     , &
      &      z_TempInpolB     , &
      &      z_TempInpolKappa , &
      &      z_TempAvrXY      , &
      &      r_TempAvrXY      , &
      &      z_siMtxG        , &
      &      zz_siMtxGCt     , &
      &      zz_siMtxW       , &
      &      zz_siMtxWH      , &
      &      zz_siMtxH       , &
      &      zz_siMtxQ      , &
      &      zz_siMtxS      , &
      &      zz_siMtxQS     , &
      &      zz_siMtxR      )


    deallocate  & 
      & ( xy_lnPs       , & ! π＝ ln Ｐs           (t)
      &   xyz_lnPsAdv   , & ! π の移流 ｖ・∇π
      &   xyz_lnPsAdvSum, & ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ
      &   xyz_DivSum    , & ! 発散の積下げ Σ[k〜K] ＤΔσ

      &   xy_DlnPsDtN     , & ! πの時間変化 dπ/dt   (Δt)
      &   w_DlnPsDtN      , & ! πの時間変化 (スペクトル) (Δt)
      &   w_lnPsA      , & ! π (t+Δt)

      &   xyr_VSigma     , & ! 鉛直σ速度(半整数レベル) (t)
      &   xyr_VSigmaNonG , & ! 鉛直σ速度 非重力波成分 (半整数レベル) (t)

      &   xyz_TempEdd    , &! Ｔ' ：温度の擾乱 (整数レベル) (t)
      &   xyz_TempVir    , &! Ｔv ：仮温度 (virtual temperature)(t)
      &   xyz_TempVirEdd , &! Ｔv'：仮温度の擾乱               (t)
      &   xyr_TempEdd, &! Ｔ' ：温度の擾乱 (半整数レベル)   (t)

      &   xyz_UAN      , & ! 東西運動量変化項ＵＡ
      &   xyz_VAN      , & ! 南北運動量変化項ＶＡ

      &   wz_DVorDtN      , & ! 渦度変化のスペクトルデータ (Δt) 
      &   wz_VorA      , & ! 渦度のスペクトルデータ    (t+Δt)

      &   xyz_KE        , & ! 運動エネルギー + 仮温度補正
                          ! ｕ**2+ｖ**2    + ΣＷ(Ｔv-Ｔ)

      &   wz_DDivDtN        , & ! 発散のスペクトルデータ (Δt)
      &   wz_DivA        , & ! 発散のスペクトルデータ (t+Δt)
      &   wz_PresTendTemp , & ! 温度による圧力傾度のスペクトルデータ
      &   wz_PresTendPs   , & ! 地表圧力による圧力傾度のスペクトルデータ

      &   xyz_DTempLocalDtN , & ! 温度の局所時間変化項  Ｈ
                              ! (水平発散 Ｔ'Ｄ + σ移流
                              !  + π変化の効果) (全て非重力波項)

      &   wz_DTempDtN       , & ! 温度の変化スペクトルデータ (Δt)
      &   wz_TempA       , & ! 温度のスペクトルデータ   (t+Δt)

      &   xyz_QVapDivVSigmaAdv, & ! Ｒ＝ｑＤ＋σ移流
      &   wz_DQVapDtN           , & ! 比湿の変化スペクトルデータ (Δt)
      &   wz_QVapA             & ! 比湿のスペクトルデータ   (t+Δt)
      & )

    deallocate  &
      & ( wz_PsiA , & ! スペクトル(流線関数)
      &   wz_ChiA   & ! スペクトル(ポテンシャル)
      & )

    call EndSub(subname)
  end subroutine dynamics_end


!--
!!  !=== Function wa_NumVis_wa : 数値粘性項の計算
!!  !
!!  !数値粘性項を計算する内部関数である。
!!  !
!!  function wa_NumVis_wa(wa_data)
!!  !==== Dependency
!!    use type_mod           , only: STRING, REKIND, DBKIND, INTKIND
!!    use constants_mod      , only: R0, VisOrder, EFoldTime
!!    use grid_3d_mod        , only: km
!!    use grid_wavenumber_mod, only: nm
!!    use dc_trace           , only: BeginSub, EndSub, DbgMessage, DataDump
!!    implicit none
!!    !==== Input
!!    !
!!    real(DBKIND), intent(in) :: wa_data((nm+1)*(nm+1), km)
!!    !
!!    !==== Output
!!    !
!!    real(DBKIND)             :: wa_NumVis_wa((nm+1)*(nm+1), km)
!!
!!    real(DBKIND)      :: VisCoef   ! 超粘性係数
!!    character(STRING),  parameter:: subname = "wa_NumVis_wa"
!!  continue
!!    call BeginSub(subname)
!!
!!    VisCoef = (   real(  ( nm*(nm+1) ), DBKIND  )  &
!!         &          / R0**2   )**(-VisOrder/2) &
!!         &     / EFoldTime
!!
!!!!$    call DbgMessage('VisCoef=<%f>', d=(/VisCoef/))
!!
!!    wa_NumVis_wa = &
!!         & wa_NumVisScaler_wa(wa_data)  &
!!         &   - VisCoef                  &
!!         &      * (  - (2.0d0/R0**2)**(VisOrder/2)  ) * wa_data
!!
!!!!$    call DbgMessage('-VisCoef*(2/R0**2)**(VisOrder/2)=<%f>', &
!!!!$         & d=(/ - VisCoef * (  - (2.0d0/R0**2)**(VisOrder/2)  ) /)   )
!!
!!
!!!!!$  hdc    = ( dble( (hdmaxn*(hdmaxn+1) ) ) / pradisq )**(-hdord/2) &
!!!!!$       / hdts
!!!!!$
!!!!!$  hdord   : 水平拡散の次数（ナブラの hdord 乗、ラプラシアンの hdord/2 乗）
!!!!!$  hdmaxn  : 打ち切り波数（ T21 なら 21、T42 なら 42、... ）
!!!!!$  pradisq : 惑星半径の二乗 [m2]
!!!!!$  hdts    : 全波数 hdmaxn における減衰の時定数 [s]
!!
!!    call EndSub(subname)
!!  end function wa_NumVis_wa
!!
!!  !=== Function wa_NumVisScaler_wa : スカラーに対する数値粘性項の計算
!!  !
!!  !スカラーに対する数値粘性項を計算する内部関数である。
!!  !
!!  function wa_NumVisScaler_wa(wa_data)
!!    !==== Dependency
!!    use type_mod           , only: STRING, REKIND, DBKIND, INTKIND
!!    use constants_mod      , only: R0, VisOrder, EFoldTime
!!    use grid_3d_mod        , only: km
!!    use grid_wavenumber_mod, only: nm
!!    use spml_mod           , only: rn
!!    use dc_trace           , only: BeginSub, EndSub, DbgMessage, DataDump
!!    implicit none
!!    !==== Input
!!    !
!!    real(DBKIND), intent(in) :: wa_data((nm+1)*(nm+1), km)
!!    !
!!    !==== Output
!!    !
!!    real(DBKIND)             :: wa_NumVisScaler_wa((nm+1)*(nm+1), km)
!!
!!    real(DBKIND)             :: VisCoef   ! 超粘性係数
!!    real(DBKIND)             :: wa_rn((nm+1)*(nm+1), km)
!!    integer(INTKIND)         :: k
!!    character(STRING),  parameter:: subname = "wa_NumVisScaler_wa"
!!    continue
!!    call BeginSub(subname)
!!
!!    do k = 1, km
!!       wa_rn(:,k) = rn(:,1)
!!    enddo
!!
!!    VisCoef = (   real(  ( nm*(nm+1) ), DBKIND  )  &
!!         &          / R0**2   )**(-VisOrder/2) &
!!         &     / EFoldTime
!!
!!!!$    call DbgMessage('VisCoef=<%f>', d=(/VisCoef/))
!!
!!    wa_NumVisScaler_wa  =   &
!!         & - VisCoef * ( ( - wa_rn / R0**2 )**(VisOrder/2) ) * wa_data
!!
!!    call EndSub(subname)
!!  end function wa_NumVisScaler_wa
!++

end module dynamics_mod
 
