!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Module dynamics_mod
!
!   * Developers: Morikawa Yasuhiro
!   * Version: $Id: dynamics.f90,v 1.21 2005/05/02 11:23:41 morikawa Exp $
!   * Tag Name: $Name: dcpam2-20050522-1 $
!   * Change History: 
!
!== 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
!
!現在、力学コアの全てがこのモジュール内にある。本来ならば、
!解く方程式や項の種類によってモジュール化がなされるべきである。
!
                                                                 !=end
module dynamics_mod
                                                                 !=begin
  !== Dependency
  use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
                                                                 !=end
  implicit none
                                                                 !=begin
  !== Public Interface
  private
  public :: dynamics_init, dynamics_leapfrog  ! subroutines
  public :: dynamics_diagnostic, dynamics_end ! subroutines
  public :: dynamics_diffusion                ! subroutines
                                                                 !=end

  !-------------------------------------------------------------------
  !   dynamics_init において値を設定するもの
  !-------------------------------------------------------------------
  real(DBKIND) :: Kappa               ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱)
  real(DBKIND), pointer, save :: &

    & xy_Coli(:,:)       =>null() , & ! 格子点データ(コリオリ力)
    & xy_UVFact(:,:)     =>null() , & ! u→U, v→V のファクター
    & z_DelSigma(:)      =>null() , & ! Δσ(整数レベル)

    & wz_DiffVorDiv(:,:) =>null() , & ! 運動量 (渦度発散) 水平拡散係数
    & wz_DiffTherm(:,:)  =>null() , & ! 熱、比湿 水平拡散係数

    & z_TempAve(:)       =>null() , & ! 平均温度 (整数レベル)
    & z_TempAveHalf(:)   =>null() , & ! 平均温度 (半整数レベル)

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

    ! semi-implicit用行列
    & zz_DivMtx(:,:)   =>null() , & ! Ｗ:発散の式での温度補正 
                                    !    (線形重力波項の効果)
    & zz_TempMtx(:,:)  =>null() , & ! ｈ: ＱＳ (σ移流) - Ｒ
    & zz_TempMtxQ(:,:) =>null() , & ! Ｑ: dＴ/dσ (線形重力波項の効果)
    & zz_TempMtxS(:,:) =>null() , & ! Ｓ: dσ/dｔ (線形重力波項の効果)
    & zz_TempMtxR(:,:) =>null()     ! Ｒ: ＲＤ =
                                    !      κＴ（∂π/∂ｔ＋ｄσ/ｄｔ／σ）
                                    !     気圧変化の効果の係数
                                    !     (線形重力波項の効果)


  !-------------------------------------------------------------------
  !   dynamics_leapfrog において利用するもの
  !-------------------------------------------------------------------
  real(DBKIND), pointer, save :: &

    & xy_lnPs(:,:)   =>null() , & ! π＝ ln Ｐs           (t)
    & xy_lnPs_T(:,:) =>null() , & ! πの時間変化 dπ/dt   (Δt)
    & w_lnPs_T(:)    =>null() , & ! πの時間変化 (スペクトル) (Δt)
    & w_lnPs_a(:)    =>null() , & ! π (t+Δt)

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

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

    & xyz_UA_T(:,:,:) =>null() , & ! 東西運動量変化項ＵＡ
    & xyz_VA_T(:,:,:) =>null() , & ! 南北運動量変化項ＶＡ

    & wz_Vor_T(:,:)   =>null() , & ! 渦度変化のスペクトルデータ (Δt)
    & wz_Vor_a(:,:)   =>null() , & ! 渦度のスペクトルデータ    (t+Δt)

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

    & wz_Div_T(:,:) =>null() , & ! 発散のスペクトルデータ (Δt)
    & wz_Div_a(:,:) =>null() , & ! 発散のスペクトルデータ (t+Δt)
    & wz_PresTendTemp(:,:) =>null() , & ! 温度による圧力傾度のスペクトルデータ
    & wz_PresTendPs(:,:)   =>null() , & ! 地表圧力による圧力傾度のスペクトルデータ

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

    & wz_Temp_T(:,:) =>null() , & ! 温度の変化スペクトルデータ (Δt)
    & wz_Temp_a(:,:) =>null() , & ! 温度のスペクトルデータ   (t+Δt)

    & xyz_QVapDivVSigmaAdv(:,:,:) => null() , & ! Ｒ＝ｑＤ＋σ移流
    & wz_QVap_T(:,:) =>null() , & ! 比湿の変化スペクトルデータ (Δt)
    & wz_QVap_a(:,:) =>null()     ! 比湿のスペクトルデータ   (t+Δt)

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


  logical, save :: dynamics_initialized = .false.
  character(STRING),parameter:: version = &
       & '$Id: dynamics.f90,v 1.21 2005/05/02 11:23:41 morikawa Exp $'
  character(STRING),parameter:: tagname = '$Name: dcpam2-20050522-1 $'

