!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
! dynamics.f90 - DCPAM Dynamical Core Module
!
!   [支配方程式系 : 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 : リープフロッグ
!
! History
!   2004/07/29 Morikawa Yasuhiro     create
!   2004/11/04 Morikawa Yasuhiro     update

module dynamics_mod

  use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
  implicit none

  private
  public :: dynamics_init, dynamics_leapfrog, dynamics_diagnostic, &
       &    dynamics_end

  !-------------------------------------------------------------------
  !   dynamics_init において値を設定するもの
  !-------------------------------------------------------------------
  real(DBKIND) :: Kappa               ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱)
  real(DBKIND), pointer, save :: &
       &
       & xy_Coli(:,:)        , & ! 格子点データ(コリオリ力)
       & xy_UVFact(:,:)      , & ! u→U, v→V のファクター
       & z_DelSigma(:)       , & ! Δσ(整数レベル)
       &
       & w_DiffV(:)          , & ! 運動量水平拡散係数
       & w_DiffTherm(:)      , & ! 熱、水水平拡散係数
       &
       & z_TempAve(:)        , & ! 平均温度 (整数レベル)
       & z_TempAveHalf(:)    , & ! 平均温度 (半整数レベル)
       &
       & z_HydroAlpha(:)     , & ! 静水圧の式の係数 α
       & z_HydroBeta(:)      , & ! 静水圧の式の係数 β
       & z_TempInpolKappa(:) , & ! 温度鉛直補間の係数κ
       & z_TempInpolA(:)     , & ! 温度鉛直補間の係数ａ
       & z_TempInpolB(:)     , & ! 温度鉛直補間の係数ｂ
       &
       ! semi-implicit用行列
       & zz_DivMtx(:,:)      , & ! Ｗ:発散の式での温度補正 (線形重力波項の効果)
       & zz_TempMtx(:,:)     , & ! ｈ: ＱＳ (σ移流) - Ｒ
       & zz_TempMtxQ(:,:)    , & ! Ｑ: dＴ/dσ (線形重力波項の効果)
       & zz_TempMtxS(:,:)    , & ! Ｓ: dσ/dｔ (線形重力波項の効果)
       & zz_TempMtxR(:,:)        ! Ｒ: ＲＤ = κＴ（∂π/∂ｔ＋ｄσ/ｄｔ／σ）
                                 !     気圧変化の効果の係数
                                 !     (線形重力波項の効果)


  !-------------------------------------------------------------------
  !   dynamics_leapfrog において利用するもの
  !-------------------------------------------------------------------
  real(DBKIND), pointer, save :: &
       &
       & xy_lnPs(:,:)          , & ! π＝ ln Ｐs           (t)
       & xy_lnPs_T(:,:)        , & ! πの時間変化 dπ/dt   (Δt)
       & w_lnPs_T(:)           , & ! πの時間変化 (スペクトル) (Δt)
       & w_lnPs_A(:)           , & ! π (t+Δt)
       &
       & xyz_lnPsAdv(:,:,:)    , & ! π の移流 ｖ・∇π
       & xyz_lnPsAdvSum(:,:,:) , & ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ
       & xyz_DivSum(:,:,:)     , & ! 発散の積下げ Σ[k〜K] ＤΔσ
       & xyz_VSigmaHalf(:,:,:) , & ! 鉛直σ速度(半整数レベル) (t)
       & xyz_VSigmaHalfNonG(:,:,:), & ! 鉛直σ速度 非重力波成分 (半整数レベル) (t)
       &
       & xyz_TempEdd(:,:,:)    , & ! Ｔ' ：温度の擾乱 (整数レベル)     (t)
       & xyz_TempEddHalf(:,:,:), & ! Ｔ' ：温度の擾乱 (半整数レベル)   (t)
       & xyz_TempVir(:,:,:)    , & ! Ｔv ：仮温度 (virtual temperature)(t)
       & xyz_TempVirEdd(:,:,:) , & ! Ｔv'：仮温度の擾乱                (t)
       &
       & xyz_UA_T(:,:,:)       , & ! 東西運動量変化項ＵＡ
       & xyz_VA_T(:,:,:)       , & ! 南北運動量変化項ＶＡ
       &
       & wz_Vor_T(:,:)         , & ! 渦度変化のスペクトルデータ (Δt)
       & wz_Vor_A(:,:)         , & ! 渦度のスペクトルデータ    (t+Δt)
       &
       & xyz_KE(:,:,:)         , & ! 運動エネルギー + 仮温度補正
       &                           ! ｕ**2+ｖ**2    + ΣＷ(Ｔv-Ｔ)
       & xy_Phi(:,:)           , & ! 地表ジオポテンシャルΦ
       &
       & 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)


  logical, save :: dynamics_initialized = .false.

