!Title(表題): Two-dimensional barotropic model on a double cyclic region
!             (2 重周期境界領域上の2次元順圧モデル)
!
!History(履歴):
!   2016/01/28   岡崎 正悟     新規作成
!   2016/07/04   岡崎 正悟     L&H と整合的な結果が得られたことを確認した版
!   2016/09/16   岡崎 正悟     ラプラシアンを改訂
!   2016/09/30   岡崎 正悟     ヤコビアンを改訂
!   2016/10/31   岡崎 正悟     上層と下層で粘性係数を変える場合に対応
!
! The governing equation(支配方程式): Larichev and Held(1995)
!  2層モデル
!    Barotropic モード:
!      d\nabla^2\psi/dt
!          + J(\psi, \nabla^2 \psi) + J(\tau, \nabla^2 \tau)
!            + U d\nabla^2\tau/dx
!              = -\kappa\nabla^2(\phi-\tau)/2-\nu\nabla^{8}\nabla^2\phi
!    Baroclinic モード:
!      d(\nabla^2-k_d^2)\tau/dt
!          + J(\psi, (\nabla^2 -k_d^2)\tau) + J(\tau, \nabla^2\phi)
!            + Ud\nabla^2\phi/dx + U k_d^2 d\phi/dx
!              = -\kappa\nabla^2(\tau-\phi)/2-\nu\nabla^{8}(\nabla^2-k_d^2)\tau
program two_layer_double_cyclic

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

