!= dcpam 主プログラム
!
!= dcpam main program
!
! Authors::   Yasuhiro Morikawa, Satoshi Noda, Yoshiyuki O. Takahashi, Shin-ichi Takehiro
! Version::   $Id: dcpam_main_hs94.f90,v 1.67 2026/08/15 00:00:00 takepiro Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
program dcpam_tem_test
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! モデルの使い方については {チュートリアル}[link:../../../doc/tutorial/rakuraku/] を
  ! 参照してください.
  !
  ! See {Tutorial}[link:../../../doc/tutorial/rakuraku/index.htm.en] for usage of the 
  ! model. 
  !

  ! モジュール引用 ; USE statements
  !
  ! 時刻管理
  ! Time control
  !
  use timeset, only: TimesetProgress, &
    & TimeB, &                ! ステップ $ t - \Delta t $ の時刻. 
                              ! Time of step $ t - \Delta t $. 
    & TimeN, &                ! ステップ $ t $ の時刻. 
                              ! Time of step $ t $. 
    & TimeA, &                ! ステップ $ t + \Delta t $ の時刻. 
                              ! Time of step $ t + \Delta t $. 
    & EndTime, &              ! 計算終了時刻. 
                              ! End time of calculation
    & DelTime                 ! $ \Delta t $ [s]

!!$  ! リスタートデータ入出力
!!$  ! Restart data input/output
!!$  !
!!$  use restart_file_io, only: RestartFileOutPut

  ! ヒストリデータ出力
  ! History data output
  !
  use gtool_historyauto, only: HistoryAutoPut, HistoryAutoAllVarFix

  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
