!--------------------------------------------------------------------------
!  Copyright (C) 2005--2008 SPMODEL Development Group. All rights reserved.
!--------------------------------------------------------------------------
!
!Title(表題): Two-dimensional barotropic model on a rotating sphere. 
!             (回転球面上の2次元順圧モデル)
!
!       The time integration is performed with
!       the 4th order Runge-Kutta scheme.
!       The viscous linear term ant the beta term are transformed 
!         by the exp., cos, sin functions.  
!       Program is organized for decaying tubulence problem. 
!
!History(履歴):  2005/06/24  S. Takehiro, created
!                2005/08/24  S. Takehiro, lon_weight, lat_weight output. 
!                2005/09/04  S. Takehiro, lon_weight, lat_weight modified.
!                2006/05/01  S. Takehiro, xtypes in HistoryCreate fixed.
!                2007/12/07  S. Takehiro, namelist file is set from 
!                                         the command line argument. 
!                2008/02/11  S. Takehiro, Freq of Rossby waves fixed.
!                2008/02/13  S. Takehiro, Energy and Enstrophy output fixed.
!                2008/02/14  S. Takehiro, Freq of Rossby waves fixed again.
!                2008/02/26  S. Takehiro, Energy, Enstrophy dissipation added.
!                2008/08/10  S. Takehiro, index changed, im -> 0:im-1
!                2018/11/24  S. Takehiro, w_HVisc(l_nm(0,0)) fixed
!
! The governing equation(支配方程式): 
!
!      d\nabla^2\psi/dt 
!          = -J(\psi,\nabla^2\psi)/a^2 - 2\Omega/a^2 d\psi/d\lambda 
!            +(-1)^{p+1}\nu_{2p}(\nabla^2 + 2/a^2)^p \nabla^2\psi
!
program spbaro_rn4expcs_freedecay

