!--------------------------------------------------------------------------
!  Copyright (C) 2002--2022 SPMODEL Development Group. All rights reserved.
!--------------------------------------------------------------------------
!
!表題  回転球殻ブシネスク流体モデル
!
!      * 球面調和函数変換 + チェビシェフ選点法
!      * 時間積分 Adams-Bashforth + Crank Nicholson スキーム
!      * ダイナモベンチマーク用設定
!
!履歴  2002/11/01  佐々木洋平  ベンチマーク Case 1 リスタート可能バージョン
!      2002/11/02  佐々木洋平  spmodel 書法に従う修正, 読み込み初期値のポテンシャル化
!      2002/11/19  佐々木洋平  初期時刻の読み込みの修正
!      2002/11/19  佐々木洋平  境界条件をタウ法から変更
!      2004/02/16  佐々木 洋平 gt4f90io の使用に伴うライブラリ名変更
!      2004/11/26  竹広 真一   Adams-Bashforth scheme 用の各時間変化量を
!                              リスタート変数に追加
!      2004/11/27  竹広 真一   内部初期条件の設定ルーチン追加.
!      2004/12/30  竹広 真一   HistoryCreate に xtype を追加
!                              (Alpha machine での出力時エラーへの暫定対処).
!                              結果を出力した時間を表示.
!      2005/01/03  竹広 真一   解像度追加, 外部ファイルで与えた初期値は出力せず
!      2005/01/04  竹広 真一   HistoryCreate のオプション xtypes typo 修正
!
!      2005/01/06  竹広 真一   ベンチマークプログラムより一般化
!                              球殻半径, 物理パラメター, 境界条件を namelist 化
!                              変数読み込みを HistoryGet で行う.
!      2005/01/07  竹広 真一   Namelist 変数等を属性として出力.
!      2005/01/19  竹広 真一   半径比で球殻内外半径を設定可能にする.
!      2005/01/20  竹広 真一   渦度, 電流, ヘリシティーの計算追加.
!      2005/01/25  竹広 真一   タイトルを Namelist で設定可能にする.
!      2005/02/04  竹広 真一   磁場拡散項の係数 1/Pm が抜けていたのを修正.
!      2005/02/11  竹広 真一   係数 1/Pm の修正再度.
!      2005/04/17  竹広 真一   解像度を Namelist で指定
!                              エネルギースペクトルの出力追加
!      2005/04/29  竹広 真一   非線形項を V・▽V から (▽xV)xV に変更
!                              渦度の計算をポテンシャルから直接行なう.
!      2005/05/09  竹広 真一   初期値のタイプ(case0/case1)のスイッチ導入
!      2005/05/19  竹広 真一   Clank Nicolson 第 1 段後に境界条件計算を追加
!      2006/02/19  竹広 真一   Crank-Nicolson 行列・計算変更
!      2006/02/21  竹広 真一   ソース名(日付)をログに出力
!      2006/02/28  竹広 真一   RDoc 用にコメント変更, ファイル改名
!      2007/06/26  竹広 真一   lon_weight, lat_weight 出力
!      2007/09/13  竹広 真一   時間積分前処理を変更
!      2007/12/19  竹広 真一   温度 Crank-Nicolson 行列を修正
!      2007/12/20  竹広真一    コマンドライン引数を解析し,
!                              NAMELIST ファイル名を取得
!      2008/01/25  竹広 真一   磁場エネルギー計算 bug fix
!      2008/06/11  竹広 真一   OPENMP スイッチ導入
!      2008/08/09  竹広 真一   経度添え字を im->0:im-1 に変更
!      2011/09/16  佐々木 洋平 gt4_hisyory -> gtool_history
!      2012/06/03  竹広 真一   Crank-Nicoslon 行列省メモリ版
!      2012/08/24  竹広 真一   wt_LaplaPol2Pol_wt を使用
!      2013/02/19  竹広 真一   wt_module_sjpack を使用
!      2013/03/09  竹広 真一   スケーリング変更
!      2013/03/16  竹広 真一   MPI 並列版
!      2013/03/17  竹広 真一   磁場なし版
!      2013/05/21  竹広 真一   コリオリ項を実空間計算へ
!      2013/07/14  竹広 真一   LU 行列大規模に, 出力スイッチ導入
!      2013/08/09  竹広 真一   スペクトル MPI 分割
!      2013/09/06  竹広 真一   超粘性導入
!      2013/09/18  竹広 真一   ランダム温度擾乱初期値導入
!      2014/03/12  竹広 真一   wt_mpi_module_nm 用
!      2014/03/21  竹広 真一   wt_mpi_module_nmomp 用
!      2014/03/22  竹広 真一   output_diag bug fix
!      2017/10/28  竹広 真一   ut_mpi_module_nmomp 用
!      2018/12/24  竹広 真一   nm_l 引数修正
!      2019/04/30  竹広 真一   水平最高波数成分除去
!      2019/11/19  竹広 真一   スペクトル診断出力
!      2020/06/03  竹広 真一   sypack/ISPACK3 用に変更
!      2020/07/07  竹広 真一   dc_messege MPI 出力コントロールスイッチ導入
!      2020/11/12  竹広 真一   セクター計算スイッチ導入
!      2022/12/22  竹広 真一   ektor_nmz, ekpol_nmz 出力復活
!                              ektoreqs_nmz ektoreqa_nmz 出力追加
!
!program bsnsq_rot_3d_shell_ut_mpi_module_sypack_abcn_hdiff_Omega_mint
program bsnsq_rot_3d_shell
  !
  !  回転球殻ブシネスク流体モデル
  !
  !   * 球面調和函数変換 + チェビシェフガラーキン法
  !   * 時間積分 Adams-Bashforth + Crank Nicholson スキーム
  !   * ダイナモベンチマーク用設定
  !
  use lumatrix
  use at_module,  only : a_Int_ag
  use ut_mpi_module
  use w_mpi_module_base, only : js => j1, je => j2
  use dc_trace,   only : SetDebug, BeginSub, EndSub, DbgMessage
  use dc_message, only : MessageNotify, MessageSuppressMPI
  use dc_types
  use dc_string,   only : StoA, Printf
  use dc_args
  use gtool_history
  use mpi
  implicit none

