!----------------------------------------------------------------------
!  Copyright (C) 2005--2010 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gt4f90io and ISPACK  2005/09/28 S.Takehiro
! 
! Solving a non-linear 2-D shallow water system on a sphere. 
!     
!     d\zeta/dt = -\Dinv{a(1-\mu^2)}\DP{[(f+\zeta)U)]}{\lamdba} 
!                 -\Dinv{a}\DP{[(f+\zeta)V)]}{\mu} - F_\zeta^{Diff},
!
!     dD/dt = \Dinv{a(1-\mu^2)}\DP{[(f+\zeta)V)]}{\lamdba} 
!             -\Dinv{a}\DP{[(f+\zeta)U)]}{\mu}
!             -\Dlapla[g(h+h_s)+E] - F_D^{Diff},
!
!     dh/dt = -\Dinv{a(1-\mu^2)}\DP{(hU)]}{\lamdba} 
!             -\Dinv{a}\DP{(hV)}{\mu} - F_h^{Diff},
!
!     \Dlapla\psi = \zeta, \Dlapla\Xi = D,
!
!     U = -\frac{(1-\mu^2)}{a}\DP{\psi}{\mu} + \Dinv{a}\DP{\Xi}{\lambda}
!     V = \Dinv{a}\DP{\psi}{\lambda} + \frac{(1-\mu^2)}{a}\DP{\Xi}{\mu}
!
! where hs is hight of topography, E is Kinetic energy, (U^2+V^2)/2(1-\mu^2), 
! F_*^{Diff} are hyper viscosity
!
! The time integration is performed by using 4th-order Runge-Kutta scheme
! except for the dissipation terms (Crank-Nicolson scheme). 
!
! The program is organized for freedecaying turbulence problem. 
! 
! History: 2005/10/09 S. Takehiro, created
!          2005/10/10 S. Takehiro, longname of div fixed.
!          2005/10/31 S. Nishizawa, define work variables for subroutines as global variable
!          2005/11/02 S. Nishizawa, change arguments of subroutine
!          2008/08/10 S. Takehiro, index changed, im -> 0:im-1
!          2010/02/07 S. Takehiro, initial perturbation can bu given 
!                                  as hsfc, div as well as vor.
!
program spshallow_zd_rn4cn_freedecay