!========= 宣言部 =========
  !--------- 変数(格子点データ) ---------
   !----- 順圧(barotropic)モード -----
    double precision, allocatable :: yx_tro_VelX(:,:)      ! 速度経度成分
    double precision, allocatable :: yx_tro_VelY(:,:)      ! 速度緯度成分
    double precision, allocatable :: yx_tro_Vor(:,:)       ! 渦度(鉛直成分)
    double precision, allocatable :: yx_tro_StrFunc(:,:)   ! 流線関数
    double precision              :: tro_ektot
    double precision              :: tro_enstot
   !----- 傾圧(baroclinic)モード -----
    double precision, allocatable :: yx_cli_VelX(:,:)      ! 速度経度成分
    double precision, allocatable :: yx_cli_VelY(:,:)      ! 速度緯度成分
    double precision, allocatable :: yx_cli_Vor(:,:)       ! 渦度(鉛直成分)
    double precision, allocatable :: yx_cli_StrFunc(:,:)   ! 流線関数
    double precision              :: cli_ektot
    double precision              :: cli_enstot
  !--------- 変数(スペクトルデータ) ---------
   !----- 順圧(barotropic)モード -----
    double precision, allocatable :: ee_tro_Vor(:,:)       ! 渦度(鉛直成分)
    double precision, allocatable :: ee_tro_StrFunc(:,:)   ! 流線関数
    double precision, allocatable :: ee_tro_Vor_tmp(:,:)   ! 渦度
    
    double precision, allocatable :: ee_tro_ESp(:,:)       ! エネルギースペクトル
    double precision, allocatable :: ee_tro_EnsSp(:,:)     ! エンストロフィースペクトル
    
    double precision, allocatable :: e_tro_ESp(:)          ! 1次元エネルギースペクトル
    double precision, allocatable :: e_tro_EnsSp(:)        ! 1次元エンストロフィースペクトル
    
    double precision, allocatable :: ee_tro_k1(:,:)        ! 拡散以外の渦度時間変化(1段目)
    double precision, allocatable :: ee_tro_k2(:,:)        ! 拡散以外の渦度時間変化(2段目)
    double precision, allocatable :: ee_tro_k3(:,:)        ! 拡散以外の渦度時間変化(3段目)
    double precision, allocatable :: ee_tro_k4(:,:)        ! 拡散以外の渦度時間変化(4段目)

    double precision, allocatable :: ee_tro_J(:,:)
    
    double precision              :: tro_ektotsp
    double precision              :: tro_enstotsp
    
    ! ----- 傾圧(baroclinic)モード -----
    double precision, allocatable :: ee_cli_Vor(:,:)       ! 渦度(鉛直成分)
    double precision, allocatable :: ee_cli_T(:,:)
    double precision, allocatable :: ee_cli_StrFunc(:,:)   ! 流線関数
    double precision, allocatable :: ee_cli_Vor_tmp(:,:)   ! 渦度
    
    double precision, allocatable :: ee_cli_ESp(:,:)       ! エネルギースペクトル
    double precision, allocatable :: ee_cli_EnsSp(:,:)     ! エンストロフィースペクトル
    
    double precision, allocatable :: e_cli_ESp(:)          ! 1次元エネルギースペクトル
    double precision, allocatable :: e_cli_EnsSp(:)        ! 1次元エンストロフィースペクトル
    
    double precision, allocatable :: ee_cli_k1(:,:)        ! 拡散以外の渦度時間変化(1段目)
    double precision, allocatable :: ee_cli_k2(:,:)        ! 拡散以外の渦度時間変化(2段目)
    double precision, allocatable :: ee_cli_k3(:,:)        ! 拡散以外の渦度時間変化(3段目)
    double precision, allocatable :: ee_cli_k4(:,:)        ! 拡散以外の渦度時間変化(4段目)

    double precision, allocatable :: ee_cli_J(:,:)
    
    double precision              :: cli_ektotsp
    double precision              :: cli_enstotsp
    
  !--------- 固定パラメータ ---------
    double precision, parameter  :: pi = acos(-1.0D0) ! 円周率
    character(len=20)   :: DbgMessageFmt='*** DbgMESSAGE ***'
    double precision, parameter  :: vmiss = -999.0           ! 欠損値

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

   !--- NAMELIST 入力用デフォルトファイル名 ---
    character(len=100)  :: nmlfile='2-layer-cyclic.nml'
    logical             :: Verbose=.false.          ! 冗長なメッセージ出力 On/Off
    logical             :: DebugOn=.false.          ! デバッグメッセージ On/Off
    namelist /message/  Verbose, DebugOn           !
   !--- 格子点・スペクトル ---
    integer :: km=21                               ! X 方向切断全波数
    integer :: lm=21                               ! Y 方向切断全波数
    integer :: im=64                               ! X 方向格子点数 (>3*km)
    integer :: jm=64                               ! Y 方向格子点数 (>3*lm)
    namelist /gridset/ km, lm, im, jm
   !--- 初期エネルギー ---
    double precision    :: Nmin=0.0D0             ! 初期エネルギー分布の全波数領域域最小値
    double precision    :: Nmax=21.0D0            ! 初期エネルギー分布の全波数領域最大値
    integer             :: Seed=0                 ! seed(1)に設定する種の値
    double precision    :: tro_Etotal=1.0D0       ! 初期平均エネルギーの値
    double precision    :: cli_Etotal=1.0D0       ! 初期平均エネルギーの値 
    namelist /initvalue/ Seed, Nmin, Nmax, tro_Etotal, cli_Etotal
   !--- 物理パラメータ ---
    double precision            :: XLength=1.0D0         ! 領域の大きさ(X 方向)
    double precision            :: YLength=1.0D0         ! 領域の大きさ(Y 方向)
    double precision            :: Beta=0.0D0            ! βパラメター
    integer                     :: HVOrder=4             ! 超粘性の次数(1 で普通の粘性,
                                                         ! 水平ラプラシアンの階数)
    double precision            :: HVisc=1.0D0           ! 超粘性係数
    double precision            :: delta_HVisc=-1.0D0    ! 超粘性係数
    double precision            :: U=5.0D-3              ! 平均流の速さ
    double precision            :: kappa=4.0D-2          ! 粘性係数
    double precision            :: kd=5.0D1              ! 内部変形半径
    namelist /physics/  XLength, YLength, Beta, HVOrder, HVisc, delta_HVisc, &
                      & U, kappa, kd
   !--- 初期値 ---
    character(len=100) :: initial_file=''          ! 初期値データファイル名
                                                   ! (空なら内部で初期値を計算)
    real               :: initial_time=0.0         ! 初期時刻
    namelist /initial/ initial_file, initial_time  
   !--- 時間積分 ---
    double precision :: delta_t=1.0e-7              ! 時間積分刻み
    integer          :: nstep=2000                  ! 時間積分ステップ数
    namelist /tint/    delta_t, nstep
   !--- ヒストリー出力 ---
    character(len=100) :: hst_file='2-layer-cyclic.nc'      ! ヒストリーファイル名
    character(len=100) :: title = &                         ! タイトル
         'Simuration of QG turbulence with two-layer model'
    integer            :: hst_intstep=200                   ! ヒストリー出力間隔ステップ数
    namelist /history/    hst_file, title, hst_intstep
   !--- リスタート出力 ---
    character(len=100) :: rst_file='2-layer-cyclic-rst_out.nc' ! リスタート出力ファイル名
    integer :: rst_intstep=200                                 ! リスタート出力間隔ステップ数
    namelist /restart/   rst_file, rst_intstep 
  !---- 作業変数 ----
    integer                       :: it=0                   ! 時間ステップ
    double precision              :: time                   ! モデル内時間
    integer                       :: k, l, n, wn            ! 波数
    integer                       :: i, j
    integer                       :: Nmax_int
    integer,allocatable           :: k_tmp(:), l_tmp(:)
    type(GT_HISTORY)              :: hst_rst                ! リスタート GT_HISTORY 変数

  !--------- NAMELIST 読み込み ---------
    call Open(arg)
    call Debug(arg) ; call Help(arg) ; call Strict(arg)
    call Get(arg, argv)
    if ( size(argv) .le. 0 ) then
       call MessageNotify('W','plbaro-beta_rn4excs_exp1',&
            'There is no argument. The defalut namlist file, '&
            //trim(nmlfile)//' is used.')
    else
       nmlfile=argv(1)
       call MessageNotify('M','plbaro-beta_rn4excs_exp1',&
            'Namelist file is '//trim(nmlfile))
    endif
    deallocate(argv)
    call Close(arg)
    
    write(6,nml=message) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=message) ; write(6,nml=message) ; close(10)
    
    if (verbose) write(6,nml=gridset) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=gridset) ; write(6,nml=gridset) ; close(10)
    
    if (verbose) write(6,nml=physics) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=physics) ; write(6,nml=physics) ; close(10)
    
    if (verbose) write(6,nml=initial) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=initial) ; write(6,nml=initial) ; close(10)
    
    if (verbose) write(6,nml=initvalue) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=initvalue) ; write(6,nml=initvalue) ; close(10)
    
    if (verbose) write(6,nml=tint) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=tint) ; write(6,nml=tint) ; close(10)
    
    if (verbose) write(6,nml=history) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=history) ; write(6,nml=history) ; close(10)
    
    if (verbose) write(6,nml=restart) 
    open(10,file=nmlfile,status='OLD')
    read(10,nml=restart) ; write(6,nml=restart) ; close(10)
  
  !--------- デバッグ出力制御設定 ---------
    if (DebugOn) then
       call SetDebug
    end if
  !--------- 変数の割り付け ---------
    wn = nint(sqrt(real(lm**2+km**2)))
    Nmax_int = idint(Nmax)
   !----- 順圧(barotropic)モード -----
    allocate(yx_tro_VelX(0:jm-1,0:im-1),yx_tro_VelY(0:jm-1,0:im-1))
    allocate(yx_tro_Vor(0:jm-1,0:im-1),yx_tro_StrFunc(0:jm-1,0:im-1))
    allocate(ee_tro_Vor(-lm:lm,-km:km),ee_tro_StrFunc(-lm:lm,-km:km))   
    allocate(ee_tro_k1(-lm:lm,-km:km),ee_tro_k2(-lm:lm,-km:km))
    allocate(ee_tro_k3(-lm:lm,-km:km),ee_tro_k4(-lm:lm,-km:km))
    allocate(ee_tro_ESp(-lm:lm,-km:km),ee_tro_EnsSp(-lm:lm,-km:km))
    allocate(ee_tro_J(-lm:lm,-km:km))
    allocate(e_tro_ESp(1:Nmax_int), e_tro_EnsSp(1:Nmax_int))
   !----- 傾圧 (baroclinic)モード -----
    allocate(yx_cli_VelX(0:jm-1,0:im-1),yx_cli_VelY(0:jm-1,0:im-1))
    allocate(yx_cli_Vor(0:jm-1,0:im-1),yx_cli_StrFunc(0:jm-1,0:im-1))
    allocate(ee_cli_Vor(-lm:lm,-km:km),ee_cli_StrFunc(-lm:lm,-km:km))
    allocate(ee_cli_T(-lm:lm,-km:km))
    allocate(ee_cli_k1(-lm:lm,-km:km),ee_cli_k2(-lm:lm,-km:km))
    allocate(ee_cli_k3(-lm:lm,-km:km),ee_cli_k4(-lm:lm,-km:km))
    allocate(ee_cli_ESp(-lm:lm,-km:km),ee_cli_EnsSp(-lm:lm,-km:km))
    allocate(ee_cli_J(-lm:lm,-km:km))
    allocate(e_cli_ESp(1:Nmax_int), e_cli_EnsSp(1:Nmax_int))
  !--------- 座標値の設定 ---------
    call DbgMessage(fmt='call %c', c1='ee_initial') 
    call ee_Initial(im,jm,km,lm,0.0D0,XLength,0.0D0,YLength)
  !--------- 初期値設定 ---------
    time = initial_time

    if ( initial_file == "") then
       ! リスタートファイルを指定しない場合内部で ee_Vor を与える. 
       call set_initial_values
    else
       call HistoryGet( trim(initial_file), 'ee_tro_vor', ee_tro_Vor)
       call HistoryGet( trim(initial_file), 'ee_cli_vor', ee_cli_Vor)
    end if
       
    ee_tro_Strfunc = ee_LaplacianInv_ee(ee_tro_Vor)
    ee_cli_Strfunc = ee_cli_Vor2Strfunc_ee(ee_cli_Vor)
    
    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.')

  !--------- 時間積分(Runge-Kutta 法) ---------
  do it=1,nstep
     time = initial_time + it * delta_t
     ee_tro_k1 = 0.0D0
     ee_tro_k2 = 0.0D0
     ee_tro_k3 = 0.0D0
     ee_tro_k4 = 0.0D0
     ee_cli_k1 = 0.0D0
     ee_cli_k2 = 0.0D0
     ee_cli_k3 = 0.0D0
     ee_cli_k4 = 0.0D0
    !---- 1 段目 [ k1 = f( x_n, t_n) ] ----
     !--- Barotropic モード
      ee_tro_k1 = ee_tro_f_ee_ee(ee_tro_Vor,ee_cli_Vor) 
     !--- Baroclinic モード
      ee_cli_k1 = ee_cli_f_ee_ee(ee_tro_Vor,ee_cli_Vor)
    !---- 2 段目 [ k2 = f( x_n+k1*dt/2, t_n+dt/2 ) ] ----
     !--- Baroclinic モード
      ee_tro_k2 = ee_tro_f_ee_ee(ee_tro_Vor + delta_t * ee_tro_k1 / 2.0D0 &
                                ,ee_cli_Vor + delta_t * ee_cli_k1 / 2.0D0)
     !--- Baroclinic モード
      ee_cli_k2 = ee_cli_f_ee_ee(ee_tro_Vor + delta_t * ee_tro_k1 / 2.0D0 &
                                ,ee_cli_Vor + delta_t * ee_cli_k1 / 2.0D0)
    !---- 3 段目 [ k3 = f( x_n+k2*dt/2, t_n+dt/2 ) ] ----
     !--- Barotropic モード
      ee_tro_k3 = ee_tro_f_ee_ee(ee_tro_Vor + delta_t * ee_tro_k2 / 2.0D0 &
                                ,ee_cli_Vor + delta_t * ee_cli_k2 / 2.0D0) 
     !--- Baroclinic モード
      ee_cli_k3 = ee_cli_f_ee_ee(ee_tro_Vor + delta_t * ee_tro_k2 / 2.0D0 &
                                ,ee_cli_Vor + delta_t * ee_cli_k2 / 2.0D0)
    !---- 4 段目 [ k4 = f( x_n+k3*dt, t_n+dt ) ] ----
     !--- Barotropic モード
      ee_tro_k4 = ee_tro_f_ee_ee(ee_tro_Vor + delta_t * ee_tro_k3 &
                                ,ee_cli_Vor + delta_t * ee_cli_k3)
     !--- Baroclinic モード
      ee_cli_k4 = ee_cli_f_ee_ee(ee_tro_Vor + delta_t * ee_tro_k3 &
                                ,ee_cli_Vor + delta_t * ee_cli_k3)
    !--------- 積分 ---------
     !--- Barotropic モード
      ee_tro_Vor = ee_tro_Vor + delta_t * (ee_tro_k1 / 6.0D0 + ee_tro_k2 / 3.0D0 &
                                         + ee_tro_k3 / 3.0D0 + ee_tro_k4 / 6.0D0 )
     !--- Baroclinic モード
      ee_cli_Vor = ee_cli_Vor + delta_t * (ee_cli_k1 / 6.0D0 + ee_cli_k2 / 3.0D0 &
                                         + ee_cli_k3 / 3.0D0 + ee_cli_k4 / 6.0D0 )
           
     if(mod(it,hst_intstep) .eq. 0)then                    ! ヒストリー出力
        call output_history
     endif

     if(mod(it,rst_intstep) .eq. 0)then                    ! リスタート出力
        call output_restart
     endif
  enddo
 
  call DbgMessage(fmt='%c %c',      &
&                 c1=DbgMessageFmt, &
&                 c2='Time integration end.')

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

  call output_restart_close
  call output_history_close
 
