!----------------------------------------------------------------------
!  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.
!          2020/04/13 R. Suzuki  mpi用の編集を行っている.
!          2020/04/15 R. Suzuki  case1用に初期値を書き換えた
program spshallow_zd_rn4cn_freedecay

!== モジュール引用宣言 ================================================
  use mpi
  use w_mpi_module ! add rsuzuki
  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, Printf ! Printf add rsuzuki
  use dc_args
  implicit none

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

 !---- 変数(格子点データ)---
  real(8), allocatable :: xv_Hsfc(:,:)          ! 表面変位擾乱 (t)
  real(8), allocatable :: xv_Vor(:,:)           ! 渦度 (t)
  real(8), allocatable :: xv_Div(:,:)           ! 発散 (t)

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

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

  real(8), allocatable :: xv_VelLon(:,:)        ! 速度経度成分(出力用変数)
  real(8), allocatable :: xv_VelLat(:,:)        ! 速度緯度成分(出力用変数)
  real(8), allocatable :: xv_StrFunc(:,:)       ! 流線関数
  real(8), allocatable :: xv_VelPot(:,:)        ! 速度ポテンシャル

  real(8), allocatable :: xv_FluxVorLon(:,:)    ! (ζ+f)U (副プログラム作業用)
  real(8), allocatable :: xv_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(:)            ! エンストロフィースペクトル

