Class dynamics_mod
In: dynamics/dynamics.f90

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

現在、力学過程の全てがこのモジュール内にある。本来ならば、 解く方程式や項の種類によってモジュール化がなされるべきである。

Methods

dynamics_diagnostic   dynamics_diffusion   dynamics_end   dynamics_init   dynamics_leapfrog  

Included Modules

type_mod grid_3d_mod grid_wavenumber_mod constants_mod time_mod spml_mod io_gt4_out_mod dc_trace dc_string

Public Instance methods

dynamics_diagnostic( x_Lon, y_Lat, z_Sigma, r_Sigma, xyz_VelLonA, xyz_VelLatA, xyz_VorA, xyz_DivA, xyz_TempA, xyz_QVapA, xy_PsA )
Subroutine :
x_Lon(:) :real(DBKIND), intent(in)
: 経度座標
y_Lat(:) :real(DBKIND), intent(in)
: 緯度座標
z_Sigma(:) :real(DBKIND), intent(in)
: σレベル(整数)座標
r_Sigma(:) :real(DBKIND), intent(in)
: σレベル(半整数)座標
xyz_VelLonA(:,:,:) :real(DBKIND), intent(out)
: 速度経度成分 (t+Δt)
xyz_VelLatA(:,:,:) :real(DBKIND), intent(out)
: 速度緯度成分 (t+Δt)
xyz_VorA(:,:,:) :real(DBKIND), intent(in)
: 渦度 (t+Δt)
xyz_DivA(:,:,:) :real(DBKIND), intent(in)
: 発散 (t+Δt)
xyz_TempA(:,:,:) :real(DBKIND), intent(in)
: 温度 (t+Δt)
xyz_QVapA(:,:,:) :real(DBKIND), intent(in)
: 比湿 (t+Δt)
xy_PsA(:,:) :real(DBKIND), intent(in)
: 地表面気圧 (t+Δt)

(in)

Calculate Diagnostic Values.

診断的に得られる量の演算を行なう。 現在は渦度発散から速度成分 (経度方向、緯度方向) の演算を行なうのみである。

[Source]

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

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

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

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

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

    real(DBKIND) :: TotalMass      ! 全質量

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

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

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

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

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

!

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

    call io_gt4_out_Put('TotalMass', TotalMass)

    call EndSub(subname)
  end subroutine dynamics_diagnostic
dynamics_diffusion( xyz_VorB, xyz_DivB, xyz_TempB, xyz_QVapB, xyz_VorA, xyz_DivA, xyz_TempA, xyz_QVapA )
Subroutine :
xyz_VorB(:,:,:) :real(DBKIND), intent(in)
: 渦度 (t-Δt)
xyz_DivB(:,:,:) :real(DBKIND), intent(in)
: 発散 (t-Δt)
xyz_TempB(:,:,:) :real(DBKIND), intent(in)
: 温度 (t-Δt)
xyz_QVapB(:,:,:) :real(DBKIND), intent(in)
: 比湿 (t-Δt)
xyz_VorA(:,:,:) :real(DBKIND), intent(inout)
: 渦度 (t+Δt)
xyz_DivA(:,:,:) :real(DBKIND), intent(inout)
: 発散 (t+Δt)
xyz_TempA(:,:,:) :real(DBKIND), intent(inout)
: 温度 (t+Δt)
xyz_QVapA(:,:,:) :real(DBKIND), intent(inout)
: 比湿 (t+Δt)

[Source]

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

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

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

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

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

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

!

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

!

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

!

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

!


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

    call EndSub(subname)
  end subroutine dynamics_diffusion
dynamics_end( )
Subroutine :

Terminate module

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

[Source]

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

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

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

  continue

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

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


    deallocate ( xy_lnPs       , xyz_lnPsAdv   , xyz_lnPsAdvSum, xyz_DivSum    , xy_DlnPsDtN     , w_DlnPsDtN      , w_lnPsA      , xyr_VSigma     , xyr_VSigmaNonG , xyz_TempEdd    , xyz_TempVir    , xyz_TempVirEdd , xyr_TempEdd, xyz_UAN      , xyz_VAN      , wz_DVorDtN      , wz_VorA      , xyz_KE        , wz_DDivDtN        , wz_DivA        , wz_PresTendTemp , wz_PresTendPs   , xyz_DTempLocalDtN , wz_DTempDtN       , wz_TempA       , xyz_QVapDivVSigmaAdv, wz_DQVapDtN           , wz_QVapA )

    deallocate ( wz_PsiA , wz_ChiA )

    call EndSub(subname)
  end subroutine dynamics_end