! 以上 メインプログラム 
!========================================================================
! 以下 サブルーチン

contains
 
  function e_tro_eksp_ee(ee_tro_Vor)
  !----------------------------------------------------------------------------
  ! Barotropic 渦度から 1 次元エネルギースペクトル計算
  !----------------------------------------------------------------------------
    double precision                :: e_tro_eksp_ee(wn)
    double precision, intent(IN)    :: ee_tro_Vor(-lm:lm,-km:km)         ! 渦度(鉛直成分)
    ! ---- エネルギースペクトル計算
    e_tro_eksp_ee = 0.0D0
    do k=-km,km
       do l=-lm,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if ((l .eq. 0) .AND. (k .eq. 0)) then
             goto 10
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 10
          else
             e_tro_eksp_ee(n) &
                  = e_tro_eksp_ee(n) &
                  + ee_tro_Vor(l,k) * ee_tro_Vor(l,k) / dble(l**2+k**2)
          endif
10        continue
       enddo
    enddo
  end function e_tro_eksp_ee

  function e_cli_eksp_ee(ee_cli_Vor)
  !----------------------------------------------------------------------------
  ! Baroclinic 渦度から 1 次元エネルギースペクトル計算
  !----------------------------------------------------------------------------
    double precision                :: e_cli_eksp_ee(wn)
    double precision, intent(IN)    :: ee_cli_Vor(-lm:lm,-km:km)         ! 渦度(鉛直成分)
    ! ---- エネルギースペクトル計算
    e_cli_eksp_ee = 0.0D0
    do k=-km,km
       do l=-lm,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if ((l .eq. 0) .AND. (k .eq. 0)) then
             goto 60
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 60
          else
             e_cli_eksp_ee(n) &
                  = e_cli_eksp_ee(n) &
                  + ee_cli_Vor(l,k) * ee_cli_Vor(l,k) / dble(l**2+k**2+kd**2)
          endif