!== モジュール引用宣言 ================================================

  use w_module 
  use gt4_history, only : HistoryCreate, HistoryPut, HistoryGet, &
                          HistoryAddVariable, HistoryAddAttr,    &
                          HistoryClose, GT_HISTORY
  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_Hsfc(:,:)          ! 表面変位擾乱 (t)
  real(8), allocatable :: xy_Vor(:,:)           ! 渦度 (t)
  real(8), allocatable :: xy_Div(:,:)           ! 発散 (t)

  real(8), allocatable :: xy_Coli(:,:)          ! 惑星渦度
  real(8), allocatable :: xy_Htopo(:,:)         ! 表面地形

  real(8), allocatable :: xy_VelLonCosLat(:,:)  ! U = xy_VelLon*Cos(xy_Lat)
  real(8), allocatable :: xy_VelLatCosLat(:,:)  ! V = xy_VelLat*Cos(xy_Lat)

  real(8), allocatable :: xy_VelLon(:,:)        ! 速度経度成分(出力用変数)
  real(8), allocatable :: xy_VelLat(:,:)        ! 速度緯度成分(出力用変数)
  real(8), allocatable :: xy_StrFunc(:,:)       ! 流線関数
  real(8), allocatable :: xy_VelPot(:,:)        ! 速度ポテンシャル

  real(8), allocatable :: xy_FluxVorLon(:,:)    ! (ζ+f)U (副プログラム作業用)
  real(8), allocatable :: xy_FluxVorLat(:,:)    ! (ζ+f)V (副プログラム作業用)

 !---- 変数(スペクトルデータ) ----
  real(8), allocatable :: w_Hsfc(:)             ! 表面変位擾乱 (t)
  real(8), allocatable :: w_Vor(:)              ! 渦度 (t)
  real(8), allocatable :: w_Div(:)              ! 発散 (t)

  real(8), allocatable :: w_HsfcTmp(:)          ! 表面変位擾乱 (Runge-Kutta 用)
  real(8), allocatable :: w_VorTmp(:)           ! 渦度 (Runge-Kutta 用)
  real(8), allocatable :: w_DivTmp(:)           ! 発散 (Runge-Kutta 用)

  real(8), allocatable :: w_StrFunc(:)          ! 流線関数
  real(8), allocatable :: w_VelPot(:)           ! 速度ポテンシャル

  real(8), allocatable :: w_Htopo(:)            ! 表面地形

  real(8), allocatable :: w_KEnergy(:)          ! 運動エネルギー

  real(8), allocatable :: w_DHsfcDtK1(:)        ! 変位時間変化(RK 第1段)
  real(8), allocatable :: w_DHsfcDtK2(:)        ! 変位時間変化(RK 第2段)
  real(8), allocatable :: w_DHsfcDtK3(:)        ! 変位時間変化(RK 第3段)
  real(8), allocatable :: w_DHsfcDtK4(:)        ! 変位時間変化(RK 第4段)

  real(8), allocatable :: w_DVorDtK1(:)         ! 渦度時間変化(RK 第1段)
  real(8), allocatable :: w_DVorDtK2(:)         ! 渦度時間変化(RK 第2段)
  real(8), allocatable :: w_DVorDtK3(:)         ! 渦度時間変化(RK 第3段)
  real(8), allocatable :: w_DVorDtK4(:)         ! 渦度時間変化(RK 第4段)

  real(8), allocatable :: w_DDivDtK1(:)         ! 発散時間変化(RK 第1段)
  real(8), allocatable :: w_DDivDtK2(:)         ! 発散時間変化(RK 第2段)
  real(8), allocatable :: w_DDivDtK3(:)         ! 発散時間変化(RK 第3段)
  real(8), allocatable :: w_DDivDtK4(:)         ! 発散時間変化(RK 第4段)

  real(8), allocatable :: nm_EvSp(:,:)          ! 渦エネルギースペクトル
  real(8), allocatable :: nm_EdSp(:,:)          ! 発散エネルギースペクトル
  real(8), allocatable :: nm_EkSp(:,:)          ! 運動エネルギースペクトル
  real(8), allocatable :: nm_EpSp(:,:)          ! 位置エネルギースペクトル
  real(8), allocatable :: nm_ESp(:,:)           ! 全エネルギースペクトル
  real(8), allocatable :: nm_EnsSp(:,:)         ! エンストロフィースペクトル

  real(8), allocatable :: n_EvSp(:)             ! 渦エネルギースペクトル
  real(8), allocatable :: n_EdSp(:)             ! 発散エネルギースペクトル
  real(8), allocatable :: n_EkSp(:)             ! 運動エネルギースペクトル
  real(8), allocatable :: n_EpSp(:)             ! 位置エネルギースペクトル
  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       ! 欠損値

 !---- 作業変数 ----
  real(8), allocatable :: w_HVisc(:)         ! 超粘性係数(運動方程式)
  real(8), allocatable :: w_HDiff(:)         ! 超拡散係数(質量保存式)

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

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

  type(ARGS) :: arg                          ! コマンド引数解析用
  character(STRING), pointer :: argv(:) => null()
  character(STRING)  :: nmlfile              ! NAMELIST 設定ファイル名

 !---- NAMELIST 変数 ----

  ! (メッセージ出力制御)
  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   = 6.37122D6     ! 惑星半径
  real(8)            :: Omega    = 7.292D-5      ! 回転角速度
  real(8)            :: Alpha    = 0.0           ! 回転軸方向(極との角度,deg.)
  real(8)            :: Grav     = 9.80616D0     ! 重力加速度
  real(8)            :: Hbar     = 4.0D3         ! 平均水深
  integer            :: HVOrder  = 2             ! 超粘性の次数
                                                 ! (水平ラプラシアンの階数)
  real(8)            :: HVisc    = 0.5D16        ! 超粘性係数
  integer            :: HDOrder  = 2             ! 超拡散の次数
                                                 ! (水平ラプラシアンの階数)
  real(8)            :: HDiff    = 0.5D16        ! 超拡散係数
  namelist /physics/ Radius, Omega, Alpha, Grav, Hbar, &
                     HVOrder, HVisc, HDOrder, HDiff

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

  ! (実験パラメター, 内部設定初期値)
  integer    :: seed=0             ! seed(1)に設定する種の値
  integer    :: nzero=10           ! 初期エネルギースペクトル分布のパラメタ
  integer    :: gamma=100          ! 初期エネルギースペクトル分布のパラメタ
  real(8)    :: Etotal=1.0D0       ! 初期全エネルギーの値
  character(len=10) :: disttype='YIHY2002' ! 初期エネルギー分布のタイプ
                                           ! (YIHY2002/YY1993)
  character(len=4)  :: fieldtype='Vor'     ! 初期場の種類(VOR/DIV/HSFC)
  namelist /initvalue/ seed, nzero, gamma, Etotal, disttype, fieldtype

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

  ! (ヒストリー出力)
  character(len=100) :: hst_file=   ''           ! ヒストリーファイル名
  character(len=100) :: source = &               ! ソースファイル名
       'spshallow_zd_rn4cn_freedecay.f90 (2010/02/07)'
  character(len=100) :: title = &                ! タイトル
       '2-dim shallow water 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 

