!== エネルギースペクトル分布導出プログラム
!
! Authors::   Takashi TSUCHIYA, Yasuhiro MORIKAWA
! Version::   $Id:  $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2007. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
! netCDF に存在する渦度の実空間データより
! エネルギースペクトル分布を求めるプログラム
!
! メモ: 時間ループを行う際のタイムステップを netCDF データから
!       取得できないだろうか？
!       → 現在は無理らしい
!
!       netCDF データの渦度は有次元化されているので,エネルギー
!       スペクトルを計算するにあたり渦度を無次元化する必要がある.
!       渦度は U_scale / Rplanet で有次元化する.
!       これも netCDF データから取得できない
!       →
!
program translate_spectrum

!== 宣言文 =============================================

  use gt4_history
  use dc_string, only: StoA, toChar
  use w_module,  only: l_nm, w_Lapla_w
  use wa_module, only: xy_w, w_xy, wa_Initial
  use dc_types, only: DP, STRING
  use dc_string, only: StoA
  use dc_args, only: ARGS, Open, HelpMsg, Option, Debug, Help, Strict, Close
  use dc_message, only: MessageNotify

  !---------------------------------------------------------
  !  実験の表題, モデルの名称, 所属機関名
  !  Title of a experiment, name of model, sub-organ
  !---------------------------------------------------------
  character(*), parameter:: title = &
    & 'energy_spectrum ver4 $Name:  $ :: '
  character(*), parameter:: source = &
    & 'dcpam4 ' // &
    & '(See http://www.gfd-dennou.org/library/dcpam)'
  character(*), parameter:: institution = &
    & 'GFD Dennou Club (See http://www.gfd-dennou.org)'

  !---------------------------------------------------------
  !  作業変数
  !  Work variables
  !---------------------------------------------------------
  type(ARGS):: arg            ! コマンドライン引数. 
                              ! Command line arguments
  logical:: OPT_input         ! -I, --input オプションの有無. 
                              ! Existence of '-I', '--input' option
  character(STRING):: VAL_input
                              ! -I, --input オプションの値. 
                              ! Value of '-I', '--input' option
  logical:: OPT_output        ! -O, --output オプションの有無. 
                              ! Existence of '-O', '--output' option
  character(STRING):: VAL_output
                              ! -O, --output オプションの値. 
                              ! Value of '-O', '--output' option

  character(STRING):: input_filename = 'Vor.nc'
  character(STRING):: input_filename2 = 'Div.nc'
  character(STRING):: input_filename3 = 'U.nc'
  character(STRING):: input_filename4 = 'V.nc'
  character(STRING):: input_filename5 = '/home/tutaka/tutaka_weight/weight_T170.nc'
  character(STRING):: input_filename6 = 'Ps.nc'
  character(STRING):: input_filename7 = 'Temp.nc'
  character(STRING):: output_filename = 'ESp.nc'

  ! 無次元化に必要な変数
  real(8) :: U_scale = 10D0
  real(8) :: RPlanet = 6.371E6

  ! 作業用変数
  integer :: i, t, k, n

  ! dimensions
  real(8), pointer  :: lon(:)  => null()
  real(8), pointer  :: lat(:)  => null()
  real(8), pointer  :: sig(:)  => null()
  integer, pointer  :: time(:) => null()
  real(8), pointer  :: lonw(:)  => null()
  real(8), pointer  :: latw(:)  => null()