60        continue
       enddo
    enddo
  end function e_cli_eksp_ee
  
  function e_enssp_ee(ee_Vor)
  !----------------------------------------------------------------------------
  ! 渦度から 1 次元エンストロフィースペクトル計算
  !----------------------------------------------------------------------------
    double precision              :: e_enssp_ee(wn)
    double precision,intent(IN)   :: ee_Vor(-lm:lm,-km:km)         ! 渦度(鉛直成分)
    e_enssp_ee = 0.0D0
    do k=-km,km
       do l=-lm,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if ((l .eq. 0) .AND. (k .eq. 0)) then
             goto 20
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 20
          else
             e_enssp_ee(n) = e_enssp_ee(n) + ee_Vor(l,k) * ee_Vor(l,k)
          endif
20        continue
       enddo
    enddo
  end function e_enssp_ee

  function ee_Laplacian_ee(ee)
  !----------------------------------------------------------------------------
  ! ラプラシアンの計算(円形切断あり)
  !----------------------------------------------------------------------------
    double precision,intent(in) :: ee(-lm:lm,-km:km)
    double precision            :: ee_Laplacian_ee(-lm:lm,-km:km)

    ee_Laplacian_ee = 0.0D0
    do l=-lm,lm
       do k=-km,km
          n = nint(sqrt(dble(l**2+k**2)))
          if (n .eq. 0) then
             goto 100
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 100
          else
             ee_Laplacian_ee(l,k) = - ee(l,k) * dble(k**2+l**2)
          endif
100       continue
       end do
    end do
  end function ee_Laplacian_ee

  function ee_LaplacianInv_ee(ee)
  !----------------------------------------------------------------------------
  ! 逆ラプラシアンの計算(円形切断あり)
  !----------------------------------------------------------------------------
    double precision,intent(in) :: ee(-lm:lm,-km:km)
    double precision            :: ee_LaplacianInv_ee(-lm:lm,-km:km)
    
    ee_LaplacianInv_ee = 0.0D0
    do l=-lm,lm
       do k=-km,km
          n = nint(sqrt(dble(l**2+k**2)))
          if (n .eq. 0) then
             goto 110
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 110
          else
             ee_LaplacianInv_ee(l,k) = - ee(l,k) / dble(k**2+l**2)
          endif
110       continue
       end do
    end do
  end function ee_LaplacianInv_ee
  
  function ee_cli_Strfunc2Vor_ee(ee_cli_Strfunc)
  !----------------------------------------------------------------------------
  ! 傾圧(Baroclinic)モードの流線函数から渦度計算
  !----------------------------------------------------------------------------
    double precision,intent(in) :: ee_cli_Strfunc(-lm:lm,-km:km)
    double precision            :: ee_cli_Strfunc2Vor_ee(-lm:lm,-km:km)
    
    ee_cli_Strfunc2Vor_ee = ee_Lapla_ee(ee_cli_StrFunc) - kd**2.0 * ee_cli_Strfunc
    
  end function ee_cli_Strfunc2Vor_ee
 
  function ee_cli_Vor2Strfunc_ee(ee_cli_Vor)
  !----------------------------------------------------------------------------
  ! 傾圧(Baroclinic)モードの渦度から流線函数計算
  !----------------------------------------------------------------------------
    double precision :: ee_cli_Vor(-lm:lm,-km:km)
    double precision :: ee_cli_Vor2Strfunc_ee(-lm:lm,-km:km)
    integer          :: k,l

    ee_cli_Vor2Strfunc_ee = 0.0D0
    do k=-km,km
       do l=-lm,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if (k==0 .AND. l==0) then
             goto 120
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 120
          else
             ee_cli_Vor2Strfunc_ee(l,k) &
                  & = - ee_cli_Vor(l,k) / dble(k**2+l**2+kd**2)
          end if
120       continue
       end do
    end do
  end function ee_cli_Vor2Strfunc_ee

  function ee_Jac_ee_ee(ee_A,ee_B)
    !----------------------------------------------------------------------------
    ! ヤコビアンの計算(円形切断あり Ver)
    !----------------------------------------------------------------------------
    double precision, intent(IN) :: ee_A(-lm:lm,-km:km)
    double precision, intent(IN) :: ee_B(-lm:lm,-km:km)
    double precision             :: ee_Jac_ee_ee(-lm:lm,-km:km)
    double precision             :: ee_DxA(-lm:lm,-km:km)
    double precision             :: ee_DyB(-lm:lm,-km:km)
    double precision             :: ee_DyA(-lm:lm,-km:km)
    double precision             :: ee_DxB(-lm:lm,-km:km)
    double precision             :: yx_DxA(0:jm-1,0:im-1)
    double precision             :: yx_DyB(0:jm-1,0:im-1)
    double precision             :: yx_DyA(0:jm-1,0:im-1)
    double precision             :: yx_DxB(0:jm-1,0:im-1)
    double precision             :: yx_Jac_tmp(0:jm-1,0:im-1)
    
    ee_DxA = 0.0D0
    ee_DyB = 0.0D0
    ee_DyA = 0.0D0
    ee_DxB = 0.0D0
    do k=-km,km
       do l=-lm,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if (n .eq. 0) then
             goto 130
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 130
          else
             ee_DxA(l,k) = -dble(k)*ee_A(-l,-k)
             ee_DyB(l,k) = -dble(l)*ee_B(-l,-k)
             ee_DyA(l,k) = -dble(l)*ee_A(-l,-k)
             ee_DxB(l,k) = -dble(k)*ee_B(-l,-k)
          end if