!!$  use composition, only: &
!!$    &                    ncmax, &
!!$                              ! 成分の数
!!$                              ! Number of composition
!!$    &                    a_QMixName, &
!!$                              ! 成分の変数名
!!$                              ! Name of variables for composition
!!$    &                    a_QMixLongName, &
!!$                              ! 成分の長い変数名
!!$                              ! Long name of variables for composition
!!$    &                    IndexH2OVap, &
!!$    &                    IndexH2OLiq, &
!!$    &                    IndexH2OSol, &
!!$    &                    IndexTKE,    &
!!$    &                    CompositionInqIndex


  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, &  ! 経度格子点数. 
                              ! Number of grid points in longitude
    &                jmax, &  ! 緯度格子点数. 
                              ! Number of grid points in latitude
    &                kmax     ! 鉛直層数. 
                              ! Number of vertical level

  ! 座標データ設定
  ! Axes data settings
  !
  use axesset, only: &
       x_Lon, &
                            ! $ \lambda $ [rad.] . 経度. Longitude
       y_Lat, &
                            ! $ \varphi $ [rad.] . 緯度. Latitude
       z_Sigma, &
                            ! $ \sigma $ レベル (整数). 
                            ! Full $ \sigma $ level
       r_Sigma
                            ! $ \sigma $ レベル (半整数). 
                            ! Half $ \sigma $ level
  !
  !
  use dc_message, only: MessageNotify

  ! 物理・数学定数
  use constants0, only : PI
                            ! $ \pi $ .
                            ! 円周率.  Circular constant


  ! 物理定数設定
  ! Physical constants settings
  !
  use constants, only: ConstantsInit, &
    & RPlanet,&
                            ! $ a $ [m]. 
                            ! 惑星半径. 
                            ! Radius of planet
    & GasRDry, &
                            ! $ R $ [J kg-1 K-1]. 
                            ! 乾燥大気の気体定数. 
                            ! Gas constant of air
    & CpDry,&
                            ! $ C_p $ [J kg-1 K-1]. 
                            ! 乾燥大気の定圧比熱. 
                            ! Specific heat of air at constant pressure
    & Omega,&
                            ! $ \Omega $ [s-1]. 
                            ! 回転角速度. 
                            ! Angular velocity
    & Grav
                            ! $ g $ [m s-2]. 
                            ! 重力加速度. 
                            ! Gravitational acceleration

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision. 
    &                 STRING, &  ! 文字列.       Strings. 
    &                 TOKEN      ! キーワード.   Keywords. 

  !
  ! Calc Mass stream function
  !
  use mass_streamfunc, only :  yr_MassStreamFunc

  !
  ! TEM analysis
  !
  use tem, only :  TEMAnalysis

  ! 宣言文 ; Declaration statements
  !
  implicit none

  character(*), parameter:: prog_name = 'dcpam_main'
                            ! 主プログラム名. 
                            ! Main program name

  ! 予報変数
  !
  real(DP), allocatable :: xyz_U (:,:,:)
                              ! $ u $ .   東西風速. Eastward wind (m s-1)
  real(DP), allocatable :: xyz_V (:,:,:)
                              ! $ v $ .   南北風速. Northward wind (m s-1)
  real(DP), allocatable :: xyr_SigDot (:,:,:)
                              ! $ \dot\sigma $ . 鉛直流
                              ! vertical flow
  real(DP), allocatable :: xyz_Temp (:,:,:)
                              ! $ T $ .   温度. Temperature (K)
  real(DP), allocatable :: xyz_Entro (:,:,:)
                              ! $ T $ .   エントロピー. Entropy (J/K)
  real(DP), allocatable :: xyz_QVap(:,:,:)
                              ! $ q $ .   混合比. Mass mixing ratio of constituents (1)
  real(DP), allocatable :: xy_Ps (:,:)
                              ! $ p_s $ . 地表面気圧. Surface pressure (Pa)

  real(DP), allocatable :: xyz_EpLat (:,:,:)
               ! $ Ep_lat $ . EP フラックス緯度成分. 
               ! latitudinal EP flux

  real(DP), allocatable :: xyr_EpSig (:,:,:)
               ! $ Ep_lat $ . EP フラックス鉛直成分. 
               ! vertical EP flux

  real(DP), allocatable :: xyz_PsVRes (:,:,:)
               ! $ vres $ . 残差南北風速. 
               ! Residual northward wind

  real(DP), allocatable :: xyr_PsSigDotRes (:,:,:)
               ! $ \dot\sigma_res $ . 残差鉛直風速. 
               ! Residual vertical wind

  real(DP), allocatable :: xyr_StrRes (:,:,:)
               ! $ \psi_res $ . 残差質量流線関数. 
               ! Residual mass stream fkuction

  real(DP), allocatable :: yr_Msf    (:,:)
               ! $ \psi $ . 質量流線関数. 
               ! Mass stream fkuction


  real(DP), allocatable :: xyz_EpLatAns (:,:,:)
               ! $ Ep_lat $ . EP フラックス緯度成分. 
               ! latitudinal EP flux

  real(DP), allocatable :: xyr_EpSigAns (:,:,:)
               ! $ Ep_lat $ . EP フラックス鉛直成分. 
               ! vertical EP flux

  real(DP), allocatable :: xyz_PsVResAns (:,:,:)
               ! $ vres $ . 残差南北風速. 
               ! Residual northward wind

  real(DP), allocatable :: xyr_PsSigDotResAns (:,:,:)
               ! $ \dot\sigma_res $ . 残差鉛直風速. 
               ! Residual vertical wind

  real(DP), allocatable :: xyr_StrResAns (:,:,:)
               ! $ \psi_res $ . 残差質量流線関数. 
               ! Residual mass stream fkuction

  character(len=1)       :: TestMode = 'A'
                            ! TEM テスト計算モード
                            !   Z: 帯状流れ
                            !   W : 波擾乱
                            !   R : 残差+波擾乱
                            !   A : すべて
    
  !
  ! データパラメター 
  !
  real(DP), parameter :: Ps = 1D5
  real(DP), parameter :: T0 = 300.0D0
  real(DP), parameter :: s0    = 5700.0D0
  real(DP), parameter :: Gamma = 1000.0D0
  real(DP), parameter :: s1 = 100.D0
  real(DP), parameter :: U0 = 30.0D0
  real(DP), parameter :: V0 = 1.0D0
  real(DP), parameter :: W0 = 1.0D0
  real(DP), parameter :: U1 = 1.0D0
  real(DP), parameter :: V1 = 1.0D-1
  real(DP), parameter :: W1 = 1.0D-1

  ! 実行文 ; Executable statement
  !

  ! 主プログラムの初期化 (内部サブルーチン)
  ! Initialization for the main program (Internal subroutine)
  !
  call MainInit
  !
  ! データ設定
  ! Set data 
  !
  call set_basicstate

  if ( TestMode == 'Z' .OR. TestMode == 'A' ) then
     call set_zonal_test
     call set_zonal_answer
  end if
  if ( TestMode == 'W' .OR. TestMode == 'A' ) then
     call set_disturbance1_test
     call set_disturbance1_answer
  endif
  if ( TestMode == 'R' .OR. TestMode == 'A' ) then
     call set_disturbance2_test
     call set_disturbance2_answer
  endif

  ! ヒストリデータ出力 (スタート時刻)
  ! History data output (Start time)
  !
  call HistoryAutoPut( TimeN, 'U', xyz_U )
  call HistoryAutoPut( TimeN, 'V', xyz_V )
  call HistoryAutoPut( TimeN, 'SigDot', xyr_SigDot )
  call HistoryAutoPut( TimeN, 'Temp', xyz_Temp )
  call HistoryAutoPut( TimeN, 'Ps', xy_Ps )

  ! 診断変数の割付
  ! Allocation of diagnostic variables
  !
  yr_MSF = yr_MassStreamFunc( xy_Ps, xyz_V ) 
  call HistoryAutoPut( TimeN, 'Msf', spread(yr_MSF,1,imax) )

  !
  ! TEM 計算
  !
  call TEMAnalysis(                                              &
       & xyz_U,  xyz_V, xyr_SigDot, xyz_QVap, xyz_Temp, xy_Ps,   & ! (in)
       & xyz_EPLat,   xyr_EPSig, xyz_PsVRes, xyr_PsSigDotRes, & ! (out)
       & xyr_StrRes                                              & ! (out)
       & )

  call HistoryAutoPut( TimeN, 'VRes', xyz_PsVRes )
  call HistoryAutoPut( TimeN, 'SigDotRes', xyr_PsSigDotRes )
  call HistoryAutoPut( TimeN, 'StrRes', xyr_StrRes )
  call HistoryAutoPut( TimeN, 'EPLat', xyz_EPLat )
  call HistoryAutoPut( TimeN, 'EPSig', xyr_EPSig )

  call HistoryAutoPut( TimeN, 'VResAns', xyz_PsVResAns )
  call HistoryAutoPut( TimeN, 'SigDotResAns', xyr_PsSigDotResAns )
  call HistoryAutoPut( TimeN, 'StrResAns', xyr_StrResAns )
  call HistoryAutoPut( TimeN, 'EPLatAns', xyz_EPLatAns )
  call HistoryAutoPut( TimeN, 'EPSigAns', xyr_EPSigAns )

  !
  ! TEM チェック
  !
  call CheckTEM
  
  ! 主プログラムの終了処理 (内部サブルーチン)
  ! Termination for the main program (Internal subroutine)
  !
  call MainTerminate