! add rsuzuki case1 解析用
  real(8) :: VelLonRef
  real(8) :: VelLatRef
  real(8) :: Lon_ref
  real(8) :: Lat_ref
  real(8), allocatable :: xv_HsfcRef(:,:)       ! 表面変位 解析解
  real(8) :: U0    = 3.861068D1                 ! 剛体回転流の速度(m/sec)
  real(8) :: H0    = 1.0D3                      ! 表面変位の振幅(m)
  real(8) :: Rh    = 2.1040667D6                ! 表面変位の半径(m)
  real(8) :: LonC  = 270.0                      ! 表面変位の位置(経度,deg.)
  real(8) :: LatC  = 0.0                        ! 表面変位の位置(緯度,deg.)
  real(8), allocatable :: r(:,:)                ! LatC, LonC からの距離
  real(8), allocatable :: xv_HsfcError(:,:)     ! 解析解との表面変位の誤差
  real(8) :: HsfcmaxDef, HsfcmaxRef             ! 全プロセスの中から表面変位の解析解との差の
                                                ! 最大値と解析解の最大値を取得するための変数

 !---- 固定パラメタ -----
  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                            ! 全波数, 帯状波数
  integer :: i, j                            ! 経度, 緯度格子インデックス add rsuzuki

  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  !

  ! (時間積分)
  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

  ! mpi 用変数 add
  integer :: np, ip, ierr ! MPI用変数
  character(len=STRING) :: cip !ip 文字変数

  ! 計算時間計測用 add
  real(8) :: TIM0, TIM1 ! start, end time

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

 ! MPI スタート ! add
  call MPI_INIT(IERR)
  call MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)

 !---------------- コマンドライン解析 ---------------------
  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 読み込み ---------------------

  if( ip==0 ) write(6,nml=message)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=message) ; if( ip==0 ) write(6,nml=message) ; close(10)

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

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

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

  call Printf(unit=cip,fmt="%04d",i=(/ip/))
  if ( initial_file /= '' ) then
     initial_file=trim(initial_file)//'-ip'//trim(adjustl(cip))//'.nc'
  endif

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

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

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

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

 ! MPI用座標値設定 add
  call DbgMessage(fmt='call %c', c1='w_mpi_initial')
  call w_mpi_Initial(nm,im,jm) !nm:全波数 im:東西格子点数 jm:南北格子点数 np_in:ダミー変数?

 !------------------ 変数の割り付け ---------------------
  allocate(xv_Vor(0:im-1,jc),xv_Div(0:im-1,jc))
  allocate(xv_Hsfc(0:im-1,jc),xv_Htopo(0:im-1,jc))
  allocate(xv_HsfcRef(0:im-1,jc),r(0:im-1,jc))
  allocate(xv_HsfcError(0:im-1,jc))

  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(xv_VelLonCosLat(0:im-1,jc),xv_VelLatCosLat(0:im-1,jc))

  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(xv_Coli(0:im-1,jc))
  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(xv_VelLon(0:im-1,jc),xv_VelLat(0:im-1,jc))
  allocate(xv_StrFunc(0:im-1,jc),xv_VelPot(0:im-1,jc))

  allocate(xv_FluxVorLon(0:im-1,jc),xv_FluxVorLat(0:im-1,jc))

  !------ 計算開始時間の取得 ------ add rsuzuki
   call APTIME(TIM0)

  w_spectrum_VMiss = vmiss

 !------------------ 物理係数の設定 -----------------------
  xv_Coli = 2 * Omega * ( -cos(xv_Lon)*cos(xv_Lat)*sin(Alpha*pi/180.0) &
                         + sin(xv_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 を与える.
     call set_initial_values

  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)  渦度, 発散

  xv_Vor = xv_w(w_Vor) ; xv_Div = xv_w(w_Div) ; xv_HSfc = xv_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

! case1 解析用
     Lon_ref = Lon_ref + delta_t*VelLonRef
     Lat_ref = Lat_ref + delta_t*VelLatRef

     r = Radius * acos(sin(Lat_ref)*sin(xv_Lat) &
        + cos(Lat_ref)*cos(xv_Lat)*cos(xv_Lon-Lon_ref))

     where ( r .le. Rh )
        xv_HsfcRef = H0/2.0*(1.0+cos(pi*r/Rh))
     elsewhere
        xv_HsfcRef = 0.0
     end where

! case1 では連続の式のみを解く

     !---- 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

     xv_Vor = xv_w(w_VorTmp) !
     xv_Div = xv_w(w_DivTmp) !
     xv_HSfc = xv_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

     xv_Vor = xv_w(w_VorTmp) !
     xv_Div = xv_w(w_DivTmp) !
     xv_HSfc = xv_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

     xv_Vor = xv_w(w_VorTmp) !
     xv_Div = xv_w(w_DivTmp) !
     xv_HSfc = xv_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

     xv_Vor = xv_w(w_Vor) ; xv_Div = xv_w(w_Div) ; xv_HSfc = xv_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

!------ 計算終了時間の取得 ------ add rsuzuki
  call APTIME(TIM1)
  print *, '計算時間[sec] : ', TIM1-TIM0

  call output_restart_close
  call output_history_close

  ! MPI 終了
  call w_mpi_Finalize()
  call MPI_FINALIZE(IERR)

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

contains

!=========================== 初期値設定 ============================
 !
 !  初期値設定(リスタートファイルない場合のデフォルト設定,case1)
 !
  subroutine set_initial_values    ! w_Vor, w_Div, w_Hsfc, w_Htopo を与える.

    real(8) :: U0    = 3.861068D1           ! 剛体回転流の速度(m/sec)
    real(8) :: H0    = 1.0D3                ! 表面変位の振幅(m)
    real(8) :: Rh    = 2.1040667D6          ! 表面変位の半径(m)
    real(8) :: LonC  = 270.0                ! 表面変位の位置(経度,deg.)
    real(8) :: LatC  = 0.0                  ! 表面変位の位置(緯度,deg.)

    namelist /case1/  U0, H0, Rh, LonC, LatC

    real(8) :: r                            ! LatC, LonC からの距離

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

    w_Vor = w_Lapla_w(w_xy( &
                 - Radius * U0 * &
                   (sin(xy_Lat)*cos(Alpha*pi/180.0) &
                 - cos(xy_Lon)*cos(xy_Lat)*sin(Alpha*pi/180.0)) &
             ))/Radius**2

    w_Div= 0.0D0

    do j=1,jc
       do i=0,im-1
          r = Radius * acos(sin(LatC*pi/180.0)*sin(xv_Lat(i,j)) &
               + cos(LatC*pi/180.0)*cos(xv_Lat(i,j)) &
                    *cos(xv_Lon(i,j)-LonC*pi/180.0))
          if ( r .le. Rh ) then
             xv_Hsfc(i,j) = H0/2.0*(1.0+cos(pi*r/Rh))
          else
             xv_Hsfc(i,j) = 0.0
          end if
       end do
    end do

    w_Hsfc = w_xv(xv_Hsfc)
    w_Htopo = 0.0

! 解析計算用

    Lon_ref = LonC*pi/180.0
    Lat_ref = LatC*pi/180.0
    xv_HsfcRef = xv_Hsfc

    VelLonRef = (U0/Radius)*(cos(Lat_ref)*cos(Alpha*pi/180.0)+ &
                 sin(Lat_ref)*cos(Lon_ref)*sin(Alpha*pi/180.0)) !(62)
    VelLatRef = -(U0/Radius)*sin(Lon_ref)*sin(Alpha*pi/180.0)

  end subroutine set_initial_values

 !---------------- ポテンシャルエネルギースペクトル計算 ----------------
  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)  渦度, 発散

     xv_FluxVorLon = (xv_Vor + xv_Coli) * xv_VelLonCosLat !
     xv_FluxVorLat = (xv_Vor + xv_Coli) * xv_VelLatCosLat !

     w_KEnergy = w_xv((xv_VelLonCosLat**2 + xv_VelLatCosLat**2) & !
                       /(2.0*cos(xv_Lat)**2))

     ! 時間変化項の見積り
     w_DVorDtNoDisp = - w_DivLambda_xv(xv_FluxVorLon)/Radius & !
                      - w_DivMu_xv(xv_FluxVorLat)/Radius !

     w_DDivDtNoDisp =   w_DivLambda_xv(xv_FluxVorLat)/Radius            & !
                    - w_DivMu_xv(xv_FluxVorLon)/Radius                  & !
                    - w_Lapla_w( Grav*(w_Hsfc+ w_Htopo)+w_KEnergy )/Radius**2

     w_DHsfcDtNoDisp = - w_DivLambda_xv(xv_Hsfc*xv_VelLonCosLat)/Radius & !
                       - w_DivMu_xv(xv_Hsfc*xv_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

    ! 流線関数・速度ポテンシャルから速度成分を計算
    xv_VelLonCosLat =  -xv_Gradmu_w(w_StrFunc)/Radius & !
                      + xv_GradLambda_w(w_VelPot)/Radius  !
    xv_VelLatCosLat =   xv_GradLambda_w(w_StrFunc)/Radius & !
                      + xv_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,jc,(nm+1)**2,0/),& ! jm => jc
           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',v_Lat, hst_rst)  ! y_Lat => v_Lat   ! 変数出力
    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',v_Lat_weight,hst_rst) ! 変数出力
                                       !y_Lat_weight => v_Lat_weight

   !---- 物理変数定義 ----
    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,&
               '  ip = ', ip ! add ip を追加
    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,jc,(nm+1)**2,nm+1,2*nm+1,0/),& ! jm => jc
           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',v_Lat/pi*180)     ! y_Lat => v_Lat  ! 変数出力
                        ! y_Lat 分割ファイル両方同じ緯度
                        ! v_Lat 赤道をまたいで分割される

    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',v_Lat_weight)         ! 変数出力
                                       !y_Lat_weight => v_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')