!== メインルーチン ==========================================================

 !---------------- コマンドライン解析 ---------------------
  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: spshallow_zd_rn4cn_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)

 !---------------- 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=initvalue)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=initvalue) ; write(6,nml=initvalue) ; 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_Vor(0:im-1,jm),xy_Div(0:im-1,jm))
  allocate(xy_Hsfc(0:im-1,jm),xy_Htopo(0:im-1,jm))

  allocate(w_Vor((nm+1)*(nm+1)),w_Div((nm+1)*(nm+1)))
  allocate(w_Hsfc((nm+1)*(nm+1)),w_Htopo((nm+1)*(nm+1)))
  allocate(w_StrFunc((nm+1)*(nm+1)),w_VelPot((nm+1)*(nm+1)))

  allocate(w_VorTmp((nm+1)*(nm+1)),w_DivTmp((nm+1)*(nm+1)))
  allocate(w_HsfcTmp((nm+1)*(nm+1)))
  allocate(w_KEnergy((nm+1)*(nm+1)))

  allocate(xy_VelLonCosLat(0:im-1,jm),xy_VelLatCosLat(0:im-1,jm))

  allocate(w_DVorDtK1((nm+1)*(nm+1)),w_DVorDtK2((nm+1)*(nm+1)),  &
           w_DVorDtK3((nm+1)*(nm+1)),w_DVorDtK4((nm+1)*(nm+1)) )
  allocate(w_DDivDtK1((nm+1)*(nm+1)),w_DDivDtK2((nm+1)*(nm+1)),  &
           w_DDivDtK3((nm+1)*(nm+1)),w_DDivDtK4((nm+1)*(nm+1)) )
  allocate(w_DHsfcDtK1((nm+1)*(nm+1)),w_DHsfcDtK2((nm+1)*(nm+1)), &
           w_DHsfcDtK3((nm+1)*(nm+1)),w_DHsfcDtK4((nm+1)*(nm+1)) )

  allocate(xy_Coli(0:im-1,jm))
  allocate(w_HVisc((nm+1)*(nm+1)),w_HDiff((nm+1)*(nm+1)))

  allocate(nm_EvSp(0:nm,-nm:nm),nm_EdSp(0:nm,-nm:nm))
  allocate(nm_EkSp(0:nm,-nm:nm),nm_EpSp(0:nm,-nm:nm))
  allocate(nm_ESp(0:nm,-nm:nm),nm_EnsSp(0:nm,-nm:nm))

  allocate(n_EvSp(0:nm),n_EdSp(0:nm),n_EkSp(0:nm),n_EpSp(0:nm))
  allocate(n_ESp(0:nm),n_EnsSp(0:nm))

  allocate(xy_VelLon(0:im-1,jm),xy_VelLat(0:im-1,jm))
  allocate(xy_StrFunc(0:im-1,jm),xy_VelPot(0:im-1,jm))
  
  allocate(xy_FluxVorLon(0:im-1,jm),xy_FluxVorLat(0:im-1,jm))

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

 !------------------ 物理係数の設定 -----------------------
  xy_Coli = 2 * Omega * ( -cos(xy_Lon)*cos(xy_Lat)*sin(Alpha*pi/180.0) &
                         + sin(xy_Lat)*cos(Alpha*pi/180.0) )

  w_HVisc = HVisc &
          *( (-rn(:,1)/Radius**2)**HVOrder &
              -(2.0D0/Radius**2)**HVOrder )

  w_HDiff = HDiff * (-rn(:,1)/Radius**2)**HDOrder

  ! rn(ln(0,0,1) は正の値なので修正しておく.
  w_HVisc(l_nm(0,0)) = 0.0D0

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

  if ( initial_file == "") then
     ! リスタートファイルを指定しない場合, 
     ! 内部で w_Vor, w_Div, w_Hsfc, w_Htopo を与える. 
     if ( fieldtype(1:3) == 'VOR' ) then
        call set_initial_values_vor
     else if ( fieldtype(1:3) == 'DIV' ) then
        call set_initial_values_div
     else if ( fieldtype(1:4) == 'HSFC' ) then
        call set_initial_values_hsfc
     else
        call MessageNotify('E','main','fieldtype in NML initivaule invalid')
     endif
  else
     !  初期値設定(リスタートファイルからの読みこみ)
     call HistoryGet( trim(initial_file), 'w_vor',   w_Vor,   time )
     call HistoryGet( trim(initial_file), 'w_div',   w_div,   time )
     call HistoryGet( trim(initial_file), 'w_hsfc',  w_Hsfc,  time )
     call HistoryGet( trim(initial_file), 'w_htopo', w_Htopo )
  endif

 ! 渦度・発散から流線関数・速度ポテンシャルならびに速度成分を計算
  call VorDiv2Velocity(                  &
       w_Vor, w_Div         )              !(in)  渦度, 発散

  xy_Vor = xy_w(w_Vor) ; xy_Div = xy_w(w_Div) ; xy_HSfc = xy_w(w_Hsfc)

 !------------- 時間積分(Runge-Kutta + Crank-Nicolson 法) --------------
  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

     !---- 1 段目 [ k1 = f( x_n, t_n) ] ----
     call TimeDerivetives_noDisppation(                 &
       w_Vor, w_Div, w_Hsfc,                            &!(in) 渦度, 発散, 変位
       w_DVorDtK1, w_DDivDtK1, w_DHsfcDtK1)              !(out) 時間変化

     !---- 2 段目 [ k2 = f( x_n+k1*dt/2, t_n+dt/2 ) ] ----
     w_VorTmp  = w_Vor  + w_DVorDtK1  * delta_t/2.0D0
     w_DivTmp  = w_Div  + w_DDivDtK1  * delta_t/2.0D0
     w_HsfcTmp = w_Hsfc + w_DHsfcDtK1 * delta_t/2.0D0

     xy_Vor = xy_w(w_VorTmp)
     xy_Div = xy_w(w_DivTmp)
     xy_HSfc = xy_w(w_HsfcTmp)

     call TimeDerivetives_noDisppation(                 &
       w_VorTmp, w_DivTmp, w_HsfcTmp,                   &!(in) 渦度, 発散, 変位
       w_DVorDtK2, w_DDivDtK2, w_DHsfcDtK2)              !(out) 時間変化

     !---- 3 段目 [ k3 = f( x_n+k2*dt/2, t_n+dt/2 ) ] ----
     w_VorTmp  = w_Vor  + w_DVorDtK2  * delta_t/2.0D0
     w_DivTmp  = w_Div  + w_DDivDtK2  * delta_t/2.0D0
     w_HsfcTmp = w_Hsfc + w_DHsfcDtK2 * delta_t/2.0D0

     xy_Vor = xy_w(w_VorTmp)
     xy_Div = xy_w(w_DivTmp)
     xy_HSfc = xy_w(w_HsfcTmp)

     call TimeDerivetives_noDisppation(                 &
       w_VorTmp, w_DivTmp, w_HsfcTmp,                   &!(in) 渦度, 発散, 変位
       w_DVorDtK3, w_DDivDtK3, w_DHsfcDtK3)              !(out) 時間変化

     !---- 4 段目 [ k4 = f( x_n+k3*dt, t_n+dt ) ] ----
     w_VorTmp  = w_Vor  + w_DVorDtK3  * delta_t
     w_DivTmp  = w_Div  + w_DDivDtK3  * delta_t
     w_HsfcTmp = w_Hsfc + w_DHsfcDtK3 * delta_t

     xy_Vor = xy_w(w_VorTmp)
     xy_Div = xy_w(w_DivTmp)
     xy_HSfc = xy_w(w_HsfcTmp)

     call TimeDerivetives_noDisppation(                 &
       w_VorTmp, w_DivTmp, w_HsfcTmp,                   &!(in) 渦度, 発散, 変位
       w_DVorDtK4, w_DDivDtK4, w_DHsfcDtK4)              !(out) 時間変化

     !---- 積分(Runge-Kutta) ----
     w_Vor  = w_Vor  + delta_t * (   w_DVorDtK1/6.0D0 + w_DVorDtK2/3.0D0 &
                                   + w_DVorDtK3/3.0D0 + w_DVorDtK4/6.0D0   )
     w_Div  = w_Div  + delta_t * (  w_DDivDtK1/6.0D0 + w_DDivDtK2/3.0D0 &
                                   + w_DDivDtK3/3.0D0 + w_DDivDtK4/6.0D0   )
     w_Hsfc = w_Hsfc + delta_t * (   w_DHsfcDtK1/6.0D0 + w_DHsfcDtK2/3.0D0 &
                                   + w_DHsfcDtK3/3.0D0 + w_DHsfcDtK4/6.0D0  )

     !---- 散逸項積分(Crank-Nicolson) ----

     w_Vor  = (1-w_HVisc*delta_t/2.0D0)/(1+w_HVisc*delta_t/2.0D0) * w_Vor
     w_Div  = (1-w_HVisc*delta_t/2.0D0)/(1+w_HVisc*delta_t/2.0D0) * w_Div
     w_Hsfc = (1-w_HDiff*delta_t/2.0D0)/(1+w_HDiff*delta_t/2.0D0) * w_Hsfc

     xy_Vor = xy_w(w_Vor) ; xy_Div = xy_w(w_Div) ; xy_HSfc = xy_w(w_Hsfc)

     if(mod(it,hst_intstep) .eq. 0)then                    ! ヒストリー出力
       ! 渦度・発散から流線関数・速度ポテンシャルならびに速度成分を計算
        call VorDiv2Velocity(                  &
          w_Vor, w_Div         )              !(in)  渦度, 発散

        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_vor    ! w_Vor 擾乱を与える. 

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

    ! 乱数設定
    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_vor',&
            '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_vor',&
            '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_vor',&
            '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)

    ! 渦度場以外は 0
    w_Div = 0.0D0
    w_Hsfc = 0.0D0
    w_Htopo = 0.0D0

    deallocate(n_EspInit,seedarray)

  end subroutine set_initial_values_vor

 !
 !  初期値設定(リスタートファイルない場合のデフォルト設定, 発散擾乱)
 !
  subroutine set_initial_values_div    ! w_Div 擾乱を与える. 

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

    ! 乱数設定
    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_div',&
            '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_div',&
            '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_div',&
            'The parameter "disttype" should be "YIHY2002" or "YY1993"')
    endif

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

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

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

    ! エネルギーの大きさを定める. 
    w_VelPot = w_VelPot * sqrt(Etotal/sum(n_ESp))
    xy_VelPot = xy_w(w_VelPot)
    w_Div = w_Lapla_w(w_VelPot)/Radius**2
    xy_Div =  xy_w(w_Div)

    ! 渦度場以外は 0
    w_Vor = 0.0D0
    w_Hsfc = 0.0D0
    w_Htopo = 0.0D0

    deallocate(n_EspInit,seedarray)

  end subroutine set_initial_values_div

 !
 !  初期値設定(リスタートファイルない場合のデフォルト設定, 表面変位場擾乱)
 !
  subroutine set_initial_values_hsfc    ! w_Hsfc 擾乱を与える. 

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

    ! 乱数設定
    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_hsfc',&
            '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_hsfc',&
            '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_hsfc',&
            'The parameter "disttype" should be "YIHY2002" or "YY1993"')
    endif

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

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

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

    ! エネルギーの大きさを定める. 
    w_Hsfc = w_Hsfc * sqrt(Etotal/sum(n_ESp))

    ! 表面変位場以外は 0
    w_Vor = 0.0D0
    w_Div = 0.0D0
    w_Htopo = 0.0D0

    deallocate(n_EspInit,seedarray)

  end subroutine set_initial_values_hsfc

 !---------------- ポテンシャルエネルギースペクトル計算 ----------------
  function nm_EnergyFromSurfaceHeight_w(w_Hsfc)
    ! 
    ! 表面変位のスペクトルデータからポテンシャルエネルギーの球面調和函数成分
    ! (スペクトル)を計算する(1 層用).
    !
    !  * 全波数 n, 帯状波数 m の表面変位のスペクトル成分Hsfc(n,m) から
    !    エネルギースペクトルは (1/2)Grav Hsfc(n,m)^2 と計算される.
    !    Rd は変形半径である
    !
    !  * 全てのエネルギースペクトル成分の和に4πをかけたものが球面上での
    !    全エネルギーに等しい.
    !
    !  * データの存在しない全波数 n, 帯状波数 m の配列には欠損値が格納される.
    !    欠損値の値はモジュール変数 w_spectrum_VMiss によって設定できる
    !    (初期値は -999.0)
    !
    real(8), intent(in)   :: w_Hsfc(:)
    !(in) 表面変位(スペクトルデータ)

    real(8), dimension(0:nm,-nm:nm) :: nm_EnergyFromSurfaceHeight_w
    !(out) エネルギースペクトル(水平全波数 n, 帯状波数 m 空間)

    integer n,m                               ! DO 変数

    nm_EnergyFromSurfaceHeight_w = w_spectrum_VMiss

    do n=0,nm
       do m=-n,n
          nm_EnergyFromSurfaceHeight_w(n,m) &
               = 0.5 * Grav* w_Hsfc(l_nm(n,m))**2
       enddo
    enddo
  end function nm_EnergyFromSurfaceHeight_w

  function n_EnergyFromSurfaceHeight_w(w_Hsfc)
    !
    ! 表面変位のスペクトルデータから各全波数のポテンシャルエネルギー成分
    ! (スペクトル)を計算する(1 層用).
    !
    !  * 全波数 n の表面変位のスペクトル成分H(n,m) から
    !    エネルギースペクトルはΣ[m=-nm]^nm(1/2)Grav*H(n,m)^2
    !    と計算される. 
    !
    !  * 全てのエネルギースペクトル成分の和に 4πをかけたものが
    !    球面上での全エネルギーに等しい.
    !

    real(8), intent(in)      :: w_Hsfc(:)
    !(in) 表面変位(スペクトルデータ)

    real(8), dimension(0:nm) :: n_EnergyFromSurfaceHeight_w
    !(out) エネルギースペクトル (水平全波数 n 空間) 

    integer n,m                                 ! DO 変数

    do n=0,nm
       n_EnergyFromSurfaceHeight_w(n)  &
            = 0.5 * Grav * sum(w_Hsfc(l_nm(n,(/(m,m=-n,n)/)))**2,1) 
    enddo

  end function n_EnergyFromSurfaceHeight_w