contains

  subroutine dynamics_init(x_Lon, y_Lat, z_Sigma, z_SigmaHalf)
    !----- 使用モジュール -----
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use grid_3d_mod,         only: im, jm, km
    use grid_wavenumber_mod, only: nm
    use physics_mod, only: physics_init, R0, Omega, Cp, RAir, Nu, &
         &                 TempAve, ViscosOrder
    use time_mod,    only: DelTime
    use spml_mod,    only: spml_init, xy_Lat
    use io_gt4_out_mod,only: io_gt4_out_init, io_gt4_out_SetVars
    use dc_trace,    only: SetDebug, Message, BeginSub, EndSub, DataDump
    use dc_string,   only: toChar
    implicit none

    !-------------------------------------------------------------------
    !   変数定義
    !-------------------------------------------------------------------
    !----- 入力変数 -----
    real(DBKIND), intent(in) :: &
         & x_Lon(:)           , & ! intent(in): 経度座標
         & y_Lat(:)           , & ! intent(in): 緯度座標
         & z_Sigma(:)         , & ! intent(in): σレベル(整数)座標
         & z_SigmaHalf(:)         ! intent(in): σレベル(半整数)座標

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

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

    !----------------------------------------------------------------
    !   初期化確認
    !----------------------------------------------------------------
    call BeginSub(subname)
    if (dynamics_initialized) then
       call EndSub( subname, '%c is already called.', c1=trim(subname) )
       return
    else
       dynamics_initialized = .true.
    endif

    !-------------------------------------------------------------------
    !   物理定数の取得
    !-------------------------------------------------------------------
    call physics_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) = z_SigmaHalf(k) - z_SigmaHalf(k+1)
    enddo

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

    !-------------------------------------------------------------------
    !   運動量拡散係数 w_DiffV, 熱拡散係数 w_DiffTherm の計算
    !     ※ 面倒なので後回し
    !-------------------------------------------------------------------
!!$    allocate(w_DiffV((nm+1)*(nm+1)), w_DiffTherm((nm+1)*(nm+1)))
!!$    
!!$    do l = 1, (nm+1)*(nm+1)
!!$       w_DiffTherm(l) = Nu * &   ! マイナスが付くのが正しいはずだけど…
!!$            & ( (- l * (l + 1.0))/(R0**2) )**(ViscosOrder/2)
!!$       w_DiffV(l)     = par%Nu * &   ! マイナスが付くのが正しいはずだけど…
!!$            & (        &
!!$            &  ( (- l * (l + 1.0))/(R0**2) )**(ViscosOrder/2) &
!!$            &  - ( 2.0/(R0**2) )**(ViscosOrder/2)             &
!!$            & )
!!$    enddo
!!$
!!$    do l = 1, (nm+1)*(nm+1)
!!$       call Message('w_DiffTherm(%d)=<%c>', &
!!$            i=(/l/), c1=trim(toChar(w_DiffTherm(l))))
!!$       call Message('w_DiffV(%d)=<%c>', &
!!$            i=(/l/), c1=trim(toChar(w_DiffV(l))))
!!$    enddo



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

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

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

