!
! L&H(1995) の支配方程式の各項の寄与を計算するためのプログラム
!
! 2016/09/30  岡崎 正悟  一応完成版
! 2016/10/31  岡崎 正悟  上層と下層で粘性係数を変える場合に対応
! 
program cal_ekterm
!========= モジュール引用宣言 =========
  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, CPrintf
  use dc_args
  implicit none

!========= 宣言部 =========
  double precision, allocatable :: ee_tro_Vor(:,:)       ! 渦度(鉛直成分)
  double precision, allocatable :: ee_tro_StrFunc(:,:)   ! 流線関数
  double precision, allocatable :: ee_cli_Vor(:,:)   ! 流線関数
  double precision, allocatable :: ee_cli_T(:,:)
  double precision, allocatable :: ee_cli_StrFunc(:,:)   ! 流線関数
  !--------- 各項のエネルギー輸送函数 ---------
   !-- barotropic mode
     double precision, allocatable :: ee_T_psi_1(:,:)
     double precision, allocatable :: e_T_psi_1(:)
     double precision, allocatable :: ee_T_psi_2(:,:)
     double precision, allocatable :: e_T_psi_2(:)
     double precision, allocatable :: ee_T_psi_3(:,:)
     double precision, allocatable :: e_T_psi_3(:)
     double precision, allocatable :: ee_T_psi_4(:,:)
     double precision, allocatable :: e_T_psi_4(:)
     double precision, allocatable :: ee_T_psi_5(:,:)
     double precision, allocatable :: e_T_psi_5(:)
   !-- baroclinic mode
     double precision, allocatable :: ee_T_tau_1(:,:)
     double precision, allocatable :: e_T_tau_1(:)
     double precision, allocatable :: ee_T_tau_2(:,:)
     double precision, allocatable :: e_T_tau_2(:)
     double precision, allocatable :: ee_T_tau_3(:,:)
     double precision, allocatable :: e_T_tau_3(:)
     double precision, allocatable :: ee_T_tau_4(:,:)
     double precision, allocatable :: e_T_tau_4(:)
     double precision, allocatable :: ee_T_tau_5(:,:)
     double precision, allocatable :: e_T_tau_5(:)
     double precision, allocatable :: ee_T_tau_6(:,:)
     double precision, allocatable :: e_T_tau_6(:)
     double precision, allocatable :: ee_T_tau_7(:,:)
     double precision, allocatable :: e_T_tau_7(:)
  !----- エネルギー散逸率
     double precision              :: epsilon_tro_1
     double precision              :: epsilon_tro_2
     double precision              :: epsilon_tro_3
     double precision              :: epsilon_tro_4
     double precision              :: epsilon_tro_5
     double precision              :: epsilon_tro_total
     double precision              :: epsilon_cli_1
     double precision              :: epsilon_cli_2
     double precision              :: epsilon_cli_3
     double precision              :: epsilon_cli_4
     double precision              :: epsilon_cli_5
     double precision              :: epsilon_cli_6
     double precision              :: epsilon_cli_7
     double precision              :: epsilon_cli_total
     
     double precision, allocatable :: epsilon_tro_1_tmp(:) 
     double precision, allocatable :: epsilon_tro_2_tmp(:) 
     double precision, allocatable :: epsilon_tro_3_tmp(:) 
     double precision, allocatable :: epsilon_tro_4_tmp(:) 
     double precision, allocatable :: epsilon_tro_5_tmp(:) 
     double precision, allocatable :: epsilon_cli_1_tmp(:) 
     double precision, allocatable :: epsilon_cli_2_tmp(:) 
     double precision, allocatable :: epsilon_cli_3_tmp(:) 
     double precision, allocatable :: epsilon_cli_4_tmp(:) 
     double precision, allocatable :: epsilon_cli_5_tmp(:) 
     double precision, allocatable :: epsilon_cli_6_tmp(:) 
     double precision, allocatable :: epsilon_cli_7_tmp(:) 
  !----- エネルギーフラックス   
     double precision, allocatable :: e_tro_ekflux_1(:)
     double precision, allocatable :: e_tro_ekflux_2(:)
     double precision, allocatable :: e_tro_ekflux_total(:)
     double precision, allocatable :: e_cli_ekflux_1(:)
     double precision, allocatable :: e_cli_ekflux_2(:)
     double precision, allocatable :: e_cli_ekflux_3(:)
     double precision, allocatable :: e_cli_ekflux_total(:)
  !--- コマンドライン引数用変数 ---
  type(ARGS) :: arg
  character(STRING), pointer :: argv(:) => null()
  logical             :: Verbose=.false.          ! 冗長なメッセージ出力 On/Off
  logical             :: DebugOn=.false.          ! デバッグメッセージ On/Off
  
  !--- 格子点・スペクトル ---
  integer :: km=341                               ! X 方向切断全波数
  integer :: lm=341                               ! Y 方向切断全波数
  integer :: im=1024                              ! X 方向格子点数 (>3*km)
  integer :: jm=1024                              ! Y 方向格子点数 (>3*lm)
  integer :: wn=341
  
  !--- 物理パラメータ ---
  double precision :: kd = 1.0D1
  double precision :: kappa = 4.0D-2
  double precision :: U = 2.5D-2
  double precision :: HVisc = 2.017D-19
  double precision :: delta_HVisc = 0.0D0
  double precision :: pi = acos(-1.0D0)            !3.1415926D0
  double precision :: XLength = 2.0D0*acos(-1.0D0) !6.2831853D0         ! 領域の大きさ(X 方向)
  double precision :: YLength = 2.0D0*acos(-1.0D0) !6.2831853D0         ! 領域の大きさ(Y 方向)
  integer :: HVOrder = 4
  
  double precision :: t,tm
  character(string) :: time
  
  double precision :: init = 1800.0D0
  integer :: timestep=200
  double precision :: tint = 1.0D0
  character(len=100) :: initial_file='../2-layer-cyclic.nc'          ! 初期値データファイル名

  !--- 一時変数
  integer :: k,l,n,it,i,j
  double precision, allocatable :: ee_tro_term1(:,:)
  double precision, allocatable :: ee_tro_term2(:,:)
  double precision, allocatable :: ee_tro_term3(:,:)
  double precision, allocatable :: ee_tro_term4(:,:)
  double precision, allocatable :: ee_cli_term1(:,:)
  double precision, allocatable :: ee_cli_term2(:,:)
  double precision, allocatable :: ee_cli_term3(:,:)
  double precision, allocatable :: ee_cli_term4(:,:)
  double precision, allocatable :: ee_cli_term5(:,:)
  double precision, allocatable :: ee_cli_term6(:,:)

  !--------- 変数割り付け ---------
  allocate(ee_tro_Vor(-lm:lm,-km:km),ee_tro_StrFunc(-lm:lm,-km:km))
  allocate(ee_cli_Vor(-lm:lm,-km:km),ee_cli_T(-lm:lm,-km:km),ee_cli_StrFunc(-lm:lm,-km:km))
  allocate(ee_T_psi_1(-lm:lm,-km:km),e_T_psi_1(1:wn))
  allocate(ee_T_psi_2(-lm:lm,-km:km),e_T_psi_2(1:wn))
  allocate(ee_T_psi_3(-lm:lm,-km:km),e_T_psi_3(1:wn))
  allocate(ee_T_psi_4(-lm:lm,-km:km),e_T_psi_4(1:wn))
  allocate(ee_T_psi_5(-lm:lm,-km:km),e_T_psi_5(1:wn))
  allocate(ee_T_tau_1(-lm:lm,-km:km),e_T_tau_1(1:wn))
  allocate(ee_T_tau_2(-lm:lm,-km:km),e_T_tau_2(1:wn))
  allocate(ee_T_tau_3(-lm:lm,-km:km),e_T_tau_3(1:wn))
  allocate(ee_T_tau_4(-lm:lm,-km:km),e_T_tau_4(1:wn))
  allocate(ee_T_tau_5(-lm:lm,-km:km),e_T_tau_5(1:wn))
  allocate(ee_T_tau_6(-lm:lm,-km:km),e_T_tau_6(1:wn))
  allocate(ee_T_tau_7(-lm:lm,-km:km),e_T_tau_7(1:wn))
  allocate(e_tro_ekflux_1(1:wn),e_tro_ekflux_2(1:wn))
  allocate(e_cli_ekflux_1(1:wn),e_cli_ekflux_2(1:wn),e_cli_ekflux_3(1:wn))
  allocate(e_tro_ekflux_total(1:wn),e_cli_ekflux_total(1:wn))
  allocate(ee_tro_term1(-lm:lm,-km:km),ee_tro_term2(-lm:lm,-km:km))
  allocate(ee_tro_term3(-lm:lm,-km:km),ee_tro_term4(-lm:lm,-km:km))
  allocate(ee_cli_term1(-lm:lm,-km:km),ee_cli_term2(-lm:lm,-km:km))
  allocate(ee_cli_term3(-lm:lm,-km:km),ee_cli_term4(-lm:lm,-km:km))
  allocate(ee_cli_term5(-lm:lm,-km:km),ee_cli_term6(-lm:lm,-km:km))
  allocate(epsilon_tro_1_tmp(0:timestep), epsilon_tro_2_tmp(0:timestep))
  allocate(epsilon_tro_3_tmp(0:timestep), epsilon_tro_4_tmp(0:timestep)) 
  allocate(epsilon_tro_5_tmp(0:timestep), epsilon_cli_1_tmp(0:timestep))
  allocate(epsilon_cli_2_tmp(0:timestep), epsilon_cli_3_tmp(0:timestep))
  allocate(epsilon_cli_4_tmp(0:timestep), epsilon_cli_5_tmp(0:timestep)) 
  allocate(epsilon_cli_6_tmp(0:timestep), epsilon_cli_7_tmp(0:timestep)) 
    
  !--------- 座標値の設定 ---------
  call DbgMessage(fmt='call %c', c1='ee_initial') 
  call ee_Initial(im,jm,km,lm,0.0D0,XLength,0.0D0,YLength)
  
  call output_history_init

  do it=0,timestep
     t = init + tint * it
     time = CPrintf("t=%f",d=(/t/))
     !----- 計算結果のデータから, 渦度を取得 -----
     call HistoryGet(   trim(initial_file), 'tro_vorsp', ee_tro_Vor, &
           & time=t)
     call HistoryGet(   trim(initial_file), 'cli_vorsp', ee_cli_Vor, &
           & time=t)
     ee_tro_Vor = transpose(ee_tro_Vor)
     ee_cli_Vor = transpose(ee_cli_Vor)

     !----- 流線函数などを計算 -----
     ee_tro_Strfunc = 0.0D0
     ee_tro_Strfunc = ee_LaplacianInv_ee(ee_tro_Vor)
     ee_cli_Strfunc = 0.0D0
     ee_cli_Strfunc = ee_cli_Vor2Strfunc_ee(ee_cli_Vor)
     ee_cli_T = 0.0D0
     ee_cli_T = ee_Laplacian_ee(ee_cli_Strfunc)
     
     !--------- エンストロフィー輸送函数計算 ---------
     ee_tro_term1 = ee_Jac_ee_ee(ee_tro_Strfunc,ee_tro_Vor)
     ee_tro_term2 = ee_Jac_ee_ee(ee_cli_Strfunc,ee_cli_T)
     ee_tro_term3 = U * ee_Dx_ee(ee_cli_T)
     ee_tro_term4 = 0.5D0 * kappa * (ee_Laplacian_ee(ee_tro_Strfunc) &
          & - ee_Laplacian_ee(ee_cli_Strfunc))
     ee_cli_term1 = ee_Jac_ee_ee(ee_cli_Strfunc,ee_tro_Vor)
     ee_cli_term2 = ee_Jac_ee_ee(ee_tro_Strfunc,ee_cli_T)
     ee_cli_term3 = ee_Jac_ee_ee(ee_tro_Strfunc, -kd**2*ee_cli_Strfunc)
     ee_cli_term4 = U * ee_Dx_ee(ee_tro_Vor)
     ee_cli_term5 = U * kd**2 * ee_Dx_ee(ee_tro_Strfunc)
     ee_cli_term6 = 0.5D0 * kappa * (ee_Laplacian_ee(ee_cli_Strfunc) &
          & - ee_Laplacian_ee(ee_tro_Strfunc))
     ee_T_psi_1 = 0.0D0
     ee_T_psi_2 = 0.0D0
     ee_T_psi_3 = 0.0D0
     ee_T_psi_4 = 0.0D0
     ee_T_psi_5 = 0.0D0
     ee_T_tau_1 = 0.0D0
     ee_T_tau_2 = 0.0D0
     ee_T_tau_3 = 0.0D0
     ee_T_tau_4 = 0.0D0
     ee_T_tau_5 = 0.0D0
     ee_T_tau_6 = 0.0D0
     do l = -lm,lm
        do k = -km,km
           n = nint(sqrt(real(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
              ee_T_psi_1(l,k) =  ee_tro_Strfunc(l,k) * ee_tro_term1(l,k) &
                   & + ee_tro_Strfunc(-l,-k) * ee_tro_term1(-l,-k)
              ee_T_psi_2(l,k) =  ee_tro_Strfunc(l,k) * ee_tro_term2(l,k) &
                   & + ee_tro_Strfunc(-l,-k) * ee_tro_term2(-l,-k)
              ee_T_psi_3(l,k) =  ee_tro_Strfunc(l,k) * ee_tro_term3(l,k) &
                   & + ee_tro_Strfunc(-l,-k) * ee_tro_term3(-l,-k)
              ee_T_psi_4(l,k) &
                   & =  ee_tro_Strfunc(l,k) * ee_tro_term4(l,k) &
                   &  + ee_tro_Strfunc(-l,-k) * ee_tro_term4(-l,-k)
              ee_T_psi_5(l,k) =  ee_tro_Strfunc(l,k) * HVisc &
                   & * ((2.0D0*pi*k/XLength)**2+(2.0D0*pi*l/YLength)**2)**HVOrder &
                   & * ee_tro_Vor(l,k) &
                   & +  ee_tro_Strfunc(-l,-k) * HVisc &
                   & * ((2.0D0*pi*(-k)/XLength)**2+(2.0D0*pi*(-l)/YLength)**2)**HVOrder &
                   & * ee_tro_Vor(-l,-k) &
                   & + ee_tro_Strfunc(l,k) * 0.5D0 * delta_HVisc &
                   & * ((2.0D0*pi*(k)/XLength)**2+(2.0D0*pi*(l)/YLength)**2)**HVOrder &
                   & * (ee_tro_Vor(l,k) - ee_cli_Vor(l,k)) &
                   & + ee_tro_Strfunc(-l,-k) * 0.5D0 * delta_HVisc &
                   & * ((2.0D0*pi*(-k)/XLength)**2+(2.0D0*pi*(-l)/YLength)**2)**HVOrder &
                   & * (ee_tro_Vor(-l,-k) - ee_cli_Vor(-l,-k))
              ee_T_tau_1(l,k) = ee_cli_Strfunc(l,k) * ee_cli_term1(l,k) &
                   & + ee_cli_Strfunc(-l,-k) * ee_cli_term1(-l,-k)
              ee_T_tau_2(l,k) = ee_cli_Strfunc(l,k) * ee_cli_term2(l,k) &
                   & + ee_cli_Strfunc(-l,-k) * ee_cli_term2(-l,-k)
              ee_T_tau_3(l,k) = ee_cli_Strfunc(l,k) * ee_cli_term3(l,k) &
                   & + ee_cli_Strfunc(-l,-k) * ee_cli_term3(-l,-k)
              ee_T_tau_4(l,k) = ee_cli_Strfunc(l,k) * ee_cli_term4(l,k) &
                   & + ee_cli_Strfunc(-l,-k) * ee_cli_term4(-l,-k)
              ee_T_tau_5(l,k) = ee_cli_Strfunc(l,k) * ee_cli_term5(l,k) &
                   & + ee_cli_Strfunc(-l,-k) * ee_cli_term5(-l,-k)
              ee_T_tau_6(l,k) &
                   & = ee_cli_Strfunc(l,k)*ee_cli_term6(l,k) &
                   &  + ee_cli_Strfunc(-l,-k)*ee_cli_term6(-l,-k)
              ee_T_tau_7(l,k) = HVisc * ee_cli_Strfunc(l,k) &
                   & * ((2.0D0*pi*k/XLength)**2+(2.0D0*pi*l/YLength)**2)**HVOrder &
                   & * ee_cli_Vor(l,k) &
                   & + HVisc * ee_cli_Strfunc(-l,-k) &
                   & * ((2.0D0*pi*(-k)/XLength)**2+(2.0D0*pi*(-l)/YLength)**2)**HVOrder &
                   & * ee_cli_Vor(-l,-k) &
                   & + ee_cli_Strfunc(l,k) * 0.5D0 * delta_HVisc &
                   & * ((2.0D0*pi*(k)/XLength)**2+(2.0D0*pi*(l)/YLength)**2)**HVOrder &
                   & * (ee_cli_Vor(l,k) - ee_tro_Vor(l,k)) &
                   & + ee_cli_Strfunc(-l,-k) * 0.5D0 * delta_HVisc &
                   & * ((2.0D0*pi*(-k)/XLength)**2+(2.0D0*pi*(-l)/YLength)**2)**HVOrder &
                   & * (ee_cli_Vor(-l,-k) - ee_tro_Vor(-l,-k))    
           end if
10        continue
        end do
     end do
     e_T_psi_1 = e_ee(ee_T_psi_1)
     e_T_psi_2 = e_ee(ee_T_psi_2)
     e_T_psi_3 = e_ee(ee_T_psi_3)
     e_T_psi_4 = e_ee(ee_T_psi_4)
     e_T_psi_5 = e_ee(ee_T_psi_5)
     e_T_tau_1 = e_ee(ee_T_tau_1)
     e_T_tau_2 = e_ee(ee_T_tau_2)
     e_T_tau_3 = e_ee(ee_T_tau_3)
     e_T_tau_4 = e_ee(ee_T_tau_4)
     e_T_tau_5 = e_ee(ee_T_tau_5)
     e_T_tau_6 = e_ee(ee_T_tau_6)
     e_T_tau_7 = e_ee(ee_T_tau_7)
          
     !--------- フラックス計算 ---------
      !-- barotropic mode
       e_tro_ekflux_1 = - calc_flux(e_T_psi_1)
       e_tro_ekflux_2 = - calc_flux(e_T_psi_2)
       e_tro_ekflux_total = e_tro_ekflux_1 + e_tro_ekflux_2
      !-- baroclinic mode
       e_cli_ekflux_1 = - calc_flux(e_T_tau_1)
       e_cli_ekflux_2 = - calc_flux(e_T_tau_2)
       e_cli_ekflux_3 = - calc_flux(e_T_tau_3)
       e_cli_ekflux_total = e_cli_ekflux_1 + e_cli_ekflux_2 + e_cli_ekflux_3
     !----- エンストロフィー散逸率・注入率を計算 -----
      !-- barotropic mode
       epsilon_tro_1_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_tro_1_tmp(it) = epsilon_tro_1_tmp(it) + e_T_psi_1(n)
       end do

       epsilon_tro_2_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_tro_2_tmp(it) = epsilon_tro_2_tmp(it) + e_T_psi_2(n)
       end do

       epsilon_tro_3_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_tro_3_tmp(it) = epsilon_tro_3_tmp(it) + e_T_psi_3(n)
       end do

       epsilon_tro_4_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_tro_4_tmp(it) = epsilon_tro_4_tmp(it) + e_T_psi_4(n)
       end do

       epsilon_tro_5_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_tro_5_tmp(it) = epsilon_tro_5_tmp(it) + e_T_psi_5(n)
       end do

       !-- baroclinic mode
       epsilon_cli_1_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_cli_1_tmp(it) = epsilon_cli_1_tmp(it) + e_T_tau_1(n)
       end do

       epsilon_cli_2_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_cli_2_tmp(it) = epsilon_cli_2_tmp(it) + e_T_tau_2(n)
       end do

       epsilon_cli_3_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_cli_3_tmp(it) = epsilon_cli_3_tmp(it) + e_T_tau_3(n)
       end do

       epsilon_cli_4_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_cli_4_tmp(it) = epsilon_cli_4_tmp(it) + e_T_tau_4(n)
       end do

       epsilon_cli_5_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_cli_5_tmp(it) = epsilon_cli_5_tmp(it) + e_T_tau_5(n)
       end do

       epsilon_cli_6_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_cli_6_tmp(it) = epsilon_cli_6_tmp(it) + e_T_tau_6(n)
       end do

       epsilon_cli_7_tmp(it) = 0.0D0
       do n=1,wn
          epsilon_cli_7_tmp(it) = epsilon_cli_7_tmp(it) + e_T_tau_7(n)
       end do
     !----- 計算結果出力 -----
     call output_history
  end do

  epsilon_tro_1 = sum(epsilon_tro_1_tmp) / size(epsilon_tro_1_tmp)
  epsilon_tro_2 = sum(epsilon_tro_2_tmp) / size(epsilon_tro_2_tmp)
  epsilon_tro_3 = sum(epsilon_tro_3_tmp) / size(epsilon_tro_3_tmp)
  epsilon_tro_4 = sum(epsilon_tro_4_tmp) / size(epsilon_tro_4_tmp)
  epsilon_tro_5 = sum(epsilon_tro_5_tmp) / size(epsilon_tro_5_tmp)
  epsilon_cli_1 = sum(epsilon_cli_1_tmp) / size(epsilon_cli_1_tmp)
  epsilon_cli_2 = sum(epsilon_cli_2_tmp) / size(epsilon_cli_2_tmp)
  epsilon_cli_3 = sum(epsilon_cli_3_tmp) / size(epsilon_cli_3_tmp)
  epsilon_cli_4 = sum(epsilon_cli_4_tmp) / size(epsilon_cli_4_tmp)
  epsilon_cli_5 = sum(epsilon_cli_5_tmp) / size(epsilon_cli_5_tmp)
  epsilon_cli_6 = sum(epsilon_cli_6_tmp) / size(epsilon_cli_6_tmp)
  epsilon_cli_7 = sum(epsilon_cli_7_tmp) / size(epsilon_cli_7_tmp)

  call HistoryPut('tro_epsilon1',epsilon_tro_1)
  call HistoryPut('tro_epsilon2',epsilon_tro_2)
  call HistoryPut('tro_epsilon3',epsilon_tro_3)
  call HistoryPut('tro_epsilon4',epsilon_tro_4)
  call HistoryPut('tro_epsilon5',epsilon_tro_5)
  call HistoryPut('cli_epsilon1',epsilon_cli_1)
  call HistoryPut('cli_epsilon2',epsilon_cli_2)
  call HistoryPut('cli_epsilon3',epsilon_cli_3)
  call HistoryPut('cli_epsilon4',epsilon_cli_4)
  call HistoryPut('cli_epsilon5',epsilon_cli_5)
  call HistoryPut('cli_epsilon6',epsilon_cli_6)
  call HistoryPut('cli_epsilon7',epsilon_cli_7)
  
  call HistoryClose
  
! 以上 メインプログラム 
!========================================================================
! 以下 サブルーチン

contains
  function e_ee(ee)
  !----------------------------------------------------------------------------
  ! 2 次元スペクトルから 1 次元スペクトルを計算
  !----------------------------------------------------------------------------
    double precision                :: e_ee(1:wn)
    double precision, intent(IN)    :: ee(-lm:lm,-km:km)

    e_ee = 0.0D0
    do k=-km,km
       do l=-lm,lm
          n = nint(sqrt(real(l**2+k**2)))
          if ((l .eq. 0) .AND. (k .eq. 0)) then
             goto 210
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 210
          else
             e_ee(n) = e_ee(n) + ee(l,k)
          endif
210       continue
       enddo
    enddo
  end function e_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(real(l**2+k**2)))
          if (n .eq. 0) then
             goto 220
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 220
          else
             ee_Laplacian_ee(l,k) = - ee(l,k) * dble((2.0D0*pi*k/XLength)**2+(2.0D0*pi*l/YLength)**2)
          endif
220       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(real(l**2+k**2)))
          if (n .eq. 0) then
             goto 230
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 230
          else
             ee_LaplacianInv_ee(l,k) = - ee(l,k) / dble((2.0D0*pi*k/XLength)**2+(2.0D0*pi*l/YLength)**2)
          endif
230       continue
       end do
    end do
  end function ee_LaplacianInv_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

    do k=-km,km
       do l=-lm,lm
          n = nint(sqrt(real(l**2+k**2)))
          if ((l .eq. 0) .AND. (k .eq. 0)) then
             ee_cli_Vor2Strfunc_ee(l,k) = 0.0D0
          else if ((n .gt. km) .OR. (n .gt. km)) then
             goto 130
          else
             ee_cli_Vor2Strfunc_ee(l,k) &
                  & = - ee_cli_Vor(l,k) / ((2.0D0*pi*k/XLength)**2+(2.0D0*pi*l/YLength)**2+kd**2)
          end if
130       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(real(l**2+k**2)))
          if (n .eq. 0) then
             goto 300
          else if ((n .gt. km) .OR. (n .gt. lm)) then
             goto 300
          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
