!----------------------------------------------------------------------
!  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. 
!     
!     (   ) はBrueshaber の強制項
!
!     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+E]
!             - F_D^{Diff},
!
!     dh/dt = -\Dinv{a(1-\mu^2)}\DP{(hU)]}{\lamdba} 
!             -\Dinv{a}\DP{(hV)}{\mu} 
!             + (\sum S_storm + S_mass + S_rad)/g'
!             - 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/11/26 R. Suzuki  緯度経度の乱数を確かめる 
!          2020/11/26 R. Suzuki  ストームの開始時刻の乱数を確かめる
!          2020/11/27 R. Suzuki  ストームの正負 を確かめる
!          2020/11/27 R. Suzuki  1つのストームの生成し，消滅後，再び生成することを確かめる
!          2020/11/28 R. Suzuki  2つのストームの生成し，消滅後，再び生成することを確かめる
!          2020/11/28 R. Suzuki  Storm_init の範囲を delta_t <=  Storm_init <= Tau_interval かつ 
!                                                 mod (Storm_init, delta_t) == 0 になるように修正
!          2020/12/01 R. Suzuki  2.2Tau_storm + Tau_interval がdelta_t の倍数でなければ計算が走らないように修正
!          2020/12/01 R. Suzuki  spml2 用に乱数を書きかえる
!          2020/12/02 R. Suzuki  spml2 用に乱数を書きかえる, 放射緩和項を追加
!          2020/12/03 R. Suzuki  乱数を修正
!          2020/12/03 R. Suzuki  ストームの発生条件を修正

program spshallow_zd_rn4cn_freedecay

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

  use mpi
  use w_mpi_module ! add
  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
  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_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 :: xv_KEnergy(:,:)       ! 運動エネルギー add rsuzuki


 !---- 変数(スペクトルデータ) ----
  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_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段)

 !---- Bruehsnaer 用変数追加 -----
  !-------   Storm lat-lon  -------------
  integer :: Storm_num
  real(8) :: Lat_north_max                  ! ストームの発生最大緯度(north)
  real(8) :: Lat_north_min                  ! ストームの発生最小緯度
  real(8) :: Lat_south_max                  ! ストームの発生最大緯度(south)
  real(8) :: Lat_south_min                  ! ストームの発生最小緯度
  real(8) :: Lon_max                        ! ストームの発生最大経度
  real(8) :: Lon_min                        ! ストームの発生最小経度
  namelist /storm_lat_lon/  Storm_num, Lat_north_max, Lat_north_min,  &
                                       Lat_south_max, Lat_south_min,  &
                                       Lon_max, Lon_min

  real(8) :: Lat_ran           ! 緯度の乱数 (degree)
  real(8) :: Lon_ran           ! 経度の乱数 (degree)
  integer :: zero_or_one       ! 0 or 1
  real(8), allocatable :: Storm_Lat(:)      ! 経度の乱数を取得する (degree)
  real(8), allocatable :: Storm_Lon(:)      ! 経度の乱数を取得する (degree)

  !-------   Storm init time  -------------


  real(8) :: Storm_init            ! ストームが最初に発生する時間の乱数
  real(8), allocatable :: Storm_time_init(:)    ! ストームが最初に発生する時間の乱数を取得する 
 

  !-------   one Storm check  -------------

  real(8)            :: Tau_storm      = 1.0D5         ! ストームの緩和時間 S07 (5)
  real(8)            :: Tau_interval   = 2.0D5         ! ストームの注入時間間隔
  real(8)            :: t0             = 1.0D5         ! 最大注入率に達する時間, Taustorm と同じにしておく
  real(8)            :: Rstorm         = 1.5625D6      ! ストームの半径 S07 (5)
  real(8)            :: Smax           = 5.0D0         ! ストームの振幅 S07 (5)
  real(8)            :: Storm_alpha    = 1.0D0         ! 1.0(高気圧だけ) or 0.5 or 0.0 (低気圧だけ)
  integer            :: int_num        = 1.0           ! Storm_time_init  の最大値が Tau_interval の何倍か
  real(8)            :: Tau_mass       = 1.0D0         ! 放射緩和項(mass) の緩和時間
  real(8)            :: Tau_ape        = 1.0D0         ! 放射緩和項(ape) の緩和時間
  namelist /storm_para/  Tau_storm, Tau_interval, t0, Rstorm, Smax, Storm_alpha, int_num, &
                         Tau_mass, Tau_ape

  real(8), allocatable :: xv_OneStorm(:,:)    ! 1つのストーム
  real(8), allocatable :: w_OneStorm(:)
  real(8), allocatable :: Storm_Smax(:)       ! ストームの振幅の配列 S07 (5)
  real(8), allocatable :: t_star(:)       ! ストームの発生してからの時間
  integer :: Storm_index         ! ストームのインデックス

 !----------- sum Storm check ---------------------
  real(8), allocatable :: xv_SumStorm(:,:)    ! ストームの和
  real(8), allocatable :: w_SumStorm(:)       

 !----------- Storm rad  ---------------------
  real(8) :: Radmass                         ! 放射緩和項 (mass)
  real(8), allocatable :: xv_Radape(:,:)     ! 放射緩和項 (ape)
  real(8), allocatable :: xv_Rad(:,:)     ! 放射緩和項 (mass + ape)
  real(8), allocatable :: w_Rad(:)           ! 放射緩和項 (mass + ape)

 !----------- Storm other ---------------------
  integer :: loop_num              ! ループ用変数
  integer :: loop_storm_num        ! ループ用変数
  integer :: loop_storm_num_max    !
  real(8) :: Start_time
  real(8) :: rand1
  integer :: seedsize, i 
  integer, allocatable :: seed(:)
  integer :: SYS_TIME

 !----------- 出力用 ---------------------
  real(8), allocatable :: xv_AllHsfc(:,:)    ! xv_hsfc + Hbar
  real(8), allocatable :: xv_qv(:,:)         ! potential vorticity
  real(8), allocatable :: xv_Absvel(:,:)     ! xv_vellon**2 + xv_Vellat**2