!!$    call Message('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) = &
            & ( &
            &    z_SigmaHalf(k)   * z_HydroAlpha(k) &
            &  + z_SigmaHalf(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

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

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

!!$    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) = z_SigmaHalf(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')


    !-------------------------------------------------------------------
    !   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)
      & )

    call EndSub(subname)
  end subroutine dynamics_init




  subroutine dynamics_leapfrog                              &
       & ( x_Lon     , y_Lat     , z_Sigma  , z_SigmaHalf , &
       &
       &   xyz_VLon_B, xyz_VLat_B, xyz_Vor_B, xyz_Div_B   , &
       &   xyz_Temp_B, xyz_QVap_B, xy_Ps_B  ,               &
       &
       &   xyz_VLon  , xyz_VLat  , xyz_Vor  , xyz_Div     , &
       &   xyz_Temp  , xyz_QVap  , xy_Ps    ,               &
       &
       &   xyz_VLon_A, xyz_VLat_A, xyz_Vor_A, xyz_Div_A   , &
       &   xyz_Temp_A, xyz_QVap_A, xy_Ps_A               )

    !----- 使用モジュール -----
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use grid_3d_mod,         only: im, jm, km
    use grid_wavenumber_mod, only: nm
    use physics_mod, only: R0, Cp, EpsVT
    use time_mod,    only: DelTime
    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: SetDebug, Message, BeginSub, EndSub, DataDump
    use dc_string,   only: toChar
    implicit none

    !-------------------------------------------------------------------
    !   変数定義
    !-------------------------------------------------------------------
    !----- 入力変数 -----
    real(DBKIND), intent(in) :: &
         & x_Lon(:)          , & ! intent(in): 経度座標
         & y_Lat(:)          , & ! intent(in): 緯度座標
         & z_Sigma(:)        , & ! intent(in): σレベル(整数)座標
         & z_SigmaHalf(:)    , & ! intent(in): σレベル(半整数)座標
         &
         & xyz_VLon_B(:,:,:) , & ! intent(in): 速度経度成分 (t-Δt)
         & xyz_VLat_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_VLon(:,:,:)   , & ! intent(in): 速度経度成分 (t)
         & xyz_VLat(:,:,:)   , & ! intent(in): 速度緯度成分 (t)
         & xyz_Vor(:,:,:)    , & ! intent(in): 渦度         (t)
         & xyz_Div(:,:,:)    , & ! intent(in): 発散         (t)
         & xyz_Temp(:,:,:)   , & ! intent(in): 温度         (t)
         & xyz_QVap(:,:,:)   , & ! intent(in): 比湿         (t)
         & xy_Ps(:,:)            ! intent(in): 地表面気圧   (t)

    !----- 出力変数 -----
    real(DBKIND), intent(out) :: &
         & xyz_VLon_A(:,:,:) , & ! intent(out): 速度経度成分 (t+Δt)
         & xyz_VLat_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)

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


    character(STRING),  parameter:: subname = "dynamics_leapfrog"
    !----------------------------------------------------------------
    !   初期化確認
    !----------------------------------------------------------------
    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(:,:) )

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

    ! 地表面気圧 π の移流 = ｖ・∇π
    do k = 1, km
       xyz_lnPsAdv(:,:,k) = &
            &  (  &
            &   xyz_VLon(:,:,k)                         &
            &      * xy_GradLon_w( w_xy(xy_lnPs(:,:)) ) &
            &
            &   + xyz_VLat(:,:,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(:,:, km ) * z_DelSigma( km )

    do k = km -1 , 1, -1
       xyz_DivSum(:,:,k) = &
            &   xyz_DivSum(:,:,k+1)           &
            & + xyz_Div(:,:,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*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)

    !----- xyz_All_A%xy_Ps に代入-----
    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](Ｄ + ｖ・∇π)Δσ )
            & z_SigmaHalf(k) &
            &     * ( xyz_lnPsAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) &
            &
            ! - Σ[k〜K] (Ｄ + ｖ・∇π) Δσ
            & - ( xyz_lnPsAdvSum(:,:,k) + xyz_DivSum(:,:,k) )

       xyz_VSigmaHalfNonG(:,:,k) =  &
            &
            ! - σ[k-1/2] × ( Σ[1〜K](ｖ・∇π)Δσ )
            & z_SigmaHalf(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

    ! 順圧渦度のための強制 (本来ここにあるべきではない)
    if (km == 2) then
       xyz_VSigmaHalf(:,:,2) = 0.0
    endif
    

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

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

       xyz_TempEdd(:,:,k)   = xyz_Temp(:,:,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(:,:,k)   &
            & + z_TempInpolB(k-1) * xyz_Temp(:,:,k-1) &
            & - z_TempAveHalf(k)
    enddo

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

    !-------------------------------------------------------------------
    !   ＵＡ＝Ｖζ＋σ移流、ＶＡ＝Ｕζ＋σ移流 (渦度と発散の式で利用)
    !-------------------------------------------------------------------

    do k = 1, km
       xyz_UA_T(:,:,k) = &
            ! (ζ + ｆ) Ｖ
            & ( xyz_Vor(:,:,k) + xy_Coli(:,:) ) &
            &   * ( xyz_VLat(:,:,k) * xy_UVFact (:,:) ) &
            ! -Ｃp κ (1/ａ) Ｔv (dπ/dλ)
            & - Cp * z_TempInpolKappa(k) * (1.0 / R0) &
            &     * xyz_TempVirEdd(:,:,k) &
            &        * xy_GradLon_w(  w_xy( xy_lnPs(:,:) )  ) / R0

       xyz_VA_T(:,:,k) = &
            ! - (ζ + ｆ) Ｕ
            & - ( xyz_Vor(:,:,k) + xy_Coli(:,:) ) &
            &     * ( xyz_VLon(:,:,k) * xy_UVFact (:,:) ) &
            ! -Ｃp κ (1/ａ) Ｔv (dπ/dμ)
            & - Cp * z_TempInpolKappa(k) * (1.0 / R0) &
            &     * xyz_TempVirEdd(:,:,k) &
            &        * xy_GradLat_w(  w_xy( xy_lnPs(:,:) )  ) / R0
    enddo

!!$    call DataDump('xyz_VLon', xyz_VLon(:,:,:), strlen=80)
!!$    call DataDump('xyz_VLat', xyz_VLat(:,:,:), 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 * z_DelSigma(k)) &
            &    * xyz_VSigmaHalf(:,:,k)  &
            &       * (xyz_VLon(:,:,k-1) - xyz_VLon(:,:,k)) &
            &         * xy_UVFact (:,:)

       xyz_VA_T(:,:,k) = xyz_VA_T(:,:,k) &
            ! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (Ｖ[k-1] - Ｖ[k])
            & - 1.0 / (2 * z_DelSigma(k)) &
            &    * xyz_VSigmaHalf(:,:,k)  &
            &       * (xyz_VLat(:,:,k-1) - xyz_VLat(:,:,k)) &
            &         * xy_UVFact (:,:)
    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 * z_DelSigma(k)) &
            &    * xyz_VSigmaHalf(:,:,k+1) &
            &       * (xyz_VLon(:,:,k) - xyz_VLon(:,:,k+1)) &
            &         * xy_UVFact(:,:)

       xyz_VA_T(:,:,k) = xyz_VA_T(:,:,k) &
            ! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (Ｖ[k] - Ｖ[k+1])
            & - 1.0 / (2 * z_DelSigma(k)) &
            &    * xyz_VSigmaHalf(:,:,k+1)  &
            &       * (xyz_VLat(:,:,k) - xyz_VLat(:,:,k+1)) &
            &         * xy_UVFact(:,:)
    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*DelTime * ( wz_Vor_T(:,:) )

    !----- xyz_All_A%xyz_Vor に代入-----
    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(:,:,1) )

    do k = 2, km
       xyz_KE(:,:,k) = xyz_KE(:,:,k-1)   &
            & + Cp * z_HydroAlpha(k) &
            &     * ( xyz_TempVir(:,:,k) - xyz_Temp(:,:,k) ) &
            & + Cp * z_HydroBeta(k-1) &
            &     * ( xyz_TempVir(:,:,k-1) - xyz_Temp(:,:,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_VLon(:,:,:)**2   &
         &     + xyz_VLat(:,:,:)**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(:,:,1) )

    do k = 2, km
       wz_PresTendTemp(:,k) = wz_PresTendTemp(:,k-1) &
            & + Cp * z_HydroAlpha(k)  * w_xy( xyz_Temp(:,:,k) )   &
            & + Cp * z_HydroBeta(k-1) * w_xy( xyz_Temp(:,:,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(:,:) )  )
    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*DelTime * ( wz_Div_T(:,:) )

    !----- xyz_All_A%xyz_Div に代入-----
    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(:,:,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 DataDump('xyz_TempLocal_T', xyz_TempLocal_T(:,:,:), strlen=80)

    !-------------------------------------------------------------------
    !   ＵＴ' 、ＶＴ'   (熱力学の式)
    !-------------------------------------------------------------------
!!$      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_VLon(:,:,k) * xy_UVFact(:,:)   &
            &                  * xyz_TempEdd(:,:,k)                  , &
            &               xyz_VLat(:,:,k) * xy_UVFact(:,:)   &
            &                  * 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(:,:,kk) )
               !&
               ! 拡散項は後回し
       enddo
    enddo


    !----- A = B + 2Δt * T-----
    wz_Temp_A(:,:) = &
         & wa_xya( xyz_Temp_B(:,:,:) ) &
         & + 2*DelTime * ( wz_Temp_T(:,:) )

    !----- xyz_All_A%xyz_Div に代入-----
    xyz_Temp_A(:,:,:) = xya_wa( wz_Temp_A(:,:) )

!!$    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(:,:,:) * xyz_Div(:,:,:)

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

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

    ! 水平移流 dＵｑ/dλ + dＶｑ/dμ と水平拡散項、水蒸気ソースを計算
    do k = 1, km
       wz_QVap_T(:,k) = &
            &
            ! dＵｑ/dλ + dＶｑ/dμ
            & - w_Div_xy_xy( &
            &               xyz_VLon(:,:,k) * xy_UVFact(:,:)   &
            &                  * xyz_QVap(:,:,k)             , &
            &               xyz_VLat(:,:,k) * xy_UVFact(:,:)   &
            &                  * xyz_QVap(:,:,k)               &
            &              ) / R0                                  &
            &
            ! Ｒ＝ｑＤ＋σ移流
            & + w_xy( xyz_QVapDivVSigmaAdv(:,:,k) )
            !&
            ! 水平拡散項、水蒸気ソースは後回し
    enddo

    !----- A = B + 2Δt * T-----
    wz_QVap_A(:,:) = &
         & wa_xya( xyz_QVap_B(:,:,:) ) &
         & + 2*DelTime * ( wz_QVap_T(:,:) )

    !----- xyz_All_A%xyz_Div に代入-----
    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

  !
  ! 本当は力学部分に入れるものではないが、とりあえず診断方程式
  ! による ｕ と ｖ の計算をおこなう。
  !
  subroutine dynamics_diagnostic                            &
       & ( x_Lon     , y_Lat     , z_Sigma  , z_SigmaHalf , &
       &
       &   xyz_VLon_A, xyz_VLat_A, xyz_Vor_A, xyz_Div_A   , &
       &   xyz_Temp_A, xyz_QVap_A, xy_Ps_A               )

    !----- 使用モジュール -----
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use grid_3d_mod,         only: im, jm, km
    use grid_wavenumber_mod, only: nm
    use physics_mod, only: R0
    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 dc_trace,    only: SetDebug, Message, BeginSub, EndSub, DataDump
    use dc_string,   only: toChar
    implicit none

    !-------------------------------------------------------------------
    !   変数定義
    !-------------------------------------------------------------------
    !----- 入力変数 -----
    real(DBKIND), intent(in) :: &
         & x_Lon(:)          , & ! intent(in): 経度座標
         & y_Lat(:)          , & ! intent(in): 緯度座標
         & z_Sigma(:)        , & ! intent(in): σレベル(整数)座標
         & z_SigmaHalf(:)    , & ! 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)

    !----- 入出力変数 -----
    real(DBKIND), intent(out) :: &
         & xyz_VLon_A(:,:,:) , & ! intent(out): 速度経度成分 (t+Δt)
         & xyz_VLat_A(:,:,:)     ! intent(out): 速度緯度成分 (t+Δt)

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

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

    !-------------------------------------------------------------------
    !   ｕ と ｖ を出力
    !-------------------------------------------------------------------
    do k = 1, km
       xyz_VLon_A(:,:,k) = &
            & xy_w( &
            &    w_LaplaInv_w( &
            &         w_Div_xy_xy( &
            &                     xyz_Div_A(:,:,k)   , &
            &                     - xyz_Vor_A(:,:,k)   &
            &                    ) / R0                    &
            &                ) * R0**2                     &
            &     ) / xy_UVFact(:,:)

       xyz_VLat_A(:,:,k) = &
            & xy_w( &
            &    w_LaplaInv_w( &
            &         w_Div_xy_xy( &
            &                     xyz_Vor_A(:,:,k) , &
            &                     xyz_Div_A(:,:,k)   &
            &                    ) / R0                  &
            &                ) * R0**2                   &
            &     ) / xy_UVFact(:,:)
    enddo

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

    call EndSub(subname)
  end subroutine dynamics_diagnostic


  subroutine dynamics_end
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use physics_mod, only: physics_end
    use dc_trace,    only: SetDebug, BeginSub, EndSub, Message
    implicit none

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

  continue

    !-----------------------------------------------------------------
    !   初期化確認
    !-----------------------------------------------------------------
    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

    !-----------------------------------------------------------------
    !   物理過程の方をとりあえず終了
    !-----------------------------------------------------------------
    call physics_end

    !-----------------------------------------------------------------
    !   allocate した変数の nullify とか
    !-----------------------------------------------------------------


    call EndSub(subname)
  end subroutine dynamics_end

end module dynamics_mod