!== 宣言部 ============================================================

 !---- 変数 ----
  real(8), allocatable :: xvh_Torvel(:,:,:)      ! トロイダル速度
  real(8), allocatable :: ut_Torvel(:,:)         ! トロイダル速度
  real(8), allocatable :: xvh_Polvel(:,:,:)      ! ポロイダル速度
  real(8), allocatable :: ut_Polvel(:,:)         ! ポロイダル速度
  real(8), allocatable :: xvh_LaplaPolvel(:,:,:) ! ポロイダル速度
  real(8), allocatable :: ut_LaplaPolvel(:,:)    ! ▽^2ポロイダル速度

  real(8), allocatable :: xvh_Temp(:,:,:)       ! 温度
  real(8), allocatable :: ut_Temp(:,:)          ! 温度

  real(8), allocatable :: xvh_Vlon(:,:,:)       ! 速度(経度)
  real(8), allocatable :: xvh_Vlat(:,:,:)       ! 速度(緯度)
  real(8), allocatable :: xvh_Vrad(:,:,:)       ! 速度(動径)

  real(8), allocatable :: xvh_Zlon(:,:,:)       ! 渦度(経度)
  real(8), allocatable :: xvh_Zlat(:,:,:)       ! 渦度(緯度)
  real(8), allocatable :: xvh_Zrad(:,:,:)       ! 渦度(動径)

  real(8), allocatable :: xva_TempBndry(:,:,:)  ! 温度境界値
  real(8), allocatable :: ua_TempBndry(:,:)     ! 温度境界値スペクトル

  real(8), allocatable :: nz_EkTor(:,:)   ! 運動エネルギー(トロイダル成分)
  real(8), allocatable :: nz_EkPol(:,:)   ! 運動エネルギー(ポロイダル成分)

  real(8), allocatable :: nz_EnsTor(:,:)  ! エンストロフィー(トロイダル成分)
  real(8), allocatable :: nz_EnsPol(:,:)  ! エンストロフィー(ポロイダル成分)

 !---- 超拡散係数 ----
  real(8), allocatable :: ut_HVisc(:,:)   ! 超粘性
  real(8), allocatable :: ut_HTDiff(:,:)  ! 超熱拡散

 !---- 時間積分用変数 ----

  real(8), allocatable :: ut_DtDTorvel1(:,:), ut_DtDTorvel0(:,:)
                          ! トロイダル速度時間変化
  real(8), allocatable :: ut_DtDLaplaPolvel1(:,:), ut_DtDLaplaPolvel0(:,:)
                          ! ▽^2ポロイダル速度時間変化
  real(8), allocatable :: ut_DtDTemp1(:,:), ut_DtDTemp0(:,:)
                          ! 温度時間変化

  real(8), allocatable :: utt_DifLUMT_Torvel(:,:,:)
                          ! Crank Nicholson 陰的計算行列(トロイダル速度場)
  integer, allocatable :: ut_kpivot_Torvel(:,:)
                          ! Crank Nicholson 行列ピボット(トロイダル速度場)

  real(8), allocatable :: utt_DifLUMT_Polvel(:,:,:)
                          ! Crank Nicholson 陰的計算行列
                          ! (ポロイダル速度場)
  integer, allocatable :: ut_kpivot_Polvel(:,:)
                          ! Crank Nicholson 行列ピボット(ポロイダル速度場)

  real(8), allocatable :: utt_DifLUMT_Temp(:,:,:)
                          ! Crank Nicholson 陰的計算行列(温度場)
  integer, allocatable :: ut_kpivot_Temp(:,:)    ! ピボット
                          ! Crank Nicholson 陰的計算行列ピボット(温度場)

 !---- その他 ----
  real(8), parameter :: pi=3.1415926535897932385D0
  character(len=20)  :: DbgMessageFmt='*** DbgMESSAGE ***  '
  real(8), parameter  :: vmiss = -999.0       ! 欠損値
  integer :: it=0
  real(8) :: time
  integer :: nmp, myid, ierr                  ! MPI 用変数
  integer :: n, m, l, nmary(2)
  character(len=STRING) :: cip                ! ip 文字変数
  type(GT_HISTORY) :: hst_rst                 ! リスタート GT_HISTORY 変数

 !---- NAMELIST 変数 ----
  ! コマンドライン引数用変数
  type(ARGS) :: arg
  character(STRING), pointer :: argv(:) => null()

  ! NAMELIST 入力用デフォルトファイル名
  character(STRING)  :: nmlfile='bsnsq_rot-3d-shell_ut_mpi_module_sypack_abcn_hdiff_Omega_mint.nml'

  logical            :: verbose=.false.          ! 出力メッセージレベル
  logical            :: DebugOn=.false.          ! デバッグ出力コントロール
  integer            :: Output_Rank=-1           ! dc_message MPI 出力rank
                                                 ! 負のときすべての rank の
                                                 ! メッセージ出力
  namelist /message/  Verbose, DebugOn, Output_Rank  !

  integer  :: im=64, jm=32, km=16                ! 格子点(経度, 緯度, 動径)
  integer  :: nm=21, lm=16                       ! 切断波数(水平, 動径)
  integer  :: omp=1                              ! OPENMP プロセス数
  integer  :: mint=1                             ! 東西波数間隔
  integer  :: npv=1                              ! 鉛直プロセス分割数
  namelist /gridset/ im, jm, km, nm, lm, omp, mint, npv !-- 格子点・スペクトル --

  real(8)            :: ri=0.0D0                 ! 内半径
  real(8)            :: ro=0.0D00                ! 外半径
  real(8)            :: eta=0.0D0                ! 内外半径比
  namelist /radius/  ri, ro, eta                 !-- 座標変数など --

  character(STRING)  :: initial_file=''          ! 初期値データファイル名
  real               :: initial_time=0.0         ! 初期時刻
  namelist /initial/ initial_file, initial_time  ! -- 初期値 --

  real(8) :: delta_t=1.0e-7                      ! 時間積分刻み
  integer :: nstep=2000                          ! 時間積分ステップ数
  namelist /tint/    delta_t, nstep              ! -- 時間積分 --

  character(len=3)  :: output_mode="ALL"         ! 出力モードスイッチ
                                                 ! ALL: 全て出力
                                                 ! HST: ヒストリーのみ出力
                                                 ! RST: リスタートのみ出力
                                                 ! NON: ファイル出力なし
  namelist /output/ output_mode

  character(STRING) :: hst_file=''               ! 出力ファイル名
  character(STRING) :: title = &                 ! タイトル
            'Thermal convection of Boussinesq fluid in a 3-dim rotationg spherical shell'!
  character(STRING) :: source = &                ! ソース名
            'bsnsq_rot_3d_shell_ut_mpi_module_sypack_abcn_hdiff_Omega_mint.f90 (2022/12/22)'!
  integer           :: hst_intstep=100           ! 出力間隔ステップ数
  namelist /history/ hst_file, title, hst_intstep   ! -- ヒストリー出力 --

  character(STRING) :: rst_file=''               ! リスタート出力ファイル名
  integer           :: rst_intstep=200           ! リスタート出力間隔ステップ数
  namelist /restart/   rst_file, rst_intstep     ! -- リスタート出力 --

  real(8)            :: Ra=100.0D0               ! レイリー数
  real(8)            :: Pr=1.0D0                 ! プランドル数
  real(8)            :: Ekman=1.0d-3             ! エクマン数
  real(8)            :: Ektherm=1.0d-3           ! 熱エクマン数
  real(8)            :: Pm=5.0D0                 ! 磁気プランドル数
  namelist /physics/  Ra, Pr, Ekman, Pm          ! -- 物理パラメター --

  character(len=2)   :: VelBC='RR'               ! 速度境界条件(RR/RF/FR/FF)
  character(len=2)   :: TempBC='DD'              ! 温度境界条件(DD/DN/ND/NN)
  real(8)            :: Temptop=0.0D0            ! 境界温度(上端)
  real(8)            :: Tempbottom=1.0D0         ! 境界温度(下端)
  namelist /boundary/ VelBC, TempBC, Temptop, Tempbottom
                                                 ! -- 境界条件 --
 !---- 超粘性・拡散 ----
  real(8)            :: epsilon = 0.0d0          ! \nu=\nu_0
  integer            :: n0      = 10             ! \nu=\nu_0((1+epsilon(n-n0)**2
  namelist /hyperdiff/  epsilon, n0              ! -- 超拡散係数設定 --

 !---- 初期値の種類 ----
  character(len=1)  :: type  = 'P'               ! 内部初期値タイプ
                                                 ! P : 一点温度擾乱
                                                 ! R : ランダム温度擾乱
                                                 ! B : ベンチマーク
  real(8)           :: amp = 1.0D-2              ! P,R での擾乱振幅
  integer           :: seed = 1                  ! R でのランダム種
  namelist /inittype/  type, amp, seed

 !---------------- MPI スタート ---------------------
  call MPI_INIT(IERR)
  call MPI_COMM_RANK(MPI_COMM_WORLD,MYID,IERR)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,NMP,IERR)
  ! call AMPINT

 !---------------- NAMELIST 読み込み ---------------------
  call Open(arg)
  call Debug(arg) ; call Help(arg) ; call Strict(arg)
  call Get(arg, argv)
  if ( size(argv) .le. 0 ) then
     call MessageNotify('W','main',&
          'Usage: bsnsq_rot_3d_shell_ut_mpi_module_extomp_abcn_hdiff_Omega_ES_mint.out [nmlfilename]')
     call MessageNotify('E','main','There is no argument. ')
  else
     nmlfile=argv(1)
!!$     call MessageNotify('M','main','Namelist file is '//trim(nmlfile))
  endif
  deallocate(argv)
  call Close(arg)

  !write(6,nml=message)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=message) ; if( myid==0 ) write(6,nml=message) ; close(10)
  call MessageSuppressMPI( output_rank )
  call MessageNotify('M','main', trim(source))


  if(verbose) write(6,nml=gridset)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=gridset)  ; if( myid==0 ) write(6,nml=gridset) ; close(10)

  if (verbose) write(6,nml=radius)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=radius)
   !
   ! ri,ro か eta のどちらかを NAMELIST radius で与える. eta の設定が優先される
   !
  if ( eta /= 0.0D0 ) then
     ri=eta/(1.0D0-eta) ; ro=1.0D0/(1.0D0-eta)
     call MessageNotify('M','main', &
             'The inner/outer radii were calculatied with their ratio, eta.')
  else if ( ro /= 0.0D0 .and. ro /= 0.0D0 ) then
     eta = ri/ro
     call MessageNotify('M','main', &
             'The ratio of inner/outer radii were calculatied with ro and ri.')
  else
     call MessageNotify('E','main', &
                 'The inner/outer radii or their ratio should be given')
  endif
  if( myid==0 ) write(6,nml=radius)  ; close(10)

  if(verbose)write(6,nml=initial)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=initial) ; if( myid==0 ) write(6,nml=initial) ; close(10)
  call Printf(unit=cip,fmt="%04d",i=(/myid/))
  if ( initial_file /= '' ) then
     initial_file=trim(initial_file)//'-ip'//trim(adjustl(cip))//'.nc'
  endif

  if(verbose)write(6,nml=tint)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=tint)     ; if( myid==0 ) write(6,nml=tint)    ; close(10)

  if(verbose)write(6,nml=output)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=output)   ; if( myid==0 ) write(6,nml=output)  ; close(10)

  if(verbose)write(6,nml=history)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=history)   ; if( myid==0 ) write(6,nml=history)  ; close(10)
  hst_file=trim(hst_file)//'-ip'//trim(adjustl(cip))//'.nc'

  if(verbose)write(6,nml=restart)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=restart)   ; if( myid==0 ) write(6,nml=restart)  ; close(10)
  rst_file=trim(rst_file)//'-ip'//trim(adjustl(cip))//'.nc'

  if(verbose)write(6,nml=physics)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=physics)  ; if( myid==0 ) write(6,nml=physics) ; close(10)
  Ektherm = Ekman/Pr

  if(verbose)write(6,nml=boundary)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=boundary) ; if( myid==0 ) write(6,nml=boundary); close(10)

  if(verbose)write(6,nml=hyperdiff)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=hyperdiff) ; if( myid==0 ) write(6,nml=hyperdiff); close(10)

  if(verbose)write(6,nml=inittype)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=inittype) ; if( myid==0 ) write(6,nml=inittype); close(10)

 !---------------- デバッグ出力制御設定 -----------------
  if (DebugOn) then
    call SetDebug
  end if

 !---------------- 座標値の設定 ---------------------
  call DbgMessage(fmt='call %c', c1='ut_Initial')

  call ut_mpi_Initial(im,jm,km,nm,lm,ri,ro,npv=npv,mint=mint)