!----------------------------------------------------------------------

 !---- 変数(保存量初期値) ---- add rsuzuki
  real(8) :: MassInit                        ! 質量初期値
  real(8) :: TEnergyInit                     ! エネルギー初期値
  real(8) :: EnstrophyInit                   ! エンストロフィー初期値

 !---- 固定パラメタ -----
  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)            :: RedGrav  = 9.8             ! 低減重力加速度
  real(8)            :: Hbar     = 4000.0D0  ! 初期の平均水深
  integer            :: HVOrder  = 2             ! 超粘性の次数
                                                 ! (水平ラプラシアンの階数)
  real(8)            :: HVisc    = 0.5D16        ! 超粘性係数
  integer            :: HDOrder  = 2             ! 超拡散の次数
                                                 ! (水平ラプラシアンの階数)
  real(8)            :: HDiff    = 0.5D16        ! 超拡散係数
  namelist /physics/ Radius, Omega, Alpha, RedGrav, 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=storm_lat_lon)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=storm_lat_lon) ; if( ip==0 ) write(6,nml=storm_lat_lon) ; close(10)

  if (verbose) write(6,nml=storm_para)
  open(10,file=nmlfile,status='OLD')
  read(10,nml=storm_para) ; if( ip==0 ) write(6,nml=storm_para) ; 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'

! Storm_time_init の乱数確認用
!  open(1, file = 'Storm_time_init.dat')

! Storm_alpha の確認用
!  open(2, file = 'Storm_alpha.dat')

! Storm_para の確認用
!  open(3, file = 'Storm_para.dat')

! Storm_lat_lon の確認用
!  open(4, file = 'Storm_lat_lon.dat')

! Storm_info の確認用
  open(5, file = 'Storm_info.dat')
  open(7, file = 'Storm_info_correct.dat')
  open(8, file = 'Storm_init.dat')

 !---------------- デバッグ出力制御設定 -----------------
  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))
  allocate(xv_KEnergy(0:im-1,jc))

  allocate(w_Vor(nc),w_Div(nc))
  allocate(w_Hsfc(nc))
  allocate(w_StrFunc(nc),w_VelPot(nc))

  allocate(w_VorTmp(nc),w_DivTmp(nc))
  allocate(w_HsfcTmp(nc))
  allocate(w_KEnergy(nc))

  allocate(xv_VelLonCosLat(0:im-1,jc),xv_VelLatCosLat(0:im-1,jc))

  allocate(w_DVorDtK1(nc),w_DVorDtK2(nc),  &
           w_DVorDtK3(nc),w_DVorDtK4(nc) )
  allocate(w_DDivDtK1(nc),w_DDivDtK2(nc),  &
           w_DDivDtK3(nc),w_DDivDtK4(nc) )
  allocate(w_DHsfcDtK1(nc),w_DHsfcDtK2(nc), &
           w_DHsfcDtK3(nc),w_DHsfcDtK4(nc) )

  allocate(xv_Coli(0:im-1,jc))
  allocate(w_HVisc(nc),w_HDiff(nc))

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