!=========================== 時間変化項計算 ============================

  !
  ! 散逸項を除いた時間変化計算
  !
  subroutine TimeDerivetives_noDisppation( &
       w_Vor, w_Div, w_Hsfc,               &             !(in) 渦度, 発散, 変位
       w_DVorDtNoDisp, w_DDivDtNoDisp, w_DHsfcDtNoDisp)  !(out) 時間変化

   !---- 変数(入力) ----
    real(8), intent(in) :: w_Vor(:)  ! 渦度
    real(8), intent(in) :: w_Div(:)  ! 発散
    real(8), intent(in) :: w_Hsfc(:) ! 変位

   !---- 変数(出力) ----
    real(8), intent(out) :: w_DHsfcDtNoDisp(:) ! 変位時間変化
    real(8), intent(out) :: w_DVorDtNoDisp(:)   ! 渦度時間変化
    real(8), intent(out) :: w_DDivDtNoDisp(:)   ! 発散時間変化

    ! 渦度・発散から流線関数・速度ポテンシャルならびに速度成分を計算
    call VorDiv2Velocity(                &
       w_Vor, w_Div         )              !(in)  渦度, 発散

     xy_FluxVorLon = (xy_Vor + xy_Coli) * xy_VelLonCosLat
     xy_FluxVorLat = (xy_Vor + xy_Coli) * xy_VelLatCosLat

     w_KEnergy = w_xy((xy_VelLonCosLat**2 + xy_VelLatCosLat**2) &
                       /(2.0*cos(xy_Lat)**2))                

     ! 時間変化項の見積り
     w_DVorDtNoDisp = - w_DivLambda_xy(xy_FluxVorLon)/Radius &
                      - w_DivMu_xy(xy_FluxVorLat)/Radius

     w_DDivDtNoDisp =   w_DivLambda_xy(xy_FluxVorLat)/Radius            &
                    - w_DivMu_xy(xy_FluxVorLon)/Radius                  &
                    - w_Lapla_w( Grav*(w_Hsfc+ w_Htopo)+w_KEnergy )/Radius**2 

     w_DHsfcDtNoDisp = - w_DivLambda_xy(xy_Hsfc*xy_VelLonCosLat)/Radius &
                       - w_DivMu_xy(xy_Hsfc*xy_VelLatCosLat)/Radius &
                       - Hbar * w_Div
    
  end subroutine TimeDerivetives_noDisppation

 !
 ! 渦度・発散から流線関数・速度ポテンシャルならびに速度成分を計算
 !
  subroutine VorDiv2Velocity(            &
       w_Vor, w_Div         )              !(in)  渦度, 発散

   !---- 変数(入力) ----
    real(8), intent(in)  :: w_Vor(:)                ! 渦度
    real(8), intent(in)  :: w_Div(:)                ! 発散

    ! 渦度・発散から流線関数・速度ポテンシャルを計算
    w_Strfunc = w_LaplaInv_w(w_Vor) * Radius**2
    w_VelPot  = w_LaplaInv_w(w_Div) * Radius**2

    ! 流線関数・速度ポテンシャルから速度成分を計算
    xy_VelLonCosLat =  -xy_Gradmu_w(w_StrFunc)/Radius &
                      + xy_GradLambda_w(w_VelPot)/Radius 
    xy_VelLatCosLat =   xy_GradLambda_w(w_StrFunc)/Radius &
                      + xy_Gradmu_w(w_VelPot)/Radius 

  end subroutine VorDiv2Velocity