dynamics_init( x_Lon, y_Lat, z_Sigma, r_Sigma )
Subroutine :
x_Lon(:) :real(DBKIND), intent(in)
: 経度座標
y_Lat(:) :real(DBKIND), intent(in)
: 緯度座標
z_Sigma(:) :real(DBKIND), intent(in)
: σレベル(整数)座標
r_Sigma(:) :real(DBKIND), intent(in)
: σレベル(半整数)座標

Initialize module and Calculate Non-Predictional Value.

以降のサブルーチンで用いる変数の allocate 、および 時間発展しない量の演算を行なう。 また、変数データ出力のための初期設定も行なう。

[Source]

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

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

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

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

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

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

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

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

!

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

!

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

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

!

    !-------------------------------------------------------------------
    !   運動量 (渦度発散) 拡散係数 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))

!

    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

!

    !-------------------------------------------------------------------
    !   温度鉛直補間の係数ａ、ｂの計算
    !-------------------------------------------------------------------
    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

!

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

    z_TempAvrXY = TempAve

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

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

!


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

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

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


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

    ! 初期化
    zz_siMtxW = 0.0

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

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

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

    ! 初期化
    zz_siMtxS = 0.0

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

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

    enddo
!

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

    ! 初期化
    zz_siMtxQ = 0.0

    do k = 1, km
      zz_siMtxQ(k,k) = ( r_TempAvrXY(k) - z_TempAvrXY(k) ) / z_DelSigma(k)
    enddo
    do k = 1, km - 1
      zz_siMtxQ(k,k+1) = ( z_TempAvrXY(k) - r_TempAvrXY(k+1) ) / z_DelSigma(k)
    enddo
!

    ! ＱＳ

    ! 初期化
    zz_siMtxQS = 0.0

    zz_siMtxQS = matmul(zz_siMtxQ, zz_siMtxS)


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

    ! 初期化
    zz_siMtxR = 0.0

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

      do kk = k+1, km
        zz_siMtxR(k,kk) = zz_siMtxR(k,kk) - z_HydroBeta(k) / z_DelSigma(k) * z_DelSigma(kk) * z_TempAvrXY(k)
      enddo
    enddo
!

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

    ! 初期化
    zz_siMtxH = 0.0

    zz_siMtxH = zz_siMtxQS - zz_siMtxR

    ! Ｗｈ:

    ! 初期化
    zz_siMtxWH = 0.0

    zz_siMtxWH = matmul(zz_siMtxW, zz_siMtxH)

!

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

    call io_gt4_out_SetVars('TotalMass')

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

    !-------------------------------------------------------------------
    !   dynamics_leapfrog で計算する変数の allocate 
    !-------------------------------------------------------------------
    allocate ( xy_lnPs (im, jm)          , xyz_lnPsAdv(im, jm, km)   , xyz_lnPsAdvSum(im, jm, km), xyz_DivSum(im, jm, km)    , xy_DlnPsDtN( im, jm )       , w_DlnPsDtN( (nm+1)*(nm+1) ) , w_lnPsN( (nm+1)*(nm+1) ) , xy_GradLonlnPsN( im, jm ), xy_GradLatlnPsN( im, jm ), wz_TempN( (nm+1)*(nm+1), km ) , w_lnPsA( (nm+1)*(nm+1) ) , xyr_VSigma(im, jm, km+1) , xyr_VSigmaNonG(im, jm, km+1) , xyz_TempEdd(im, jm, km)   , xyz_TempVir(im, jm, km)   , xyz_TempVirEdd(im, jm, km ), xyr_TempEdd(im, jm, km+1), xyz_UAN(im, jm, km)      , xyz_VAN(im, jm, km)      , wz_DVorDtN( (nm+1)*(nm+1), km ) , wz_VorA( (nm+1)*(nm+1), km ) , xyz_KE(im, jm, km)      , wz_DDivDtN( (nm+1)*(nm+1), km) , wz_DivA( (nm+1)*(nm+1), km) , wz_PresTendTemp( (nm+1)*(nm+1), km) , wz_PresTendPs( (nm+1)*(nm+1), km) , xyz_DTempLocalDtN(im, jm, km)  , wz_DTempDtN( (nm+1)*(nm+1), km ) , wz_TempA( (nm+1)*(nm+1), km ) , xyz_QVapDivVSigmaAdv( im, jm, km), wz_DQVapDtN( (nm+1)*(nm+1), km) , wz_QVapA( (nm+1)*(nm+1), km) )


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


    call EndSub(subname)
  end subroutine dynamics_init