contains
                                                                 !=begin
  !== Procedure Interface
  !
  !=== Initialize module and Calculate Non-Predictional Value.
  !
  !以降のサブルーチンで用いる変数の allocate 、および
  !時間発展しない量の演算を行なう。
  !また、変数データ出力のための初期設定も行なう。
  !
  subroutine dynamics_init(x_Lon, y_Lat, z_Sigma, r_Sigma)
  !==== Dependency
    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
                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & x_Lon(:)           , & ! intent(in): 経度座標
         & y_Lat(:)           , & ! intent(in): 緯度座標
         & z_Sigma(:)         , & ! intent(in): σレベル(整数)座標
         & r_Sigma(:)             ! intent(in): σレベル(半整数)座標
                                                                 !=end

    !----- 拡散係数計算用変数 -----
    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)
    if (dynamics_initialized) then
       call EndSub( subname, '%c is already called.', c1=trim(subname) )
       return
    else
       dynamics_initialized = .true.
    endif

    !----------------------------------------------------------------
    !   Version identifier
    !----------------------------------------------------------------
    call DbgMessage('%c :: %c', c1=trim(version), c2=trim(tagname))

    !-------------------------------------------------------------------
    !   物理定数の取得
    !-------------------------------------------------------------------
    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_TempAve(km) )
    allocate( z_TempAveHalf(km+1) )

    z_TempAve(:) = TempAve

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

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

!!$    call DataDump('z_TempAve', z_TempAve(:), strlen=80)
!!$    call DataDump('z_TempAveHalf', z_TempAveHalf(:), strlen=80)


    !-------------------------------------------------------------------
    !   semi-implicit 用行列の計算
    !-------------------------------------------------------------------
    allocate( zz_DivMtx(km,km) )
    allocate( zz_TempMtx(km,km) )
    allocate( zz_TempMtxQ(km,km) )
    allocate( zz_TempMtxS(km,km) )
    allocate( zz_TempMtxR(km,km) )

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

    ! 初期化
    zz_DivMtx(:,:) = 0.0

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

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

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

    ! 初期化
    zz_TempMtxS(:,:) = 0.0

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

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

    enddo
!!$    call DataDump('zz_TempMtxS', zz_TempMtxS(:,:), strlen=80)

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

    ! 初期化
    zz_TempMtxQ(:,:) = 0.0

    do k = 1, km
       zz_TempMtxQ(k,k) = ( z_TempAveHalf(k) - z_TempAve(k) ) / z_DelSigma(k)
    enddo
    do k = 1, km - 1
       zz_TempMtxQ(k,k+1) = ( z_TempAve(k) - z_TempAveHalf(k+1) ) / z_DelSigma(k)
    enddo
!!$    call DataDump('zz_TempMtxQ', zz_TempMtxQ(:,:), strlen=80)

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

    ! 初期化
    zz_TempMtxR(:,:) = 0.0

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

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

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

    ! 初期化
    zz_TempMtx(:,:) = 0.0

    zz_TempMtx(:,:) = &
         & matmul(zz_TempMtxQ(:,:), zz_TempMtxS(:,:)) - zz_TempMtxR(:,:)