contains

  !-------------------------------------------------------------------

  subroutine MainInit
    !
    ! 主プログラムの初期化手続き. 
    !
    ! Initialization procedure for the main program. 
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperInit

    !
    ! MPI 計算での Message 制御
    ! Message control for MPI calculation
    !
    use mpi_messagecntl, only : MPIMessageCntlInit

    use dc_message, only: MessageNotify

    ! コマンドライン引数処理
    ! Command line option parser
    !
    use option_parser, only: OptParseInit

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: NmlutilInit, NmlutilMsg

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetInit, TimesetDelTimeHalf, &
      & TimeN                 ! ステップ $ t $ の時刻. Time of step $ t $. 

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    use fileset, only: FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: GridsetInit, &
      &                  imax, & ! 経度格子点数. 
                                 ! Number of grid points in longitude
      &                  jmax, & ! 緯度格子点数. 
                                 ! Number of grid points in latitude
      &                  kmax, & ! 鉛直層数. 
                                 ! Number of vertical level
      &                kslmax, & ! 地下の鉛直層数. 
                                 ! Number of subsurface vertical level
      &                ksimax    ! 海氷の鉛直層数. 
                                 ! Number of sea ice vertical level

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only: CompositionInit


    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetInit

    ! ヒストリデータ出力
    ! History data output
    !
    use history_file_io, only: HistoryFileOpen
    use gtool_historyauto, only: &
      & HistoryAutoAddVariable, HistoryAutoAddAttr, HistoryAutoPut

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: &
      & STDOUT               ! 標準出力の装置番号. Unit number of standard output

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 時系列データの読み込み
    ! Reading time series
    !
    use read_time_series, only : ReadTimeSeriesInit

    !
    ! Calc Mass stream function
    !
    use mass_streamfunc, only :  yr_MassStreamFunc

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level,
    ! and calculate pressure and height
    !
    use auxiliary, only: AuxVarsInit

    ! メッセージ制御
    ! Message control
    !
    use mpi_messagecntl, only : DoesOutputMPIMessage

    ! 宣言文 ; Declaration statements
    !
    implicit none

    character(*), parameter:: version = &
      & '$Name:  $' // &
      & '$Id: dcpam_main_hs94.f90,v 1.67 2026/08/15 00:00:00 takepiro Exp $'
                              ! 主プログラムのバージョン
                              ! Main program version

    character(STRING)      :: namelist_filename
                              ! NAMELIST ファイルの名称. 
                              ! NAMELIST file name

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
                              ! IOSTAT of NAMELIST read

    integer:: j,k             ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituents
    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /tem_test_nml/ &
      & TestMode
          !
          ! デフォルト値については初期化手続 "main/dcpam_main.F90#MainInit" 
          ! のソースコードを参照のこと. 
          !
          ! Refer to source codes in the initialization procedure
          ! "main/dcpam_main.F90#MainInit" for the default values. 
          !

    ! 実行文 ; Executable statement
    !

    ! Initialize MPI
    !
    call MPIWrapperInit

    ! コマンドライン引数処理
    ! Command line option parser
    !
    call OptParseInit(       &
      & namelist_filename,    & ! (out)
      & prog_name            & ! (in )
      & )

    ! NAMELIST ファイル名入力
    ! Input NAMELIST file name
    !
    call NmlutilInit( &
      & namelist_filename  & ! (in)
      & )
    !
    ! MPI 計算での Message 制御
    ! Message control for MPI calculation
    !
    call MPIMEssageCntlInit

    ! 計算モードの設定
    ! Configure calculation mode
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, &            ! (in)
        & nml = tem_test_nml, &    ! (out)
        & iostat = iostat_nml )    ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, prog_name ) ! (in)
      if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = tem_test_nml )
    end if

    ! Initialization of modules used in this module
    !


    ! 時刻管理
    ! Time control
    !
    call TimesetInit

    ! 出力ファイルの基本情報管理
    ! Management basic information for output files
    ! 
    call FilesetInit

    ! 格子点設定
    ! Grid points settings
    !
    call GridsetInit

    ! 物理定数設定
    ! Physical constants settings
    !
    call ConstantsInit

    ! 座標データ設定
    ! Axes data settings
    !
    call AxessetInit

    ! 時系列データの読み込み
    ! Reading time series
    !
    call ReadTimeSeriesInit

    ! 温度の半整数σレベルの補間, 気圧と高度の算出
    ! Interpolate temperature on half sigma level,
    ! and calculate pressure and height
    !
    call AuxVarsInit

    
    ! 予報変数の割付
    ! Allocation of prediction variables
    !
    allocate( xyz_U      (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_V      (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyr_SigDot (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyz_Temp   (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_Entro  (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyz_QVap   (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xy_Ps      (0:imax-1, 1:jmax) )

    allocate( xyz_EPLat       (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyr_EPSig       (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyz_PsVRes      (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyr_PsSigDotRes (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyr_StrRes      (0:imax-1, 1:jmax, 0:kmax) )

    xyz_EPLat       = 0.0D0
    xyr_EPSig       = 0.0D0
    xyz_PsVRes      = 0.0D0
    xyr_PsSigDotRes = 0.0D0
    xyr_StrRes      = 0.0D0

    allocate( xyz_EPLatAns       (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyr_EPSigAns       (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyz_PsVResAns      (0:imax-1, 1:jmax, 1:kmax) )
    allocate( xyr_PsSigDotResAns (0:imax-1, 1:jmax, 0:kmax) )
    allocate( xyr_StrResAns      (0:imax-1, 1:jmax, 0:kmax) )

    xyz_EPLatAns       = 0.0D0
    xyr_EPSigAns       = 0.0D0
    xyz_PsVResAns      = 0.0D0
    xyr_PsSigDotResAns = 0.0D0
    xyr_StrResAns      = 0.0D0

    allocate( yr_Msf    (1:jmax,0:kmax) )
    !
    ! ヒストリデータファイルの初期化
    ! Initialization of history data files
    !
    call HistoryFileOpen
    !
    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'U' , &         ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'eastward wind', 'm s-1' )               ! (in)

    call HistoryAutoAddVariable( 'V' , &         ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'northward wind', 'm s-1' )              ! (in)

    call HistoryAutoAddVariable( 'SigDot' , &    ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &  ! (in)
      & 'vertical wind', 'm s-1' )               ! (in)

    call HistoryAutoAddVariable( 'Temp' , &      ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'temperature', 'K' )                     ! (in)

    call HistoryAutoAddVariable( 'Ps' , &        ! (in)
      & (/ 'lon ', 'lat ', 'time' /), &          ! (in)
      & 'surface pressure', 'Pa' )               ! (in)

    call HistoryAutoAddVariable( 'Msf' , &       ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &  ! (in)
      & 'mass stream function', 'kg s-1' )       ! (in)

    call HistoryAutoAddVariable( 'VRes' , &      ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'residual northward wind', 'm s-1' )     ! (in)

    call HistoryAutoAddVariable( 'SigDotRes', &  ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &  ! (in)
      & 'residual vertical wind', 'm s-1' )      ! (in)

    call HistoryAutoAddVariable( 'StrRes' , &       ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &     ! (in)
      & 'residual mass stream function', 'kg s-1' ) ! (in)

    call HistoryAutoAddVariable( 'EPLat' , &      ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'northward EP flux', 'm s-1' )     ! (in)

    call HistoryAutoAddVariable( 'EPSig', &  ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &  ! (in)
      & 'vertical EP flux', 'm s-1' )      ! (in)

    call HistoryAutoAddVariable( 'VResAns' , &      ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'residual northward wind', 'm s-1' )     ! (in)

    call HistoryAutoAddVariable( 'SigDotResAns', &  ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &  ! (in)
      & 'residual vertical wind', 'm s-1' )      ! (in)

    call HistoryAutoAddVariable( 'StrResAns' , &       ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &     ! (in)
      & 'residual mass stream function', 'kg s-1' ) ! (in)

    call HistoryAutoAddVariable( 'EPLatAns' , &      ! (in)
      & (/ 'lon ', 'lat ', 'sig ', 'time' /), &  ! (in)
      & 'northward EP flux', 'm s-1' )     ! (in)

    call HistoryAutoAddVariable( 'EPSigAns', &  ! (in)
      & (/ 'lon ', 'lat ', 'sigm', 'time' /), &  ! (in)
      & 'vertical EP flux', 'm s-1' )      ! (in)


  end subroutine MainInit

  !-------------------------------------------------------------------

  subroutine MainTerminate
    !
    ! 主プログラムの終了処理手続き. 
    !
    ! Termination procedure for the main program. 
    !

    ! モジュール引用 ; USE statements
    !

    ! MPI
    !
    use mpi_wrapper, only : MPIWrapperFinalize

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: AxessetFinalize

    ! 時刻管理
    ! Time control
    !
    use timeset, only: TimesetClose

    ! ヒストリデータ出力
    ! History data output
    !
    use history_file_io, only: HistoryFileClose

    ! 宣言文 ; Declaration statements
    !
    implicit none

    ! 実行文 ; Executable statement
    !

    ! ヒストリデータファイルクローズ
    ! Close history data files
    !
    call HistoryFileClose

    ! 予報変数の割付解除
    ! Deallocation of prediction variables
    !
    deallocate( xyz_U     )
    deallocate( xyz_V     )
    deallocate( xyz_Temp  )
    deallocate( xyz_QVap )
    deallocate( xy_Ps     )

    deallocate( xyz_EPLat       )
    deallocate( xyr_EPSig       )
    deallocate( xyz_PsVRes      )
    deallocate( xyr_PsSigDotRes )
    deallocate( xyr_StrRes      )

    deallocate( xyz_EPLatAns       )
    deallocate( xyr_EPSigAns      )
    deallocate( xyz_PsVResAns      )
    deallocate( xyr_PsSigDotResAns )
    deallocate( xyr_StrResAns      )

    ! 診断変数の割付解除
    ! Dellocation of diagnostic variables
    !
    deallocate( yr_MSF )

    ! 各モジュール内の変数の割付解除
    ! Dellocation of variables in modules
    !
    call AxessetFinalize

    ! 時刻管理終了処理
    ! Termination of time control
    !
    call TimesetClose

    ! Finalize MPI
    !
    call MPIWrapperFinalize

  end subroutine MainTerminate
  !
  !----------------------------------------------------------------------
  !
  subroutine set_basicstate
    !
    ! データ入力(基本場)
    ! Restart data input
    !
    integer :: j, k
    
    xy_PS    = Ps

    do k=1,kmax
       do j=1,jmax
          xyz_U(:,j,k) = U0*cos(y_Lat(j))*z_Sigma(k)
       end do
    end do

    do k=1,kmax
       xyz_Entro(:,:,k) = CpDry * log(T0) + Gamma * (z_Sigma(k)-1.0D0)
       xyz_Temp(:,:,k)  = z_Sigma(k)**(GasRDry/CpDry) * exp(xyz_Entro(:,:,k)/CpDry)
    end do

    xyz_V = 0.0D0
    xyr_SigDot = 0.0D0

  end subroutine set_basicstate
  !
  !----------------------------------------------------------------------
  !
  subroutine set_zonal_test
    !
    ! データ入力(帯状) : 残差循環 1 項目チェック
    ! Restart data input
    !
    integer j,k

    do k=1,kmax
       do j=1,jmax
          xyz_V(:,j,k) = xyz_V(:,j,k) &
               + V0 * sin(2*y_Lat(j))*cos(PI*z_Sigma(k))
       end do
    end do
    do k=0,kmax
       do j=1,jmax
          xyr_SigDot(:,j,k) = xyr_SigDot(:,j,k) &
               + W0 * (1.0D0 - 3.0D0*sin(y_Lat(j))**2)*cos(PI*r_Sigma(k))
       end do
    end do
  end subroutine set_zonal_test
  !
  !----------------------------------------------------------------------
  !
  subroutine set_disturbance1_test
    !
    ! データ入力(擾乱), EP フラックス 1 項目チェック
    ! Restart data input
    !
    integer j,k

    do k=1,kmax
       do j=1,jmax
          xyz_U(:,j,k) = xyz_U(:,j,k) &
               + U1*sin(2*y_Lat(j))**2 * cos(PI*z_Sigma(k)) * cos(x_Lon)
          xyz_V(:,j,k) = xyz_V(:,j,k) &
               + V1*sin(2*y_Lat(j))*cos(PI*z_Sigma(k)) * cos(x_Lon)
       end do
    end do
    do k=0,kmax
       do j=1,jmax
          xyr_SigDot(:,j,k) = xyr_SigDot(:,j,k) &
               + W1 * cos(y_Lat(j))*sin(PI*r_Sigma(k)) * cos(x_Lon)
       end do
    end do
  end subroutine set_disturbance1_test
  !
  !----------------------------------------------------------------------
  !
  subroutine set_disturbance2_test
    !
    ! データ入力(擾乱) , 残差循環, EP フラックス 2 項目チェック
    ! Restart data input
    !
    integer j,k

    do k=1,kmax
       do j=1,jmax
          xyz_Entro(:,j,k) = xyz_Entro(:,j,k) &
               + s1*cos(y_Lat(j))*cos(2.0D0*x_Lon)
          xyz_Temp(:,j,k) = z_Sigma(k)**(GasRDry/CpDry) * exp(xyz_Entro(:,j,k)/CpDry)
       end do
    end do
    
    do k=1,kmax
       do j=1,jmax
          xyz_V(:,j,k) = xyz_V(:,j,k) &
               + V1*sin(2*y_Lat(j))*cos(PI*z_Sigma(k))* cos(2.0D0*x_Lon)
       end do
    end do
    do k=0,kmax
       do j=1,jmax
          xyr_SigDot(:,j,k) = xyr_SigDot(:,j,k) &
               + W1 * (1.0D0 - 3.0D0*sin(y_Lat(j))**2)*cos(PI*r_Sigma(k)) &
                    * cos(2.0D0*x_Lon)
       end do
    end do
  end subroutine set_disturbance2_test
  !
  !============================== 正解 ==============================
  !
  subroutine set_zonal_answer
    !
    ! 正解(擾乱) : 残差循環 1 項目チェック
    ! Ansser data input
    !
    integer j,k

    do k=1,kmax
       do j=1,jmax
          xyz_PsVResAns(:,j,k) = xyz_PsVResAns(:,j,k) &
               + Ps * V0 * sin(2*y_Lat(j))*cos(PI*z_Sigma(k))
       end do
    end do

    do k=0,kmax
       do j=1,jmax
          xyr_PsSigDotResAns(:,j,k) = xyr_PsSigDotResAns(:,j,k) &
               + Ps * W0 * (1.0D0 - 3.0D0*sin(y_Lat(j))**2)*cos(PI*r_Sigma(k))
       end do
    end do
    
    do k=0,kmax
       do j=1,jmax
          xyr_StrResAns(:,j,k) = xyr_StrResAns(:,j,k) &
               - 2 * Rplanet * Ps * V0/Grav &
                  * sin(2*y_Lat(j))*cos(y_Lat(j))*sin(PI*r_Sigma(k))
       end do
    end do

  end subroutine set_zonal_answer
  !
  !----------------------------------------------------------------------
  !
  subroutine set_disturbance1_answer
    !
    ! 正解(擾乱) EP フラックス 1 項目チェック
    ! Ansser data input
    !
    integer j,k

    do k=1,kmax
       do j=1,jmax
          xyz_EPLatAns(:,j,k) = xyz_EPLatAns(:,j,k) &
               - Ps*U1*Rplanet*V1/2.0D0 &
                  * sin(2*y_Lat(j))**3 * cos(y_Lat(j)) * cos(PI*z_Sigma(k))**2
       end do
    end do

    do k=0,kmax
       do j=1,jmax
          xyr_EPSigAns(:,j,k) = xyr_EPSigAns(:,j,k) &
               - Ps*U1*Rplanet*W1/4.0D0 &
                  * sin(2*y_Lat(j))**2 * cos(y_Lat(j))**2 * sin(2*PI*r_Sigma(k))
       end do
    end do
  end subroutine set_disturbance1_answer
  !
  !----------------------------------------------------------------------
  !
  subroutine set_disturbance2_answer
    !
    ! 正解(擾乱) : 残差循環, EP フラックス 2 項目チェック
    ! Ansser data input
    !
    integer j,k

    do k=1,kmax
       do j=1,jmax
          xyz_PsVResAns(:,j,k) = xyz_PsVResAns(:,j,k) &
               + PI*Ps*s1*V1/Gamma &
                   * sin(y_Lat(j)) * cos(y_Lat(j))**2*sin(PI*z_Sigma(k))
       end do
    end do

    do k=0,kmax
       do j=1,jmax
          xyr_PsSigDotResAns(:,j,k) = xyr_PsSigDotResAns(:,j,k) &
               + Ps*s1*V1/(Gamma*Rplanet) &
                   * cos(y_Lat(j)) * (1.0D0-4.0D0*sin(y_Lat(j))**2)*cos(PI*r_Sigma(k))
       end do
    end do

    do k=0,kmax
       do j=1,jmax
          xyr_StrResAns(:,j,k) = xyr_StrResAns(:,j,k)&
               + 2*PI*Rplanet*Ps*s1*V1/(Gamma*Grav) &
                   * sin(y_Lat(j)) * cos(y_Lat(j))**3*(cos(PI*r_Sigma(k))-1.0D0)
       end do
    end do
    
    do k=1,kmax
       do j=1,jmax
          xyz_EPLatAns(:,j,k) = xyz_EPLatAns(:,j,k)&
               + U0*Rplanet*Ps*s1*V1/Gamma &
                   * sin(y_Lat(j))* cos(y_Lat(j))**4 * cos(PI*z_Sigma(k))
       end do
    end do

    do k=0,kmax
       do j=1,jmax
          xyr_EPSigAns(:,j,k) = xyr_EPSigAns(:,j,k) &
               + 2*Ps*s1*V1/Gamma &
                   * sin(y_Lat(j))**2 * cos(y_Lat(j))**3 &
                   * (U0*r_Sigma(k) + Omega*Rplanet) * cos(PI*r_Sigma(k))
       end do
    end do
  end subroutine set_disturbance2_answer
  !
  !============================== チェック ==============================
  !
  subroutine CheckTEM
    use dc_test, only : AssertEqual

    ! 判定誤差設定
    integer, parameter :: check_digits = 2
    integer :: ignore = -3

    ignore = ignore_xya(xyz_PsVResAns,check_digits)
    call AssertEqual(&
         message='PsVRes',                                          &
         answer = xyz_PsVRes,                                       &
         check = xyz_PsVResAns,                                     &
         significant_digits = check_digits, ignore_digits = ignore  &
         )

    ignore = ignore_xya(xyr_PsSigDotResAns,check_digits)
    call AssertEqual(&
         message='PsSigDotRes',                                     &
         answer = xyr_PsSigDotRes,                                  &
         check = xyr_PsSigDotResAns,                                &
         significant_digits = check_digits, ignore_digits = ignore  &
         )

    ignore = ignore_xya(xyr_StrResAns,check_digits)
    call AssertEqual(&
         message='StrRes',                                          &
         answer = xyr_StrRes,                                       &
         check = xyr_StrResAns,                                     &
         significant_digits = check_digits, ignore_digits = ignore  &
         )

    ignore = ignore_xya(xyz_EPLatAns,check_digits)
    call AssertEqual(&
         message='EPLat',                                          &
         answer = xyz_EPLat,                                       &
         check = xyz_EPLatAns,                                     &
         significant_digits = check_digits, ignore_digits = ignore  &
         )

    ignore = ignore_xya(xyr_EPSigAns,check_digits)
    call AssertEqual(&
         message='EPSig',                                          &
         answer = xyr_EPSig,                                       &
         check = xyr_EPSigAns,                                     &
         significant_digits = check_digits, ignore_digits = ignore  &
         )

  end subroutine CheckTEM

  function ignore_xya(xya_Data,check) result(ignore)
    real(DP), intent(IN) :: xya_Data(:,:,:)
    integer,  intent(IN) :: check
    integer              :: ignore
    

    real(DP), parameter :: Epsilon = 1.0D-10  ! 微小パラメター
    integer, parameter  :: ignore_zero = -6   ! 0 field での判定基準

    if ( maxval(abs(xya_Data)) < Epsilon ) then
       ignore = ignore_zero
    else
       ignore = int(log10(maxval(abs(xya_Data))))-check
    end if
    
  end function ignore_xya
  
end program dcpam_tem_test