300       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

  function calc_flux(e)
  !----------------------------------------------------------------------------
  ! 輸送函数からフラックス計算
  !----------------------------------------------------------------------------
    double precision, intent(IN) :: e(1:wn)
    double precision             :: calc_flux(0:wn)
    ! ----- エンストロフィー輸送関数の計算 -----
    calc_flux(0) = 0.0D0
    ! ----- エンストロフィーフラックス計算 -----
    do n=1,wn
       calc_flux(n) = calc_flux(n-1) + e(n)
    enddo
  end function calc_flux
  
!--------- ヒストリー出力 ---------
  subroutine output_history_init
    call HistoryCreate( &
         file="ek-analysis.nc", &
         title="two-layer turbulence energy analysis", &
         source='2d-layer-cyclic.f90 (2016/01/05)', &
         institution='OKAZAKI Shogo',     &
         dims=(/'k','l','n','t'/), &
         dimsizes=(/2*km+1,2*lm+1,wn,0/),&
         longnames=(/'X-wavenumber','Y-wavenumber',&
                     'wavenumber  ','time        '/),&
         units=(/'1','1','1','1'/),   &
         origin=real(init), interval=real(tint) )

    !---- 座標変数定義, 出力 ----
    call HistoryPut('k',(/(dble(k),k=-km,km)/) )        ! 変数出力
    call HistoryPut('l',(/(dble(l),l=-lm,lm)/) )        ! 変数出力    
    call HistoryPut('n',(/(dble(n),n=1,wn)/) )

    !---- パラメター定義, 出力 ----
     !-- エンストロフィー輸送函数
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_T1', dims=(/'n','t'/), & 
           longname='barotropic energy transfer function 1', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_T2', dims=(/'n','t'/), & 
           longname='barotropic energy transfer function 2', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_T3', dims=(/'n','t'/), & 
           longname='barotropic energy transfer function 3', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_T4', dims=(/'n','t'/), & 
           longname='barotropic energy transfer function 4', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_T5', dims=(/'n','t'/), & 
           longname='barotropic energy transfer function 5', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_T1', dims=(/'n','t'/), & 
           longname='baroclinic energy transfer function 1', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_T2', dims=(/'n','t'/), & 
           longname='baroclinic energy transfer function 2', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_T3', dims=(/'n','t'/), & 
           longname='baroclinic energy transfer function 3', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_T4', dims=(/'n','t'/), & 
           longname='baroclinic energy transfer function 4', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_T5', dims=(/'n','t'/), & 
           longname='baroclinic energy transfer function 5', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_T6', dims=(/'n','t'/), & 
           longname='baroclinic energy transfer function 6', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_T7', dims=(/'n','t'/), & 
           longname='baroclinic energy transfer function 7', units='1', xtype='double')
     !-- エンストロフィー散逸率
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_epsilon1', dims=(/''/), & 
           longname='barotropic energy dissipation rate 1', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_epsilon2', dims=(/''/), & 
           longname='barotropic energy dissipation rate 2', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_epsilon3', dims=(/''/), & 
           longname='barotropic energy dissipation rate 3', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_epsilon4', dims=(/''/), & 
           longname='barotropic energy dissipation rate 4', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_epsilon5', dims=(/''/), & 
           longname='barotropic energy dissipation rate 5', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_epsilon1', dims=(/''/), & 
           longname='baroclinic energy dissipation rate 1', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_epsilon2', dims=(/''/), & 
           longname='baroclinic energy dissipation rate 2', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_epsilon3', dims=(/''/), & 
           longname='baroclinic energy dissipation rate 3', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_epsilon4', dims=(/''/), & 
           longname='baroclinic energy dissipation rate 4', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_epsilon5', dims=(/''/), & 
           longname='baroclinic energy dissipation rate 5', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_epsilon6', dims=(/''/), & 
           longname='baroclinic energy dissipation rate 6', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_epsilon7', dims=(/''/), & 
           longname='baroclinic energy dissipation rate 7', units='1', xtype='double')
     !-- エンストロフィーフラックス
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_flux1', dims=(/'n','t'/), & 
           longname='barotropic energy flux 1', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_flux2', dims=(/'n','t'/), & 
           longname='barotropic energy flux 2', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='tro_flux_total', dims=(/'n','t'/), & 
           longname='barotropic total energy flux', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_flux1', dims=(/'n','t'/), & 
           longname='baroclinic energy flux 1', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_flux2', dims=(/'n','t'/), & 
           longname='baroclinic energy flux 2', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_flux3', dims=(/'n','t'/), & 
           longname='baroclinic energy flux 3', units='1', xtype='double')
      call HistoryAddVariable( &                                ! 変数定義
           varname='cli_flux_total', dims=(/'n','t'/), & 
           longname='baroclinic total energy flux', units='1', xtype='double')
      
  end subroutine output_history_init

  subroutine output_history
  !------------------------------------------------------
  ! ヒストリー出力
  !------------------------------------------------------
    write(6,*) ' History file output at it = ',it, '  time = ', t
    call HistoryPut('t',real(t))
    call HistoryPut('tro_T1',e_T_psi_1)
    call HistoryPut('tro_T2',e_T_psi_2)
    call HistoryPut('tro_T3',e_T_psi_3)
    call HistoryPut('tro_T4',e_T_psi_4)
    call HistoryPut('tro_T5',e_T_psi_5)
    call HistoryPut('cli_T1',e_T_tau_1)
    call HistoryPut('cli_T2',e_T_tau_2)
    call HistoryPut('cli_T3',e_T_tau_3)
    call HistoryPut('cli_T4',e_T_tau_4)
    call HistoryPut('cli_T5',e_T_tau_5)
    call HistoryPut('cli_T6',e_T_tau_6)
    call HistoryPut('cli_T7',e_T_tau_7)
    call HistoryPut('tro_flux1',e_tro_ekflux_1)
    call HistoryPut('tro_flux2',e_tro_ekflux_2)
    call HistoryPut('tro_flux_total',e_tro_ekflux_total)
    call HistoryPut('cli_flux1',e_cli_ekflux_1)
    call HistoryPut('cli_flux2',e_cli_ekflux_2)
    call HistoryPut('cli_flux3',e_cli_ekflux_3)
    call HistoryPut('cli_flux_total',e_cli_ekflux_total)
  end subroutine output_history

end program cal_ekterm