! Brueshaber 用
  !---------  Storm  lat-lon --------------
  allocate(Storm_Lat(1:Storm_num), Storm_Lon(1:Storm_num))
  !---------  Storm  init t  --------------
  allocate(Storm_time_init(1:Storm_num))
  !---------  one Storm check  ------------
  allocate(xv_OneStorm(0:im-1,jc))       
  allocate(w_OneStorm(nc))       
  allocate(Storm_Smax(1:Storm_num))
  allocate(w_SumStorm(nc))   
  allocate(t_star(1:Storm_num))
  !---------  two Storm check  ------------
  allocate(xv_SumStorm(0:im-1,jc))
  !---------  rad Storm check  ------------
  allocate(xv_Radape(0:im-1,jc))
  allocate(xv_Rad(0:im-1,jc))
  allocate(w_Rad((nm+1)*(nm+1)))
  !---------  rad Storm check  ------------
  allocate(xv_AllHsfc(0:im-1,jc))
  allocate(xv_qv(0:im-1,jc))
  allocate(xv_Absvel(0:im-1,jc))
  
  
  !------ 2.2Tau_Storm + Tau_interval がdelta_t の倍数かを判断 -----
!  if (ip == 0) then
!   write(*,*) 'dt, 2.2T_s + T_int, dt/(2.2T_s + T_int)', &
!               int(delta_t), int(2.2*Tau_storm + Tau_interval) , (2.2*Tau_Storm + Tau_interval)/delta_t
! 
!   if (mod(int(2.2*Tau_storm + Tau_interval), int(delta_t)) /= 0) then
!   write(*,*) '2.2Tau_storm + Tau_interval がdelta_t の倍数ではありません'
!
!   call MPI_ABORT(MPI_COMM_WORLD, 1, IERR)
!   else
!   write(*,*) 'delta_t OK!'
!   end if
!  end if

  ! ----- 計算開始時間の取得 ------ add rsuzuki
  call MXTIME(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 ) ! ok

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

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


  xv_KEnergy = (xv_VelLonCosLat**2 + xv_VelLatCosLat**2) & ! add rsuzuki
                       /(2.0*cos(xv_Lat)**2)




 !------------------ 緯度経度の乱数を作成 ----------------

    call random_seed(size=seedsize) !初期値のサイズを取得
    allocate(seed(seedsize)) !配列の割り当て

!    do i = 1, seedsize
!    call system_clock(SYS_TIME) !時間を取得
    call MXTIME(SYS_TIME)
    seed = SYS_TIME
    call random_seed(put = seed)
!    end do

     if (ip == 0) then
     do loop_num = 1, Storm_num
!   緯度方向の乱数
!   北,南半球のどちらかを選択する
      do i = 1, 10
      call random_number(rand1)    
      end do
      zero_or_one = int(rand1*2) 
    
      if (zero_or_one == 0) then
      do i = 1, 10
      call random_number(rand1)    
      end do
      Lat_ran = int(rand1*((Lat_north_max-Lat_north_min) +1)) + Lat_north_min  
                ! (Lat_north_max - Lat_north_min)  > int(R*( (Lat_north_max-Lat_north_min) +1)) >= 0
                !       Lat_north_max        >  Lat_ran   >=  Lat_north_min
      else
      do i = 1, 10
      call random_number(rand1)    
      end do
      Lat_ran = int(rand1*((Lat_south_max-Lat_south_min) +1)) + Lat_south_min  
      end if
  
!     経度方向の乱数
      do i = 1, 10
      call random_number(rand1)    
      end do
      Lon_ran = int(rand1*((Lon_max-Lon_min) +1)) + Lon_min 
  
      Storm_Lat(loop_num) = Lat_ran
      Storm_Lon(loop_num) = Lon_ran

    end do  