130       continue
       end do
    end do
    
    yx_DxA = yx_ee(ee_DxA)
    yx_DyB = yx_ee(ee_DyB)
    yx_DyA = yx_ee(ee_DyA)
    yx_DxB = yx_ee(ee_DxB)
    
    yx_Jac_tmp = yx_DxA * yx_DyB - yx_DyA * yx_DxB
    ee_Jac_ee_ee = ee_yx(yx_Jac_tmp)
    
  end function ee_Jac_ee_ee
  
 !========================== 時間変化計算 ==========================
  !---------------------------
  ! Barotrobic mode
  !--------------------------- 
  function ee_tro_f_ee_ee(ee_tro_Vor,ee_cli_Vor)
    double precision :: ee_tro_Vor(-lm:lm,-km:km)    ! 順圧モードの渦度
    double precision :: ee_cli_Vor(-lm:lm,-km:km)    ! 傾圧モードの渦度
    double precision :: ee_tro_f_ee_ee(-lm:lm,-km:km)
    
    double precision :: ee_tro_Strfunc(-lm:lm,-km:km)
    double precision :: ee_cli_Strfunc(-lm:lm,-km:km)
    double precision :: ee_cli_T(-lm:lm,-km:km)
    
    double precision :: ee_tro_J1(-lm:lm,-km:km)   ! 非線型項
    double precision :: ee_tro_J2(-lm:lm,-km:km)   ! 非線型項
    double precision :: ee_tro_lin(-lm:lm,-km:km)  ! 線型項
    double precision :: ee_tro_dump(-lm:lm,-km:km) ! 拡散項
    double precision :: ee_tro_vis(-lm:lm,-km:km)  ! 粘性項
    
    ee_tro_J1      = 0.0D0
    ee_tro_J2      = 0.0D0
    ee_tro_lin     = 0.0D0
    ee_tro_dump    = 0.0D0 
    ee_tro_vis     = 0.0D0
    ee_tro_f_ee_ee = 0.0D0
    
    ee_tro_Strfunc = ee_LaplacianInv_ee(ee_tro_Vor)
    ee_cli_Strfunc = ee_cli_Vor2Strfunc_ee(ee_cli_Vor)
    ee_cli_T       = ee_Laplacian_ee(ee_cli_Strfunc)
    
    ! 非線型項の計算
    ee_tro_J1 = ee_Jac_ee_ee(ee_tro_Strfunc,ee_tro_Vor)
    ee_tro_J2 = ee_Jac_ee_ee(ee_cli_Strfunc,ee_cli_T)
    
    ! 線型項の計算
    ee_tro_lin = U * ee_Dx_ee(ee_cli_T)
    
    ! エクマンダンピングの項の計算
    ee_tro_dump = 0.5D0 * kappa  &
         & * (ee_Laplacian_ee(ee_tro_Strfunc) - ee_Laplacian_ee(ee_cli_Strfunc))
    
    ! 粘性項の計算
    do k=-km,km
       do l=-lm,lm
!          n = nint(sqrt(dble(l**2+k**2)))
!          if (n .eq. 0) then
!             goto 200
!          else if ((n .gt. km) .OR. (n .gt. lm)) then
!             goto 200
!          else
             ee_tro_vis(l,k) = HVisc * dble(k**2+l**2)**HVOrder * ee_tro_Vor(l,k) &
                  & + 0.5D0 * delta_HVisc  &
                  &   * dble(k**2+l**2)**HVOrder &
                  &   * (ee_tro_Vor(l,k) - ee_cli_Vor(l,k))
!          end if
!200    continue
       end do
    end do
    
    !----- 時間変化計算 -----
    ee_tro_f_ee_ee = - ee_tro_J1 - ee_tro_J2 - ee_tro_lin - ee_tro_dump - ee_tro_vis
  end function ee_tro_f_ee_ee
  
 !---------------------------
 ! Baroclinic mode
 !--------------------------- 
  function ee_cli_f_ee_ee(ee_tro_Vor,ee_cli_Vor)
    double precision :: ee_tro_Vor(-lm:lm,-km:km)    ! 順圧モードの渦度
    double precision :: ee_cli_Vor(-lm:lm,-km:km)    ! 傾圧モードの渦度
    double precision :: ee_cli_f_ee_ee(-lm:lm,-km:km)
 
    double precision :: ee_tro_Strfunc(-lm:lm,-km:km)
    double precision :: ee_cli_Strfunc(-lm:lm,-km:km)
    
    double precision :: ee_cli_J1(-lm:lm,-km:km)   ! 非線型項
    double precision :: ee_cli_J2(-lm:lm,-km:km)   ! 非線型項
    double precision :: ee_cli_lin(-lm:lm,-km:km)  ! 線型項
    double precision :: ee_cli_dump(-lm:lm,-km:km) ! 拡散項
    double precision :: ee_cli_vis(-lm:lm,-km:km)  ! 粘性項
    
    ee_cli_J1         = 0.0D0
    ee_cli_J2         = 0.0D0
    ee_cli_lin        = 0.0D0
    ee_cli_dump       = 0.0D0 
    ee_cli_vis        = 0.0D0
    ee_cli_f_ee_ee    = 0.0D0
    
    ee_tro_Strfunc = ee_LaplacianInv_ee(ee_tro_Vor)
    ee_cli_Strfunc = ee_cli_Vor2Strfunc_ee(ee_cli_Vor)
    
    ! 非線型項の計算
    ee_cli_J1 = ee_Jac_ee_ee(ee_tro_Strfunc,ee_cli_Vor)
    ee_cli_J2 = ee_Jac_ee_ee(ee_cli_Strfunc,ee_tro_Vor)
    
    ! 線型項の計算
    ee_cli_lin = U * (ee_Dx_ee(ee_tro_Vor) + kd**2 * ee_Dx_ee(ee_tro_Strfunc))
    
    ! エクマンダンピングの項の計算
    ee_cli_dump = 0.5D0 * kappa  &
         & * (ee_Laplacian_ee(ee_cli_Strfunc) - ee_Laplacian_ee(ee_tro_Strfunc))
    
    ! 粘性項の計算
    do k=-km,km
       do l=-lm,lm
