!--------------------------------------------------------------------------
!  Copyright (C) 2002--2005 SPMODEL Development Group. All rights reserved.
!--------------------------------------------------------------------------
!
!表題  回転球殻ナビエーストークス流体モデル : 
!      2nd order Adams-Bashforth + Crank Nicholson scheme
!      球殻回転粘着境界問題
!
!履歴  2005/02/11  竹広 真一   mhd-boussinesq モデル元に簡略化
!      2005/02/19  竹広 真一   球殻回転粘着境界問題用に改造
!      2008/08/10  竹広 真一   配列添え字変更 im -> 0:im-1
!
program NS_rot_3d_shell_abcn_restart_namelist_rotbndry

  use lumatrix
  use wt_module
  use dc_message, only : MessageNotify
  use gt4_history
  implicit none

 !---- 空間解像度設定 ----
  ! T10 N8
  integer,parameter  :: im=32, jm=16, km=8  ! 格子点の設定(経度, 緯度, 動径)
  integer,parameter  :: nm=10, lm=8         ! 切断波数の設定(水平, 動径)
  ! T10 N16
  !integer,parameter  :: im=32, jm=16, km=16  ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter  :: nm=10, lm=16         ! 切断波数の設定(水平, 動径)
  !! T21 N16
  !integer,parameter  :: im=64, jm=32, km=16  ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter  :: nm=21, lm=16         ! 切断波数の設定(水平, 動径)
  !! T21 N24
  !integer,parameter  :: im=64, jm=32, km=24  ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter  :: nm=21, lm=24         ! 切断波数の設定(水平, 動径)
  !! T21 N32
  !integer,parameter  :: im=64, jm=32, km=32  ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter  :: nm=21, lm=32         ! 切断波数の設定(水平, 動径)
  !! T42 N32
  !integer,parameter  :: im=128, jm=64, km=32 ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter  :: nm=42, lm=32         ! 切断波数の設定(水平, 動径)
  !! T42 N48
  !integer,parameter  :: im=128, jm=64, km=48 ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter  :: nm=42, lm=48         ! 切断波数の設定(水平, 動径)
  !  T85 N32
  !integer,parameter :: im=256, jm=128, km=32 ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter :: nm=85, lm=32          ! 切断波数の設定(水平, 動径)
  !  T85 N48
  !integer,parameter :: im=256, jm=128, km=48 ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter :: nm=85, lm=48          ! 切断波数の設定(水平, 動径)
  !  T85 N64
  !integer,parameter :: im=256, jm=128, km=64 ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter :: nm=85, lm=64          ! 切断波数の設定(水平, 動径)
  !  T170 N48
  !integer,parameter :: im=512, jm=256, km=48 ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter :: nm=170, lm=48         ! 切断波数の設定(水平, 動径)
  !  T170 N64
  !integer,parameter :: im=512, jm=256, km=64 ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter :: nm=170, lm=64         ! 切断波数の設定(水平, 動径)
  !  T341 N80
  !integer,parameter :: im=1024, jm=512, km=80 ! 格子点の設定(経度, 緯度, 動径)
  !integer,parameter :: nm=341, lm=80          ! 切断波数の設定(水平, 動径)

 !---- 変数 ----
  real(8), dimension(0:im-1,jm,0:km) :: xyz_Torvel      ! トロイダル速度
  real(8), dimension((nm+1)**2,0:lm) :: wt_Torvel       ! トロイダル速度
  real(8), dimension(0:im-1,jm,0:km) :: xyz_Polvel      ! ポロイダル速度
  real(8), dimension((nm+1)**2,0:lm) :: wt_Polvel       ! ポロイダル速度
  real(8), dimension(0:im-1,jm,0:km) :: xyz_LaplaPolvel ! ポロイダル速度
  real(8), dimension((nm+1)**2,0:lm) :: wt_LaplaPolvel  ! ▽^2ポロイダル速度

  real(8), dimension(0:im-1,jm,0:km) :: xyz_Vlon       ! 速度(経度)
  real(8), dimension(0:im-1,jm,0:km) :: xyz_Vlat       ! 速度(緯度)
  real(8), dimension(0:im-1,jm,0:km) :: xyz_Vrad       ! 速度(動径)

  real(8), dimension(0:im-1,jm,0:km) :: xyz_Zlon       ! 渦度(経度)
  real(8), dimension(0:im-1,jm,0:km) :: xyz_Zlat       ! 渦度(緯度)
  real(8), dimension(0:im-1,jm,0:km) :: xyz_Zrad       ! 渦度(動径)

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

 ! トロイダル速度時間変化
  real(8), dimension((nm+1)**2,0:lm) :: wt_DtDTorvel1, wt_DtDTorvel0
 ! ▽^2ポロイダル速度時間変化
  real(8), dimension((nm+1)**2,0:lm) :: wt_DtDLaplaPolvel1, wt_DtDLaplaPolvel0

 ! Crank Nicholson 陰的計算用拡散行列(速度場)
  real(8), dimension((nm+1)*(nm+1),0:lm,0:lm)   :: DifLUMT_Vel
  integer, dimension((nm+1)*(nm+1),0:lm)        :: kpivot_Vel     ! ピボット

 ! 上下トロイダルポテンシャル境界値
  real(8), dimension(0:im-1,jm)        :: xy_TorvelBoundaryTop
  real(8), dimension(0:im-1,jm)        :: xy_TorvelBoundaryBottom
  real(8), dimension((nm+1)*(nm+1))    :: w_TorvelBoundaryTop
  real(8), dimension((nm+1)*(nm+1))    :: w_TorvelBoundaryBottom

 !---- その他 ----
  real(8), parameter :: pi=3.1415926535897932385D0
  integer :: it=0
  real(8) :: time

 !---- NAMELIST 変数 ----
  logical            :: verbose=.false.          ! 出力メッセージレベル
  namelist /message/  verbose                    !
                                                 !-- 座標変数など --
  real(8)            :: ri=0.0D0                 ! 内半径
  real(8)            :: ro=0.0D00                ! 外半径
  real(8)            :: eta=0.0D0                ! 内外半径比
  namelist /radius/  ri, ro, eta
                                                 ! -- 初期値 --
  character(len=100) :: 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=100) :: output_file=''           ! 出力ファイル名
  character(len=100) :: title = &                ! タイトル
            'Navier-Stokes equation in a 3-dim spherical shell'!
  integer :: ndisp=200                           ! 出力間隔ステップ数
  namelist /output/    output_file, title, ndisp ! 
                                                 ! -- 物理パラメター --
  real(8)            :: Rs=0.0D0                 ! ロスビー数
  real(8)            :: Re=1.0D0                 ! レイノルズ数
  real(8)            :: RsInv=0.0D0              ! 1/ロスビー数
  real(8)            :: ReInv=0.0D0              ! 1/レイノルズ数
  namelist /physics/  Rs, Re                     ! 
                                                 ! -- 境界条件 --
  character(len=2)   :: VelBC='RR'               ! 速度境界条件(RR/RF/FR/FF)
  real(8)            :: OmegaTop=1.0D0           ! 上境界回転角速度
  real(8)            :: LonOmegaTop=0.0D0        ! 上境界回転軸経度
  real(8)            :: LatOmegaTop=90.0D0       ! 上境界回転軸緯度
  real(8)            :: OmegaBottom=1.0D0        ! 下境界回転角速度
  real(8)            :: LonOmegaBottom=0.0D0     ! 下境界回転軸経度
  real(8)            :: LatOmegaBottom=90.0D0    ! 下境界回転軸緯度
  namelist /boundary/ VelBC, &
                      OmegaTop,    LonOmegaTop,    LatOmegaTop,   &
                      OmegaBottom, LonOmegaBottom, LatOmegaBottom

 !------------------ NAMELIST 読み込み ---------------------
  write(6,nml=message)  ; read(5,nml=message) ; write(6,nml=message)

  if (verbose) write(6,nml=radius) ; read(5,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
  write(6,nml=radius)

  if(verbose)write(6,nml=initial) ; read(5,nml=initial) ;write(6,nml=initial)
  if(verbose)write(6,nml=tint)    ; read(5,nml=tint)    ;write(6,nml=tint)
  if(verbose)write(6,nml=output)  ; read(5,nml=output)  ;write(6,nml=output)

  if(verbose)write(6,nml=physics) ; read(5,nml=physics)
  if ( Rs == 0.0D0 ) then
     RsInv = 0.0D0
     call MessageNotify('M','main', &
             'Inverse of Rossby Number is set to zero.')
  else
     RsInv = 1.0D0/Rs
  endif
  if ( Re == 0.0D0 ) then
     ReInv = 0.0D0
     call MessageNotify('M','main', &
             'Inverse of Reynolds Number is set to zero.')
  else
     ReInv = 1.0D0/Re
  endif
  write(6,nml=physics)

  if(verbose)write(6,nml=boundary); read(5,nml=boundary);write(6,nml=boundary)

 !---------------- 座標値の設定 ---------------------
  call wt_Initial(im,jm,km,nm,lm,ri,ro)

 !------------------- 初期値設定 ----------------------
  if ( initial_file == "") then
     call set_initial_values( xyz_Torvel, xyz_Polvel )

  else
     call initial_read_gtool4(initial_file, initial_time, &
          xyz_Torvel, xyz_Polvel, wt_DtDTorvel1, wt_DtDLaplaPolvel1 )
  endif

  time = initial_time
  
  wt_Torvel = wt_xyz(xyz_Torvel) 
  wt_Polvel = wt_xyz(xyz_Polvel)
  wt_LaplaPolvel = wt_Lapla_wt(wt_Polvel)
  xyz_LaplaPolvel = xyz_wt(wt_LaplaPolvel)

  call wt_Potential2Vector(xyz_Vlon, xyz_Vlat, xyz_Vrad, wt_Torvel, wt_Polvel)

  if ( initial_file == "") then
     call TimeDerivetives_noDiffusion( &
          wt_Torvel, wt_Polvel, wt_LaplaPolvel, &
          xyz_Vlon, xyz_Vlat, xyz_Vrad, &
          wt_DtDTorvel1, wt_DtDLaplaPolvel1 )
  endif

  call CNDiffusionMatrix( ReInv, delta_t, DifLUMT_Vel, kpivot_Vel )

 !------------------ 境界条件設定 -------------------
  xy_TorvelBoundaryTop = OmegaTop * Ro &
                          * (  cos(LatOmegaTop/180.0*pi)*cos(xy_Lat)     &
                                 * cos(LonOmegaTop/180.0D0*pi-xy_Lon)    &
                              + sin(LatOmegaTop/180.0D0*pi)*sin(xy_Lat) )
  w_TorvelBoundaryTop=w_xy(xy_TorvelBoundaryTop)

  xy_TorvelBoundaryBottom = OmegaBottom * Ri &
                          * (  cos(LatOmegaBottom/180.0*pi)*cos(xy_Lat)   &
                                 * cos(LonOmegaBottom/180.0D0*pi-xy_Lon)  &
                              + sin(LatOmegaBottom/180.0D0*pi)*sin(xy_Lat) )
  w_TorvelBoundaryBottom=w_xy(xy_TorvelBoundaryBottom)

 !------------------- 時間積分(Adams-Bashforth 法) --------------------
  call output_gtool4_init
  if ( initial_file == '' ) call output_gtool4    ! 内部で与えた初期値は出力

  do it=1,nstep
     time = initial_time + it * delta_t

     ! 昔の時間変化項を保存
     wt_DtDTorvel0 = wt_DtDTorvel1
     wt_DtDLaplaPolvel0 = wt_DtDLaplaPolvel1

     call TimeDerivetives_noDiffusion( &
       wt_Torvel, wt_Polvel, wt_LaplaPolvel, &
       xyz_Vlon, xyz_Vlat, xyz_Vrad, &
       wt_DtDTorvel1, wt_DtDLaplaPolvel1 )

     ! トロイダル速度場
     wt_Torvel = wt_Torvel + delta_t *( 3/2.0*wt_DtDTorvel1 - 1/2.0*wt_DtDTorvel0 ) &
                           + delta_t/2.0 * ReInv * wt_Lapla_wt(wt_Torvel) 
     wt_Torvel = LUSolve(DifLUMT_Vel,kpivot_Vel,wt_Torvel)

     ! ポロイダル速度場
     wt_LaplaPolvel = wt_LaplaPolvel &
           + delta_t *(3/2.0*wt_DtDLaplaPolvel1 - 1/2.0*wt_DtDLaplaPolvel0 ) &
           + delta_t/2.0 * ReInv * wt_Lapla_wt(wt_LaplaPolvel) 
     wt_LaplaPolvel = LUSolve(DifLUMT_Vel,kpivot_Vel,wt_LaplaPolvel)

     ! 境界条件
     call wt_TorBoundariesGrid(wt_Torvel,cond=velBC, &
          values=(/w_TorvelBoundaryTop,w_TorvelBoundaryBottom/))
     wt_Polvel = wt_wz(wz_LaplaPol2Pol_wz(wz_wt(wt_LaplaPolvel),cond=velBC))
     wt_LaplaPolvel = wt_Lapla_wt(wt_Polvel)

     ! グリッド値計算
     xyz_Torvel = xyz_wt(wt_Torvel)
     xyz_Polvel = xyz_wt(wt_Polvel)
     xyz_LaplaPolvel = xyz_wt(wt_LaplaPolvel)

     call wt_Potential2Vector(xyz_Vlon,xyz_Vlat,xyz_Vrad,wt_Torvel,wt_Polvel)

     if(mod(it,ndisp) .eq. 0)then                    ! 出力
        call output_gtool4
     endif
  enddo
  call output_gtool4_close

contains

! 以上 メインプログラム 
!--------------------------------------------------------------------------------
! 以下 サブルーチン
  !------------------- 時間変化項(拡散項以外) ----------------------
  subroutine TimeDerivetives_noDiffusion( &
       wt_Torvel, wt_Polvel, wt_LaplaPolvel, &
       xyz_Vlon, xyz_Vlat, xyz_Vrad, &
       wt_DtDTorvel, wt_DtDLaplaPolvel )

   ! トロイダル速度
    real(8), dimension((nm+1)**2,0:lm), intent(IN) :: wt_Torvel
   ! ポロイダル速度
    real(8), dimension((nm+1)**2,0:lm), intent(IN) :: wt_Polvel
   ! ▽^2 ポロイダル速度
    real(8), dimension((nm+1)**2,0:lm), intent(IN) :: wt_LaplaPolvel
   ! 温度

    real(8), dimension(0:im-1,jm,0:km), intent(IN) :: xyz_Vlon    ! 速度(経度)
    real(8), dimension(0:im-1,jm,0:km), intent(IN) :: xyz_Vlat    ! 速度(緯度)
    real(8), dimension(0:im-1,jm,0:km), intent(IN) :: xyz_Vrad    ! 速度(動径)

   ! トロイダル速度時間変化項
    real(8), dimension((nm+1)**2,0:lm), intent(OUT) :: wt_DtDTorvel
   ! ▽^2 ポロイダル速度時間変化項
    real(8), dimension((nm+1)**2,0:lm), intent(OUT) :: wt_DtDLaplaPolvel

    real(8), dimension(0:im-1,jm,0:km) :: xyz_VGradxyz_Vlon   ! v・▽v(経度)
    real(8), dimension(0:im-1,jm,0:km) :: xyz_VGradxyz_Vlat   ! v・▽v(緯度)
    real(8), dimension(0:im-1,jm,0:km) :: xyz_VGradxyz_Vrad   ! v・▽v(動径)

   ! v・▽v 非線形項計算
    call wt_VGradV(xyz_VGradxyz_Vlon,xyz_VGradxyz_Vlat,xyz_VGradxyz_Vrad, &
                   xyz_Vlon, xyz_Vlat, xyz_Vrad )

  ! トロイダル時間変化
    wt_DtDTorvel = &
           wt_L2Inv_wt(  RsInv * wt_KxRGrad_wt(wt_Torvel)   &
                       - RsInv * wt_Qoperator_wt(wt_Polvel) ) &
         - wt_L2Inv_wt(wt_RadRot_xyz_xyz(xyz_VGradxyz_Vlon,xyz_VGradxyz_Vlat )) 
  ! ポロイダル時間変化
    wt_DtDLaplaPolvel = &
           wt_L2Inv_wt(  RsInv * wt_KxRGrad_wt(wt_LaplaPolvel) &
                       + RsInv * wt_Qoperator_wt(wt_Torvel) )  &
         + wt_L2Inv_wt(wt_RadRotRot_xyz_xyz_xyz(                 &
               xyz_VGradxyz_Vlon, xyz_VGradxyz_Vlat, xyz_VGradxyz_Vrad ))

  end subroutine TimeDerivetives_noDiffusion

 !------------------- 拡散項 ----------------------
  subroutine CNDiffusionMatrix( Diffc, dt, DiffLUMatrix, kpivot )

    real(8), intent(IN)          :: Diffc    ! 拡散係数
    real(8), intent(IN)          :: dt       ! 時間刻

    ! Crank Nicholson 拡散陰的計算用行列(1-D dt/2▽^2, LU 分解)
    real(8), dimension((nm+1)*(nm+1),0:lm,0:lm), intent(OUT)  :: DiffLUMatrix
    ! ピボット情報
    integer, dimension((nm+1)*(nm+1),0:lm), intent(OUT)       :: kpivot

    ! 作業用変数
    real(8), dimension((nm+1)*(nm+1),0:lm)           :: wt_I

    integer :: l

    DiffLUMatrix = 0.0
    do l=0,lm
       wt_I = 0.0 ; wt_I(:,l) = 1.0             ! 各波数成分独立
       DiffLUMatrix(:,:,l) = wt_I - Diffc * dt/2.0 * wt_Lapla_wt(wt_I) 
    enddo

    call ludecomp(DiffLUMatrix,kpivot)

  end subroutine CNDiffusionMatrix

 !------------------- 初期値設定(ファイルなし) ----------------------
  subroutine set_initial_values(xyz_Torvel, xyz_Polvel )

    real(8), dimension(0:im-1,jm,0:km), intent(out)  :: xyz_Torvel ! トロイダル速度
    real(8), dimension(0:im-1,jm,0:km), intent(out)  :: xyz_Polvel ! ポロイダル速度

   ! 静止速度場
    xyz_Torvel =  0.0D0
    xyz_Polvel =  0.0D0

  end subroutine set_initial_values

 !------------------- 初期値設定(ファイルから) ----------------------  
  subroutine initial_read_gtool4( filename, time, &
       xyz_Torvel, xyz_Polvel, wt_DtDTorvel, wt_DtDLaplaPolvel )


    character(len=*), intent(in)   :: filename
    real, intent(in)               :: time
    real(8), dimension(0:im-1,jm,0:km), intent(out)  :: xyz_Torvel  ! トロイダル速度
    real(8), dimension(0:im-1,jm,0:km), intent(out)  :: xyz_Polvel  ! ポロイダル速度
 ! トロイダル速度時間変化 
    real(8), dimension((nm+1)**2,0:lm), intent(out) :: wt_DtDTorvel
 ! ▽^2ポロイダル速度時間変化
    real(8), dimension((nm+1)**2,0:lm), intent(out) :: wt_DtDLaplaPolvel

 ! 読み込み用配列
    real(8), dimension(0:im-1,jm,0:km) :: xyz_DtDVar

    call HistoryGet( trim(filename), 'torvel', xyz_Torvel, time )
    call HistoryGet( trim(filename), 'polvel', xyz_Polvel, time )

    call HistoryGet( trim(filename), 'DtDTorvel', xyz_DtDVar, time )
      wt_DtDTorvel = wt_xyz(xyz_DtDVar)
    call HistoryGet( trim(filename), 'DtDLaplaPolvel', xyz_DtDVar, time )
      wt_DtDLaplaPolvel = wt_xyz(xyz_DtDVar)

  end subroutine initial_read_gtool4

 !------------------- 出力 ----------------------
  subroutine output_gtool4_init

   !---- ヒストリーファイル作成 ----
    call HistoryCreate( &
           file=trim(output_file), &
           title=trim(title), &
           source='NS_rot-3d-shell_abcn_restart_namelist.f90 (2008/08/10)', &
           institution='GFD_Dennou Club SPMODEL project',     &
           dims=(/'lon','lat','rad','t  '/), dimsizes=(/im,jm,km+1,0/),&
           longnames=(/'Longitude','Latitude ','Radius   ','time     '/),&
           units=(/'1','1','1','1'/),   &
           origin=real(time), interval=real(ndisp*delta_t), &
           xtypes=(/'double','double','double','double'/))

   !---- 座標変数定義, 出力 ----
    call HistoryPut('lon',x_Lon/pi*180)                       ! 変数出力
    call HistoryAddattr('lon','topology','circular')          ! 周期属性
    call HistoryAddattr('lon','modulo',360.0)                 ! 周期属性
    call HistoryPut('lat',y_Lat/pi*180)                       ! 変数出力
    call HistoryPut('rad',z_Rad)                              ! 変数出力

   !---- 物理変数定義 ----
    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='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='DtDTorvel', dims=(/'lon','lat','rad','t  '/), & 
           longname='time variation of velocity toroidal potential', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='DtDLaplaPolvel', dims=(/'lon','lat','rad','t  '/), & 
           longname='time variation of laplacian of velocity poloidal potential', &
           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='um', dims=(/'lat','rad','t  '/), & 
           longname='zonal flow', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='hel', dims=(/'lon','lat','rad','t  '/), & 
           longname='helicity', units='1', xtype='double')

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

    call HistoryAddAttr('lon', '+Rossby',   Rs    )
    call HistoryAddAttr('lon', '+Reynolds', Re    )

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

   end subroutine output_gtool4_init

  subroutine output_gtool4
    write(6,*) 'it = ',it, '  time = ', time
    call HistoryPut('t',real(time))

   !---- 物理変数出力 ----
    call HistoryPut('torvel',xyz_Torvel)
    call HistoryPut('polvel',xyz_Polvel)
    call HistoryPut('laplapolvel',xyz_LaplaPolvel)
    call HistoryPut('vlon',xyz_Vlon)
    call HistoryPut('vlat',xyz_Vlat)
    call HistoryPut('vrad',xyz_Vrad)

   !---- 時間変化項出力 ----
    call HistoryPut('DtDTorvel',xyz_wt(wt_DtDTorvel1))
    call HistoryPut('DtDLaplaPolvel',xyz_wt(wt_DtDLaplaPolvel1))

   !---- 診断量出力 ----
    xyz_Zlon = xyz_RotLon_wt_wt(wt_xyz(xyz_Vrad),wt_xyz(xyz_Vlat))
    xyz_Zlat = xyz_RotLat_wt_wt(wt_xyz(xyz_Vlon),wt_xyz(xyz_Vrad))
    xyz_Zrad = xyz_wt(wt_RotRad_xyz_xyz(xyz_Vlat,xyz_Vlon))

    call HistoryPut('zlon',xyz_Zlon)
    call HistoryPut('zlat',xyz_Zlat)
    call HistoryPut('zrad',xyz_Zrad)

    call HistoryPut('ek',AvrLonLatRad_xyz((xyz_Vlon**2+xyz_Vlat**2+xyz_Vrad**2)/2))
    call HistoryPut('um',yz_avrLon_xyz(xyz_Vlon))
    call HistoryPut('hel',xyz_Vlon*xyz_Zlon + xyz_Vlat*xyz_Zlat + xyz_Vrad*xyz_Zrad)
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program NS_rot_3d_shell_abcn_restart_namelist_rotbndry