dynamics_leapfrog( x_Lon, y_Lat, z_Sigma, r_Sigma, xyz_VelLonB, xyz_VelLatB, xyz_VorB, xyz_DivB, xyz_TempB, xyz_QVapB, xy_PsB, xyz_VelLonN, xyz_VelLatN, xyz_VorN, xyz_DivN, xyz_TempN, xyz_QVapN, xy_PsN, xyz_VelLonA, xyz_VelLatA, xyz_VorA, xyz_DivA, xyz_TempA, xyz_QVapA, xy_PsA )
Subroutine :
x_Lon(:) :real(DBKIND), intent(in)
: 経度座標
y_Lat(:) :real(DBKIND), intent(in)
: 緯度座標
z_Sigma(:) :real(DBKIND), intent(in)
: σレベル(整数)座標
r_Sigma(:) :real(DBKIND), intent(in)
: σレベル(半整数)座標
xyz_VelLonB(:,:,:) :real(DBKIND), intent(in)
: 速度経度成分 (t-Δt)
xyz_VelLatB(:,:,:) :real(DBKIND), intent(in)
: 速度緯度成分 (t-Δt)
xyz_VorB(:,:,:) :real(DBKIND), intent(in)
: 渦度 (t-Δt)
xyz_DivB(:,:,:) :real(DBKIND), intent(in)
: 発散 (t-Δt)
xyz_TempB(:,:,:) :real(DBKIND), intent(in)
: 温度 (t-Δt)
xyz_QVapB(:,:,:) :real(DBKIND), intent(in)
: 比湿 (t-Δt)
xy_PsB(:,:) :real(DBKIND), intent(in)
: 地表面気圧 (t-Δt)
xyz_VelLonN(:,:,:) :real(DBKIND), intent(in)
: 速度経度成分 (t)
xyz_VelLatN(:,:,:) :real(DBKIND), intent(in)
: 速度緯度成分 (t)
xyz_VorN(:,:,:) :real(DBKIND), intent(in)
: 渦度 (t)
xyz_DivN(:,:,:) :real(DBKIND), intent(in)
: 発散 (t)
xyz_TempN(:,:,:) :real(DBKIND), intent(in)
: 温度 (t)
xyz_QVapN(:,:,:) :real(DBKIND), intent(in)
: 比湿 (t)
xy_PsN(:,:) :real(DBKIND), intent(in)
: 地表面気圧 (t)
xyz_VelLonA(:,:,:) :real(DBKIND), intent(out)
: 速度経度成分 (t+Δt)
xyz_VelLatA(:,:,:) :real(DBKIND), intent(out)
: 速度緯度成分 (t+Δt)
xyz_VorA(:,:,:) :real(DBKIND), intent(out)
: 渦度 (t+Δt)
xyz_DivA(:,:,:) :real(DBKIND), intent(out)
: 発散 (t+Δt)
xyz_TempA(:,:,:) :real(DBKIND), intent(out)
: 温度 (t+Δt)
xyz_QVapA(:,:,:) :real(DBKIND), intent(out)
: 比湿 (t+Δt)
xy_PsA(:,:) :real(DBKIND), intent(out)
: 地表面気圧 (t+Δt)

(out)

Calculate Predictional Values.

時間発展する量の演算を行なう。 演算するデータの出力も行なう。