!          n = nint(sqrt(dble(l**2+k**2)))
!          if (n .eq. 0) then
!             goto 210
!          else if ((n .gt. km) .OR. (n .gt. lm)) then
!             goto 210
!          else
             ee_cli_vis(l,k) = HVisc * dble(k**2+l**2)**HVOrder * ee_cli_Vor(l,k) &
                  & + 0.5D0 * delta_HVisc  &
                  &   * dble(k**2+l**2)**HVOrder &
                  &   * (ee_cli_Vor(l,k) - ee_tro_Vor(l,k))
!          end if
!210       continue
       enddo
    enddo
    
    !----- 時間変化計算 -----
    ee_cli_f_ee_ee = - ee_cli_J1 - ee_cli_J2 - ee_cli_lin - ee_cli_dump - ee_cli_vis
  end function ee_cli_f_ee_ee
  
!=========================== 初期値設定 ============================
 !------------------------------------------------------
 !  初期値設定(リスタートファイルない場合のデフォルト設定)
 !------------------------------------------------------
  subroutine set_initial_values
    integer              :: random_seed_size   ! 乱数の種の長さ
    integer, allocatable :: seedarray(:)       ! 乱数の種
    real                 :: harvest            ! 乱数の値
    real                 :: random_phase(-lm:lm,-km:km)
    integer              :: i_count(1:Nmax_int)
    
   !---- 初期の1次元エネルギースペクトルを与える----
    e_tro_ESp = 6.0D-8
    e_cli_ESp = 6.0D-8
   
   !---- カウンタの設定 ----
    i_count = 0
    k=0
    do l=1,Nmax_int
       i_count(l) = i_count(l)+1
    enddo
 
    do k = 1,km
       do l = -lm,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if (n .EQ. 0) then 
             goto 30
          endif
          if (n .gt. km) then
             goto 30
          endif
          i_count(n) = i_count(n) + 1
30        continue
       enddo
    enddo
    
   !---- 乱数設定
    call random_seed(size=random_seed_size)
    allocate(seedarray(random_seed_size))
    call random_seed(get=seedarray)
    seedarray(1)=Seed
    call random_seed(put=seedarray)
    
    do k=-km,km
       do l=-lm,lm
             call random_number(harvest)
             random_phase(l,k)= harvest            ! [0,1] の一様乱数
       enddo
    enddo
 
  !---- 1次元エネルギースペクトルから流線函数を計算
   !---------- barotropic ----------
    ee_tro_StrFunc = 0.0D0
    k=0
       do l=1,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if ( (n .ge. Nmin) .AND. (n .le. Nmax) ) then
             ee_tro_StrFunc(l,k) &
                  & = 2.0D+0 * sqrt(e_tro_ESp(n) / dble(k**2+l**2) / dble(i_count(n))) &
                  &   * cos(2.0*pi*random_phase(l,k))
             ee_tro_StrFunc(-l,k) &
                  & = 2.0D+0 * sqrt(e_tro_ESp(n) / dble(k**2+l**2) / dble(i_count(n))) &
                  &   * sin(2.0*pi*random_phase(l,k))
          endif
       enddo
    
    do k=1,km
       do l=-lm,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if ( (n .ge. Nmin) .AND. (n .le. Nmax) ) then
             ee_tro_StrFunc(l,k) &
                  & = 2.0D+0 * sqrt(e_tro_ESp(n) / dble(k**2+l**2) / dble(i_count(n))) &
                  &   * cos(2.0*pi*random_phase(l,k))
             ee_tro_StrFunc(l,-k) &
                  & = 2.0D+0 * sqrt(e_tro_ESp(n) / dble(k**2+l**2) / dble(i_count(n))) &
                  &   * sin(2.0*pi*random_phase(l,k))
          endif
       enddo
    enddo
   !---------- baroclinic ----------
    ee_cli_StrFunc = 0.0D0
    k=0
       do l=1,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if ( (n .ge. Nmin) .AND. (n .le. Nmax) ) then
             ee_cli_StrFunc(l,k) &
                  & = 2.0D+0 * sqrt(e_cli_ESp(n) &
                  & / (dble(k**2+l**2)+kd**2) / dble(i_count(n))) &
                  & * cos(2.0*pi*random_phase(l,k))
             ee_cli_StrFunc(-l,k) &
                  & = 2.0D+0 * sqrt(e_cli_ESp(n) &
                  & / (dble(k**2+l**2)+kd**2) / dble(i_count(n))) &
                  & * sin(2.0*pi*random_phase(l,k))
          endif
       enddo
   
    do k=1,km
       do l=-lm,lm
          n = nint(sqrt(dble(l**2+k**2)))
          if ( (n .ge. Nmin) .AND. (n .le. Nmax) ) then
             ee_cli_StrFunc(l,k) &
                  & = 2.0D+0 * sqrt(e_cli_ESp(n) &
                  & / (dble(k**2+l**2)+kd**2) / dble(i_count(n))) &
                  & * cos(2.0*pi*random_phase(l,k))
             ee_cli_StrFunc(l,-k) &
                  & = 2.0D+0 * sqrt(e_cli_ESp(n) &
                  & / (dble(k**2+l**2)+kd**2) / dble(i_count(n))) &
                  & * sin(2.0*pi*random_phase(l,k))
          endif
       enddo
    enddo
 
   !---- 全エネルギーの計算
    !---------- barotropic ----------
     tro_Etotal=0.0
     do k=1,Nmax_int
        tro_Etotal = tro_Etotal + e_tro_ESp(k)
     enddo
    !---------- baroclinic ----------
     cli_Etotal=0.0
     do k=1,Nmax_int
        cli_Etotal = cli_Etotal + e_cli_ESp(k)
     enddo

   !---- (平均)エネルギーを規格化
    !---------- barotropic ----------     
     ee_tro_StrFunc = ee_tro_StrFunc / sqrt(tro_Etotal*1.0D5)
     ee_tro_Vor = ee_Laplacian_ee(ee_tro_StrFunc)
    !---------- baroclinic ----------     
     ee_cli_StrFunc = ee_cli_StrFunc / sqrt(cli_Etotal*1.0D5)
     ee_cli_Vor = ee_cli_Strfunc2Vor_ee(ee_cli_Strfunc)

   end subroutine set_initial_values