! test
!     Storm_Lat(1) = 0.0D0
!     Storm_Lon(1) = 100.0D0
!     Storm_Lat(2) = 0.0D0
!     Storm_Lon(2) = 300.0D0
     end if


     call MPI_Bcast(Storm_Lat(1), Storm_num, MPI_REAL8, 0, MPI_COMM_WORLD, IERR)
     call MPI_Bcast(Storm_Lon(1), Storm_num, MPI_REAL8, 0, MPI_COMM_WORLD, IERR)

 !------------------ ストームの最初に発生する時刻 : Tau_init の乱数を作成 ----------------
    if (ip == 0) then
     do loop_num = 1, Storm_num
!      m = int(Tau_interval / delta_t)
!      write(*,*) 'tau_interval / delta_t = ', m
      do i = 1, 10
      call random_number(rand1)    
      end do
      Storm_init = int(rand1*(Tau_interval*int_num - delta_t + 1)) + delta_t
                 ! Tau_interval*int_num - dt > int(rand1*(Tau_interval*int_num - dt + 1)) >= 0
                 ! Tau_interval*int_num  > int(rand1*(Tau_interval*int_num - dt + 1)) + dt  >= dt

 !     Storm_init = (int(rand1*m*int_num) + 1) * delta_t
                 ! 2*m-1 >= int(rand*2*m) >= 0
                 ! 2*m  >= int(rand*2*m) + 1 >= 1
      Storm_time_init(loop_num) = Storm_init
      write(*,*) 'time_init =', Storm_time_init(loop_num)
     end do
!     close(1)
    end if
! test
!     Storm_time_init(1) = 10.0D0
!     Storm_time_init(2) = 20.0D0

     call MPI_Bcast(Storm_time_init(1), Storm_num, MPI_REAL8, 0, MPI_COMM_WORLD, IERR)
!     write(*,*) "Storm_time_init", Storm_time_init
 !------------------ ストームの正負をStorm_alpha によって振り分ける ----------------
     if (ip == 0) then
     if (Storm_alpha == 1.0D0) then
       do loop_num = 1, Storm_num
       Storm_Smax(loop_num) = Smax
!       write(2,*) Storm_Smax(loop_num)
       end do
     else if (Storm_alpha == 0.5D0) then
       do loop_num = 1, Storm_num/2  ! Storm_num が奇数なら切り捨てられる
       Storm_Smax(loop_num) = Smax
       end do
       do loop_num = Storm_num/2 + 1, Storm_num  
       Storm_Smax(loop_num) = - Smax
       end do
     else
       do loop_num = 1, Storm_num
       Storm_Smax(loop_num) = - Smax
!       write(2,*) Storm_Smax(loop_num)
       end do 
     end if
!     close(2)
     end if

     call MPI_Bcast(Storm_Smax(1), Storm_num, MPI_REAL8, 0, MPI_COMM_WORLD, IERR)

 !------------------ ストームの情報出力 ---------------


!----------------------------- Storm parameter -------------------------------------
     100 FORMAT('delta_t =',es8.1, ',     Storm_num =',i3, ',     Rstorm =', es8.1)
     101 FORMAT('Tau_storm =',es8.1, ',   Tau_max =', es8.1, ',  Tau_interval =', es8.1)
     104 FORMAT('Tau_mass =',es8.1, ',   Tau_ape =', es8.1/)
     write(5,100) delta_t, Storm_num, Rstorm
     write(5,101) Tau_storm, 1.1*Tau_storm, Tau_interval
     write(5,104) Tau_mass, Tau_ape
!-----------------------------------------------------------------------------------


!----------------------------- ストームの初期の位置と，時間-------------------------
     102 FORMAT(i7,f7.3,i12,i12,i15,i14,i14)
     write(5,*) 'loop_num   Smax   Storm_Lat   Storm_Lon   Storm_init_t   Storm_max_t   Storm_disap' 
     do loop_num = 1, Storm_num
     write(5,102) loop_num, Storm_Smax(loop_num),  &
                  int(Storm_Lat(loop_num)), int(Storm_Lon(loop_num)), & 
                  int(storm_time_init(loop_num)), &
                  int(storm_time_init(loop_num)+1.1*tau_storm), &
                  int(storm_time_init(loop_num)+2.2*tau_storm)
     write(8,102) loop_num, Storm_Smax(loop_num),  &
                  int(Storm_Lat(loop_num)), int(Storm_Lon(loop_num)), & 
                  int(storm_time_init(loop_num)), &
                  int(storm_time_init(loop_num)+1.1*tau_storm), &
                  int(storm_time_init(loop_num)+2.2*tau_storm)
     end do
     close(8)