!  real(8), pointer  :: lon_weight(:)  => null()
!  real(8), pointer  :: coslat_lat_weight(:)  => null()
  ! variables
  real(8), pointer  :: xy_Vor(:,:) => null()
  real(8), pointer  :: xy_Div(:,:) => null()
  real(8), pointer  :: xy_U(:,:) => null()
  real(8), pointer  :: xy_V(:,:) => null()
  real(8), pointer  :: xy_Ps(:,:) => null()
  real(8), pointer  :: xy_Temp(:,:) => null()


  ! 条件
  integer :: nmax ! 最大全波数
  integer :: imax ! 経度格子点数
  integer :: jmax ! 緯度格子点数
  integer :: kmax ! 鉛直層数

  real(8), allocatable  :: xy_Vor_bt(:,:)    ! 順圧渦度(2次元)
  real(8), allocatable  :: xy_Vor_bc(:,:)    ! 傾圧渦度(2次元)
  real(8), allocatable  :: xy_Vor_bt_n(:,:)  ! 無次元順圧渦度(2次元)
  real(8), allocatable  :: xy_Vor_bc_n(:,:)  ! 無次元傾圧渦度(2次元)
  real(8), allocatable  :: xy_Div_bc(:,:)    ! 傾圧発散(2次元)
  real(8), allocatable  :: xy_Div_bt(:,:)    ! 順圧発散(2次元)
  real(8), allocatable  :: xyz_Vor(:,:,:)    ! 渦度(3次元)
  real(8), allocatable  :: xyz_Vor_bt(:,:,:) ! 順圧渦度(3次元)
  real(8), allocatable  :: xyz_Vor_bc(:,:,:) ! 傾圧渦度(3次元)
  real(8), allocatable  :: xy_Div_bc_n(:,:)  ! 無次元傾圧発散(2次元)
  real(8), allocatable  :: xy_Div_bt_n(:,:)  ! 無次元順圧発散(2次元)
  real(8), allocatable  :: xyz_Div(:,:,:)    ! 発散(3次元)
  real(8), allocatable  :: xyz_Div_bc(:,:,:) ! 傾圧発散(3次元)
  real(8), allocatable  :: xyz_Div_bt(:,:,:) ! 順圧発散(3次元)
  real(8), allocatable  :: xy_U_bt(:,:)      ! 東西速度(2次元)
  real(8), allocatable  :: xy_U_bc(:,:)      ! 東西速度(2次元)
  real(8), allocatable  :: xyz_U(:,:,:)      ! 東西速度(3次元)
  real(8), allocatable  :: xyz_V(:,:,:)      ! 南北速度(3次元)
  real(8), allocatable  :: xyz_E(:,:,:)      ! エネルギー(3次元)
  real(8), allocatable  :: we_xyz_E(:,:,:)    ! 重みをつけたエネルギー(3次元)
  real(8), allocatable  :: we_xyz_U(:,:,:)    ! 重みをつけた東西速度(3次元)
  real(8), allocatable  :: we_xyz_V(:,:,:)    ! 重みをつけた南北速度(3次元)
  real(8), allocatable  :: we_xyz_Vor(:,:,:)  ! 重みをつけた渦度(3次元)
  real(8), allocatable  :: w_Vor(:)          ! 渦度 (スペクトルデータ)
  real(8), allocatable  :: w_Vor_bt(:)       ! 順圧渦度 (スペクトルデータ)
  real(8), allocatable  :: w_Vor_bc(:)       ! 傾圧渦度 (スペクトルデータ)
  real(8), allocatable  :: w_Div_bc(:)       ! 傾圧発散 (スペクトルデータ)
  real(8), allocatable  :: w_Div_bt(:)       ! 順圧発散 (スペクトルデータ)
  real(8), allocatable  :: ESp_btv(:)        ! 順圧渦度のエネルギースペクトル
  real(8), allocatable  :: ESp_bcv(:)        ! 傾圧渦度のエネルギースペクトル
  real(8), allocatable  :: ESp_bcd(:)        ! 傾圧発散のエネルギースペクトル
  real(8), allocatable  :: ESp_btd(:)        ! 順圧発散のエネルギースペクトル
  integer, allocatable  :: nn(:)             ! 全波数
  real(8), allocatable  :: xyz_Ens(:,:,:)    ! エンストロフィー(3次元)
  real(8), allocatable  :: xy_Vor_up_n(:,:)  ! 無次元上層渦度
  real(8), allocatable  :: xy_Vor_low_n(:,:) ! 無次元下層渦度
  real(8), allocatable  :: w_Vor_up(:)       ! 上層渦度 (スペクトルデータ)
  real(8), allocatable  :: w_Vor_low(:)      ! 下層渦度 (スペクトルデータ)
  real(8), allocatable  :: ESp_vup(:)        ! 傾圧発散のエネルギースペクトル
  real(8), allocatable  :: ESp_vlow(:)       ! 傾圧発散のエネルギースペクトル
  real(8), allocatable  :: K_Energy(:)       ! 運動エネルギー
  real(8), allocatable  :: we_xy_Ps(:,:)      ! 重みをつけた地表面気圧(2次元)
  real(8), allocatable  :: xyz_Temp(:,:,:)   ! 温度(3次元)
  real(8), allocatable  :: we_xyz_Temp(:,:,:) ! 重みをつけた温度(3次元)
  real(8), allocatable  :: Temp_mean(:)      ! 温度の全球平均
  real(8), allocatable  :: Sk_n(:)           ! 渦度歪度
  real(8), allocatable  :: Sk_s(:)           ! 渦度歪度
  real(8), allocatable  :: Sk_n_a(:)         ! 渦度歪度
  real(8), allocatable  :: Sk_n_b(:)         ! 渦度歪度
  real(8), allocatable  :: Sk_s_a(:)         ! 渦度歪度
  real(8), allocatable  :: Sk_s_b(:)         ! 渦度歪度


  real(8)  :: Energy            ! エネルギー
  real(8)  :: Energy_btv        ! 順圧渦度エネルギー
  real(8)  :: Energy_bcv        ! 傾圧渦度エネルギー
  real(8)  :: Energy_bcd        ! 傾圧発散エネルギー
  real(8)  :: Energy_btd        ! 順圧発散エネルギー
  real(8)  :: Bcd_vup           ! 傾圧発散と上層渦度の比
  real(8)  :: Bcd_vlow          ! 傾圧発散と下層渦度の比
  real(8)  :: Ps_mean           ! 地表面気圧の全球平均

  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: d $'
  character(*), parameter:: subname = 'energy_spectrum'