!=========================== リスタート出力 ============================
 !
 ! リスタート出力初期化
 !
  subroutine output_restart_init
    call HistoryCreate( &
           file=trim(rst_file), title=trim(title), source=trim(source), &
           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     ','sec   '/),   &
           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/sec', xtype='double', history=hst_rst)
    call HistoryAddVariable( &                                ! 変数定義
           varname='w_div', dims=(/'nm','t '/), & 
           longname='Divergence', &
           units='1/sec', xtype='double', history=hst_rst)
    call HistoryAddVariable( &                                ! 変数定義
           varname='w_hsfc', dims=(/'nm','t '/), & 
           longname='Suface height', &
           units='m', xtype='double', history=hst_rst)
    call HistoryAddVariable( &                                ! 変数定義
           varname='w_htopo', dims=(/'nm'/), & 
           longname='Bottom topography', &
           units='m', xtype='double', history=hst_rst)

   !---- 物理変数出力 ----
    call HistoryPut('w_htopo', w_Htopo, hst_rst)

   !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ----
    call HistoryAddAttr('lon','+Radius',  Radius,  hst_rst)
    call HistoryAddAttr('lon','+Omega',   Omega,   hst_rst)
    call HistoryAddAttr('lon','+Alpha',   Alpha,   hst_rst)
    call HistoryAddAttr('lon','+Grav',    Grav,    hst_rst)
    call HistoryAddAttr('lon','+Hbar',    Hbar,    hst_rst)
    call HistoryAddAttr('lon','+HVOrder', HVOrder, hst_rst)
    call HistoryAddAttr('lon','+HVisc',   HVisc,   hst_rst)
    call HistoryAddAttr('lon','+HDOrder', HDOrder, hst_rst)
    call HistoryAddAttr('lon','+HDiff',   HDiff,   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)
    call HistoryPut('w_div', w_Div, hst_rst)
    call HistoryPut('w_hsfc',w_Hsfc,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','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     ','sec   '/),   &
           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='vor', dims=(/'lon','lat','t  '/), & 
           longname='Vorticity', units='1/sec', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='div', dims=(/'lon','lat','t  '/), & 
           longname='Divergence', units='1/sec', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='hsfc', dims=(/'lon','lat','t  '/), & 
           longname='Surface height', units='m', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='strfunc', dims=(/'lon','lat','t  '/), & 
           longname='Stream function', units='m2/sec', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='velpot', dims=(/'lon','lat','t  '/), & 
           longname='Velocity potential', units='m2/sec', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='vellon', dims=(/'lon','lat','t  '/), & 
           longname='lon-velocity', units='m/sec', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='vellat', dims=(/'lon','lat','t  '/), & 
           longname='lat-velocity', units='m/sec', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='am', dims=(/'lat','t  '/), & 
           longname='Zonal mean angular momentum (U cos\phi)',  &
           units='m^2/sec)', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='qv', dims=(/'lon','lat','t  '/), & 
           longname='Potential vorticity', units='1/(m sec)', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='htopo', dims=(/'lon','lat'/), & 
           longname='Bottom topography', units='m', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='ek', dims=(/'t'/), & 
           longname='Kinetic energy', units='m^3/sec^2', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ep', dims=(/'t'/), & 
           longname='Potential energy', units='m^3/sec^2', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='en', dims=(/'t'/), & 
           longname='Total energy', units='m^3/sec^2', xtype='double')

   !---- エネルギースペクトル診断変数定義 ----
    call HistoryAddVariable( &                                ! 変数定義
           varname='evspnm', dims=(/'n','m','t'/), & 
           longname='Vortical kinetic energy spectrum (n-m space)', &
           units='1', xtype='double')
    call HistoryAddAttr('evspnm','missing_value', vmiss  )
    call HistoryAddVariable( &                                ! 変数定義
           varname='evspn', dims=(/'n','t'/), & 
           longname='Vortical kinetic energy spectrum', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='evsp', dims=(/'t'/), & 
           longname='Mean vortical kinetic energy (by spectrum)', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='evsptot', dims=(/'t'/), & 
           longname='Total vortical kinetic energy (by spectrum)', &
           units='1', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='edspnm', dims=(/'n','m','t'/), & 
           longname='Divergent kinetic energy spectrum (n-m space)', &
           units='1', xtype='double')
    call HistoryAddAttr('edspnm','missing_value', vmiss  )
    call HistoryAddVariable( &                                ! 変数定義
           varname='edspn', dims=(/'n','t'/), & 
           longname='Divergent kinetic energy spectrum', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='edsp', dims=(/'t'/), & 
           longname='Mean divergent kinetic energy (by spectrum)', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='edsptot', dims=(/'t'/), & 
           longname='Total divergent kinetic energy (by spectrum)', &
           units='1', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='ekspnm', dims=(/'n','m','t'/), & 
           longname='Kinetic energy spectrum (n-m space)', &
           units='1', xtype='double')
    call HistoryAddAttr('ekspnm','missing_value', vmiss  )
    call HistoryAddVariable( &                                ! 変数定義
           varname='ekspn', dims=(/'n','t'/), & 
           longname='Kinetic energy spectrum', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='eksp', dims=(/'t'/), & 
           longname='Mean kinetic energy (by spectrum)', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='eksptot', dims=(/'t'/), & 
           longname='Total kinetic energy (by spectrum)', &
           units='1', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='epspnm', dims=(/'n','m','t'/), & 
           longname='potential energy spectrum (n-m space)', &
           units='1', xtype='double')
    call HistoryAddAttr('epspnm','missing_value', vmiss  )
    call HistoryAddVariable( &                                ! 変数定義
           varname='epspn', dims=(/'n','t'/), & 
           longname='potential energy spectrum', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='epsp', dims=(/'t'/), & 
           longname='Mean potential energy (by spectrum)', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='epsptot', dims=(/'t'/), & 
           longname='Total potential energy (by spectrum)', &
           units='1', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='enspnm', dims=(/'n','m','t'/), & 
           longname='Total energy spectrum (n-m space)', &
           units='1', xtype='double')
    call HistoryAddAttr('enspnm','missing_value', vmiss  )

    call HistoryAddVariable( &                                ! 変数定義
           varname='enspn', dims=(/'n','t'/), & 
           longname='Total energy spectrum', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ensp', dims=(/'t'/), & 
           longname='Mean total energy (by spectrum)', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='ensptot', dims=(/'t'/), & 
           longname='Total energy (by spectrum)', &
           units='1', xtype='double')

   !---- 物理変数出力 ----
    call HistoryPut('htopo', xy_Htopo)

   !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ----
    call HistoryAddAttr('lon','+Radius',  Radius )
    call HistoryAddAttr('lon','+Omega',   Omega  )
    call HistoryAddAttr('lon','+Alpha',   Alpha  )
    call HistoryAddAttr('lon','+Grav',    Grav   )
    call HistoryAddAttr('lon','+Hbar',    Hbar   )
    call HistoryAddAttr('lon','+HVOrder', HVOrder)
    call HistoryAddAttr('lon','+HVisc',   HVisc  )
    call HistoryAddAttr('lon','+HDOrder', HDOrder)
    call HistoryAddAttr('lon','+HDiff',   HDiff  )
    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))

   !---- 物理変数出力 ----
    call HistoryPut('vor', xy_Vor)
    call HistoryPut('div', xy_Div)
    call HistoryPut('hsfc',xy_Hsfc)

    call HistoryPut('strfunc', xy_w(w_StrFunc))
    call HistoryPut('velpot',  xy_w(w_VelPot))

    xy_VelLon = xy_VelLonCosLat/cos(xy_Lat)
    xy_VelLat = xy_VelLatCosLat/cos(xy_Lat)
    call HistoryPut('vellon',xy_VelLon)
    call HistoryPut('vellat',xy_VelLat)

    call HistoryPut('am', y_AvrLon_xy(xy_VelLonCosLat)*Radius)
    call HistoryPut('qv', (xy_Coli + xy_Vor)/(Hbar + xy_Hsfc))

    call HistoryPut('ek',0.5*AvrLonLat_xy((Hbar+xy_hsfc) &
                              *(xy_Vellon**2+xy_Vellat**2)))
    call HistoryPut('ep',0.5*Grav*AvrLonLat_xy(xy_Hsfc**2))

    call HistoryPut('en',0.5*AvrLonLat_xy(&
            (Hbar+xy_hsfc)*(xy_Vellon**2+xy_Vellat**2)+Grav*xy_Hsfc**2))

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

    call HistoryPut('evspnm',nm_EvSp/Radius**2)
    call HistoryPut('evspn',n_EvSp/Radius**2)
    call HistoryPut('evsptot',(4*pi)*sum(n_Evsp))
    call HistoryPut('evsp',sum(n_Evsp)/Radius**2)

    nm_EdSp = nm_EnergyFromStreamfunc_w(w_VelPot)
    n_EdSp = n_EnergyFromStreamfunc_w(w_VelPot)

    call HistoryPut('edspnm',nm_EdSp/Radius**2)
    call HistoryPut('edspn',n_EdSp/Radius**2)
    call HistoryPut('edsptot',(4*pi)*sum(n_Edsp))
    call HistoryPut('edsp',sum(n_Edsp)/Radius**2)

    nm_EkSp = w_spectrum_VMiss
    do n=0,nm
       do m=-n,n
          nm_EkSp(n,m) = nm_EvSp(n,m) + nm_EdSp(n,m)
       enddo
    enddo
    n_EkSp = n_EvSp + n_EdSp

    call HistoryPut('ekspnm',nm_EkSp/Radius**2)
    call HistoryPut('ekspn',n_EkSp/Radius**2)
    call HistoryPut('eksptot',(4*pi)*sum(n_EkSp))
    call HistoryPut('eksp',sum(n_EkSp)/Radius**2)

    nm_EpSp = nm_EnergyFromSurfaceHeight_w(w_Hsfc)
    n_EpSp = n_EnergyFromSurfaceHeight_w(w_Hsfc)

    call HistoryPut('epspnm',nm_EpSp/Radius**2)
    call HistoryPut('epspn',n_EpSp/Radius**2)
    call HistoryPut('epsptot',(4*pi)*sum(n_Epsp))
    call HistoryPut('epsp',sum(n_Epsp)/Radius**2)

    nm_ESp = w_spectrum_VMiss
    do n=0,nm
       do m=-n,n
          nm_ESp(n,m) = nm_EkSp(n,m) + nm_EpSp(n,m)
       enddo
    enddo
    n_ESp = n_EkSp + n_EpSp

    call HistoryPut('enspnm',nm_ESp/Radius**2)
    call HistoryPut('enspn',n_ESp/Radius**2)
    call HistoryPut('ensptot',(4*pi)*sum(n_ESp))
    call HistoryPut('ensp',sum(n_ESp)/Radius**2)

  end subroutine output_history

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

end program spshallow_zd_rn4cn_freedecay
