!
! 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=170                                ! X æ–¹å‘åˆ‡æ–å…¨æ³¢æ•°
  integer :: lm=170                                ! Y æ–¹å‘åˆ‡æ–å…¨æ³¢æ•°
  integer :: im=512                               ! X æ–¹å‘æ ¼åç‚¹æ•° (>3*km)
  integer :: jm=512                               ! Y æ–¹å‘æ ¼åç‚¹æ•° (>3*lm)
  integer :: wn=170

  !--- ç‰©ç†ãƒ‘ãƒ©ãƒ¡ãƒ¼ã‚¿ ---
  double precision :: kd = 1.0D1
  double precision :: kappa = 4.0D-2
  double precision :: U = 2.5D-2
  double precision :: HVisc = 5.280D-17
  double precision :: delta_HVisc = -5.280D-17
  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 = 1300.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