continue

!== 実行文 =============================================

  !---------------------------------------------------------
  !  コマンドライン引数の処理
  !  Command line arguments handling
  !---------------------------------------------------------
  call Open( arg )
  call HelpMsg( arg, 'Title', title )
  call HelpMsg( arg, 'Usage', &
    & './energy_spectrum [Options]' )
  call HelpMsg( arg, 'Source', source )
  call HelpMsg( arg, 'Institution', institution )
  call Option( arg, StoA('-I', '--input'), &
    & OPT_input, VAL_input, &
    & help = 'Input filename (default is ' // trim(input_filename) // ') ')
  call Option( arg, StoA('-O', '--output'), &
    & OPT_output, VAL_output, &
    & help = 'Output filename (default is ' // trim(output_filename) // ')' )
  call Debug( arg ) ; call Help( arg ) ; call Strict( arg, severe = .true. )
  call Close( arg )

  ! ファイル名の設定
  if ( trim(VAL_input) /= '' ) input_filename = VAL_input
  if ( trim(VAL_output) /= '' ) output_filename = VAL_output

  ! 初期の時間, タイムステップ (暫定)
  tinit = 0.0
  history_interval_time = 1.0

  ! HistoryGetPointer
  call HistoryGetPointer(input_filename, 'lon', lon)
  call HistoryGetPointer(input_filename, 'lat', lat)
  call HistoryGetPointer(input_filename, 'sig', sig)
  call HistoryGetPointer(input_filename, 'time', time)
  call HistoryGetPointer(input_filename5, 'lon_weight', lonw)
  call HistoryGetPointer(input_filename5, 'coslat_lat_weight', latw)

  ! 格子点数, 層数, 波数を netCDF データから取得
  imax = size(lon)    ! 経度格子点数
  jmax = size(lat)    ! 緯度格子点数
  kmax = size(sig)    ! 鉛直層数
  nmax = imax / 3     ! 最大全波数 (小数点切り捨て)

  call MessageNotify('M', subname, &
    & 'Grid points are (%d, %d, %d)', &
    & i = (/imax, jmax, kmax/) )
  call MessageNotify('M', subname, &
    & 'Maximum truncated wavenumber is (%d)', &
    & i = (/nmax/) )

  allocate( w_Vor((nmax+1)*(nmax+1)) )
  allocate( w_Vor_bt((nmax+1)*(nmax+1)) )
  allocate( w_Vor_bc((nmax+1)*(nmax+1)) )
  allocate( w_Div_bc((nmax+1)*(nmax+1)) )
  allocate( w_Div_bt((nmax+1)*(nmax+1)) )
  allocate( xy_Vor_bt(1:imax,1:jmax) )
  allocate( xy_Vor_bc(1:imax,1:jmax) )
  allocate( xy_Vor_bt_n(1:imax,1:jmax) )
  allocate( xy_Vor_bc_n(1:imax,1:jmax) )
  allocate( xyz_Vor(1:imax,1:jmax,1:kmax) )
  allocate( xyz_Vor_bt(1:imax,1:jmax,1:kmax) )
  allocate( xyz_Vor_bc(1:imax,1:jmax,1:kmax) )
  allocate( xy_Div_bc(1:imax,1:jmax) )
  allocate( xy_Div_bc_n(1:imax,1:jmax) )
  allocate( xy_Div_bt(1:imax,1:jmax) )
  allocate( xy_Div_bt_n(1:imax,1:jmax) )
  allocate( xyz_Div(1:imax,1:jmax,1:kmax) )
  allocate( xyz_Div_bc(1:imax,1:jmax,1:kmax) )
  allocate( xyz_Div_bt(1:imax,1:jmax,1:kmax) )
  allocate( xy_U(1:imax,1:jmax) )
  allocate( xy_U_bt(1:imax,1:jmax) )
  allocate( xy_U_bc(1:imax,1:jmax) )
  allocate( xyz_U(1:imax,1:jmax,1:kmax) )
  allocate( xyz_V(1:imax,1:jmax,1:kmax) )
  allocate( xyz_E(1:imax,1:jmax,1:kmax) )
  allocate( we_xyz_E(1:imax,1:jmax,1:kmax) )
  allocate( we_xyz_U(1:imax,1:jmax,1:kmax) )
  allocate( we_xyz_V(1:imax,1:jmax,1:kmax) )
  allocate( we_xyz_Vor(1:imax,1:jmax,1:kmax) )
  allocate( ESp_btv(1:nmax) )
  allocate( ESp_bcv(1:nmax) )
  allocate( ESp_bcd(1:nmax) )
  allocate( ESp_btd(1:nmax) )
  allocate( nn(1:nmax) )
  allocate( xyz_Ens(1:imax,1:jmax,1:kmax) )
  allocate( xy_Vor_up_n(1:imax,1:jmax) )
  allocate( xy_Vor_low_n(1:imax,1:jmax) )
  allocate( w_Vor_up((nmax+1)*(nmax+1)) )
  allocate( w_Vor_low((nmax+1)*(nmax+1)) )
  allocate( ESp_vup(1:nmax) )
  allocate( ESp_vlow(1:nmax) )
  allocate( K_Energy(1:kmax) )
  allocate( we_xy_Ps(1:imax,1:jmax) )
  allocate( xyz_Temp(1:imax,1:jmax,1:kmax) )
  allocate( we_xyz_Temp(1:imax,1:jmax,1:kmax) )
  allocate( Temp_mean(1:kmax) )
  allocate( Sk_n(1:kmax) )
  allocate( Sk_s(1:kmax) )
  allocate( Sk_n_a(1:kmax) )
  allocate( Sk_n_b(1:kmax) )
  allocate( Sk_s_a(1:kmax) )
  allocate( Sk_s_b(1:kmax) )

  nn = (/(i,i=1,nmax)/)

! HistoryCreate
   call HistoryCreate( &
    & file=output_filename, title='energy_spectrum_and_etc', &
    & source='source', institution='institution', &
    & dims=(/'nn  ','lon ','lat ','sig ','time','lonw','latw','slat'/), &
    & dimsizes=(/nmax,imax,jmax,kmax,0,imax,jmax,jmax/), &
    & longnames=(/'total_wavunumber', &
    &             'longitude       ', &
    &             'latitude        ', &
    &             'sigma           ', &
    &             'time            ', &
    &             'lonw            ', &
    &             'latw            ', &
    &             'sine_latitude   '/), &
    & units=StoA('1','degree_east','degree_north','1', &
    &            'history_interval_unit','degree_east-w', &
    &            'degree_north-w','degree_north'), &
    & origin=real(tinit + history_interval_time), &
    & interval=real(history_interval_time) )

! 次元変数出力
  call HistoryPut('nn', nn)
  call HistoryPut('lon', lon)
  call HistoryPut('lat', lat)
  call HistoryPut('slat', sin(lat/180*3.14159265))
  call HistoryPut('sig', sig)
  call HistoryPut('lonw', lonw)
  call HistoryPut('latw', latw)

! 変数定義
  call HistoryAddVariable( &
    & varname = 'Vor_bt', &
    & dims = (/'lon ', 'lat ', 'time'/), &
    & longname = 'barotoropic_vorticity', &
    & units = '1/s', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Vor_bc', &
    & dims = (/'lon ', 'lat ', 'time'/), &
    & longname = 'baroclinic_vorticity', &
    & units = '1/s', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Div_bc', &
    & dims = (/'lon ', 'lat ', 'time'/), &
    & longname = 'baroclinic_diffusion', &
    & units = '1/s', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Div_bt', &
    & dims = (/'lon ', 'lat ', 'time'/), &
    & longname = 'barotropic_diffusion', &
    & units = '1/s', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'U_bt', &
    & dims = (/'lon ', 'slat', 'time'/), &
    & longname = 'barotropic_mean_lon_velocity', &
    & units = 'm/s', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'U_bc', &
    & dims = (/'lon ', 'slat', 'time'/), &
    & longname = 'baroclinic_mean_lon_velocity', &
    & units = 'm/s', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'ESp_btv', &
    & dims = (/'nn  ', 'time'/), &
    & longname = 'barotoropic_vorticity_energy_spectrum', &
    & units = '1/s^2', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'ESp_bcv', &
    & dims = (/'nn  ', 'time'/), &
    & longname = 'baroclinic_vorticity_energy_spectrum', &
    & units = '1/s^2', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'ESp_bcd', &
    & dims = (/'nn  ', 'time'/), &
    & longname = 'baroclinic_divergence_energy_spectrum', &
    & units = '1/s^2', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'K_Energy', &
    & dims = (/'sig','time'/), &
    & longname = 'kinematic_energy', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Energy', &
    & dims = (/'time'/), &
    & longname = 'energy', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Energy_btv', &
    & dims = (/'time'/), &
    & longname = 'barotoropic_vorticity_energy', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Energy_bcv', &
    & dims = (/'time'/), &
    & longname = 'baroclinic_vorticity_energy', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Energy_bcd', &
    & dims = (/'time'/), &
    & longname = 'baroclinic_divergence_energy', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Energy_btd', &
    & dims = (/'time'/), &
    & longname = 'barotropic_divergence_energy', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Bcd_vup', &
    & dims = (/'time'/), &
    & longname = 'divergene/vorticity', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Bcd_vlow', &
    & dims = (/'time'/), &
    & longname = 'divergence/vorticity', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Ps_mean', &
    & dims = (/'time'/), &
    & longname = 'geopotential_mean', &
    & units = 'Ps', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Temp_mean', &
    & dims = (/'sig','time'/), &
    & longname = 'temperature_mean', &
    & units = 'K', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Sk_n', &
    & dims = (/'sig','time'/), &
    & longname = 'north_vorticity_skewness', &
    & units = '1', xtype = 'double' )
  call HistoryAddVariable( &
    & varname = 'Sk_s', &
    & dims = (/'sig','time'/), &
    & longname = 'south_vorticity_skewness', &
    & units = '1', xtype = 'double' )

!  call HistoryAddVariable( &
!    & varname = 'Ens', &
!    & dims = (/'lon ', 'lat ', 'sig ', 'time'/), &
!    & longname = 'enstrophy', &
!    & units = '1/s^2', xtype = 'double' )

  call wa_Initial(nmax, imax, jmax, kmax)

  end_time = time(1) / history_interval_time

  do t=1,int(end_time)

    z=history_interval_time*t

    call HistorySetTime(real(z))

    do k=1,kmax

      ! HistoryGetPointer
      call HistoryGetPointer( &
        & input_filename, 'Vor', &                    ! (in)
        & xy_Vor, &                                   ! (out)
        & range = 'sig=' // trim(toChar(sig(k))) // &
        &         ',time=' // trim(toChar(z)), &      ! (in)
        & quiet = .true. )                            ! (in)
      ! HistoryGetPointer
      call HistoryGetPointer( &
        & input_filename2, 'Div', &                   ! (in)
        & xy_Div, &                                   ! (out)
        & range = 'sig=' // trim(toChar(sig(k))) // &
        &         ',time=' // trim(toChar(z)), &      ! (in)
        & quiet = .true. )                            ! (in)
      ! HistoryGetPointer
      call HistoryGetPointer( &
        & input_filename3, 'U', &                     ! (in)
        & xy_U, &                                     ! (out)
        & range = 'sig=' // trim(toChar(sig(k))) // &
        &         ',time=' // trim(toChar(z)), &      ! (in)
        & quiet = .true. )                            ! (in)
      ! HistoryGetPointer
      call HistoryGetPointer( &
        & input_filename4, 'V', &                     ! (in)
        & xy_V, &                                     ! (out)
        & range = 'sig=' // trim(toChar(sig(k))) // &
        &         ',time=' // trim(toChar(z)), &      ! (in)
        & quiet = .true. )                            ! (in)
      call HistoryGetPointer( &
        & input_filename7, 'Temp', &                     ! (in)
        & xy_Temp, &                                     ! (out)
        & range = 'sig=' // trim(toChar(sig(k))) // &
        &         ',time=' // trim(toChar(z)), &      ! (in)
        & quiet = .true. )                            ! (in)

      xyz_Vor(:,:,k) = xy_Vor
      xyz_Div(:,:,k) = xy_Div
      xyz_U(:,:,k) = xy_U
      xyz_V(:,:,k) = xy_V
      xyz_Temp(:,:,k) = xy_Temp

      ! エネルギーを計算する
      do j=1,jmax
        do i=1,imax
          xyz_E(i,j,k) = 0.5 * ( xyz_U(i,j,k)**2 + xyz_V(i,j,k)**2 )
        enddo
      enddo

      ! エネルギー, 温度, 渦度に重みをつける
      do j=1,jmax
        do i=1,imax
          we_xyz_E(i,j,k) = xyz_E(i,j,k) * lonw(i) * latw(j)
          we_xyz_Temp(i,j,k) = xyz_Temp(i,j,k) * lonw(i) * latw(j)
          we_xyz_Vor(i,j,k) = xyz_Vor(i,j,k) * lonw(i) * latw(j)
        enddo
      enddo
      we_xyz_E(:,:,:) = we_xyz_E(:,:,:) / (2*2*3.1415936)
      we_xyz_Temp(:,:,:) = we_xyz_Temp(:,:,:) / (2*2*3.1415936)
        ! w_xyz_Vor を元の値の桁まで戻す
        we_xyz_Vor = we_xyz_Vor * imax * jmax / (2*2*3.1415936)

      ! 重みがついたエネルギー, 温度を足し上げる
      K_Energy(k) = 0
      Temp_mean(k) = 0
      do j=1,jmax
        do i=1,imax
          K_Energy(k) = K_Energy(k) + we_xyz_E(i,j,k)
          Temp_mean(k) = Temp_mean(k) + we_xyz_Temp(i,j,k)
        enddo
      enddo
      K_Energy(k) = K_Energy(K)/(U_scale**2)

     ! 渦度歪度を求める (North, South)
     Sk_n(k) = 0
     Sk_s(k) = 0
     Sk_n_a(k) = 0
     Sk_s_a(k) = 0
     Sk_n_b(k) = 0
     Sk_s_b(k) = 0
     do i=1,imax
       do j=1,int(jmax/2)
         Sk_s_a(k) = Sk_s_a(k) + (we_xyz_Vor(i,j,k)*RPlanet/U_scale)**3
         Sk_s_b(k) = Sk_s_b(k) + (we_xyz_Vor(i,j,k)*RPlanet/U_scale)**2
       enddo
       do j=int(jmax/2+1),jmax
         Sk_n_a(k) = Sk_n_a(k) + (we_xyz_Vor(i,j,k)*RPlanet/U_scale)**3
         Sk_n_b(k) = Sk_n_b(k) + (we_xyz_Vor(i,j,k)*RPlanet/U_scale)**2
       enddo
     enddo
     Sk_n(k) = Sk_n_a(k) / Sk_n_b(k)**(3/2)
     Sk_s(k) = Sk_s_a(k) / Sk_s_b(k)**(3/2)

!        enddo
       ! エンストロフィーを計算
!      do i=1,imax
!        do j=1,jmax
!          xyz_Ens(i,j,k) = xyz_Vor(i,j,k)*xyz_Vor(i,j,k)/2
!        enddo
!      enddo

    ! k loop 終了
    enddo

    ! HistoryGetPointer
    call HistoryGetPointer( &
      & input_filename6, 'Ps', &                     ! (in)
      & xy_Ps, &                                     ! (out)
      & range = ',time=' // trim(toChar(z)), &      ! (in)
      & quiet = .true. )                            ! (in)

      ! Ps に重みをつける
      do j=1,jmax
        do i=1,imax
          we_xy_Ps(i,j) = xy_Ps(i,j) * lonw(i) * latw(j)
        enddo
      enddo

      Ps_mean = 0
      do j=1,jmax
        do i=1,imax
          Ps_mean = Ps_mean + we_xy_Ps(i,j)
        enddo
      enddo
      Ps_mean = Ps_mean / (2*2*3.1415936)

    ! 順圧渦度、傾圧渦度、順圧発散、傾圧発散を求める
    ! sig(1)=0.81 sig(2)=0.26
    xy_Vor_bt(:,:) = ( xyz_Vor(:,:,2) + xyz_Vor(:,:,1) )/2
    xy_Vor_bc(:,:) = ( xyz_Vor(:,:,2) - xyz_Vor(:,:,1) )/2
    xy_Div_bt(:,:) = ( xyz_Div(:,:,2) + xyz_Div(:,:,1) )/2
    xy_Div_bc(:,:) = ( xyz_Div(:,:,2) - xyz_Div(:,:,1) )/2
    xy_U_bt(:,:) = ( xyz_U(:,:,2) + xyz_U(:,:,1) )/2
    xy_U_bc(:,:) = ( xyz_U(:,:,2) - xyz_U(:,:,1) )/2

    ! 順圧渦度、傾圧渦度、順圧発散、傾圧発散、上層渦度、下層渦度を無次元化する
    xy_Vor_bt_n(:,:) = xy_Vor_bt * RPlanet / U_scale
    xy_Vor_bc_n(:,:) = xy_Vor_bc * RPlanet / U_scale
    xy_Div_bt_n(:,:) = xy_Div_bt * RPlanet / U_scale
    xy_Div_bc_n(:,:) = xy_Div_bc * RPlanet / U_scale
    xy_Vor_up_n(:,:) = xyz_Vor(:,:,2) * RPlanet / U_scale
    xy_Vor_low_n(:,:) = xyz_Vor(:,:,1) * RPlanet / U_scale

    ! 無次元の順圧渦度、傾圧渦度、順圧発散、傾圧発散、上層渦度、下層渦度を
    ! 格子データからスペクトルデータへ(正)変換する
    w_Vor_bt = w_xy(xy_Vor_bt_n)
    w_Vor_bc = w_xy(xy_Vor_bc_n)
    w_Div_bt = w_xy(xy_Div_bt_n)
    w_Div_bc = w_xy(xy_Div_bc_n)
    w_Vor_up = w_xy(xy_Vor_up_n)
    w_Vor_low = w_xy(xy_Vor_low_n)

    ! エネルギースペクトルを計算する
    do n=1,nmax
      ESp_btv(n) = 0.5 / (n*(n+1)) * sum( w_Vor_bt(l_nm(n,(/(m,m=-n,n)/)))**2 )
      ESp_bcv(n) = 0.5 / (n*(n+1)) * sum( w_Vor_bc(l_nm(n,(/(m,m=-n,n)/)))**2 )
      ESp_btd(n) = 0.5 / (n*(n+1)) * sum( w_Div_bt(l_nm(n,(/(m,m=-n,n)/)))**2 )
      ESp_bcd(n) = 0.5 / (n*(n+1)) * sum( w_Div_bc(l_nm(n,(/(m,m=-n,n)/)))**2 )
      ESp_vup(n) = 0.5 / (n*(n+1)) * sum( w_Vor_up(l_nm(n,(/(m,m=-n,n)/)))**2 )
      ESp_vlow(n) = 0.5 / (n*(n+1)) * sum( w_Vor_low(l_nm(n,(/(m,m=-n,n)/)))**2 )
    enddo

    ! sig(1)=0.81 sig(2)=0.26
    xy_Vor_bt(:,:) = ( xyz_Vor(:,:,2) + xyz_Vor(:,:,1) )/2
    xy_Vor_bc(:,:) = ( xyz_Vor(:,:,2) - xyz_Vor(:,:,1) )/2
    xy_Div_bt(:,:) = ( xyz_Div(:,:,2) + xyz_Div(:,:,1) )/2
    xy_Div_bc(:,:) = ( xyz_Div(:,:,2) - xyz_Div(:,:,1) )/2

    ! 変数出力
    call HistoryPut( &
      & 'Vor_bt', xy_Vor_bt)                   ! (in)
    call HistoryPut( &
      & 'Vor_bc', xy_Vor_bc)                   ! (in)
    call HistoryPut( &
      & 'Div_bt', xy_Div_bt)                   ! (in)
    call HistoryPut( &
      & 'Div_bc', xy_Div_bc)                   ! (in)
    call HistoryPut( &
      & 'U_bt', xy_U_bt)                       ! (in)
    call HistoryPut( &
      & 'U_bc', xy_U_bc)                       ! (in)
    call HistoryPut( &
      & 'ESp_btv', ESp_btv)                    ! (in)
    call HistoryPut( &
      & 'ESp_bcv', ESp_bcv)                    ! (in)
    call HistoryPut( &
      & 'ESp_bcd', ESp_bcd)                    ! (in)
    call HistoryPut( &
      & 'Energy', sum(ESp_btv)+sum(ESp_bcv)+sum(ESp_bcd)+sum(ESp_btd))   ! (in)
    call HistoryPut( &
      & 'Energy_btv', sum(ESp_btv)/(sum(ESp_btv)+sum(ESp_bcv)+sum(ESp_bcd)+sum(ESp_btd)))   ! (in)
    call HistoryPut( &
      & 'Energy_bcv', sum(ESp_bcv)/(sum(ESp_btv)+sum(ESp_bcv)+sum(ESp_bcd)+sum(ESp_btd)))   ! (in)
    call HistoryPut( &
      & 'Energy_btd', sum(ESp_btd)/(sum(ESp_btv)+sum(ESp_bcv)+sum(ESp_bcd)+sum(ESp_btd)))   ! (in)
    call HistoryPut( &
      & 'Energy_bcd', sum(ESp_bcd)/(sum(ESp_btv)+sum(ESp_bcv)+sum(ESp_bcd)+sum(ESp_btd)))   ! (in)
    call HistoryPut( &
      & 'Bcd_vup', sum(ESp_bcd)/sum(ESp_vup) )   ! (in)
    call HistoryPut( &
      & 'Bcd_vlow', sum(ESp_bcd)/sum(ESp_vlow) )   ! (in)
    call HistoryPut( &
      & 'K_Energy', K_Energy )         ! (in)
    call HistoryPut( &
      & 'Ps_mean', Ps_mean )         ! (in)
    call HistoryPut( &
      & 'Temp_mean', Temp_mean )         ! (in)
    call HistoryPut( &
      & 'Sk_n', Sk_n )         ! (in)
    call HistoryPut( &
      & 'Sk_s', Sk_s )         ! (in)
!    call HistoryPut( &
!      & 'Ens', xyz_Ens)                                    ! (in)

  call MessageNotify('M', subname, &
    & 'time = %d', &
    & i = (/int(z)/) )

  enddo

  call HistoryClose

end program translate_spectrum