!!$  ut_VMiss = vmiss

 !------------------ 変数の割り付け ---------------------
  allocate(xvh_Torvel(0:im-1,jc,kc),ut_Torvel(nmc,0:lm))
  allocate(xvh_Polvel(0:im-1,jc,kc),ut_Polvel(nmc,0:lm))
  allocate(xvh_LaplaPolvel(0:im-1,jc,kc),ut_LaplaPolvel(nmc,0:lm))

  allocate(xvh_Temp(0:im-1,jc,kc),ut_Temp(nmc,0:lm))

  allocate(xvh_Vlon(0:im-1,jc,kc))        ! 速度(経度)
  allocate(xvh_Vlat(0:im-1,jc,kc))        ! 速度(緯度)
  allocate(xvh_Vrad(0:im-1,jc,kc))        ! 速度(動径)

  allocate(xvh_Zlon(0:im-1,jc,kc))        ! 渦度(経度)
  allocate(xvh_Zlat(0:im-1,jc,kc))        ! 渦度(緯度)
  allocate(xvh_Zrad(0:im-1,jc,kc))        ! 渦度(動径)

  allocate(xva_TempBndry(0:im-1,jc,2))      ! 温度境界値
  allocate(ua_TempBndry(nmc,2))        ! 温度境界値スペクトル

  allocate(nz_EkTor(0:nm,0:km))         ! 運動エネルギー(トロイダル成分)
  allocate(nz_EkPol(0:nm,0:km))         ! 運動エネルギー(ポロイダル成分)

  allocate(nz_EnsTor(0:nm,0:km))        ! エンストロフィー(トロイダル成分)
  allocate(nz_EnsPol(0:nm,0:km))        ! エンストロフィー(ポロイダル成分)

  allocate(ut_HVisc(nmc,0:lm))       ! 超粘性
  allocate(ut_HTDiff(nmc,0:lm))      ! 超熱拡散

 ! トロイダル速度時間変化
  allocate(ut_DtDTorvel1(nmc,0:lm),ut_DtDTorvel0(nmc,0:lm))

 ! ▽^2ポロイダル速度時間変化
  allocate(ut_DtDLaplaPolvel1(nmc,0:lm),ut_DtDLaplaPolvel0(nmc,0:lm))

 ! 温度時間変化
  allocate(ut_DtDTemp1(nmc,0:lm),ut_DtDTemp0(nmc,0:lm))

 ! Crank Nicholson 陰的計算用拡散行列(トロイダル速度場)
  allocate(utt_DifLUMT_Torvel(nmc,0:lm,0:lm))
  allocate(ut_kpivot_Torvel(nmc,0:lm))          ! ピボット

 ! Crank Nicholson 陰的計算用拡散行列(ポロイダル速度場)
  allocate(utt_DifLUMT_Polvel(nmc,0:lm,0:lm))
  allocate(ut_kpivot_Polvel(nmc,0:lm))          ! ピボット

 ! Crank Nicholson 陰的計算用拡散行列(温度場)
  allocate(utt_DifLUMT_Temp(nmc,0:lm,0:lm))
  allocate(ut_kpivot_Temp(nmc,0:lm))          ! ピボット

 !------------------- 初期値設定 ----------------------
  if ( initial_file == "") then
     select case( type )
     case ('B')
        call set_initial_bench
        ! ut_Torvel, ut_Polvel, ut_LaplaPolvel, ut_Temp を設定する(ベンチマーク)
     case ('P')
        call set_initial_pointwise
        ! ut_Torvel, ut_Polvel, ut_LaplaPolvel, ut_Temp を設定する(温度1点擾乱)
     case ('R')
        call set_initial_random
        ! ut_Torvel, ut_Polvel, ut_LaplaPolvel, ut_Temp を設定する(ランダム温度擾乱)
     case default
        call MessageNotify('E','main','type should be B/P/R.')
     end select
  else
     call initial_read_restart
       ! ut_Torvel, ut_Polvel, ut_LaplaPolvel, ut_Temp,
       ! ut_DtDTorvel1, ut_DtDLaplaPolvel1, ut_DtDTemp1 を読み込む
  endif

 !------------------- 境界値設定 ----------------------
  xva_Tempbndry(:,:,1) = Temptop
  xva_Tempbndry(:,:,2) = Tempbottom

  ua_TempBndry(:,1) = u_xv(xva_TempBndry(:,:,1))
  ua_TempBndry(:,2) = u_xv(xva_TempBndry(:,:,2))

  !------------------- 超拡散係数設定 ----------------------
  do l=1,nmc
     nmary = nm_l(l+nms-1)
     if ( nmary(1) .ge. n0 ) then
        ut_HVisc(l,:)  = Ekman * ( 1.0D0 + epsilon*(nmary(1)-n0)**2 )
        ut_HTDiff(l,:) = Ektherm * ( 1.0D0 + epsilon*(nmary(1)-n0)**2 )
     else if ( nmary(1) .ge. 0 ) then
        ut_HVisc(l,:)  = Ekman
        ut_HTDiff(l,:) = Ektherm
     endif
  end do

 !------------------- Crank-Nicoslon 用行列設定 ----------------------
  call CNDiffusionMatrixTorvel( ut_HVisc, delta_t, utt_DifLUMT_Torvel, ut_kpivot_Torvel )
  call CNDiffusionMatrixPolvel( ut_HVisc, delta_t, utt_DifLUMT_Polvel, ut_kpivot_Polvel )
  call CNDiffusionMatrixTemp( ut_HTDiff, delta_t, utt_DifLUMT_Temp,   ut_kpivot_Temp )

 !------------------- 時間積分開始の前処理 --------------------
  time = initial_time

  if ( initial_file == "") then
     call ut_TorBoundaries(ut_Torvel,cond=velBC)
     ut_Polvel = ut_LaplaPol2Pol_ut(ut_LaplaPolvel,cond=velBC)
     ut_LaplaPolvel = ut_Lapla_ut(ut_Polvel)
     call ut_Boundaries(ut_Temp,ua_TempBndry,cond=TempBC)
  endif

  xvh_temp = xvh_ut(ut_temp)
  call ut_Potential2VectorMPI(xvh_Vlon,xvh_Vlat,xvh_Vrad,ut_Torvel,ut_Polvel)
  call ut_Potential2VectorMPI &
       (xvh_ZLon,xvh_ZLat,xvh_ZRad,-ut_LaplaPolvel,ut_Torvel)

  if ( initial_file == "") then
     call TimeDerivetives_noDiffusion( &
          ut_Torvel, ut_Polvel, ut_LaplaPolvel, ut_Temp, &
          xvh_Vlon, xvh_Vlat, xvh_Vrad, &
          xvh_Zlon, xvh_Zlat, xvh_Zrad, &
          xvh_Temp, &
          ut_DtDTorvel1, ut_DtDLaplaPolvel1, ut_DtDTemp1 )
  endif

  if ( output_mode == "ALL" ) then
     call output_restart_init
     call output_history_init
  else if ( output_mode == "HST" ) then
     call output_history_init
  else if ( output_mode == "RST" ) then
     call output_restart_init
  endif

  if ( initial_file == '' ) then    ! 内部で与えた初期値は出力
     if ( output_mode=="ALL" .OR. output_mode=="HST" ) then
        call output_history
     else
        call output_diag
     endif
  endif

  !------------------- 時間積分(Adams-Bashforth 法) --------------------
  do it=1,nstep
     time = initial_time + it * delta_t

     ! 昔の時間変化項を保存
     ut_DtDTorvel0 = ut_DtDTorvel1
     ut_DtDLaplaPolvel0 = ut_DtDLaplaPolvel1
     ut_DtDTemp0   = ut_DtDTemp1


     call TimeDerivetives_noDiffusion( &
          ut_Torvel, ut_Polvel, ut_LaplaPolvel, ut_Temp, &
          xvh_Vlon, xvh_Vlat, xvh_Vrad, &
          xvh_Zlon, xvh_Zlat, xvh_Zrad, &
          xvh_Temp, &
          ut_DtDTorvel1, ut_DtDLaplaPolvel1, ut_DtDTemp1 )

     ! トロイダル速度場
     ut_Torvel = ut_Torvel &
          + delta_t *( 3/2.0*ut_DtDTorvel1 - 1/2.0*ut_DtDTorvel0 ) &
          + delta_t/2.0 * ut_HVisc * ut_Lapla_ut(ut_Torvel)
     call ut_TorBoundaries(ut_Torvel,cond=velBC)
     ut_Torvel = LUSolve(utt_DifLUMT_Torvel,ut_kpivot_Torvel,ut_Torvel)

     ! ポロイダル速度場
     ut_LaplaPolvel = ut_LaplaPolvel &
           + delta_t *(3/2.0*ut_DtDLaplaPolvel1 - 1/2.0*ut_DtDLaplaPolvel0 )&
           + delta_t/2.0 * ut_HVisc * ut_Lapla_ut(ut_LaplaPolvel)
     ut_Polvel = ut_LaplaPol2Pol_ut(ut_LaplaPolvel,cond=velBC)
     ut_Polvel = LUSolve(utt_DifLUMT_Polvel,ut_kpivot_Polvel,ut_Polvel)
     ut_LaplaPolvel = ut_Lapla_ut(ut_Polvel)

     ! 温度場
     ut_Temp = ut_Temp &
          + delta_t *( 3/2.0*ut_DtDTemp1 - 1/2.0*ut_DtDTemp0 )&
          + delta_t/2.0 * ut_HTdiff * ut_Lapla_ut(ut_Temp)
     call ut_Boundaries(ut_Temp,ua_TempBndry,cond=TempBC)
     ut_Temp = LUSolve(utt_DifLUMT_Temp,ut_kpivot_Temp,ut_Temp)

     ! 高波数除去
     do l=1,nmc
        nmary = nm_l(l+nms-1)
        if ( nmary(1) == nm ) then
           ut_Torvel(l,:) = 0.0D0
           ut_Polvel(l,:) = 0.0D0
           ut_LaplaPolvel(l,:) = 0.0D0
           ut_Temp(l,:) = 0.0D0
        endif
     end do

     ! グリッド値計算
     xvh_temp = xvh_ut(ut_temp)
     call ut_Potential2VectorMPI(xvh_Vlon,xvh_Vlat,xvh_Vrad,ut_Torvel,ut_Polvel)
     call ut_Potential2RotationMPI &
          (xvh_ZLon,xvh_ZLat,xvh_ZRad,ut_Torvel,ut_Polvel)

     if(mod(it,hst_intstep) .eq. 0 ) then ! ヒストリー出力
        if ( output_mode=="ALL" .OR. output_mode=="HST" ) then
           call output_history
        else
           call output_diag
        endif
     endif

     if(mod(it,rst_intstep) .eq. 0 ) then  ! リスタート出力
        if ( output_mode=="ALL" .OR. output_mode=="RST" ) then
           call output_restart
        endif
     endif

  enddo

  if ( output_mode=="ALL" ) then
     call output_history_close
     call output_restart_close
  else if ( output_mode=="HST" ) then
     call output_history_close
  else if ( output_mode=="RST" ) then
     call output_restart_close
  endif

 !------ MPIの終了 ------

  call MPI_FINALIZE(IERR)
  !call AMPFNL