!!$    call DataDump('zz_TempMtx', zz_TempMtx(:,:), 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_UA_T')
    call io_gt4_out_SetVars('xyz_VA_T')
    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('xyz_VSigmaHalf')
    call io_gt4_out_SetVars('xyz_VSigmaHalfNonG')
    call io_gt4_out_SetVars('xyz_TempLocal_T')
    call io_gt4_out_SetVars('xyz_TempAdv')
    call io_gt4_out_SetVars('xyz_TempEddHalf')
    call io_gt4_out_SetVars('xyz_Temp_T')
    call io_gt4_out_SetVars('xyz_TempLocal_T_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_lnPs_T( im, jm )       , & ! πの時間変化 dπ/dt   (Δt)
      &   w_lnPs_T( (nm+1)*(nm+1) ) , & ! πの時間変化 (スペクトル) (Δt)
      &   w_lnPs_a( (nm+1)*(nm+1) ) , & ! π (t+Δt)

      &   xyz_VSigmaHalf(im, jm, km+1) , & ! 鉛直σ速度(半整数レベル) (t)
      &   xyz_VSigmaHalfNonG(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)
      &   xyz_TempEddHalf(im, jm, km+1), &! Ｔ' ：温度の擾乱 (半整数レベル)   (t)

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

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

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

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

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

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

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


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


    call EndSub(subname)
  end subroutine dynamics_init


                                                                 !=begin
  !=== Calculate Predictional Values.
  !
  !時間発展する量の演算を行なう。
  !演算するデータの出力も行なう。
  !
  subroutine dynamics_leapfrog                                &
       & ( x_Lon       , y_Lat       , z_Sigma  , r_Sigma   , &

       &   xyz_VelLon_b, xyz_VelLat_b, xyz_Vor_b, xyz_Div_b , &
       &   xyz_Temp_b  , xyz_QVap_b  , xy_Ps_b  ,             &

       &   xyz_VelLon_n, xyz_VelLat_n, xyz_Vor_n, xyz_Div_n , &
       &   xyz_Temp_n  , xyz_QVap_n  , xy_Ps_n  ,             &

       &   xyz_VelLon_a, xyz_VelLat_a, xyz_Vor_a, xyz_Div_a , &
       &   xyz_Temp_a  , xyz_QVap_a  , xy_Ps_a             )

  !==== Dependency
    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
                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & x_Lon(:)          , & ! intent(in): 経度座標
         & y_Lat(:)          , & ! intent(in): 緯度座標
         & z_Sigma(:)        , & ! intent(in): σレベル(整数)座標
         & r_Sigma(:)        , & ! intent(in): σレベル(半整数)座標

         & xyz_VelLon_b(:,:,:) , & ! intent(in): 速度経度成分 (t-Δt)
         & xyz_VelLat_b(:,:,:) , & ! intent(in): 速度緯度成分 (t-Δt)
         & xyz_Vor_b(:,:,:)    , & ! intent(in): 渦度         (t-Δt)
         & xyz_Div_b(:,:,:)    , & ! intent(in): 発散         (t-Δt)
         & xyz_Temp_b(:,:,:)   , & ! intent(in): 温度         (t-Δt)
         & xyz_QVap_b(:,:,:)   , & ! intent(in): 比湿         (t-Δt)
         & xy_Ps_b(:,:)        , & ! intent(in): 地表面気圧   (t-Δt)

         & xyz_VelLon_n(:,:,:) , & ! intent(in): 速度経度成分 (t)
         & xyz_VelLat_n(:,:,:) , & ! intent(in): 速度緯度成分 (t)
         & xyz_Vor_n(:,:,:)    , & ! intent(in): 渦度         (t)
         & xyz_Div_n(:,:,:)    , & ! intent(in): 発散         (t)
         & xyz_Temp_n(:,:,:)   , & ! intent(in): 温度         (t)
         & xyz_QVap_n(:,:,:)   , & ! intent(in): 比湿         (t)
         & xy_Ps_n(:,:)            ! intent(in): 地表面気圧   (t)

    !==== Output
    !
    real(DBKIND), intent(out) :: &
         & xyz_VelLon_a(:,:,:) , & ! intent(out): 速度経度成分 (t+Δt)
         & xyz_VelLat_a(:,:,:) , & ! intent(out): 速度緯度成分 (t+Δt)
         & xyz_Vor_a(:,:,:)  , & ! intent(out): 渦度         (t+Δt)
         & xyz_Div_a(:,:,:)  , & ! intent(out): 発散         (t+Δt)
         & xyz_Temp_a(:,:,:) , & ! intent(out): 温度         (t+Δt)
         & xyz_QVap_a(:,:,:) , & ! intent(out): 比湿         (t+Δt)
         & xy_Ps_a(:,:)          ! intent(out): 地表面気圧   (t+Δt)
                                                                 !=end

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


    character(STRING),  parameter:: subname = "dynamics_leapfrog"

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

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

    xy_lnPs(:,:) = log( xy_Ps_n(:,:) )

    ! 地表面気圧 π の移流 = ｖ・∇π
    do k = 1, km
       xyz_lnPsAdv(:,:,k) = &
            &  (  &
            &   xyz_VelLon_n(:,:,k)                         &
            &      * xy_GradLon_w( w_xy(xy_lnPs(:,:)) ) &

            &   + xyz_VelLat_n(:,:,k)                       &
            &      * xy_GradLat_w( w_xy(xy_lnPs(:,:)) ) &

            &  ) / 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_Div_n(:,:, km ) * z_DelSigma( km )

    do k = km -1 , 1, -1
       xyz_DivSum(:,:,k) = &
            &   xyz_DivSum(:,:,k+1)           &
            & + xyz_Div_n(:,:,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_lnPs_T(:,:) =                 &
         & - xyz_lnPsAdvSum(:,:,1)   ! - Σ[1〜K](ｖ・∇π)Δσ

!!$    call DataDump('xy_lnPs_T', xy_lnPs_T(:,:), strlen=80)


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

    !----- A = B + 2Δt * T-----
    w_lnPs_a(:) = &
         & w_xy(  log( xy_Ps_b(:,:) )  ) &
         & + 2.0*DelTime * ( w_lnPs_T(:) )

!!$    call DataDump('xy_Ps_b', xy_Ps_b(:,:), strlen=80)
!!$    call DataDump('w_lnPs_b', w_xy(  log( xy_Ps_b(:,:) )  ), strlen=80)
!!$    call DataDump('w_lnPs_a', w_lnPs_a(:), strlen=80)

    !----- xy_Ps_a に代入-----
    xy_Ps_a(:,:) = exp(  xy_w( w_lnPs_a(:) )  )

!!$    call DataDump('xy_Ps_a', xy_Ps_a(:,:), strlen=80)

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

    ! σ速度 (3/2 〜 K-1/2)
    do k = 2, km+1 - 1
       xyz_VSigmaHalf(:,:,k) =  &

            ! - σ[k-1/2] × ( Σ[1〜K](Ｄ + ｖ・∇π)Δσ )
            & r_Sigma(k) &
            &     * ( xyz_lnPsAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) &

            ! - Σ[k〜K] (Ｄ + ｖ・∇π) Δσ
            & - ( xyz_lnPsAdvSum(:,:,k) + xyz_DivSum(:,:,k) )

       xyz_VSigmaHalfNonG(:,:,k) =  &

            ! - σ[k-1/2] × ( Σ[1〜K](ｖ・∇π)Δσ )
            & r_Sigma(k) * xyz_lnPsAdvSum(:,:,1) &

            ! - Σ[k〜K] (ｖ・∇π) Δσ
            & - xyz_lnPsAdvSum(:,:,k)
    enddo

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

    call io_gt4_out_Put('xyz_VSigmaHalf', xyz_VSigmaHalf(:,:,:))
    call io_gt4_out_Put('xyz_VSigmaHalfNonG', xyz_VSigmaHalfNonG(:,:,:))
!!$    call DataDump('xyz_VSigmaHalf', xyz_VSigmaHalf(:,:,:), strlen=80)

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

       xyz_TempEdd(:,:,k)   = xyz_Temp_n(:,:,k) - z_TempAve(k)
       xyz_TempVirEdd(:,:,k) = xyz_TempVir(:,:,k) - z_TempAve(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
       xyz_TempEddHalf(:,:,k) = &
            &   z_TempInpolA(k)   * xyz_Temp_n(:,:,k)   &
            & + z_TempInpolB(k-1) * xyz_Temp_n(:,:,k-1) &
            & - z_TempAveHalf(k)
    enddo

    call io_gt4_out_Put('xyz_TempEddHalf', xyz_TempEddHalf(:,:,:))
!!$    call DataDump('xyz_TempEddHalf', xyz_TempEddHalf(:,:,:), strlen=80)

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

    do k = 1, km
       xyz_UA_T(:,:,k) = &
            ! (ζ + ｆ) ｖ
            & ( xyz_Vor_n(:,:,k) + xy_Coli(:,:) ) * xyz_VelLat_n(:,:,k) &

            ! - (Ｃp κ (1/ａ) Ｔ'v (dπ/dλ) ) / cosφ
            & - Cp * z_TempInpolKappa(k) &
            &           * xyz_TempVirEdd(:,:,k)  &
            &             * xy_GradLon_w(  w_xy( xy_lnPs(:,:) )  ) &
            &                / R0

       xyz_VA_T(:,:,k) = &
            ! - (ζ + ｆ) ｕ
            & - ( xyz_Vor_n(:,:,k) + xy_Coli(:,:) ) * xyz_VelLon_n(:,:,k) &

            ! - (Ｃp κ (1/ａ) Ｔ'v (dπ/dμ) ) / cosφ
            & - Cp * z_TempInpolKappa(k) &
            &           * xyz_TempVirEdd(:,:,k) &
            &             * xy_GradLat_w(  w_xy( xy_lnPs(:,:) )  ) &
            &               / 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_VA_T', xyz_VA_T(1,1,:), strlen=80, multi=(/-1/))
!!$    call DataDump('xyz_UA_T', xyz_UA_T(1,1,:), strlen=80, multi=(/-1/))
!!$    call DataDump('xyz_VA_T', xyz_VA_T(1,1,:), strlen=80, multi=(/-1/))

    do k = 2, km

       xyz_UA_T(:,:,k) = xyz_UA_T(:,:,k) &
            ! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (ｕ[k-1] - ｕ[k])
            & - 1.0 / (2.0 * z_DelSigma(k)) &
            &    * xyz_VSigmaHalf(:,:,k)  &
            &       * ( xyz_VelLon_n(:,:,k-1) - xyz_VelLon_n(:,:,k) )

       xyz_VA_T(:,:,k) = xyz_VA_T(:,:,k) &
            ! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (ｖ[k-1] - ｖ[k])
            & - 1.0 / (2.0 * z_DelSigma(k)) &
            &    * xyz_VSigmaHalf(:,:,k)  &
            &       * ( xyz_VelLat_n(:,:,k-1) - xyz_VelLat_n(:,:,k) )
    enddo

!!$    call DataDump('xyz_UA_T', xyz_UA_T(1,1,:), strlen=80, multi=(/-2/))
!!$    call DataDump('xyz_VA_T', xyz_VA_T(1,1,:), strlen=80, multi=(/-2/))

    do k = 1, km - 1

       xyz_UA_T(:,:,k) = xyz_UA_T(:,:,k) &
            ! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (ｕ[k] - ｕ[k+1])
            & - 1.0 / (2.0 * z_DelSigma(k)) &
            &    * xyz_VSigmaHalf(:,:,k+1) &
            &       * ( xyz_VelLon_n(:,:,k) - xyz_VelLon_n(:,:,k+1) )

       xyz_VA_T(:,:,k) = xyz_VA_T(:,:,k) &
            ! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (ｖ[k] - ｖ[k+1])
            & - 1.0 / (2.0 * z_DelSigma(k)) &
            &    * xyz_VSigmaHalf(:,:,k+1)  &
            &       * ( xyz_VelLat_n(:,:,k) - xyz_VelLat_n(:,:,k+1) )
    enddo

!!$    call DataDump('xyz_UA_T', xyz_UA_T(1,1,:), strlen=80, multi=(/-3/))
!!$    call DataDump('xyz_VA_T', xyz_VA_T(1,1,:), strlen=80, multi=(/-3/))

!!$    call DataDump('xyz_UA_T', xyz_UA_T(:,:,:), strlen=80)
!!$    call DataDump('xyz_VA_T', xyz_VA_T(:,:,:), strlen=80)

    call io_gt4_out_Put('xyz_UA_T', xyz_UA_T(:,:,:))
    call io_gt4_out_Put('xyz_VA_T', xyz_VA_T(:,:,:))

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

    ! dζ/dt の計算
    wz_Vor_T(:,:) = &
         & wa_Div_xya_xya( xyz_VA_T(:,:,:), - xyz_UA_T(:,:,:) ) / R0
         !& + w_DiffV(k) * wa_xya(xyz_Vor(:,:,k))
         !  拡散係数の導出が厄介なので後回し

    !----- A = B + 2Δt * T-----
    wz_Vor_a(:,:) = &
         & wa_xya( xyz_Vor_b(:,:,:) )              &
         & + 2.0d0*DelTime * ( wz_Vor_T(:,:) )
!!!         & + 2.0d0*DelTime * wa_NumVis_wa(  wa_xya( xyz_Vor_b(:,:,:) )  )


    !----- xyz_Vor_a に代入-----
    xyz_Vor_a(:,:,:) = xya_wa( wz_Vor_a(:,:) )

!!$    call DataDump('wz_Vor_T', wz_Vor_T(:,:), strlen=70)
!!$    call DataDump('xyz_Vor_b', xyz_Vor_b(:,:,:), strlen=80)
!!$    call DataDump('xyz_Vor', xyz_Vor(:,:,:), strlen=80)
!!$    call DataDump('xyz_Vor_T', xya_wa(wz_Vor_T(:,:)), strlen=80)
!!$    call DataDump('xyz_Vor_a', xyz_Vor_a(:,:,:), strlen=80)

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

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

    ! 初期化
    xyz_KE(:,:,:) = 0.0

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

    do k = 2, km
       xyz_KE(:,:,k) = xyz_KE(:,:,k-1)   &
            & + Cp * z_HydroAlpha(k) &
            &     * ( xyz_TempVir(:,:,k) - xyz_Temp_n(:,:,k) ) &
            & + Cp * z_HydroBeta(k-1) &
            &     * ( xyz_TempVir(:,:,k-1) - xyz_Temp_n(:,:,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_VelLon_n(:,:,:)**2   &
         &     + xyz_VelLat_n(:,:,:)**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_PresTendTemp(:,1) = Cp * z_HydroAlpha(1) * w_xy( xyz_Temp_n(:,:,1) )

    do k = 2, km
       wz_PresTendTemp(:,k) = wz_PresTendTemp(:,k-1) &
            & + Cp * z_HydroAlpha(k)  * w_xy( xyz_Temp_n(:,:,k) )   &
            & + Cp * z_HydroBeta(k-1) * w_xy( xyz_Temp_n(:,:,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_TempAve(k) &
            &   * w_xy(  log( xy_Ps_n(:,:) )  )
    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_UA_T - wz_VA_T)', &
!!$         & wa_Div_xya_xya( xyz_UA_T(:,:,:), - xyz_VA_T(:,:,:) ), &
!!$         & strlen=70)
!!$    call DataDump('Div(wz_UA_T - wz_VA_T)/R0', &
!!$         & wa_Div_xya_xya( xyz_UA_T(:,:,:), - xyz_VA_T(:,:,:) ) /R0, &
!!$         & strlen=70)
!!$    call DataDump('xyz_Div(wz_UA_T - wz_VA_T)/R0', &
!!$         & xya_wa(  wa_Div_xya_xya( xyz_UA_T(:,:,:), - xyz_VA_T(:,:,:) )  ) /R0, &
!!$         & strlen=70)
!!$    call DataDump('Lapla(wz_KE)', &
!!$         & wa_Lapla_wa( wa_xya(xyz_KE(:,:,:)) ), strlen=70)

    ! dＤ/dt の計算
    wz_Div_T(:,:) = &
         & wa_Div_xya_xya( xyz_UA_T(:,:,:), xyz_VA_T(:,:,:) ) / R0 &
         & - wa_Lapla_wa( wa_xya(xyz_KE(:,:,:)) ) / R0**2 &
         !& + w_DiffV(k) * wa_xya(xyz_Div(:,:,k))
         !  拡散係数の導出が厄介なので後回し

         ! semi-implicitだとここには無い項達
         !  - ∇**2 (Φ + ＷＴ + Ｇπ)
         & - wa_Lapla_wa( &
         !&    w_xy( xy_Phi(:,:) )  &  ! 地表ジオポテンシャルΦは後回し
         &    wz_PresTendTemp(:,:)   & ! ＷＴ
         &    + wz_PresTendPs(:,:)   & ! Ｇπ = (Ｃp κ Ｔ)π
         &   ) / R0**2

    !----- A = B + 2Δt * T-----
    wz_Div_a(:,:) = &
         & wa_xya( xyz_Div_b(:,:,:) )          &
         & + 2.0d0*DelTime * ( wz_Div_T(:,:) )
!!!         & + 2.0d0*DelTime * wa_NumVis_wa(  wa_xya( xyz_Div_b(:,:,:) )  )

    !----- xyz_Div_a に代入-----
    xyz_Div_a(:,:,:) = xya_wa( wz_Div_a(:,:) )

!!$    call DataDump('xyz_Div_b', xyz_Div_b(:,:,:), strlen=80)
!!$    call DataDump('xyz_Div', xyz_Div(:,:,:), strlen=80)
!!$    call DataDump('xyz_Div_T', xya_wa(wz_Div_T(:,:)), strlen=80)
!!$    call DataDump('xyz_Div_a', xyz_Div_a(:,:,:), strlen=80)


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

    ! 温度の局所時間変化項  Ｈ の計算
    do k = 1, km
       xyz_TempLocal_T(:,:,k) = &

            ! Ｔ' Ｄ
            & xyz_TempEdd(:,:,k) * xyz_Div_n(:,:,k) &

            ! κＴv (ｖ・∇π)  : πの移流の効果
            & + z_TempInpolKappa(k) * xyz_TempVir(:,:,k)  &
            &    * xyz_lnPsAdv(:,:,k)                   &

            ! - (α/Δσ) {Ｔv Σ[k〜K] (ｖ・∇π) + Ｔv' Σ[k〜K] (ＤΔσ)}
            & - z_HydroAlpha(k) / z_DelSigma(k)                 &
            &    * (                                            &
            &       xyz_TempVir(:,:,k) * xyz_lnPsAdvSum(:,:,k)  &
            &       + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) &
            &      )
    enddo

    do k = 2, km
       xyz_TempLocal_T(:,:,k) = xyz_TempLocal_T(:,:,k) &

            ! - (dσ/dt) * (１/Δσ) * (Ｔ'[k-1/2] − Ｔ'[k])
            & - xyz_VSigmaHalf(:,:,k) / z_DelSigma(k)               &
            &     * ( xyz_TempEddHalf(:,:,k) - xyz_TempEdd(:,:,k) ) &

            ! - (dσ/dt)^NG * (１/Δσ) * (Ｔ￣[k-1/2] − Ｔ￣[k])
            & - xyz_VSigmaHalfNonG(:,:,k) / z_DelSigma(k)           &
            &     * ( z_TempAveHalf(k) - z_TempAve(k) )
    enddo

    do k = 1, km - 1
       xyz_TempLocal_T(:,:,k) = xyz_TempLocal_T(:,:,k) &

            ! - (dσ/dt) * (１/Δσ) * (Ｔ'[k] − Ｔ'[k+1/2])
            & - xyz_VSigmaHalf(:,:,k+1) / z_DelSigma(k)               &
            &     * ( xyz_TempEdd(:,:,k) - xyz_TempEddHalf(:,:,k+1) ) &

            ! - (dσ/dt)^NG * (１/Δσ) * (Ｔ￣[k] − Ｔ￣[k+1/2])
            & - xyz_VSigmaHalfNonG(:,:,k+1) / z_DelSigma(k)           &
            &     * ( z_TempAve(k) - z_TempAveHalf(k+1) )             &

            ! - (β/Δσ) {   Ｔv  Σ[k+1〜K] (ｖ・∇π)
            !               + Ｔv' Σ[k+1〜K] (ＤΔσ)}
            & - 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_TempLocal_T_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_TempLocal_T', xyz_TempLocal_T(:,:,:), strlen=80)
    call io_gt4_out_Put('xyz_TempLocal_T',  xyz_TempLocal_T(:,:,:))


    !-------------------------------------------------------------------
    !   ＵＴ' 、ＶＴ'   (熱力学の式)
    !-------------------------------------------------------------------
!!$      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

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

    ! 非重力波項成分 (非線形項) の部分に関する温度変化を計算
    do k = 1, km
       wz_Temp_T(:,k) = &

            ! dＵＴ'/dλ + dＶＴ'/dμ
            & - w_Div_xy_xy( &
            &               xyz_VelLon_n(:,:,k) * xyz_TempEdd(:,:,k), &
            &               xyz_VelLat_n(:,:,k) * xyz_TempEdd(:,:,k)  &
            &              ) / R0                                 &

            ! 温度の局所時間変化項 Ｈ
            & + w_xy( xyz_TempLocal_T(:,:,k) )
            !&
            ! 非断熱加熱 Ｑ もまだ入ってないよ
            !&
            ! 摩擦熱および熱拡散は後回し
    enddo

!!$    call DataDump('xyz_Temp_T_NonGrav', xya_wa(wz_Temp_T(:,:)), strlen=70)


    ! 重力波項成分 (線形項) の部分に関する温度変化を計算
    do k = 1, km
       do kk = 1, km
          wz_Temp_T(:,k) = wz_Temp_T(:,k) &

               ! ｈＤ (重力波項成分によるπ変化の効果)
               & - zz_TempMtx(k,kk) * w_xy( xyz_Div_n(:,:,kk) )
               !&
               ! 拡散項は後回し
       enddo
    enddo


    !----- A = B + 2Δt * T-----
    wz_Temp_a(:,:) = &
         & wa_xya( xyz_Temp_b(:,:,:) )          &
         & + 2.0d0*DelTime * ( wz_Temp_T(:,:) )
!!!         & + 2.0d0*DelTime * wa_NumVisScaler_wa(  wa_xya( xyz_Temp_b(:,:,:) )  )

    !----- xyz_Temp_a に代入-----
    xyz_Temp_a(:,:,:) = xya_wa( wz_Temp_a(:,:) )


    call io_gt4_out_Put(  'xyz_Temp_T', xya_wa( wz_Temp_T(:,:) )  )
!!$    call DataDump('xyz_Temp_b', xyz_Temp_b(:,:,:), strlen=80)
!!$    call DataDump('xyz_Temp', xyz_Temp(:,:,:), strlen=80)
!!$    call DataDump('xyz_Temp_T', xya_wa(wz_Temp_T(:,:)), strlen=80)
!!$    call DataDump('xyz_Temp_a', xyz_Temp_a(:,:,:), strlen=80)

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

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

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

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

    ! 水平移流 dＵｑ/dλ + dＶｑ/dμ と水平拡散項、水蒸気ソースを計算
    do k = 1, km
       wz_QVap_T(:,k) = &

            ! dＵｑ/dλ + dＶｑ/dμ
            & - w_Div_xy_xy( &
            &               xyz_VelLon_n(:,:,k) * xyz_QVap_n(:,:,k), &
            &               xyz_VelLat_n(:,:,k) * xyz_QVap_n(:,:,k)  &
            &              ) / R0                              &

            ! Ｒ＝ｑＤ＋σ移流
            & + w_xy( xyz_QVapDivVSigmaAdv(:,:,k) )
            !&
            ! 水平拡散項、水蒸気ソースは後回し
    enddo

    !----- A = B + 2Δt * T-----
    wz_QVap_a(:,:) = &
         & wa_xya( xyz_QVap_b(:,:,:) )          &
         & + 2.0d0*DelTime * ( wz_QVap_T(:,:) )
!!!         & + 2.0d0*DelTime * wa_NumVisScaler_wa(  wa_xya( xyz_QVap_b(:,:,:) )  )

    !----- xyz_QVap_a に代入-----
    xyz_QVap_a(:,:,:) = xya_wa( wz_QVap_a(:,:) )

!!$    call DataDump('xyz_QVap_b', xyz_QVap_b(:,:,:), strlen=80)
!!$    call DataDump('xyz_QVap', xyz_QVap(:,:,:), strlen=80)
!!$    call DataDump('xyz_QVap_T', xya_wa(wz_QVap_T(:,:)), strlen=80)
!!$    call DataDump('xyz_QVap_a', xyz_QVap_a(:,:,:), strlen=80)

    call EndSub(subname)
  end subroutine dynamics_leapfrog


                                                                 !=begin
  !=== Calculate Diffusion Term
  !
  !t-Δt の値から水平拡散項を求め、それを t+Δt の値に加える。
  !
  subroutine dynamics_diffusion(                         &
    & xyz_Vor_b , xyz_Div_b , xyz_Temp_b , xyz_QVap_b  , &
    & xyz_Vor_a , xyz_Div_a , xyz_Temp_a , xyz_QVap_a        )

  !==== Dependency
    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
                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & xyz_Vor_b(:,:,:)  , & ! 渦度 (t-Δt)
         & xyz_Div_b(:,:,:)  , & ! 発散 (t-Δt)
         & xyz_Temp_b(:,:,:) , & ! 温度 (t-Δt)
         & xyz_QVap_b(:,:,:)     ! 比湿 (t-Δt)

    !==== In/Out
    !
    real(DBKIND), intent(inout) :: &
         & xyz_Vor_a(:,:,:)  , & ! 渦度 (t+Δt)
         & xyz_Div_a(:,:,:)  , & ! 発散 (t+Δt)
         & xyz_Temp_a(:,:,:) , & ! 温度 (t+Δt)
         & xyz_QVap_a(:,:,:)     ! 比湿 (t+Δt)
                                                                 !=end

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

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

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

    xyz_Div_a  = &
         & xya_wa(  &
         &    wa_xya(xyz_Div_a) &
         &    + 2.0d0*DelTime   &
         &      * wz_DiffVorDiv * wa_xya( xyz_Div_b ) &
         & )

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

    xyz_Temp_a  = &
         & xya_wa(  &
         &    wa_xya(xyz_Temp_a) &
         &    + 2.0d0*DelTime    &
         &      * wz_DiffTherm * wa_xya( xyz_Temp_b )  &
         & )

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

    xyz_QVap_a  = &
         & xya_wa(  &
         &    wa_xya(xyz_QVap_a) &
         &    + 2.0d0*DelTime    &
         &      * wz_DiffTherm * wa_xya( xyz_QVap_b )  &
         & )

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


    ! スペクトルデータとして見て拡散が効いているか確認するための
    ! データ出力を記述したソース。本計算にはまったく必要でない。
    !!wz_Vor_a  = wa_xya(xyz_Vor_a )
    !!wz_Div_a  = wa_xya(xyz_Div_a )
    !!wz_Temp_a = wa_xya(xyz_Temp_a)
    !!wz_QVap_a = wa_xya(xyz_QVap_a)
    !!
    !!write(80, *) CurrentTime, &
    !!     &       wz_Vor_a( l_nm(21,0), 1) , wz_Div_a( l_nm(21,0), 1) , &
    !!     &       wz_Temp_a( l_nm(21,0), 1) , wz_QVap_a( l_nm(21,0), 1)
    !!call flush(80)

    call EndSub(subname)
  end subroutine dynamics_diffusion


                                                                 !=begin
  !=== Calculate Diagnostic Values.
  !
  !診断的に得られる量の演算を行なう。
  !現在は渦度発散から速度成分 (経度方向、緯度方向) の演算を行なうのみである。
  !
  subroutine dynamics_diagnostic                              &
       & ( x_Lon       , y_Lat       , z_Sigma  , r_Sigma   , &

       &   xyz_VelLon_a, xyz_VelLat_a, xyz_Vor_a, xyz_Div_a , &
       &   xyz_Temp_a  , xyz_QVap_a  , xy_Ps_a               )

  !==== Dependency
    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
                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) :: &
         & x_Lon(:)          , & ! intent(in): 経度座標
         & y_Lat(:)          , & ! intent(in): 緯度座標
         & z_Sigma(:)        , & ! intent(in): σレベル(整数)座標
         & r_Sigma(:)    , & ! intent(in): σレベル(半整数)座標

         & xyz_Vor_a(:,:,:)  , & ! intent(in): 渦度         (t+Δt)
         & xyz_Div_a(:,:,:)  , & ! intent(in): 発散         (t+Δt)
         & xyz_Temp_a(:,:,:) , & ! intent(in): 温度         (t+Δt)
         & xyz_QVap_a(:,:,:) , & ! intent(in): 比湿         (t+Δt)
         & xy_Ps_a(:,:)          ! intent(in): 地表面気圧   (t+Δt)

    !==== In/Output
    !
    real(DBKIND), intent(out) :: &
         & xyz_VelLon_a(:,:,:) , & ! intent(out): 速度経度成分 (t+Δt)
         & xyz_VelLat_a(:,:,:)     ! intent(out): 速度緯度成分 (t+Δt)
                                                                 !=end

    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)
    if (.not. dynamics_initialized) then
       call EndSub( subname, 'Call dynamics_init before call %c',  &
            &       c1=trim(subname) )
       return
    endif

    !-------------------------------------------------------------------
    !   ｕ と ｖ を出力
    !-------------------------------------------------------------------
    wz_Psi_a(:,:) = wa_LaplaInv_wa(  wa_xya( xyz_Vor_a(:,:,:) )  ) * R0**2
    wz_Chi_a(:,:) = wa_LaplaInv_wa(  wa_xya( xyz_Div_a(:,:,:) )  ) * R0**2

    xyz_VelLon_a(:,:,:) = (  xya_GradLon_wa( wz_Chi_a(:,:) ) &
         &                - xya_GradLat_wa( wz_Psi_a(:,:) )  ) / R0

    xyz_VelLat_a(:,:,:) = (  xya_GradLon_wa( wz_Psi_a(:,:) ) &
         &                + xya_GradLat_wa( wz_Chi_a(:,:) )  ) / R0

!!$    call DataDump('xyz_VelLon_a', xyz_VelLon_a(:,:,:), strlen=80)
!!$    call DataDump('xyz_VelLat_a', xyz_VelLat_a(:,:,:), strlen=80)

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

    call io_gt4_out_Put('TotalMass', TotalMass)

    call EndSub(subname)
  end subroutine dynamics_diagnostic


                                                                 !=begin
  !=== Terminate module
  !
  !((<dynamics_init>)) で allocate した変数を deallocate し、
  !演算した値も全て破棄する。
  !
  subroutine dynamics_end
  !==== Dependency
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use dc_trace,    only: BeginSub, EndSub, DbgMessage
                                                                 !=end
    implicit none

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

  continue

    !-----------------------------------------------------------------
    !   Check Initialization
    !-----------------------------------------------------------------
    call BeginSub(subname)
    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_TempAve        , &
         &      z_TempAveHalf    , &
         &      zz_DivMtx        , &
         &      zz_TempMtx       , &
         &      zz_TempMtxQ      , &
         &      zz_TempMtxS      , &
         &      zz_TempMtxR      )


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

      &   xy_lnPs_T     , & ! πの時間変化 dπ/dt   (Δt)
      &   w_lnPs_T      , & ! πの時間変化 (スペクトル) (Δt)
      &   w_lnPs_a      , & ! π (t+Δt)

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

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

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

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

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

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

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

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

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

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

    call EndSub(subname)
  end subroutine dynamics_end


!!                                                                 !=begin
!!  !=== 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
!!                                                                 !=end
!!    implicit none
!!                                                                 !=begin
!!    !==== Input
!!    !
!!    real(DBKIND), intent(in) :: wa_data((nm+1)*(nm+1), km)
!!    !
!!    !==== Output
!!    !
!!    real(DBKIND)             :: wa_NumVis_wa((nm+1)*(nm+1), km)
!!                                                                 !=end
!!
!!    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
!!                                                                 !=end
!!    implicit none
!!                                                                 !=begin
!!    !==== Input
!!    !
!!    real(DBKIND), intent(in) :: wa_data((nm+1)*(nm+1), km)
!!    !
!!    !==== Output
!!    !
!!    real(DBKIND)             :: wa_NumVisScaler_wa((nm+1)*(nm+1), km)
!!                                                                 !=end
!!
!!    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