!=========================== リスタート出力 ============================
   subroutine output_restart_init
   !------------------------------------------------------
   ! リスタート出力初期化
   !------------------------------------------------------  
     call HistoryCreate( &
          file=trim(rst_file), &
          title=trim(title), &
          source='2d-layer-cyclic.f90 (2016/01/05)', &
          institution='OKAZAKI Shogo',     &
          dims=(/'x','y','k','l','n','t'/), &
          dimsizes=(/im,jm,2*km+1,2*lm+1,Nmax_int,0/),&
          longnames=(/'X           ','Y           ',&
                      'X-wavenumber','Y-wavenumber',&
                      'wavenumber  ','time        '/),&
          units=(/'1','1','1','1','1','1'/),   &
          origin=real(time), interval=real(rst_intstep*delta_t), &
          history=hst_rst )
     
     !---- 座標変数定義, 出力 ----
     call HistoryPut('x',x_X, hst_rst)                       ! 変数出力
     call HistoryAddattr('x','topology','circular', hst_rst) ! 周期属性
     call HistoryAddattr('x','modulo',XLength, hst_rst)      ! 周期属性
     call HistoryPut('y',y_Y, hst_rst)                       ! 変数出力
     call HistoryAddattr('y','topology','circular', hst_rst) ! 周期属性
     call HistoryAddattr('y','modulo',YLength, hst_rst)      ! 周期属性
     
     call HistoryPut('k',(/(dble(k),k=-km,km)/), hst_rst)    ! 変数出力
     call HistoryPut('l',(/(dble(l),l=-lm,lm)/), hst_rst)    ! 変数出力
 
     call HistoryAddVariable( &                              ! 変数定義
          varname='x_weight', dims=(/'x'/), & 
          longname='Weight for integration in X', &
          units='1', xtype='double',history=hst_rst)
     call HistoryPut('x_weight',x_X_Weight,hst_rst)          ! 変数出力
     
     call HistoryAddVariable( &                              ! 変数定義
          varname='y_weight', dims=(/'y'/), & 
          longname='Weight for integration in Y', &
          units='1', xtype='double',history=hst_rst)
     call HistoryPut('y_weight',y_Y_Weight,hst_rst)          ! 変数出力
     
     call HistoryPut('n',(/(dble(n),n=1,Nmax_int)/), hst_rst)
     
     !---- 物理変数定義 ----
     call HistoryAddVariable( &
          varname='ee_tro_vor',dims=(/'k','l','t'/), &
          longname='Barotropic Vorticiy(spectra)',units='1',xtype='double', &
          history=hst_rst)
     call HistoryAddVariable( &
          varname='ee_cli_vor',dims=(/'k','l','t'/), &
          longname='Baroclinic Vorticiy(spectra)',units='1',xtype='double', &
          history=hst_rst)
     
     !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ----
     call HistoryAddAttr('x','+delta_t', delta_t ,hst_rst)
     call HistoryAddAttr('x','+Beta',    Beta    ,hst_rst)
     call HistoryAddAttr('x','+HVOrder', HVOrder ,hst_rst)
     call HistoryAddAttr('x','+HVisc',   HVisc   ,hst_rst)
     call HistoryAddAttr('x','+U'    ,   U       ,hst_rst)
     call HistoryAddAttr('x','+kappa',   kappa   ,hst_rst)
     call HistoryAddAttr('x','+kd'   ,   kd      ,hst_rst)
     
  end subroutine output_restart_init
  
  subroutine output_restart
   !------------------------------------------------------
   ! リスタート出力
   !------------------------------------------------------
    write(6,*) ' Restart file output at it = ',it, '  time = ', time
    call HistoryPut('t',real(time),hst_rst)
 
   !---- 物理変数出力 ----
    call HistoryPut('ee_tro_vor', ee_tro_Vor, hst_rst)
    call HistoryPut('ee_cli_vor', ee_cli_Vor, hst_rst)
  end subroutine output_restart

  subroutine output_restart_close
  !------------------------------------------------------
  ! リスタート出力終了
  !------------------------------------------------------
    call HistoryClose(hst_rst)
  end subroutine output_restart_close

