!--------------------------------------------------------------------------
!  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 1st order Euler scheme. 
!       Program is organized for stability problems of zonal flows. 
!
!History(履歴):
!   2005/07/22  竹広真一, 新規作成
!   2005/08/09  竹広真一, 渦度振幅を Bains(1976) にあわせるための修正.
!                         ISPACK は 2 で規格化されたルジャンドル陪函数.
!                         Bains(1976) の定式化は 1 で規格化されたもの. 
!   2005/08/10  竹広真一, ポテンシャル渦度, およびその緯度経度出力設定
!   2005/08/11  竹広真一, ポテンシャル渦度緯度傾度 変数名変更
!   2005/08/13  竹広真一, ポテンシャル渦度緯度傾度 半径分の1 追加
!                         座標重み出力追加
!   2005/09/05  柿並義宏, 診断量 (mvlon) 出力不具合修正．
!   2005/09/10  竹広真一, dcmodel 規約に従い大幅書き換え
!   2005/10/02  竹広真一, ファイル名変更
!   2005/10/05  柿並義宏, m_Esp, m_EnsSp 追加
!   2008/08/10  竹広真一, 配列添え字変更, im -> 0:im-1
!
! 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_euler_zonalflow

!== モジュール引用宣言 ================================================
  use w_module 
  use gt4_history, only : GT_HISTORY, HistoryCreate, HistoryPut, HistoryGet, &
                          HistoryAddVariable, HistoryClose, HistoryAddAttr
  use dc_trace,    only : SetDebug, BeginSub, EndSub, DbgMessage
  use dc_message,  only : MessageNotify
  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_DVorDt(:)         ! 渦度時間変化(t)

  real(8), allocatable :: nm_ESp(:,:)         ! エネルギースペクトル
  real(8), allocatable :: nm_EnsSp(:,:)       ! エンストロフィースペクトル
  real(8), allocatable :: n_ESp(:)            ! エネルギースペクトル
  real(8), allocatable :: m_ESp(:)            ! エネルギースペクトル
  real(8), allocatable :: n_EnsSp(:)          ! エンストロフィースペクトル
  real(8), allocatable :: m_EnsSp(:)          ! エンストロフィースペクトル

 !---- 固定パラメタ -----
  real(8), parameter  :: pi = 3.141592653589793D0 ! 円周率
  character(len=20)   :: DbgMessageFmt='*** DbgMESSAGE ***'
  real(8), parameter  :: vmiss = -999.0       ! 欠損値

 !---- NAMELIST 変数 ----
  ! NAMELIST 入力用テンポラリファイル名(名前変更禁止)
  character(len=30),parameter  :: nmlfile='spbaro_euler_zonalflow.nml'

  logical            :: Verbose=.false.          ! 冗長なメッセージ出力 On/Off
  logical            :: DebugOn=.false.          ! デバッグメッセージ On/Off
  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(len=100) :: hst_file=   ''           ! ヒストリーファイル名
  character(len=100) :: title = &                ! タイトル
       'Stability of zonal flows of 2-dim barotropic fluid on a rotating sphere'
  integer :: hst_intstep=200                     ! ヒストリー出力間隔ステップ数
  namelist /history/    hst_file, title, hst_intstep

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

 !---- 作業変数 ----
  real(8) :: Beta=0.0d0                      ! βパラメター(2*Omega/radius**2)
  real(8), allocatable :: w_HVisc(:)         ! 超粘性係数
  real(8), allocatable :: n_Damptime(:)      ! 減衰時間(e-folding time)

  integer :: it=0                            ! 時間ステップ
  real(8) :: time                            ! モデル内時間
  integer :: n, m                            ! 全波数, 帯状波数

  type(GT_HISTORY) :: hst_rst                ! リスタート GT_HISTORY 変数

 !---------------- NAMELIST 読み込み ---------------------

  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=history) 
  open(10,file=nmlfile,status='OLD')
  read(10,nml=history) ; write(6,nml=history) ; close(10)

  if (verbose) write(6,nml=restart) 
  open(10,file=nmlfile,status='OLD')
  read(10,nml=restart) ; write(6,nml=restart) ; 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_DVorDt((nm+1)*(nm+1)))
  allocate(w_HVisc((nm+1)*(nm+1)))

  allocate(nm_ESp(0:nm,-nm:nm),nm_EnsSp(0:nm,-nm:nm))
  allocate(m_ESp(-nm:nm),m_EnsSp(-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**2
  w_HVisc = (-1)**HVOrder * HVisc*((rn(:,1)+2.0D0)/Radius**2)**HVOrder
  w_HVisc(l_nm(0,0)) = vmiss   ! 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

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

  if ( initial_file == "") then
     ! リスタートファイルを指定しない場合内部で w_Vor を与える. 
     call set_initial_values
  else
     !  初期値設定(リスタートファイルからの読みこみ)
     call HistoryGet( trim(initial_file), 'w_vor', w_Vor,  time )
  endif

  w_Strfunc = w_LaplaInv_w(w_Vor) * Radius**2


 !------------------- 時間積分(Euler 法) --------------------
  call output_restart_init
  call output_history_init

  if ( initial_file == '' ) call output_history    ! 内部で与えた初期値は出力

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

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

     ! 時間変化項の見積り
     w_DVorDt = - w_Jacobian_w_w(w_StrFunc,w_Vor)/Radius**2 &
                - Beta * w_DLon_w(w_StrFunc)                &
                - w_HVisc * w_Vor

     ! 時間積分計算(Euler scheme)
     w_Vor = w_Vor + delta_t * w_DVorDt

     ! 渦度から流線関数を計算
     w_StrFunc = w_LaplaInv_w(w_Vor) * Radius**2

     if(mod(it,hst_intstep) .eq. 0)then                    ! ヒストリー出力
        call output_history
     endif

     if(mod(it,rst_intstep) .eq. 0)then                    ! リスタート出力
        call output_restart
     endif
  enddo

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

  if(.not. mod(it-1,rst_intstep) .eq. 0)then                    ! 最終出力
     call output_restart
  endif

  call output_restart_close
  call output_history_close

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

contains

!=========================== 初期値設定 ============================
 !
 !  初期値設定(リスタートファイルない場合のデフォルト設定)
 !
  subroutine set_initial_values    ! w_Vor を定める

   !---- 初期値帯状流 + ランダムホワイトノイズ
    integer :: n=4                      ! 基本流の水平全波数
    real(8) :: vor_amplitude =1.0D0     ! 渦度振幅
    real(8) :: noise_amplitude=1.0D-6   ! 微小ホワイトノイズの振幅
    integer :: seed=0                   ! seed(1)に設定するランダム種の値

    namelist /initvor/ n,vor_amplitude, noise_amplitude, seed

    integer              :: random_seed_size   ! 乱数の種の長さ
    integer, allocatable :: seedarray(:)       ! 乱数の種
    real, allocatable    :: harvest(:)         ! 乱数の値

    if (verbose) write(6,nml=initvor) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=initvor) ; write(6,nml=initvor) ; 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(harvest((nm+1)*(nm+1)))
    call random_number(harvest)

    w_Vor = 0.0d0
    w_Vor(l_nm(n,0))=vor_amplitude/sqrt(2.0d0)             ! 帯状流設定
    w_Vor = w_Vor + harvest * noise_amplitude/sqrt(2.0d0)  ! ノイズ設定
    w_Vor(l_nm(0,0)) = 0.0d0                   ! 意味のない n,m=(0,0) の成分は0

    deallocate(seedarray,harvest)

  end subroutine set_initial_values

!=========================== リスタート出力 ============================
 !
 ! リスタート出力初期化
 !
  subroutine output_restart_init
    call HistoryCreate( &
           file=trim(rst_file), &
           title=trim(title), &
           source='spbaro_euler_zonalflow.f90 (2008/08/10)', &
           institution='GFD_Dennou Club SPMODEL project',     &
           dims=(/'lon','lat','nm ','t  '/), &
           dimsizes=(/im,jm,(nm+1)**2,0/),&
           longnames=(/'Longitude            ','Latitude             ',&
                       'Hor.wave number index','time                 '/),&
           units=(/'radian','radian','1     ','1     '/),   &
           origin=real(time), interval=real(rst_intstep*delta_t), &
           xtypes=(/'real'/), history=hst_rst)

   !---- 座標変数定義, 出力 ----
    call HistoryPut('lon',x_Lon, hst_rst)                     ! 変数出力
    call HistoryAddattr('lon','topology','circular', hst_rst) ! 周期属性
    call HistoryAddattr('lon','modulo',2*pi, hst_rst)         ! 周期属性
    call HistoryPut('lat',y_Lat, hst_rst)                     ! 変数出力
    call HistoryPut('nm',(/(dble(n),n=0,(nm+1)**2)/), hst_rst)! 変数出力

    call HistoryAddVariable( &                                ! 変数定義
           varname='lon_weight', dims=(/'lon'/), & 
           longname='weight for integration in longitude', &
           units='radian', xtype='double',history=hst_rst)
    call HistoryAddVariable( &                                ! 変数定義
           varname='coslat_lat_weight', dims=(/'lat'/), & 
           longname='cos(lat) weight for integration in latitide', &
           units='1', xtype='double',history=hst_rst)
    call HistoryPut('lon_weight',x_Lon_weight,hst_rst)        ! 変数出力
    call HistoryPut('coslat_lat_weight',y_Lat_weight,hst_rst) ! 変数出力

   !---- 物理変数定義 ----
    call HistoryAddVariable( &                                ! 変数定義
           varname='w_vor', dims=(/'nm','t '/), & 
           longname='Vorticity', &
           units='1', xtype='double', history=hst_rst)

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

  end subroutine output_restart_init

 !
 ! リスタート出力
 !
  subroutine output_restart
    write(6,*) ' Restart file output at it = ',it, '  time = ', time
    call HistoryPut('t',real(time),hst_rst)

   !---- 物理変数出力 ----
    call HistoryPut('w_vor', w_Vor, 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='spbaro_euler_zonalflow.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=(/'degree','degree','1     ','1     ','1     ','1     '/),   &
           origin=real(time), interval=real(hst_intstep*delta_t), &
           xtypes=(/'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 for integration in longitude', &
           units='radian', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='coslat_lat_weight', dims=(/'lat'/), & 
           longname='cos(lat) weight for integration in latitide', &
           units='1', xtype='double')
    call HistoryPut('lon_weight',x_Lon_weight)                ! 変数出力
    call HistoryPut('coslat_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='pv', dims=(/'lon','lat','t  '/), & 
           longname='Potential vorticity', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='pvgradlat', dims=(/'lon','lat','t  '/), & 
           longname='Latitudinal gradient of PV', 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='ekspm', dims=(/'m','t'/), & 
           longname='Energy spectrum (m space)', units='1', xtype='double')
    call HistoryAddAttr('ekspm','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='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='ensspm', dims=(/'m','t'/), & 
           longname='Enstrophy spectrum (m-space)', units='1', xtype='double')
    call HistoryAddAttr('ensspm','missing_value', vmiss  )

    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')

   !---- 実験パラメターを属性として定義, 出力(全て 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_history_init

 !
 ! ヒストリー出力
 !
  subroutine output_history
    write(6,*) ' History file output at 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)*Radius**2)
    call HistoryPut('ens',&
         AvrLonLat_xy(xy_Vor**2/2.0d0)*Radius**2)
    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))
    call HistoryPut('pv',2*Omega*sin(xy_Lat)+xy_Vor)
    call HistoryPut('pvgradlat',&
         (2*Omega*cos(xy_Lat)+xy_GradLat_w(w_Vor))/Radius)

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

    do m=-nm,nm
       m_ESp(m) = sum(nm_ESp(abs(m):nm,m))
    end do
    do m=1,nm
       m_ESp(m)  = m_ESp(m) + m_ESp(-m)     ! m >= 0 に集める
       m_ESp(-m) = vmiss                    ! m < 0 は欠損値
    end do

    call HistoryPut('ekspnm',nm_ESp)
    call HistoryPut('ekspn',n_ESp)
    call HistoryPut('ekspm',m_Esp)
    call HistoryPut('eksptot',(4*pi)*sum(n_Esp))
    call HistoryPut('eksp',sum(n_Esp)/Radius**2)

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

    do m=-nm,nm
       m_EnsSp(m) = sum(nm_ESp(abs(m):nm,m))
    end do
    do m=1,nm
       m_EnsSp(m)  = m_EnsSp(m) + m_EnsSp(-m)     ! m >= 0 に集める
       m_EnsSp(-m) = vmiss                        ! m < 0 は欠損値
    end do

    call HistoryPut('ensspn',n_EnsSp)
    call HistoryPut('ensspm',m_EnsSp)
    call HistoryPut('ensspnm',nm_EnsSp)
    call HistoryPut('enssptot',(4*pi/Radius**2)*sum(n_Enssp))
    call HistoryPut('enssp',sum(n_Enssp))

  end subroutine output_history

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

end program spbaro_euler_zonalflow