!== モジュール引用宣言 ================================================
  use w_module 
  use gt4_history, only : HistoryCreate, HistoryPut, HistoryGet, &
                          HistoryAddVariable, HistoryClose, HistoryAddAttr
  use dc_trace,    only : SetDebug, BeginSub, EndSub, DbgMessage
  use dc_message,  only : MessageNotify
  use dc_types
  use dc_string,   only : StoA
  use dc_args
  implicit none

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

 !---- 変数(格子点データ) ----
  real(8), allocatable :: xy_VelLon(:,:)      ! 速度経度成分
  real(8), allocatable :: xy_VelLat(:,:)      ! 速度緯度成分
  real(8), allocatable :: xy_Vor(:,:)         ! 渦度(鉛直成分)
  real(8), allocatable :: xy_StrFunc(:,:)     ! 流線関数

 !---- 変数(スペクトルデータ) ----
  real(8), allocatable :: w_Vor(:)            ! 渦度(鉛直成分)
  real(8), allocatable :: w_StrFunc(:)        ! 流線関数

  real(8), allocatable :: w_Xi(:)             ! ξ : ζ=ξexp[-(iω+ν)t]
  real(8), allocatable :: w_Vortmp(:)         ! 渦度(ルンゲクッタ計算作業用)
  real(8), allocatable :: w_k1(:)             ! 拡散以外の渦度時間変化(1段目)
  real(8), allocatable :: w_k2(:)             ! 拡散以外の渦度時間変化(2段目)
  real(8), allocatable :: w_k3(:)             ! 拡散以外の渦度時間変化(3段目)
  real(8), allocatable :: w_k4(:)             ! 拡散以外の渦度時間変化(4段目)

  real(8), allocatable :: nm_ESp(:,:)         ! エネルギースペクトル
  real(8), allocatable :: nm_EnsSp(:,:)       ! エンストロフィースペクトル
  real(8), allocatable :: n_ESp(:)            ! エネルギースペクトル
  real(8), allocatable :: n_EnsSp(:)          ! エンストロフィースペクトル
 
 !---- 固定パラメタ -----
  real(8), parameter  :: pi = 3.141592653589793D0 ! 円周率
  character(len=20)   :: DbgMessageFmt='*** DbgMESSAGE ***'
  real(8), parameter  :: vmiss = -999.0       ! 欠損値

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

  ! NAMELIST 入力用ファイル名
  character(STRING)  :: nmlfile

  logical            :: Verbose=.false.          ! 出力メッセージレベル
  logical            :: DebugOn=.false.          ! デバッグ出力コントロール
  namelist /message/  Verbose, DebugOn           !
                                                 !-- 格子点・スペクトル --
  integer :: nm=21                               ! 切断全波数
  integer :: im=64                               ! 経度方向格子点数 (>3*nm+1)
  integer :: jm=32                               ! 緯度方向格子点数 (>3*nm/2)
  namelist /gridset/ nm, im, jm
                                                 !-- 物理パラメター --
  real(8)            :: Radius=1.0D0             ! 球の半径
  real(8)            :: Omega=0.0D0              ! 回転角速度
  integer            :: HVOrder=1                ! 超粘性の次数(1 で普通の拡散)
  real(8)            :: HVisc=1.0D0              ! 超粘性係数
  namelist /physics/  Radius, Omega, HVOrder, HVisc
                                                 ! -- 初期値 --
  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(STRING) :: output_file=''            ! 出力ファイル名
  character(STRING) :: title = &                 ! タイトル
       'Decaying turbulence of 2-dim barotropic fluid on a rotating sphere'
  character(len=100) :: source = &               ! ソース名
       'spbaro_rn4expcs_freedecay.f90 (2018/11/24)'!
  integer :: ndisp=200                           ! 出力間隔ステップ数
  namelist /output/    output_file, title, ndisp ! 

 !---- 作業変数 ----
  real(8) :: Beta=0.0d0                      ! βパラメター(2*Omega/radius)
  real(8), allocatable :: w_HVisc(:)         ! 超粘性係数
  real(8), allocatable :: n_Damptime(:)      ! 減衰時間(e-folding time)
  real(8), allocatable :: w_OmegaR(:)        ! ロスビー波の振動数
                                             ! (- Omega * abs(m)/(n*(n+1))

  integer :: it=0                            ! 時間ステップ
  real(8) :: time                            ! モデル内時間
  integer :: n, m                            ! 全波数, 帯状波数
  real(8) :: EkDispRate, EkDisp              ! エネルギー散逸率, 全散逸
  real(8) :: EnsDispRate, EnsDisp            ! エンストロフィー散逸率, 全散逸

 !---------------- NAMELIST 読み込み ---------------------
  call MessageNotify('M','main', trim(source))

  call Open(arg)
  call Debug(arg) ; call Help(arg) ; call Strict(arg)
  call Get(arg, argv)
  if ( size(argv) .le. 0 ) then
     call MessageNotify('M','main',&
          'Usage: spbaro_rn4excs_freedecay.out [namelist file]')
     call MessageNotify('E','main','There is no argument. &
          & Please set the namelist file name after the execution.')
  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) ; write(6,nml=message) ; close(10)

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

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

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

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

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

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

 !------------------ 変数の割り付け ---------------------
   allocate(xy_VelLon(0:im-1,jm),xy_VelLat(0:im-1,jm))
   allocate(xy_Vor(0:im-1,jm),xy_StrFunc(0:im-1,jm))

   allocate(w_Vor((nm+1)*(nm+1)),w_StrFunc((nm+1)*(nm+1)))
   allocate(w_Xi((nm+1)*(nm+1)))
   allocate(w_Vortmp((nm+1)*(nm+1)))
   allocate(w_k1((nm+1)*(nm+1)),w_k2((nm+1)*(nm+1)))
   allocate(w_k3((nm+1)*(nm+1)),w_k4((nm+1)*(nm+1)))
   allocate(w_HVisc((nm+1)*(nm+1)))
   allocate(w_OmegaR((nm+1)*(nm+1)))

   allocate(nm_ESp(0:nm,-nm:nm),nm_EnsSp(0:nm,-nm:nm))
   allocate(n_ESp(0:nm),n_EnsSp(0:nm))
   allocate(n_Damptime(0:nm))

 !------------------ 座標値の設定 -----------------------
  call DbgMessage(fmt='call %c', c1='w_initial') 
  call w_Initial(nm,im,jm)
  w_spectrum_VMiss = vmiss

 !------------------ 物理係数の設定 -----------------------
  Beta  = 2 *Omega/Radius
  w_HVisc = (-1)**(HVOrder)*HVisc*((rn(:,1)+2.0D0)/Radius**2)**HVOrder
  w_HVisc(l_nm(0,0)) = 0.0D0   ! rn(ln(0,0,1) は正の値なので修正しておく.

  n_Damptime(0)=vmiss                           ! 欠損値
  do n=1,nm
     if ( w_HVisc(l_nm(n,0)) .NE. 0.0D0 ) then
        n_Damptime(n)=1.0D0/w_HVisc(l_nm(n,0))
     else
        n_Damptime(n)=vmiss                     ! 欠損値
     endif
  enddo

  w_OmegaR = vmiss
  do n=1,nm
     do m=-n,n
        w_OmegaR(l_nm(n,m)) = - 2.0D0 *Omega * abs(m)/(n*(n+1))
     enddo
  enddo

 !------------------- 初期値設定 ----------------------
  time = initial_time

  if ( initial_file == "") then
     call set_initial_values       ! xy_Vor, xy_Strfunc を定める
  else
     call initial_read_gtool4      ! xy_Vor, xy_Strfunc を読みこむ
  endif

  w_Vor = w_xy(xy_Vor) ; w_Strfunc = w_xy(xy_Strfunc) 
  EkDisp = 0.0 ; EnsDisp = 0.0

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

  call DbgMessage(fmt='%c %c', &
&                 c1=DbgMessageFmt, &
&                 c2='Time integration starts.') 

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

     !---- 1 段目 [ k1 = f( x_n, t_n) ] ----
     w_k1 = w_DXiDt_w_w(w_Vor,w_Strfunc,0.0D0)

     !---- 2 段目 [ k2 = f( x_n+k1*dt/2, t_n+dt/2 ) ] ----
     w_Xi = w_Vor + w_k1*delta_t/2.0D0
     w_Vortmp = w_Xi2Vor_w(w_Xi,delta_t/2.0D0)
     w_StrFunc = w_LaplaInv_w(w_Vortmp) * Radius**2
     w_k2 = w_DXiDt_w_w(w_Vortmp,w_Strfunc,delta_t/2.0D0)

     !---- 3 段目 [ k3 = f( x_n+k2*dt/2, t_n+dt/2 ) ] ----
     w_Xi = w_Vor + w_k2*delta_t/2.0D0
     w_Vortmp = w_Xi2Vor_w(w_Xi,delta_t/2.0D0)
     w_StrFunc = w_LaplaInv_w(w_Vortmp) * Radius**2
     w_k3 = w_DXiDt_w_w(w_Vortmp,w_Strfunc,delta_t/2.0D0)

     !---- 4 段目 [ k4 = f( x_n+k3*dt, t_n+dt ) ] ----
     w_Xi = w_Vor + w_k3*delta_t
     w_Vortmp = w_Xi2Vor_w(w_Xi,delta_t)
     w_StrFunc = w_LaplaInv_w(w_Vortmp) * Radius**2
     w_k4 = w_DXiDt_w_w(w_Vortmp,w_Strfunc,delta_t)

     !---- 積分 ----
     w_Xi =w_Vor + delta_t * (   w_k1/6.0D0 + w_k2/3.0D0 &
                               + w_k3/3.0D0 + w_k4/6.0D0    )
     w_Vor = w_Xi2Vor_w(w_Xi,delta_t)
     w_StrFunc = w_LaplaInv_w(w_Vor) * Radius**2

     !---- エネルギー・エンストロフィー散逸計算 ----
     EkDispRate=0.0
     do n=0,nm
        do m=-n,n
           EkDispRate = EkDispRate &
              + HVisc * n*(n+1)/Radius**2 * ((n*(n+1)-2)/Radius**2)**HVorder &
                      * w_Strfunc(l_nm(n,m))**2
        enddo
     enddo

     EnsDispRate=0.0
     do n=0,nm
        do m=-n,n
           EnsDispRate = EnsDispRate &
                + HVisc* n**2 * (n+1)**2/Radius**4 &
                      * ((n*(n+1)-2)/Radius**2)**HVorder &
                      * w_Strfunc(l_nm(n,m))**2
        enddo
     enddo

     EkDisp = Ekdisp + EkDispRate * delta_t
     EnsDisp = Ensdisp + EnsDispRate * delta_t

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

  call DbgMessage(fmt='%c %c',      &
&                 c1=DbgMessageFmt, &
&                 c2='Time integration end.') 

  call output_gtool4_close

! 以上 メインプログラム 
!-----------------------------------------------------------------------------
! 以下 サブルーチン

contains

 !------------------- 時間変化項(非線形項のみ) ----------------------
  function w_DtDVorNonlinear_w_w(w_Vor,w_Strfunc)
    real(8)  :: w_Vor(:)                           ! 渦度(鉛直成分)
    real(8)  :: w_StrFunc(:)                       ! 流線関数
    real(8)  :: w_DtDVorNonlinear_w_w(size(w_Vor)) ! 渦度時間変化(非線形項のみ)

    w_DtDVorNonlinear_w_w = - w_Jacobian_w_w(w_StrFunc,w_Vor)/Radius**2

  end function w_DtDVorNonlinear_w_w

 !------------------- ξ から ζ を求める ----------------------
  function w_Xi2Vor_w(w_Xi,dt)
    real(8), intent(IN)  :: w_Xi(:)                ! ξ : ζ=ξexp[-(iω+ν)t]
    real(8), intent(IN)  :: dt                     ! 時間刻
    real(8)              :: w_Xi2Vor_w(size(w_Xi)) ! 渦度ζ

    w_Xi2Vor_w(l_nm(0,0)) = 0.0D0
    do n=1,nm
       do m=0,n
          w_Xi2Vor_w(l_nm(n,m)) &
               =   w_Xi(l_nm(n,m)) &
                 * exp(-w_HVisc(l_nm(n,m))*dt)*cos(w_OmegaR(l_nm(n,m))*dt) &
                 + w_Xi(l_nm(n,-m)) &
                 * exp(-w_HVisc(l_nm(n,m))*dt)*sin(w_OmegaR(l_nm(n,m))*dt)
       enddo
       do m=1,n
          w_Xi2Vor_w(l_nm(n,-m)) &
               = w_Xi(l_nm(n,-m)) &
               * exp(-w_HVisc(l_nm(n,-m))*dt)*cos(w_OmegaR(l_nm(n,-m))*dt) &
               - w_Xi(l_nm(n,m)) &
               * exp(-w_HVisc(l_nm(n,-m))*dt)*sin(w_OmegaR(l_nm(n,-m))*dt)
       enddo
    enddo

  end function w_Xi2Vor_w

 !---------------- ξ(ζ=ξexp[-(iω+ν)t])の時間変化 -------------------
  function w_DXiDt_w_w(w_Vor,w_Strfunc,dt)
    real(8), intent(IN)  :: w_Vor(:)                 ! 渦度ζ
    real(8), intent(IN)  :: w_Strfunc(:)             ! 流線関数
    real(8), intent(IN)  :: dt                       ! 時間刻
    real(8)              :: w_DXiDt_w_w(size(w_Vor)) ! ξ時間変化

    real(8)              :: w_G(size(w_Vor))         ! 作業領域
    integer n,m 

    w_G = w_DtDVorNonlinear_w_w(w_Vor,w_Strfunc)

    w_DXiDt_w_w(l_nm(0,0)) = 0.0D0
    do n=1,nm
       do m=0,n
          w_DXiDt_w_w(l_nm(n,m))   &
            =   exp(w_HVisc(l_nm(n,m))*dt)*cos(w_OmegaR(l_nm(n,m))*dt) &
              * w_G(l_nm(n,m)) &
              - exp(w_HVisc(l_nm(n,m))*dt) * sin(w_OmegaR(l_nm(n,m))*dt) &
              * w_G(l_nm(n,-m))
       enddo
       do m=1,n
          w_DXiDt_w_w(l_nm(n,-m)) &
            =   exp(w_HVisc(l_nm(n,-m))*dt) * cos(w_OmegaR(l_nm(n,-m))*dt) &
              * w_G(l_nm(n,-m))            &
              + exp(w_HVisc(l_nm(n,-m))*dt) * sin(w_OmegaR(l_nm(n,-m))*dt) &
              * w_G(l_nm(n,m))
       enddo
    enddo

  end function w_DXiDt_w_w

 !------------------- 初期値設定(ファイルなし) ----------------------
  subroutine set_initial_values    ! xy_Vor, xy_Strfunc を定める

   !---- 初期値剛体回転場
!!$    real(8), parameter :: Omega0   = 1.0D0
!!$    real(8), parameter :: LatOmega = 90.0D0/180.0D0*pi
!!$    real(8), parameter :: LonOmega = 0.0D0/180.0D0*pi
!!$
!!$    xy_StrFunc = Omega0 * Radius &
!!$                   * ( cos(LatOmega)*cos(xy_Lat)*cos(LonOmega-xy_Lon) &
!!$                       + sin(LatOmega)*sin(xy_Lat) )
!!$    w_StrFunc = w_xy(xy_StrFunc)
!!$    w_Vor = w_Lapla_w(w_StrFunc)/Radius**2
!!$    xy_Vor =  xy_w(w_Vor)
!!$
!!$   !---- 初期値 1 波数成分
!!$    integer, parameter :: ninit=4, minit=3
!!$
!!$    w_StrFunc = 0.0d0 ; w_StrFunc(l_nm(ninit,minit))=1.0d0
!!$    xy_StrFunc = xy_w(w_StrFunc)
!!$    w_Vor = w_Lapla_w(w_StrFunc)/Radius**2
!!$    xy_Vor =  xy_w(w_Vor)

   !---- 初期値乱数分布
    integer              :: random_seed_size   ! 乱数の種の長さ
    integer, allocatable :: seedarray(:)       ! 乱数の種
    real, allocatable    :: harvest(:)         ! 乱数の値
    real(8), allocatable :: n_EspInit(:)       ! 初期エネルギースペクトル分布

    integer    :: seed=0             ! seed(1)に設定する種の値
    integer    :: nzero=10           ! 初期エネルギースペクトル分布のパラメタ
    integer    :: gamma=100          ! 初期エネルギースペクトル分布のパラメタ
    real(8)    :: Etotal=1.0D0       ! 初期全エネルギーの値
    character(len=10) :: disttype='YIHY2002' ! 初期エネルギー分布のタイプ
                                            ! (YIHY2002/YY1993)

    namelist /initvalue/ seed, nzero, gamma, Etotal, disttype

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

    ! 乱数設定
    call random_seed(size=random_seed_size)
    allocate(seedarray(random_seed_size))
    call random_seed(get=seedarray)
    seedarray(1)=seed
    call random_seed(put=seedarray)

    allocate(n_EspInit(0:nm))

    ! 初期のエネルギースペクトル分布の形状を与える
    if ( trim(disttype) .eq. 'YIHY2002' ) then
       call MessageNotify('M','set_initial_values',&
            'Selected initial energy spectrum is "YIHY2002" type. ')
       n_EspInit = (/(dble(n)**(gamma/2)/dble(n+nzero)**gamma,n=0,nm)/)
    elseif ( trim(disttype) .eq. 'YY1993' ) then
       call MessageNotify('M','set_initial_values',&
            'Selected initial energy spectrum is "YY1993" type. ')
       n_EspInit = (/(dble(n)**5/exp(-dble(n)/2),n=0,nm)/)
    else
       call MessageNotify('E','set_initial_values',&
            'The parameter "disttype" should be "YIHY2002" or "YY1993"')
    endif

    w_StrFunc = 0.0d0
    do n=1,nm
       allocate(harvest(-n:n))
       call random_number(harvest)
       w_StrFunc(l_nm(n,(/(m,m=-n,n)/)))=2.0 * harvest - 1 ! [-1,1] の一様乱数
       deallocate(harvest)
    enddo
    n_ESp = n_EnergyFromStreamfunc_w(w_StrFunc)

    ! スペクトル分布の形を決める
    do n=0,nm
       do m=-n,n
          if ( n_ESp(n) .ne. 0.0d0 ) then
             w_StrFunc(l_nm(n,m)) = w_StrFunc(l_nm(n,m)) &
                   * sqrt( n_EspInit(n)/n_ESp(n) )
          endif
       enddo
    enddo

    ! エネルギースペクトル再計算
    n_ESp = n_EnergyFromStreamfunc_w(w_StrFunc)

    ! エネルギーの大きさを定める. 
    w_StrFunc = w_StrFunc * sqrt(Etotal/sum(n_ESp))
    xy_StrFunc = xy_w(w_StrFunc)
    w_Vor = w_Lapla_w(w_StrFunc)/Radius**2
    xy_Vor =  xy_w(w_Vor)

    deallocate(n_EspInit,seedarray)

  end subroutine set_initial_values

 !------------------- 初期値設定(ファイルから) ----------------------  
  subroutine initial_read_gtool4    ! xy_Vor, xy_Strfunc を読みこむ

    call HistoryGet( trim(initial_file), 'vor',    xy_Vor,     time )
    call HistoryGet( trim(initial_file), 'strfunc',xy_StrFunc, time )

    w_Vor = w_xy(xy_Vor)
    w_StrFunc = w_xy(xy_StrFunc)

  end subroutine initial_read_gtool4

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

   !---- ヒストリーファイル作成 ----
    call HistoryCreate( &
           file=trim(output_file), &
           title=trim(title), &
           source='spbaro_rn4expcs_freedecay.f90 (2008/08/10)', &
           institution='GFD_Dennou Club SPMODEL project',     &
           dims=(/'lon','lat','nm ','n  ','m  ','t  '/), &
           dimsizes=(/im,jm,(nm+1)**2,nm+1,2*nm+1,0/),&
           longnames=(/'Longitude            ','Latitude             ',&
                       'Hor.wave number index','Hor.total wave number',&
                       'zonal wave number    ','time                 '/),&
           units=(/'1','1','1','1','1','1'/),   &
           origin=real(time), interval=real(ndisp*delta_t), &
           xtypes=StoA('real','real','int','int','int','real'))

   !---- 座標変数定義, 出力 ----
    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('nm',(/(dble(n),n=0,(nm+1)**2)/))         ! 変数出力
    call HistoryPut('n',(/(dble(n),n=0,nm)/))                 ! 変数出力
    call HistoryPut('m',(/(dble(m),m=-nm,nm)/))               ! 変数出力
    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',y_Lat_weight)                ! 変数出力

   !---- パラメター定義, 出力 ----
    call HistoryAddVariable( &                                ! 変数定義
           varname='hvisc', dims=(/'n'/), & 
           longname='hyper diffusivity', units='1', xtype='double')
    call HistoryAddAttr('hvisc','missing_value', vmiss  )

    call HistoryAddVariable( &                                ! 変数定義
           varname='damptime', dims=(/'n'/), & 
           longname='dissipative damping time (e-folding)', &
           units='1', xtype='double')
    call HistoryAddAttr('damptime','missing_value', vmiss  )

    call HistoryPut('hvisc',w_HVisc(l_nm((/(n,n=0,nm)/),0)))  ! 変数出力
    call HistoryPut('damptime',n_Damptime)                    ! 変数出力

   !---- 物理変数定義 ----
    call HistoryAddVariable( &                                ! 変数定義
           varname='vor', dims=(/'lon','lat','t  '/), & 
           longname='Vorticity', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='strfunc', dims=(/'lon','lat','t  '/), & 
           longname='Stream function', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='vellon', dims=(/'lon','lat','t  '/), & 
           longname='lon-velocity', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='vellat', dims=(/'lon','lat','t  '/), & 
           longname='lat-velocity', units='1', xtype='double')

   !---- 診断量定義 ----
    call HistoryAddVariable( &                                ! 変数定義
           varname='ek', dims=(/'t  '/), & 
           longname='mean kinetic energy', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ens', dims=(/'t  '/), & 
           longname='mean enstrophy', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ektot', dims=(/'t  '/), & 
           longname='Total kinetic energy', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='enstot', dims=(/'t  '/), & 
           longname='Total enstrophy', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='mvlon', dims=(/'lat','t  '/), & 
           longname='Mean zonal flow', units='1', xtype='double')

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

    call HistoryAddVariable( &                                ! 変数定義
           varname='ekspn', dims=(/'n','t'/), & 
           longname='Energy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='eksp', dims=(/'t'/), & 
           longname='Mean energy (by spectrum)', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='eksptot', dims=(/'t'/), & 
           longname='Total energy (by spectrum)', units='1', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='ekdisprate', dims=(/'t'/), & 
           longname='Energy dissipation rage', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ekdisp', dims=(/'t'/), & 
           longname='Accumulated Energy dissipation', &
           units='1', xtype='double')


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

    call HistoryAddVariable( &                                ! 変数定義
           varname='ensspn', dims=(/'n','t'/), & 
           longname='Enstrophy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='enssp', dims=(/'t'/), & 
           longname='Mean enstrophy (by spectrum)', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='enssptot', dims=(/'t'/), & 
           longname='Total enstrophy (by spectrum)', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ensdisprate', dims=(/'t'/), & 
           longname='Enstrophy dissipation rage', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ensdisp', dims=(/'t'/), & 
           longname='Accumulated Enstrophy dissipation', &
           units='1', xtype='double')

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

  end subroutine output_gtool4_init

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

   !---- 物理変数出力 ----
    xy_Vor = xy_w(w_Vor)
    xy_StrFunc = xy_w(w_StrFunc)
    xy_VelLon = -xy_GradLat_w(w_StrFunc)/Radius
    xy_VelLat = xy_GradLon_w(w_StrFunc)/Radius

    call HistoryPut('vellon',xy_VelLon)
    call HistoryPut('vellat',xy_VelLat)
    call HistoryPut('vor',xy_Vor)
    call HistoryPut('strfunc',xy_StrFunc)

   !---- 診断量出力 ----
    call HistoryPut('ek', &
         AvrLonLat_xy((xy_VelLon**2+xy_VelLat**2)/2.0d0))
    call HistoryPut('ens',&
         AvrLonLat_xy(xy_Vor**2/2.0d0))
    call HistoryPut('ektot', &
         IntLonLat_xy((xy_VelLon**2+xy_VelLat**2)/2.0d0)*Radius**2)
    call HistoryPut('enstot',&
         IntLonLat_xy(xy_Vor**2/2.0d0)*Radius**2)
    call HistoryPut('mvlon',y_AvrLon_xy(xy_VelLon))

   !---- エネルギースペクトル診断量出力 ----
    nm_ESp = nm_EnergyFromStreamfunc_w(w_StrFunc)
    n_ESp = n_EnergyFromStreamfunc_w(w_StrFunc)

    call HistoryPut('ekspnm',nm_ESp/Radius**2)
    call HistoryPut('ekspn',n_ESp/Radius**2)
    call HistoryPut('eksptot',(4*pi)*sum(n_Esp))
    call HistoryPut('eksp',sum(n_Esp)/Radius**2)
    call HistoryPut('ekdisprate',EkDispRate)
    call HistoryPut('ekdisp',EkDisp)

   !---- エンストロフィースペクトル診断量出力 ----
    nm_EnsSp = nm_EnstrophyFromStreamfunc_w(w_StrFunc)
    n_EnsSp = n_EnstrophyFromStreamfunc_w(w_StrFunc)

    call HistoryPut('ensspn',n_EnsSp/Radius**4)
    call HistoryPut('ensspnm',nm_EnsSp/Radius**4)
    call HistoryPut('enssptot',(4*pi/Radius**2)*sum(n_Enssp))
    call HistoryPut('enssp',sum(n_Enssp)/Radius**4)
    call HistoryPut('ensdisprate',EnsDispRate)
    call HistoryPut('ensdisp',EnsDisp)

  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program spbaro_rn4expcs_freedecay