contains

! 以上 メインプログラム
!-----------------------------------------------------------------------------
! 以下 サブルーチン
 !------------------- 時間変化項(拡散項以外) ----------------------
  subroutine TimeDerivetives_noDiffusion( &
          ut_Torvel, ut_Polvel, ut_LaplaPolvel, ut_Temp, &
          xvh_Vlon, xvh_Vlat, xvh_Vrad, &
          xvh_Zlon, xvh_Zlat, xvh_Zrad, &
          xvh_Temp, &
          ut_DtDTorvel, ut_DtDLaplaPolvel, ut_DtDTemp )
    !
    ! 時間変化項(拡散項以外)の計算を行う.
    !
    real(8), intent(IN) :: ut_Torvel(:,:)        ! トロイダル速度
    real(8), intent(IN) :: ut_Polvel(:,:)        ! ポロイダル速度
    real(8), intent(IN) :: ut_LaplaPolvel(:,:)   ! ▽^2 ポロイダル速度
    real(8), intent(IN) :: ut_Temp(:,:)          ! 温度

    real(8), intent(IN)     :: xvh_Vlon(:,:,:)   ! 速度(経度)
    real(8), intent(IN)     :: xvh_Vlat(:,:,:)   ! 速度(緯度)
    real(8), intent(IN)     :: xvh_Vrad(:,:,:)   ! 速度(動径)

    real(8), intent(IN)     :: xvh_Zlon(:,:,:)   ! 渦度(経度)
    real(8), intent(IN)     :: xvh_Zlat(:,:,:)   ! 渦度(緯度)
    real(8), intent(IN)     :: xvh_Zrad(:,:,:)   ! 渦度(動径)

    real(8), intent(IN)     :: xvh_Temp(:,:,:)   ! 温度

    real(8), intent(OUT) :: ut_DtDTorvel(:,:)   ! トロイダル速度時間変化項
    real(8), intent(OUT) :: ut_DtDLaplaPolvel(:,:)
                                                ! ▽^2 ポロイダル速度時間変化項
    real(8), intent(OUT) :: ut_DtDTemp(:,:)     ! 温度時間変化項

    !--- 作業変数 ---
    real(8), dimension(size(xvh_Vlon,1),size(xvh_Vlon,2),size(xvh_Vlon,3)) &
         :: xvh_KxVlon   ! コリオリ項(経度)
    real(8), dimension(size(xvh_Vlat,1),size(xvh_Vlat,2),size(xvh_Vlat,3)) &
         :: xvh_KxVLat   ! コリオリ項(緯度)
    real(8), dimension(size(xvh_Vrad,1),size(xvh_Vrad,2),size(xvh_Vrad,3)) &
         :: xvh_KxVRad   ! コリオリ項(動径)
 
    real(8), dimension(size(xvh_Vlon,1),size(xvh_Vlon,2),size(xvh_Vlon,3)) &
         :: xvh_ZxVLon   ! (▽xV)xV (経度)
    real(8), dimension(size(xvh_Vlat,1),size(xvh_Vlat,2),size(xvh_Vlat,3)) &
         :: xvh_ZxVLat   ! (▽xV)xV(緯度)
    real(8), dimension(size(xvh_Vrad,1),size(xvh_Vrad,2),size(xvh_Vrad,3)) &
         :: xvh_ZxVRad   ! (▽xV)xV(動径)

    ! コリオリ項計算
    xvh_KxVLon = - xvh_VLat*sin(xvh_Lat) + xvh_VRad*cos(xvh_Lat)
    xvh_KxVLat =   xvh_VLon*sin(xvh_Lat)
    xvh_KxVRad = - xvh_VLon*cos(xvh_Lat)
 
    ! (▽xv)xv 非線形項計算
    xvh_ZxVLon = xvh_ZLat * xvh_VRad - xvh_ZRad * xvh_VLat
    xvh_ZxVLat = xvh_ZRad * xvh_VLon - xvh_ZLon * xvh_VRad
    xvh_ZxVRad = xvh_ZLon * xvh_VLat - xvh_ZLat * xvh_VLon

    ! トロイダル時間変化
    ut_DtDTorvel = &
         - ut_L2Inv_ut(ut_RadRot_xvh_xvh( &
                          2*xvh_KxVLon + xvh_ZxVLon, &
                          2*xvh_KxVLat + xvh_ZxVLat ))
 
    ! ポロイダル時間変化
    ut_DtDLaplaPolvel = &
         - (Ra/Ro) * ut_Temp                    &
          + ut_L2Inv_ut(ut_RadRotRot_xvh_xvh_xvh(&
                            2*xvh_KxVLon + xvh_ZxVLon, &
                            2*xvh_KxVLat + xvh_ZxVLat, &
                            2*xvh_KxVRad + xvh_ZxVRad ))
 
    ! 温度時間変化
    ut_DtDTemp = &
         - ut_Div_xvh_xvh_xvh(  &
               xvh_Vlon*xvh_Temp, xvh_Vlat*xvh_Temp, xvh_Vrad*xvh_Temp)

  end subroutine TimeDerivetives_noDiffusion

  !------------------- 拡散項(トロイダル速度場) ----------------------
  subroutine CNDiffusionMatrixTorvel( ut_Diffc, dt, DiffLUMatrix, kpivot )
    !
    ! Crank Nicolson 用拡散項陰的行列設定(トロイダル速度場)
    !
    real(8), intent(IN)          :: ut_Diffc(:,:)   ! 拡散係数
    real(8), intent(IN)          :: dt              ! 時間刻

    real(8), intent(OUT)   :: DiffLUMatrix(:,0:,0:)
                                ! Crank Nicholson 用行列(1-D dt/2▽^2, LU 分解)
    integer, intent(OUT)   :: kpivot(:,0:)        ! ピボット
                                ! Crank Nicholson 行列ピボット情報

    real(8), dimension(nmc,0:lm) :: ut_I
    ! 作業用変数
    real(8), dimension(nmc,0:lm) :: ut_DI
    ! 作業用変数

    integer :: l

    DiffLUMatrix = 0.0
    do l=0,lm
       ut_I = 0.0 ; ut_I(:,l) = 1.0             ! 各波数成分独立
       ut_DI = - ut_Diffc * dt/2.0 * ut_Lapla_ut(ut_I)
       call ut_TorBoundaries(ut_DI,cond=velBC)
       DiffLUMatrix(:,:,l) = ut_I + ut_DI
    enddo

    call ludecomp(DiffLUMatrix,kpivot)

  end subroutine CNDiffusionMatrixTorvel

  !------------------- 拡散項(ポロイダル速度場) ----------------------
  subroutine CNDiffusionMatrixPolvel( ut_Diffc, dt, DiffLUMatrix, kpivot )
    !
    ! Crank Nicolson 用拡散項陰的行列設定(ポロイダル速度場)
    !
    real(8), intent(IN)          :: ut_Diffc(:,:)   ! 拡散係数
    real(8), intent(IN)          :: dt              ! 時間刻

    real(8), intent(OUT)   :: DiffLUMatrix(:,0:,0:)
                                ! Crank Nicholson 用行列(1-D dt/2▽^2, LU 分解)
    integer, intent(OUT)   :: kpivot(:,0:)        ! ピボット
                                ! Crank Nicholson 行列ピボット情報

    real(8), dimension(nmc,0:lm) :: ut_I
    ! 作業用変数
    real(8), dimension(nmc,0:lm) :: ut_DI
    ! 作業用変数

    integer :: l

    DiffLUMatrix = 0.0
    do l=0,lm
       ut_I = 0.0 ; ut_I(:,l) = 1.0             ! 各波数成分独立
       ut_DI = - ut_Diffc * dt/2.0 * ut_Lapla_ut(ut_Lapla_ut(ut_I))
       ut_DI = ut_LaplaPol2Pol_ut(ut_DI,cond=velBC)
       DiffLUMatrix(:,:,l) = ut_I + ut_DI
    enddo

    call ludecomp(DiffLUMatrix,kpivot)

  end subroutine CNDiffusionMatrixPolvel

  !------------------- 拡散項(温度場) ----------------------
  subroutine CNDiffusionMatrixTemp( ut_Diffc, dt, DiffLUMatrix, kpivot )
    !
    ! Crank Nicolson 用拡散項陰的行列設定(温度場)
    !
    real(8), intent(IN)          :: ut_Diffc(:,:)   ! 拡散係数
    real(8), intent(IN)          :: dt              ! 時間刻

    real(8), intent(OUT)   :: DiffLUMatrix(:,0:,0:)
                                ! Crank Nicholson 用行列(1-D dt/2▽^2, LU 分解)
    integer, intent(OUT)   :: kpivot(:,0:)        ! ピボット
                                ! Crank Nicholson 行列ピボット情報

    real(8), dimension(nmc,0:lm) :: ut_I
    ! 作業用変数
    real(8), dimension(nmc,0:lm) :: ut_DI
    ! 作業用変数
    real(8), dimension(nmc,2)    :: wa_Bndry
    ! 作業用変数

    integer :: l

    DiffLUMatrix = 0.0
    wa_Bndry=0.0D0                              ! 境界の時間変化は 0
    do l=0,lm
       ut_I = 0.0 ; ut_I(:,l) = 1.0             ! 各波数成分独立
       ut_DI = - ut_Diffc * dt/2.0 * ut_Lapla_ut(ut_I)
       call ut_Boundaries(ut_DI,wa_Bndry,cond=TempBC)
       DiffLUMatrix(:,:,l) = ut_I + ut_DI
     enddo

    call ludecomp(DiffLUMatrix,kpivot)

  end subroutine CNDiffusionMatrixTemp

  !------------------- 初期値設定(一点擾乱) ----------------------
  subroutine set_initial_pointwise      ! xvh_* を設定する
    !
    ! 初期値設定ルーチン
    !
    !     初期値ファイルが指定されなかった場合に内部で初期値設定する.
    !     xvh_Torvel,xvh_Polvel,xvh_Temp を設定する
    !
    ! 擾乱温度場
    xvh_Temp =  Ro*Ri/xvh_Rad - Ri
    if ( (js .le. jm/2 .AND. jm/2 .le. je) .AND. (ks .le. km/2 .AND. km/2 .le. ke) )then
       xvh_Temp(im/2,jm/2-js+1,km/2-ks) = xvh_Temp(im/2,jm/2-js+1,km/2-ks)+Amp
    endif

    ! 速度場 ; 静止状態
    xvh_Torvel =  0.0
    xvh_Polvel =  0.0

    ut_Temp = ut_xvh(xvh_Temp)
    ut_Torvel = ut_xvh(xvh_Torvel)
    ut_Polvel = ut_xvh(xvh_Polvel)
    ut_LaplaPolvel = ut_Lapla_ut(ut_Polvel)

  end subroutine set_initial_pointwise

  !------------------- 初期値設定(ランダム) ----------------------
  subroutine set_initial_random      ! xvh_* を設定する
    !
    ! 初期値設定ルーチン
    !
    !     初期値ファイルが指定されなかった場合に内部で初期値設定する.
    !     xvh_Torvel,xvh_Polvel,xvh_Temp を設定する
    !
    !---- 初期値乱数分布
    integer              :: random_seed_size   ! 乱数の種の長さ
    integer, allocatable :: seedarray(:)       ! 乱数の種
    real, allocatable    :: harvest(:)         ! 乱数の値
    integer              :: j, k

    ! 乱数設定
    call random_seed(size=random_seed_size)
    allocate(seedarray(random_seed_size))
    call random_seed(get=seedarray)
    if ( random_seed_size > 1 ) then
       seedarray(1) = seed 
       seedarray(2) = myid
    else
       seedarray(1) = seed + myid
    endif
    call random_seed(put=seedarray)
 
    allocate(harvest(0:im-1))

    do k=1,kc
       do j=1,jc
          call random_number(harvest)
          xvh_Temp(:,j,k) = 2.0 * harvest - 1    ! [-1,1] の一様乱数
       enddo
    enddo
    deallocate(harvest,seedarray)

    ! 振幅設定
    xvh_Temp = Amp * xvh_Temp

    ! 基本場設定
    xvh_Temp =  xvh_Temp + Ro*Ri/xvh_Rad - Ri
 
    ! 速度場 ; 静止状態
    xvh_Torvel =  0.0
    xvh_Polvel =  0.0

    ut_Temp = ut_xvh(xvh_Temp)
    ut_Torvel = ut_xvh(xvh_Torvel)
    ut_Polvel = ut_xvh(xvh_Polvel)
    ut_LaplaPolvel = ut_Lapla_ut(ut_Polvel)

  end subroutine set_initial_random

  !------------------- 初期値設定(ベンチマーク) ----------------------
  subroutine set_initial_bench         ! xyz_* を設定する
    !
    ! 初期値設定ルーチン
    !
    !     初期値ファイルが指定されなかった場合に内部で初期値設定する.
    !     xyz_Torvel,xyz_Polvel,xyz_Tormag,xyz_Polmag,xyz_Temp を設定する
    !
    real(8), dimension(size(xvh_Vlon,1),size(xvh_Vlon,2),size(xvh_Vlon,3)) &
         :: xvh_X         ! x=2r-r_i-r_o
    real(8), dimension(size(xvh_Vlon,1),size(xvh_Vlon,2),size(xvh_Vlon,3)) &
         :: xvh_theta     ! 余緯度

    real(8), parameter  :: A=0.1D0                     ! 初期温度擾乱振幅
    real(8), parameter  :: pi=3.1415926535897932385D0

    xvh_X = 2*xvh_Rad - Ri - Ro

    ! 温度場
    xvh_Temp =  Ro*Ri/xvh_Rad - Ri &
       + 210*A/sqrt(17920*pi) * (1 - 3*xvh_X**2 + 3*xvh_X**4 -xvh_X**6) &
           * sin(pi/2-xvh_Lat)**4 * cos(4*xvh_Lon)

    ! 速度場 ; 静止状態
    xvh_Torvel =  0 ; xvh_Polvel =  0

    ut_Temp = ut_xvh(xvh_Temp)
    ut_Torvel = ut_xvh(xvh_Torvel)
    ut_Polvel = ut_xvh(xvh_Polvel)
    ut_LaplaPolvel = ut_Lapla_ut(ut_Polvel)

  end subroutine set_initial_bench

  !------------------- 初期値設定(ファイルから) ----------------------
  subroutine initial_read_restart
    !
    ! 初期値設定ルーチン
    !
    !     初期値ファイルから変数
    !     xyz_Torvel,xyz_Polvel,xyz_Tormag,xyz_Polmag,xyz_Temp,
    !     ut_DtDTorvel1, ut_DtDLaplaPolvel1,
    !     ut_DtDTormag1, ut_DtDPolmag1, ut_DtDTemp1 を読み込む

    call HistoryGet( trim(initial_file), 'torvel', ut_Torvel, initial_time )
    call HistoryGet( trim(initial_file), 'polvel', ut_Polvel, initial_time )
    call HistoryGet( trim(initial_file), 'laplapolvel', ut_LaplaPolvel, initial_time )
    call HistoryGet( trim(initial_file), 'temp',   ut_Temp,   initial_time )

    call HistoryGet( trim(initial_file), 'DtDTorvel', ut_DtDTorvel1, initial_time )
    call HistoryGet( trim(initial_file), 'DtDLaplaPolvel', ut_DtDLaplaPolvel1, initial_time )
    call HistoryGet( trim(initial_file), 'DtDTemp', ut_DtDTemp1, initial_time )

  end subroutine initial_read_restart
  !    
  !------------------- リスタート出力ルーチン ----------------------
  !
  subroutine output_restart_init
    !
    ! リスタート出力初期化
    !
    integer :: n, l

    call HistoryCreate( &
           file=trim(rst_file), title=trim(title), source=trim(source), &
           institution='GFD_Dennou Club SPMODEL project',     &
           dims=(/'nm','l ','t '/), &
           dimsizes=(/nmc,lm+1,0/),&
           longnames=(/'Hor.wave number index ',  &
                       'Radial chebyshev order',  &
                       'time                  '/),&
           units=(/'1','1','1'/),   &
           origin=real(time), interval=real(rst_intstep*delta_t), &
           xtypes=(/'real'/), history=hst_rst)

    !---- 座標変数定義, 出力 ----
    call HistoryPut('nm',(/(dble(n),n=1,nmc)/), history=hst_rst) ! 変数出力
    call HistoryPut('l',(/(dble(l),l=0,lm)/), history=hst_rst) ! 変数出力

    !---- 物理変数定義 ----
    ! 速度ポテンシャル 
    call HistoryAddVariable( &                                ! 変数定義
           varname='torvel', dims=(/'nm','l ','t '/), &
           longname='toroidal velocity field', &
           units='1', xtype='double', history=hst_rst)
    call HistoryAddVariable( &                                ! 変数定義
           varname='polvel', dims=(/'nm','l ','t '/), &
           longname='poloidal velocity field', &
           units='1', xtype='double', history=hst_rst)
    call HistoryAddVariable( &                                ! 変数定義
           varname='laplapolvel', dims=(/'nm','l ','t '/), &
           longname='Laplacian poloidal velocity field', &
           units='1', xtype='double', history=hst_rst)
    ! エントロピー
    call HistoryAddVariable( &                                ! 変数定義
           varname='temp', dims=(/'nm','l ','t '/), &
           longname='temperature', &
           units='1', xtype='double', history=hst_rst)

    !---- 時間変化項定義 ----
    call HistoryAddVariable( &                                ! 変数定義
           varname='DtDTorvel', dims=(/'nm','l ','t '/), &
           longname='time variation of velocity toroidal potential', &
           units='1', xtype='double', history=hst_rst)
    call HistoryAddVariable( &                                ! 変数定義
           varname='DtDLaplaPolvel', dims=(/'nm','l ','t '/), &
           longname='time variation of laplacian of velocity poloidal potential',&
           units='1', xtype='double', history=hst_rst)
    call HistoryAddVariable( &                                ! 変数定義
           varname='DtDTemp', dims=(/'nm','l ','t '/), &
           longname='time variation of temperature', &
           units='1', xtype='double', history=hst_rst)

    !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ----
    call HistoryAddAttr('nm','+Radius', (/ri,ro,eta/), history=hst_rst)
    call HistoryAddAttr('nm','+delta_t', delta_t, history=hst_rst)

    call HistoryAddAttr('nm', '+Rayleigh',   Ra, history=hst_rst)
    call HistoryAddAttr('nm', '+Prandtl',    Pr, history=hst_rst)
    call HistoryAddAttr('nm', '+Ekman',      Ekman, history=hst_rst)

    call HistoryAddAttr('nm', '+VelBoundary',  VelBC,  history=hst_rst)
    call HistoryAddAttr('nm', '+TempBoundary', TempBC, history=hst_rst)

    call HistoryAddAttr('nm', '+epsilon', epsilon, history=hst_rst)
    call HistoryAddAttr('nm', '+n0', n0, history=hst_rst)

  end subroutine output_restart_init

  subroutine output_restart
    !
    ! リスタート出力
    !
    if ( myid==0 ) &
         write(6,*) 'Output restart, it = ',it, ' time = ', time

    call HistoryPut('t',real(time), history=hst_rst)
    call HistoryPut('torvel',ut_Torvel, history=hst_rst)
    call HistoryPut('polvel',ut_Polvel, history=hst_rst)
    call HistoryPut('laplapolvel',ut_LaplaPolvel, history=hst_rst)
    call HistoryPut('temp',ut_Temp, history=hst_rst)

    call HistoryPut('DtDTorvel',ut_DtDTorvel1, history=hst_rst)
    call HistoryPut('DtDLaplaPolvel',ut_DtDLaplaPolvel1, history=hst_rst)
    call HistoryPut('DtDTemp',ut_DtDTemp1, history=hst_rst)

  end subroutine output_restart

  subroutine output_restart_close
    !
    ! リスタート出力終了処理
    !
    call HistoryClose(hst_rst)
  end subroutine output_restart_close
  !    
  !------------------- ヒストリー出力ルーチン ----------------------
  !
  subroutine output_history_init
    !
    ! ヒストリー出力初期化ルーチン
    !
    !   ファイル作成
    !   変数定義
    !
    !---- ヒストリーファイル作成 ----
    call HistoryCreate( &
           file=trim(hst_file), &
           title=trim(title), source=trim(source), &
           institution='GFD_Dennou Club SPMODEL project',     &
           dims=(/'lon','lat','rad','n  ','m  ','l  ','t  '/), &
           dimsizes=(/im,jc,kc,nm+1,2*nm+1,lm+1,0/),&
           longnames=(/'Longitude             ','Latitude              ',&
                       'Radius                ',&
                       'Hor.total wave number ',&
                       'zonal wave number     ','Radial chebyshev order',&
                       'time                  '/),&
           units=(/'1','1','1','1','1','1','1'/),   &
           origin=real(time), interval=real(hst_intstep*delta_t))

   !---- 座標変数定義, 出力 ----
    call HistoryPut('lon',x_Lon/pi*180)                       ! 変数出力
    call HistoryAddattr('lon','topology','circular')          ! 周期属性
    call HistoryAddattr('lon','modulo',360.0/mint)            ! 周期属性
    call HistoryPut('lat',v_Lat/pi*180)                       ! 変数出力
    call HistoryPut('rad',h_Rad)                              ! 変数出力
    call HistoryPut('n',(/(dble(n),n=0,nm)/))                 ! 変数出力
    call HistoryPut('m',(/(dble(m),m=-nm,nm)/))               ! 変数出力
    call HistoryPut('l',(/(dble(l),l=0,lm)/))                 ! 変数出力

    call HistoryAddVariable( &                                ! 変数定義
           varname='lon_weight', dims=(/'lon'/), &
           longname='weight function in longitude', &
           units='1', xtype='double')
    call HistoryPut('lon_weight',x_Lon_weight)                ! 変数出力
    call HistoryAddVariable( &                                ! 変数定義
           varname='lat_weight', dims=(/'lat'/), &
           longname='weight function in latitude', &
           units='1', xtype='double')
    call HistoryPut('lat_weight',v_Lat_weight)                ! 変数出力
    call HistoryAddVariable( &                                ! 変数定義
           varname='rad_weight', dims=(/'rad'/), &
           longname='weight function in radial', &
           units='1', xtype='double')
    call HistoryPut('rad_weight',h_Rad_weight)                ! 変数出力

   !---- 物理変数定義 ----
    !-- 速度ポテンシャル --
    call HistoryAddVariable( &                                ! 変数定義
           varname='torvel', dims=(/'lon','lat','rad','t  '/), &
           longname='toroidal velocity potential', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='polvel', dims=(/'lon','lat','rad','t  '/), &
           longname='poloidal velocity potential', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='laplapolvel', dims=(/'lon','lat','rad','t  '/), &
           longname='laplacian poloidal velocity potential', &
           units='1', xtype='double')

    !-- 温度 --
    call HistoryAddVariable( &                                ! 変数定義
           varname='temp', dims=(/'lon','lat','rad','t  '/), &
           longname='temperature', units='1', xtype='double')
    !-- 速度 --
    call HistoryAddVariable( &                                ! 変数定義
           varname='vlon', dims=(/'lon','lat','rad','t  '/), &
           longname='lon-velocity', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='vlat', dims=(/'lon','lat','rad','t  '/), &
           longname='lat-velocity', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='vrad', dims=(/'lon','lat','rad','t  '/), &
           longname='rad-velocity', units='1', xtype='double')
    !-- 渦度 --
    call HistoryAddVariable( &                                ! 変数定義
           varname='zlon', dims=(/'lon','lat','rad','t  '/), &
           longname='lon-vorticity', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='zlat', dims=(/'lon','lat','rad','t  '/), &
           longname='lat-vorticity', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='zrad', dims=(/'lon','lat','rad','t  '/), &
           longname='rad-voticity', units='1', xtype='double')

   !---- その他診断量定義 ----
    !-- 全エネルギー --
    call HistoryAddVariable( &                                ! 変数定義
           varname='ek', dims=(/'t  '/), &
           longname='mean kinetic energy', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ektot', dims=(/'t  '/), &
           longname='total kinetic energy', units='1', xtype='double')

    !-- 平均帯状流 --
    call HistoryAddVariable( &                                ! 変数定義
           varname='am', dims=(/'lat','rad','t  '/), &
           longname='Angular momentum', units='1', xtype='double')
    !-- ヘリシティー --
    call HistoryAddVariable( &                                ! 変数定義
           varname='hel', dims=(/'lon','lat','rad','t  '/), &
           longname='helicity', units='1', xtype='double')

    !-- 運動エネルギースペクトル診断量定義 --
    call HistoryAddVariable( &                                ! 変数定義
           varname='ektor_nmz', dims=(/'n  ','m  ','rad','t  '/), &
           longname='Toroidal kinetic energy spectrum', units='1', xtype='double')
    call HistoryAddAttr('ektor_nmz','missing_value', vmiss  )

    call HistoryAddVariable( &                                ! 変数定義
           varname='ektoreqs_nmz', dims=(/'n  ','m  ','rad','t  '/), &
           longname='Toroidal kinetic energy spectrum (eq-sym)', units='1', xtype='double')
    call HistoryAddAttr('ektoreqs_nmz','missing_value', vmiss  )

    call HistoryAddVariable( &                                ! 変数定義
           varname='ektoreqa_nmz', dims=(/'n  ','m  ','rad','t  '/), &
           longname='Toroidal kinetic energy spectrum (eq-antisym)', units='1', xtype='double')
    call HistoryAddAttr('ektoreqa_nmz','missing_value', vmiss  )

    call HistoryAddVariable( &                                ! 変数定義
           varname='ekpol_nmz', dims=(/'n  ','m  ','rad','t  '/), &
           longname='Poloidal kinetic energy spectrum', units='1', xtype='double')
    call HistoryAddAttr('ekpol_nmz','missing_value', vmiss  )

    call HistoryAddVariable( &                                ! 変数定義
           varname='ektor_nz', dims=(/'n  ','rad','t  '/), &
           longname='Toroidal kinetic energy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ekpol_nz', dims=(/'n  ','rad','t  '/), &
           longname='Poloidal kinetic energy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ektor_n', dims=(/'n  ','t  '/), &
           longname='Toroidal kinetic energy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ekpol_n', dims=(/'n  ','t  '/), &
           longname='Poloidal kinetic energy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ektor', dims=(/'t  '/), &
           longname='Toroidal kinetic energy', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ekpol', dims=(/'t  '/), &
           longname='Poloidal kinetic energy', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ekpot', dims=(/'t'/), &
           longname='Kinetic energy (by potential)', units='1', xtype='double')

   !---- エンストロフィースペクトル診断量定義 ----
    call HistoryAddVariable( &                                ! 変数定義
           varname='enstor_nmz', dims=(/'n  ','m  ','rad','t  '/), &
           longname='Toroidal enstrophy spectrum', units='1', xtype='double')
    call HistoryAddAttr('enstor_nmz','missing_value', vmiss  )

    call HistoryAddVariable( &                                ! 変数定義
           varname='enspol_nmz', dims=(/'n  ','m  ','rad','t  '/), &
           longname='Poloidal enstrophy spectrum', units='1', xtype='double')
    call HistoryAddAttr('enspol_nmz','missing_value', vmiss  )

    call HistoryAddVariable( &                                ! 変数定義
           varname='enstor_nz', dims=(/'n  ','rad','t  '/), &
           longname='Toroidal enstrophy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='enspol_nz', dims=(/'n  ','rad','t  '/), &
           longname='Poloidal enstrophy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='enstor_n', dims=(/'n  ','t  '/), &
           longname='Toroidal enstrophy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='enspol_n', dims=(/'n  ','t  '/), &
           longname='Poloidal enstrophy spectrum', units='1', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='enstor', dims=(/'t  '/), &
           longname='Toroidal enstrophy', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='enspol', dims=(/'t  '/), &
           longname='Poloidal enstrophy', units='1', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='enspot', dims=(/'t'/), &
           longname='enstrophy (by potential)', units='1', xtype='double')

    !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ----
    call HistoryAddAttr('lon','+Radius', (/ri,ro,eta/) )
    call HistoryAddAttr('lon','+delta_t', delta_t  )

    call HistoryAddAttr('lon', '+Rayleigh',   Ra    )
    call HistoryAddAttr('lon', '+Prandtl',    Pr    )
    call HistoryAddAttr('lon', '+Ekman',      Ekman )

    call HistoryAddAttr('lon', '+VelBoundary',  VelBC  )
    call HistoryAddAttr('lon', '+TempBoundary', TempBC )

    call HistoryAddAttr('lon', '+epsilon', epsilon )
    call HistoryAddAttr('lon', '+n0', n0  )

  end subroutine output_history_init

  subroutine output_history
    !
    ! ヒストリー出力
    !
    !   各物理量の gtool4 ファイルへの出力を行う
    !
    real(8) :: ek
    real(8) :: ut_TorvelEqS(nmc,0:lm)
    real(8) :: ut_TorvelEqA(nmc,0:lm)    
    real(8) :: nmz_Work(-nm:nm,0:nm,0:km)
    
    ek = AvrLonLatRad_xvh((xvh_Vlon**2+xvh_Vlat**2+xvh_Vrad**2)/2)
    
    if ( myid==0 ) &
         write(6,*) 'Output history, it = ',it, ' time = ', time, ' ek = ', ek

    call HistoryPut('t',real(time))

    !---- 物理変数出力 ----
    call HistoryPut('torvel',xvh_ut(ut_Torvel))
    call HistoryPut('polvel',xvh_ut(ut_Polvel))
    call HistoryPut('laplapolvel',xvh_ut(ut_LaplaPolvel))
    call HistoryPut('temp',xvh_ut(ut_Temp))

    call HistoryPut('vlon',xvh_Vlon)
    call HistoryPut('vlat',xvh_Vlat)
    call HistoryPut('vrad',xvh_Vrad)

    call HistoryPut('zlon',xvh_Zlon)
    call HistoryPut('zlat',xvh_Zlat)
    call HistoryPut('zrad',xvh_Zrad)

    !---- 診断量出力 ----
    !-- 平均エネルギー --
    call HistoryPut('ek',ek)

    !-- 全エネルギー --
    call HistoryPut('ektot',IntLonLatRad_xvh((xvh_Vlon**2+xvh_Vlat**2+xvh_Vrad**2)/2))

    !-- 帯状流 --
    call HistoryPut('am',vh_avrLon_xvh(xvh_Vlon*xvh_Rad*cos(xvh_Lat)))
    call HistoryPut('hel',xvh_Vlon*xvh_Zlon + xvh_Vlat*xvh_Zlat + xvh_Vrad*xvh_Zrad)

    !---- 運動エネルギースペクトル診断量出力 ----
    ut_TorvelEqS = 0.0D0 ; ut_TorvelEqA = 0.0D0

    do l=1,nmc
       nmary = nm_l(l+nms-1)
       if ( mod(nmary(1)+nmary(2),2 ) .eq. 0 ) then
          ut_TorvelEqS(l,:) = ut_Torvel(l,:)
       else
          ut_TorvelEqA(l,:) = ut_Torvel(l,:)
       end if
    end do

    nz_EkTor = nz_ToroidalEnergySpectrum_ut(ut_Torvel)
    nz_EkPol = nz_PoloidalEnergySpectrum_ut(ut_Polvel)

    nmz_Work = nmz_ToroidalEnergySpectrum_ut(ut_Torvel)
    call HistoryPut('ektor_nmz',nmz_work(:,:,ks-1:ke-1))
    nmz_Work = nmz_ToroidalEnergySpectrum_ut(ut_TorvelEqS)
    call HistoryPut('ektoreqs_nmz',nmz_work(:,:,ks-1:ke-1))
    nmz_Work = nmz_ToroidalEnergySpectrum_ut(ut_TorvelEqA)
    call HistoryPut('ektoreqa_nmz',nmz_work(:,:,ks-1:ke-1))
    nmz_Work = nmz_ToroidalEnergySpectrum_ut(ut_Polvel)
    call HistoryPut('ekpol_nmz',nmz_work(:,:,ks-1:ke-1))
    call HistoryPut('ektor_nz',nz_EkTor(:,ks-1:ke-1))
    call HistoryPut('ekpol_nz',nz_EkPol(:,ks-1:ke-1))
    call HistoryPut('ektor_n',a_Int_ag(nz_EkTor/(spread(z_Rad,1,nm+1)**2)))
    call HistoryPut('ekpol_n',a_Int_ag(nz_EkPol/(spread(z_Rad,1,nm+1)**2)))
    call HistoryPut('ektor',IntRad_z(sum(nz_EkTor,1)/z_Rad**2))
    call HistoryPut('ekpol',IntRad_z(sum(nz_EkPol,1)/z_Rad**2))
    call HistoryPut('ekpot',IntRad_z(sum(nz_EkTor,1)/z_Rad**2)+IntRad_z(sum(nz_EkPol,1)/z_Rad**2))

   !---- エンストロフィースペクトル診断量出力 ----
    nz_EnsTor = nz_ToroidalEnergySpectrum_ut(-ut_LaplaPolvel)
    nz_EnsPol = nz_PoloidalEnergySpectrum_ut(ut_Torvel)

    nmz_Work = nmz_ToroidalEnergySpectrum_ut(-ut_LaplaPolvel)
    call HistoryPut('enstor_nmz',nmz_work(:,:,ks-1:ke-1))
    nmz_Work = nmz_ToroidalEnergySpectrum_ut(ut_Torvel)
    call HistoryPut('enspol_nmz',nmz_work(:,:,ks-1:ke-1))
    call HistoryPut('enstor_nz',nz_EnsTor(:,ks-1:ke-1))
    call HistoryPut('enspol_nz',nz_EnsPol(:,ks-1:ke-1))
    call HistoryPut('enstor_n',a_Int_ag(nz_EnsTor/(spread(z_Rad,1,nm+1)**2)))
    call HistoryPut('enspol_n',a_Int_ag(nz_EnsPol/(spread(z_Rad,1,nm+1)**2)))
    call HistoryPut('enstor',IntRad_z(sum(nz_EnsTor,1)/z_Rad**2))
    call HistoryPut('enspol',IntRad_z(sum(nz_EnsPol,1)/z_Rad**2))
    call HistoryPut('enspot',IntRad_z(sum(nz_EnsTor,1)/z_Rad**2)+IntRad_z(sum(nz_EnsPol,1)/z_Rad**2))

  end subroutine output_history

  subroutine output_history_close
    !
    ! ヒストリー出力終了処理
    !
    call HistoryClose

  end subroutine output_history_close
  !    
  !------------------- 診断出力ルーチン ----------------------
  !
  subroutine output_diag

    real(8) :: ut_Null(nmc,0:lm)
    real(8) :: xvh_VLon(0:im-1,jc,kc)
    real(8) :: xvh_VLat(0:im-1,jc,kc)
    real(8) :: xvh_VRad(0:im-1,jc,kc)
    real(8) :: EkTor, EkPol

!!$       nz_EkTor = nz_ToroidalEnergySpectrum_ut(ut_Torvel)
!!$       nz_EkPol = nz_PoloidalEnergySpectrum_ut(ut_Polvel)
    ut_Null = 0.0D0

    call ut_Potential2VectorMPI &
         ( xvh_VLon, xvh_VLat, xvh_VRad, ut_Torvel, ut_Null)

    EkTor = IntLonLatRad_xvh(0.5D0*(xvh_VLon**2 + xvh_VLat**2 + xvh_VRad**2))

    call ut_Potential2VectorMPI &
         ( xvh_VLon, xvh_VLat, xvh_VRad, ut_Null, ut_Polvel)

    EkPol = IntLonLatRad_xvh(0.5D0*(xvh_VLon**2 + xvh_VLat**2 + xvh_VRad**2))

    if ( myid==0 ) then
       write(6,*) 'it = ',it, ' time = ', time,  &
                 ' EkTor = ', EkTor, ' EkPol = ', EkPol
    endif

  end subroutine output_diag

end program bsnsq_rot_3d_shell

