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               )

  !==== Dependency

                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) ::  x_Lon(:)          ,  y_Lat(:)          ,  z_Sigma(:)        ,  r_Sigma(:)    ,  xyz_VorA(:,:,:)  ,  xyz_DivA(:,:,:)  ,  xyz_TempA(:,:,:) ,  xyz_QVapA(:,:,:) ,  xy_PsA(:,:)          ! intent(in): 地表面気圧   (t+Δt)

    !==== In/Output
    !
    real(DBKIND), intent(out) ::  xyz_VelLonA(:,:,:) ,  xyz_VelLatA(:,:,:)     ! 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_PsiA = wa_LaplaInv_wa(  wa_xya( xyz_VorA )  ) * R0**2
    wz_ChiA = wa_LaplaInv_wa(  wa_xya( xyz_DivA )  ) * R0**2

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

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

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

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

    call io_gt4_out_Put('TotalMass', TotalMass)

    call EndSub(subname)
  end subroutine dynamics_diagnostic

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

  !==== Dependency

                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) ::  xyz_VorB(:,:,:)  ,  xyz_DivB(:,:,:)  ,  xyz_TempB(:,:,:) ,  xyz_QVapB(:,:,:)     ! 比湿 (t-Δt)

    !==== In/Out
    !
    real(DBKIND), intent(inout) ::  xyz_VorA(:,:,:)  ,  xyz_DivA(:,:,:)  ,  xyz_TempA(:,:,:) ,  xyz_QVapA(:,:,:)     ! 比湿 (t+Δt)
                                                                 !=end

    real(DBKIND)    ::  wz_VorA((nm+1)**2, km)  ,  wz_DivA((nm+1)**2, km)  ,  wz_TempA((nm+1)**2, km) ,  wz_QVapA((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_VorA  =  xya_wa(      wa_xya(xyz_VorA)     + 2.0d0*DelTime         * wz_DiffVorDiv * wa_xya( xyz_VorB )  )

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

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

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

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

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

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

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


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

    call EndSub(subname)
  end subroutine dynamics_diffusion

  subroutine dynamics_init(x_Lon, y_Lat, z_Sigma, r_Sigma)
  !==== Dependency

                                                                 !=end
    implicit none
                                                                 !=begin
    !==== Input
    !
    real(DBKIND), intent(in) ::  x_Lon(:)           ,  y_Lat(:)           ,  z_Sigma(:)         ,  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_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('xyz_VSigmaHalf')
    call io_gt4_out_SetVars('xyz_VSigmaHalfNonG')
    call io_gt4_out_SetVars('xyz_DTempLocalDtN')
    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_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_lnPsA( (nm+1)*(nm+1) ) ,    xyz_VSigmaHalf(im, jm, km+1) ,    xyz_VSigmaHalfNonG(im, jm, km+1) ,    xyz_TempEdd(im, jm, km)   ,    xyz_TempVir(im, jm, km)   ,    xyz_TempVirEdd(im, jm, km ),    xyz_TempEddHalf(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