!-----------------------------------------------------------------------------------

!    correct 用
     write(7,100) delta_t, Storm_num, Rstorm
     write(7,101) Tau_storm, 1.1*Tau_storm, Tau_interval
     write(7,104) Tau_mass, Tau_ape

     loop_storm_num_max = int((delta_t*nstep - Storm_time_init(loop_num))/(2.2*tau_storm + Tau_interval ))

     202 FORMAT(i7,f7.3,i15,i14,i14)
     write(7,*) 'loop_num   Smax   Storm_init_t   Storm_max_t   Storm_disap' 
     do loop_num = 1, Storm_num
     write(7,202) loop_num, Storm_Smax(loop_num),  &
                  int(storm_time_init(loop_num)), &
                  int(storm_time_init(loop_num)+1.1*tau_storm), &
                  int(storm_time_init(loop_num)+2.2*tau_storm)
     end do

     do loop_storm_num = 1, loop_storm_num_max
       do loop_num = 1, Storm_num
       Start_time =  NINT(Storm_time_init(loop_num) + &
                     loop_storm_num*2.2*tau_storm + loop_storm_num*Tau_interval )
        if(Start_time <= initial_time + nstep*delta_t) then
        write(7,202) loop_storm_num, &
                    Storm_Smax(loop_num), &
                    int(Start_time), & 
                    int(Start_time + 1.1*tau_storm), &
                    int(Start_time + 2.2*tau_storm)
        end if
       end do
     end do
     close(7)
 !------------------- 初期値設定 ----------------------
  time = initial_time

  if ( initial_file == "") then
     ! リスタートファイルを指定しない場合, 
     ! 内部で w_Vor, w_Div, w_Hsfc を与える. 
     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 )
  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) ! 

 !------------------- 保存量初期値 -------------------- add rsuzuki

  MassInit = IntLonLat_xv(xv_Hsfc + Hbar)
  TEnergyInit = IntLonLat_xv( &
                         (xv_Hsfc + Hbar)*xv_KEnergy &
                       + 0.5*RedGrav*((xv_Hsfc + Hbar)**2))

  EnstrophyInit = IntLonLat_xv( 0.5*(xv_Vor+xv_Coli)**2/(xv_Hsfc + Hbar))

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

     xv_SumStorm = 0.0D0
     w_SumStorm  = 0.0D0
     do Storm_index = 1, Storm_num
       

       t_star(Storm_index) = mod(int(time - Storm_time_init(Storm_index) ), int(2.2*tau_storm + tau_interval ))

       if (t_star(Storm_index) > 2.2*Tau_Storm &
               .and. ip == 0 .and. time >= Storm_time_init(Storm_index)) then

      
!   緯度方向の乱数
!   北,南半球のどちらかを選択する
         call random_number(rand1)    
         zero_or_one = int(rand1*2) 
     
         if (zero_or_one == 0) then
         call random_number(rand1)    
         Lat_ran = int(rand1*((Lat_north_max-Lat_north_min) +1)) + Lat_north_min  
                   ! (Lat_north_max - Lat_north_min)  > int(R*( (Lat_north_max-Lat_north_min) +1)) >= 0
                   !       Lat_north_max       > Lat_ran   >=  Lat_north_min
         else
         call random_number(rand1)    
         Lat_ran = int(rand1*((Lat_south_max-Lat_south_min) +1)) + Lat_south_min  
         end if
   
!   経度方向の乱数
         call random_number(rand1)    
         Lon_ran = int(rand1*((Lon_max-Lon_min) +1)) + Lon_min 
     
         Storm_Lat(loop_num) = Lat_ran
         Storm_Lon(loop_num) = Lon_ran
         Storm_Lat(Storm_index) = Lat_ran
         Storm_Lon(Storm_index) = Lon_ran

       end if


       call MPI_Bcast(t_star(1), Storm_num, MPI_REAL8, 0, MPI_COMM_WORLD, IERR)
       call MPI_Bcast(Storm_Lat(1), Storm_num, MPI_REAL8, 0, MPI_COMM_WORLD, IERR)
       call MPI_Bcast(Storm_Lon(1), Storm_num, MPI_REAL8, 0, MPI_COMM_WORLD, IERR)

       if ( t_star(Storm_index)  == 0.0D0 .and. time /= Storm_time_init(Storm_index)) then  
         write(5,102) Storm_index, Storm_Smax(Storm_index),  &
                  int(Storm_Lat(Storm_index)), int(Storm_Lon(Storm_index)), & 
                  int(time), &
                  int(time+1.1*tau_storm), &
                  int(time+2.2*tau_storm)
       end if

       if ( t_star(Storm_index)  <=  2.2*Tau_storm .and. time >= Storm_time_init(Storm_index)) then  
       call Storm( Storm_index, t_star, & 
                 Storm_Smax, Rstorm, t0, Tau_storm, Tau_interval,    & ! (in)
                 Storm_Lat, Storm_Lon,                               & ! (in)
                 xv_OneStorm)       ! (out)

       else 
       xv_OneStorm = 0.0D0
       end if


     xv_SumStorm = xv_SumStorm + xv_OneStorm
     end do

     w_SumStorm = w_xv(xv_SumStorm)