[Source]

  subroutine dynamics_leapfrog ( x_Lon      , y_Lat      , z_Sigma , r_Sigma  , xyz_VelLonB, xyz_VelLatB, xyz_VorB, xyz_DivB , xyz_TempB  , xyz_QVapB  , xy_PsB  , xyz_VelLonN, xyz_VelLatN, xyz_VorN, xyz_DivN , xyz_TempN  , xyz_QVapN  , xy_PsN  , xyz_VelLonA, xyz_VelLatA, xyz_VorA, xyz_DivA , xyz_TempA  , xyz_QVapA  , xy_PsA             )   !(out)
    !
    ! Calculate Predictional Values.
    !
    ! 時間発展する量の演算を行なう。
    ! 演算するデータの出力も行なう。
    !
    use type_mod,    only: STRING, REKIND, DBKIND, INTKIND
    use grid_3d_mod,         only: im, jm, km
    use grid_wavenumber_mod, only: nm
    use constants_mod, only: R0, Cp, EpsVT
    use time_mod,    only: DelTime, CurrentTime
    use spml_mod,    only: w_xy, xy_w , xy_GradLon_w, xy_GradLat_w, w_Div_xy_xy, w_LaplaInv_w, wa_xya, xya_wa, wa_Div_xya_xya, wa_Lapla_wa
    use io_gt4_out_mod,only: io_gt4_out_Put
    use dc_trace,    only: DbgMessage, BeginSub, EndSub, DataDump
    use dc_string,   only: toChar
    implicit none

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

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

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


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

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

    xy_lnPs = log( xy_PsN )

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

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


    ! π移流の積み下げ Σ[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)
!


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

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

    call io_gt4_out_Put('xyz_DivSum', xyz_DivSum)
!

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

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

!


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

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

!

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

!

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

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

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

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

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

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

    call io_gt4_out_Put('xyr_VSigma', xyr_VSigma)
    call io_gt4_out_Put('xyr_VSigmaNonG', xyr_VSigmaNonG)
!

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

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

!

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

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

    call io_gt4_out_Put('xyr_TempEdd', xyr_TempEdd)
!

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

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

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

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

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

!

!

    do k = 2, km

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

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

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

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

!

    do k = 1, km - 1

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

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

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

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

!

!

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

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

    ! dζ/dt の計算

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

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


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

!

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

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

    ! 初期化
    xyz_KE = 0.0

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

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

!

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

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

!
    call io_gt4_out_Put('xyz_KE', xyz_KE)

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

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

    ! 初期化
    wz_PresTendTemp = 0.0

    wz_TempN = wa_xya(xyz_TempN)

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

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

!


!


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

!

!

!

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

    wz_DDivDtN = wa_Div_xya_xya( xyz_UAN, xyz_VAN ) / R0 - wa_Lapla_wa( wa_xya(xyz_KE) ) / R0**2 - wa_Lapla_wa( wz_PresTendTemp + wz_PresTendPs ) / R0**2

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

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

!


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

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

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

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

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

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

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

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

!


!

!
    call io_gt4_out_Put('xyz_DTempLocalDtN',  xyz_DTempLocalDtN)


!

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

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

    do k = 1, km
      wz_DTempDtN(:,k) = - w_Div_xy_xy( xyz_VelLonN(:,:,k) * xyz_TempEdd(:,:,k), xyz_VelLatN(:,:,k) * xyz_TempEdd(:,:,k) ) / R0 + w_xy( xyz_DTempLocalDtN(:,:,k) )
      !&
      ! 非断熱加熱 Ｑ もまだ入ってないよ
      !&
      ! 摩擦熱および熱拡散は後回し
    enddo

!


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

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


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

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


    call io_gt4_out_Put(  'xyz_Temp_T', xya_wa( wz_DTempDtN )  )
!

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

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

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

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

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

    do k = 1, km
      wz_DQVapDtN(:,k) = - w_Div_xy_xy( xyz_VelLonN(:,:,k) * xyz_QVapN(:,:,k), xyz_VelLatN(:,:,k) * xyz_QVapN(:,:,k) ) / R0 + w_xy( xyz_QVapDivVSigmaAdv(:,:,k) )
        !&
        ! 水平拡散項、水蒸気ソースは後回し
    enddo

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

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

!

    call EndSub(subname)
  end subroutine dynamics_leapfrog

[Validate]