!=========================== ヒストリー出力 ============================
  subroutine output_history_init
  !------------------------------------------------------
  ! ヒストリー出力初期化
  !------------------------------------------------------  
    call HistoryCreate( &
           file=trim(hst_file), &
           title=trim(title), &
           source='2-layer-cyclic.f90 (2016/01/05)', &
           institution='OKAZAKI Shogo',     &
           dims=(/'x','y','k','l','n','t'/), &
           dimsizes=(/im,jm,2*km+1,2*lm+1,Nmax_int,0/),&
           longnames=(/'X           ','Y           ',&
                       'X-wavenumber','Y-wavenumber',&
                       'wavenumber  ','time        '/),&
           units=(/'1','1','1','1','1','1'/),   &
           origin=real(time), interval=real(rst_intstep*delta_t) )

   !---- 座標変数定義, 出力 ----
    call HistoryPut('x',x_X )                               ! 変数出力
    call HistoryAddattr('x','topology','circular')          ! 周期属性
    call HistoryAddattr('x','modulo',XLength )              ! 周期属性
    call HistoryPut('y',y_Y )                               ! 変数出力
    call HistoryAddattr('y','topology','circular')          ! 周期属性
    call HistoryAddattr('y','modulo',YLength )              ! 周期属性

    call HistoryPut('k',(/(dble(k),k=-km,km)/) )        ! 変数出力
    call HistoryPut('l',(/(dble(l),l=-lm,lm)/) )        ! 変数出力

    call HistoryAddVariable( &                              ! 変数定義
           varname='x_weight', dims=(/'x'/), & 
           longname='Weight for integration in X', &
           units='1', xtype='double')
    call HistoryPut('x_weight',x_X_Weight)                  ! 変数出力

    call HistoryAddVariable( &                              ! 変数定義
           varname='y_weight', dims=(/'y'/), & 
           longname='Weight for integration in Y', &
           units='1', xtype='double')
    call HistoryPut('y_weight',y_Y_Weight)                  ! 変数出力

    call HistoryPut('n',(/(dble(n),n=1,Nmax_int)/) )

   !---- パラメター定義, 出力 ----
    !----- barotropic mode -----
     !---- 物理変数定義 ----
      call HistoryAddVariable( &                                ! 変数定義
             varname='tro_vor', dims=(/'x','y','t'/), & 
             longname='Barotropic Vorticity', units='1', xtype='double')
      call HistoryAddVariable( &
           varname='tro_vorsp',dims=(/'k','l','t'/), &
           longname='Barotropic Vorticiy(spectra)',units='1',xtype='double')

     !---- 診断量定義 ----
      call HistoryAddVariable( &                                ! 変数定義
             varname='tro_ektot', dims=(/'t'/), & 
             longname='Barotropic total kinetic energy', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
             varname='tro_enstot', dims=(/'t'/), & 
             longname='Barotropic total enstrophy', units='1', xtype='double')
      
      call HistoryAddVariable( &
             varname='tro_1deksp', dims=(/'n','t'/), &
             longname='Barotropic 1d Energy spectrum', &
             units='1', xtype='double')
      call HistoryAddVariable( &
             varname='tro_1denssp', dims=(/'n','t'/), &
             longname='Barotropic 1d Enstrophy spectrum', &
             units='1', xtype='double')
      
   !----- baroclinic mode -----
    !---- 物理変数定義 ----
     call HistoryAddVariable( &                                ! 変数定義
            varname='cli_vor', dims=(/'x','y','t'/), & 
            longname='Baroclinic Vorticity', units='1', xtype='double')
     call HistoryAddVariable( &
          varname='cli_vorsp',dims=(/'k','l','t'/), &
          longname='Baroclinic Vorticiy(spectra)',units='1',xtype='double')
     
    !---- 診断量定義 ----
     call HistoryAddVariable( &                                ! 変数定義
            varname='cli_ektot', dims=(/'t'/), & 
            longname='Baroclinic total kinetic energy', units='1', xtype='double')
     call HistoryAddVariable( &                                ! 変数定義
            varname='cli_enstot', dims=(/'t'/), & 
            longname='Baroclinic total enstrophy', units='1', xtype='double')
     
     call HistoryAddVariable( &
            varname='cli_1deksp', dims=(/'n','t'/), &
            longname='Baroclinic 1d Energy spectrum', &
            units='1', xtype='double')
     call HistoryAddVariable( &
            varname='cli_1denssp', dims=(/'n','t'/), &
            longname='Baroclinic 1d Enstrophy spectrum', &
            units='1', xtype='double')
          
  !---- 実験パラメターを属性として定義, 出力(全て Global 属性) ----
    call HistoryAddAttr('x','+delta_t', delta_t )
    call HistoryAddAttr('x','+Beta',    Beta    )
    call HistoryAddAttr('x','+HVOrder', HVOrder )
    call HistoryAddAttr('x','+HVisc',   HVisc   )
    call HistoryAddAttr('x','+U'    ,   U       )
    call HistoryAddAttr('x','+kappa',   kappa   )
    call HistoryAddAttr('x','+kd'   ,   kd   )

  end subroutine output_history_init

  subroutine output_history
  !------------------------------------------------------
  ! ヒストリー出力
  !------------------------------------------------------
    write(6,*) ' History file output at it = ',it, '  time = ', time
    call HistoryPut('t',real(time))

   !---- 物理変数出力 ----
    !----- barotropic mode ----- 
      yx_tro_Vor = yx_ee(ee_tro_Vor)
      ee_tro_StrFunc = ee_LaplacianInv_ee(ee_tro_Vor)
      yx_tro_StrFunc = yx_ee(ee_tro_StrFunc)

      call HistoryPut('tro_vor',transpose(yx_tro_Vor))
      call HistoryPut('tro_vorsp', transpose(ee_tro_Vor))
      
     !---- 診断量出力 ----
      call HistoryPut('tro_ektot', IntYX_yx(-yx_tro_StrFunc*yx_tro_Vor)/2.0D0)
      call HistoryPut('tro_enstot', IntYX_yx(yx_tro_Vor**2/2.0d0))
   
      ! 1次元エネルギー・エンストロフィースペクトル計算
      e_tro_ESp = e_tro_eksp_ee(ee_tro_Vor)
      e_tro_EnsSp = e_enssp_ee(ee_tro_Vor)
      call HistoryPut('tro_1deksp', e_tro_ESp )
      call HistoryPut('tro_1denssp', e_tro_EnsSp )
   
   !----- baroclinic mode ----- 
      yx_cli_Vor = yx_ee(ee_cli_Vor)
      ee_cli_StrFunc = ee_cli_Vor2Strfunc_ee(ee_cli_Vor)
      yx_cli_StrFunc = yx_ee(ee_cli_StrFunc)

      call HistoryPut('cli_vor',transpose(yx_cli_Vor))
      call HistoryPut('cli_vorsp', transpose(ee_cli_Vor))
      
     !---- 診断量出力 ----
      call HistoryPut('cli_ektot', IntYX_yx(-yx_cli_StrFunc*yx_cli_Vor)/2.0D0)
      call HistoryPut('cli_enstot', IntYX_yx(yx_cli_Vor**2/2.0d0))
   
      ! 1次元エネルギー・エンストロフィースペクトル計算
      e_cli_ESp = e_cli_eksp_ee(ee_cli_Vor)
      e_cli_EnsSp = e_enssp_ee(ee_cli_Vor)
      call HistoryPut('cli_1deksp', e_cli_ESp )
      call HistoryPut('cli_1denssp', e_cli_EnsSp )
   
    end subroutine output_history

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