! --------Storm_rad ------------------
     Radmass = - AvrLonLat_xv(RedGrav*xv_Hsfc/Tau_mass)
     xv_Radape  = - (RedGrav*xv_Hsfc - AvrLonLat_xv(RedGrav*xv_Hsfc))/Tau_ape
     xv_Rad = Radmass + xv_Radape
     w_Rad = w_xv(xv_Rad)


     !---- 1 段目 [ k1 = f( x_n, t_n) ] ----
     call TimeDerivetives_noDisppation(                 &
       w_Vor, w_Div, w_Hsfc,                            &  !(in) 渦度, 発散, 変位
       w_SumStorm, w_Rad,                               &  !(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_SumStorm, w_Rad,                               &  !(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_SumStorm, w_Rad,                               &  !(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_SumStorm, w_Rad,                               &  !(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 MXTIME(TIM1)
  print *, '計算時間[sec] : ', TIM1-TIM0

  call output_restart_close
  call output_history_close

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

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

contains

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

  w_Vor  = 0.0D0
  w_Div  = 0.0D0
  w_Hsfc = 0.0D0
  w_SumStorm = 0.0D0
  xv_OneStorm = 0.0D0
  xv_SumStorm = 0.0D0
  xv_Radape = 0.0D0
  xv_Rad  = 0.0D0
  Radmass  = 0.0D0

  t_star = 1.0D30


  end subroutine set_initial_values

 !----------- one Storm を計算する----------
  subroutine Storm( Storm_index, t_star, & 
                    Storm_Smax, Rstorm, t0, Tau_storm, Tau_interval,    & ! (in)
                    Storm_Lat, Storm_Lon,                               & ! (in)
                    xv_OneStorm)       ! (out)

   !---- 変数(入力) ----
    integer, intent(in) :: Storm_index   ! ストームインデックス
    real(8), intent(in) :: t_star(:)    ! ストーム が発生する時間
    real(8), intent(in) :: Storm_Smax(:) ! Smax の配列
    real(8), intent(in) :: Rstorm        ! ストームの半径
    real(8), intent(in) :: t0            ! 最大注入率に達する時間, Taustorm と同じにしておく
    real(8), intent(in) :: Tau_storm     ! ストームの緩和時間
    real(8), intent(in) :: Tau_interval  ! ストームの注入時間間隔
    real(8), intent(in) :: Storm_Lat(:)  ! ストームの中心緯度
    real(8), intent(in) :: Storm_Lon(:)  ! ストームの中心経度

   !---- 変数(出力) ----
    real(8), intent(out) :: xv_OneStorm(:,:)    ! 1つのストーム

   !--------------------
    integer :: i, j
    real(8) :: r
    real(8) :: stime

    do j = 1, jc ! williamson testcase1 を参考にした
       do i = 0, im-1
        r = Radius * acos(sin(Storm_Lat(Storm_index)*pi/180.0D0)*sin(xv_Lat(i,j)) &
               + cos(Storm_Lat(Storm_index)*pi/180.0D0)*cos(xv_Lat(i,j)) &
                    *cos(xv_Lon(i,j)-Storm_Lon(Storm_index)*pi/180.0D0))

          if ( r <  (2.2 * Rstorm) ) then
          xv_OneStorm(i,j) = Storm_Smax(Storm_index)*   &
                             exp( - r**2/Rstorm**2 - (t_star(Storm_index) - t0 )**2/Tau_storm**2)
          else
          xv_OneStorm(i,j) = 0.0D0
          end if
        end do
    end do

  end subroutine Storm


 !---------------- ポテンシャルエネルギースペクトル計算 ----------------
  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 * RedGrav* 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 * RedGrav * 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_SumStorm, w_Rad,                  &             !(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(in) :: w_SumStorm(:) ! ストーム強制
    real(8), intent(in) :: w_Rad(:) ! 放射強制

   !---- 変数(出力) ----
    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_KEnergy)

     ! 時間変化項の見積り
     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( RedGrav*w_Hsfc + 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 &
                       + w_SumStorm / RedGrav + w_Rad / RedGrav
    
  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  ! 

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

  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,nc,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)


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


   !----- Brueshaber 用出力 ----------
   
    call HistoryAddVariable( &                                ! 変数定義
           varname='phi', dims=(/'lon','lat','t  '/), & 
           longname='Geopotential', units='m2/s2', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='OneStorm', dims=(/'lon','lat','t  '/), & 
           longname='OneStorm', units='m', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='SumStorm', dims=(/'lon','lat','t  '/), & 
           longname='SumStorm', units='m', xtype='double')

    call HistoryAddVariable( &                                ! 変数定義
           varname='Rrad', dims=(/'lon','lat','t  '/), &
           longname='Storm Rad', units='m2/s3', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='radape', dims=(/'lon','lat','t  '/), &
           longname='Storm Tau_ape', units='m2/s3', xtype='double')

   !---- 診断変数定義 ---- add rsuzuki
    call HistoryAddVariable( &                                ! 変数定義
           varname='mass', dims=(/'t'/), &
           longname='Mass error', &
           units='1', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='toeng', dims=(/'t'/), &
           longname='Total energy error', &
           units='1', xtype='double')

   !---- Brueshaber用 ---- add rsuzuki
    call HistoryAddVariable( &                                ! 変数定義
           varname='massstorm', dims=(/'t'/), &
           longname='Mass Storm', units='m3', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='radmass', dims=(/'t  '/), &
           longname='Storm Tau_mass', units='m^2/sec^3', xtype='double')

    call historyaddvariable( &                                ! 変数定義
           varname='ke', dims=(/'t'/), & 
           longname='kinetic energy (S07)', units='m^6/sec^4', xtype='double')
    call historyaddvariable( &                                ! 変数定義
           varname='ape', dims=(/'t'/), & 
           longname='available potential energy', units='m^6/sec^4', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='nonqv', dims=(/'lon','lat','t  '/), & 
           longname='Nondimensional Potential vorticity', units='1', xtype='double')


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

    xv_Absvel = xv_VelLon**2 + xv_VelLat**2
    xv_AllHsfc = xv_Hsfc + Hbar
    xv_qv = (xv_Coli + xv_Vor) / xv_AllHsfc

    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_qv) !  (lon,lat,t)
                        ! potential vorticity 
   !---- brueshaber 用変数 ----
    call HistoryPut('phi', xv_AllHsfc*RedGrav)
    call HistoryPut('OneStorm', xv_OneStorm)
    call HistoryPut('SumStorm', xv_SumStorm)
    call HistoryPut('Rrad', xv_Rad)
    call HistoryPut('radape', xv_Radape)

   !---- 診断変数出力 ---- add rsuzuki
    call HistoryPut('mass', &
         (IntLonLat_xv(xv_AllHsfc)-MassInit)/MassInit)    

    call HistoryPut('toeng', &
         ( IntLonLat_xv((xv_AllHsfc)*xv_KEnergy &
          +0.5*RedGrav*((xv_AllHsfc)**2)) & 
          -TEnergyInit &
          )/TEnergyInit) ! ok ?

   !---- brueshaber 用変数 ----
    call HistoryPut('massstorm', IntLonLat_xv(xv_SumStorm))
    call HistoryPut('radmass', Radmass)
    call HistoryPut('ke', 0.5*IntLonLat_xv(RedGrav*xv_AllHsfc*xv_Absvel))
    call HistoryPut('ape', 0.5*IntLonLat_xv((RedGrav*xv_AllHsfc)**2 - &
                           AvrLonLat_xv(RedGrav*xv_AllHsfc)**2))
    call HistoryPut('nonqv', (xv_qv - xv_Coli/AvrLonLat_xv(xv_AllHsfc)) & 
                             *AvrLonLat_xv(xv_AllHsfc)/(2*Omega))

  end subroutine output_history

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

end program spshallow_zd_rn4cn_freedecay