! 誤差診断量 add
    call HistoryAddVariable( &                                ! 変数定義
           varname='hsfcRef', dims=(/'lon','lat','t  '/), &
           longname='Surface height(Analytical)', units='m', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='hsfc_error', dims=(/'lon','lat','t  '/), &
           longname='Surface height(error)', units='m', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='L1norm', dims=(/'t'/), &
           longname='L1norm', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='L2norm', dims=(/'t'/), &
           longname='L2norm', units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='Linfnorm', dims=(/'t'/), &
           longname='Linfnorm', units='1', 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', xv_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,&
               ' ip = ', ip ! add ipを追加
    call HistoryPut('t',real(time))

   !---- 物理変数出力 ----
    call HistoryPut('vor', xv_Vor) !  (lon,lat,t)
    call HistoryPut('div', xv_Div) !  (lon,lat,t)
    call HistoryPut('hsfc',xv_Hsfc) !  (lon,lat,t)

    call HistoryPut('hsfcRef',xv_HsfcRef) !  (lon,lat,t)
    xv_HsfcError = xv_Hsfc-xv_HsfcRef
    call HistoryPut('hsfc_error',xv_HsfcError)
    call HistoryPut('L1norm',AvrLonLat_xv(abs(xv_HsfcError))/ &
                             AvrLonLat_xv(abs(xv_HsfcRef)))
    call HistoryPut('L2norm',sqrt(AvrLonLat_xv(xv_HsfcError**2))/ &
                             sqrt(AvrLonLat_xv(xv_HsfcRef**2)))

    call MPI_REDUCE(maxval(abs(xv_HsfcError)), & ! MPI_REDUCE ですべてのプロセスの中から
                                                 ! 演算し，ip==0 プロセスへ値をかえす
    HsfcmaxDef, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, MPI_COMM_WORLD, IERR)

    call MPI_REDUCE(maxval(abs(xv_HsfcRef)), &
    HsfcmaxRef, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, MPI_COMM_WORLD, IERR)
    if (ip==0) then
    call HistoryPut('Linfnorm', HsfcmaxDef/HsfcmaxRef)
    end if

    call HistoryPut('strfunc', xv_w(w_StrFunc)) !  (lon,lat,t)
    call HistoryPut('velpot',  xv_w(w_VelPot)) !  (lon,lat,t)

    xv_VelLon = xv_VelLonCosLat/cos(xv_Lat) !
    xv_VelLat = xv_VelLatCosLat/cos(xv_Lat) !
    call HistoryPut('vellon',xv_VelLon) !  (lon,lat,t)
    call HistoryPut('vellat',xv_VelLat) !  (lon,lat,t)

    call HistoryPut('am', v_AvrLon_xv(xv_VelLonCosLat)*Radius) ! (lat,t)
                        ! y_AvrLon_xv => v_AvrLon_xv
                        ! am = zonal mean angular momentum
    call HistoryPut('qv', (xv_Coli + xv_Vor)/(Hbar + xv_Hsfc)) !  (lon,lat,t)
                        ! potential vorticity
    call HistoryPut('ek',0.5*AvrLonLat_xv((Hbar+xv_hsfc) & !  (t)
                              *(xv_Vellon**2+xv_Vellat**2)))
                        ! kinetic energy
                        ! AvrLonLat_xv 2次元緯度経度格子点データの全領域平均
    call HistoryPut('ep',0.5*Grav*AvrLonLat_xv(xv_Hsfc**2)) !  (t)
                        ! potential energy
    call HistoryPut('en',0.5*AvrLonLat_xv(& !  (t)
            (Hbar+xv_hsfc)*(xv_Vellon**2+xv_Vellat**2)+Grav*xv_Hsfc**2))
                        ! Total energy

   !---- エネルギースペクトル診断出力 ----
!    nm_EvSp = nm_EnergyFromStreamfunc_w(w_StrFunc) ! 渦エネルギースペクトル
!    n_EvSp = n_EnergyFromStreamfunc_w(w_StrFunc) ! 発散エネルギースぺクトル
!
!    call HistoryPut('evspnm',nm_EvSp/Radius**2) ! (n,m,t)
!                   ! vortical kinetic energy spectrum (n,m space)
!    call HistoryPut('evspn',n_EvSp/Radius**2) ! (n,t)
!                   ! vortical kinetic energy spectrum
!    call HistoryPut('evsptot',(4*pi)*sum(n_Evsp)) ! (t)
!                   ! total vortical kinetic energy
!    call HistoryPut('evsp',sum(n_Evsp)/Radius**2) ! (t)
!                   ! mean vortical kinetic energy
!    nm_EdSp = nm_EnergyFromStreamfunc_w(w_VelPot)
!    n_EdSp = n_EnergyFromStreamfunc_w(w_VelPot)
!
!    call HistoryPut('edspnm',nm_EdSp/Radius**2) ! (n,m,t)
!                   ! divergent kinetic energy spectrum (n,m space)
!    call HistoryPut('edspn',n_EdSp/Radius**2) ! (n,t)
!                   ! divergent kinetic energy spectrum
!    call HistoryPut('edsptot',(4*pi)*sum(n_Edsp)) ! (t)
!                   ! total divergent kinetic energy
!    call HistoryPut('edsp',sum(n_Edsp)/Radius**2) ! (t)
!                   ! mean divergent kinetic energy
!    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) ! (n,m,t)
!                   ! kinetic energy spectrum (n,m space)
!    call HistoryPut('ekspn',n_EkSp/Radius**2) ! (t)
!                   ! kinetic energy spectrum
!    call HistoryPut('eksptot',(4*pi)*sum(n_EkSp)) ! (t)
!                   ! total kinetic energy
!    call HistoryPut('eksp',sum(n_EkSp)/Radius**2) ! (t)
!                   ! mean kinetic energy
!    nm_EpSp = nm_EnergyFromSurfaceHeight_w(w_Hsfc)
!    n_EpSp = n_EnergyFromSurfaceHeight_w(w_Hsfc)
!
!    call HistoryPut('epspnm',nm_EpSp/Radius**2) ! (n,m,t)
!                   ! potential energy spectrum (n,m space)
!    call HistoryPut('epspn',n_EpSp/Radius**2) ! (n,t)
!                   ! potential energy spectrum
!    call HistoryPut('epsptot',(4*pi)*sum(n_Epsp)) ! (t)
!                   ! total potential energy
!    call HistoryPut('epsp',sum(n_Epsp)/Radius**2) ! (t)
!                   ! mean potential energy
!    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) ! (n,m,t)
!                   ! total energy spectrum (n,m space)
!    call HistoryPut('enspn',n_ESp/Radius**2) ! (n,t)
                   ! total energy spectrum
!    call HistoryPut('ensptot',(4*pi)*sum(n_ESp)) ! (t)
!                   ! mean total energy
!    call HistoryPut('ensp',sum(n_ESp)/Radius**2) ! (t)
!                   ! total energy

  end subroutine output_history

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

end program spshallow_zd_rn4cn_freedecay