!  鉛直方向の運動量フラックスを求めるプログラム

! * 連続時間データを扱う場合に使用できる
! * 実行時は変数定義の箇所を確認する

program flux

  use gtool_history ! モジュール指定

  use fft8_module

  implicit none

  !%%%%%%%%%%%%%%%%%%%%%%%%%%
  !- 実行前に確認する変数
  !--- 入出力に関連すること
  ! A1
!!$  character(128), parameter :: InFileName_1 = "/home/shiraham/flux_ver2/output_data/lee-wave_test_a1_VelX_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_2 = "/home/shiraham/flux_ver2/output_data/lee-wave_test_a1_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_3 = "/home/shiraham/flux_ver2/output_data/DensBZ_a1_zcoordinate.nc"   ! 入力ファイル名
!!$  real(8) :: xrange = 200d3
!!$  real(8) :: factor = 1.0d0
!!$  integer :: out_kindex = 6
  ! A2
!!$  character(128), parameter :: InFileName_1 = "/home/shiraham/flux_ver2/output_data/lee-wave_test_a2_VelX_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_2 = "/home/shiraham/flux_ver2/output_data/lee-wave_test_a2_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_3 = "/home/shiraham/flux_ver2/output_data/DensBZ_a2_zcoordinate.nc"   ! 入力ファイル名
!!$  real(8) :: xrange = 200d3
!!$  real(8) :: factor = 1.0d0
!!$  integer :: out_kindex = 6
  ! A3
!!$  character(128), parameter :: InFileName_1 = "../a3_h100m/data/lee-wave_test_a3_VelX_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_2 = "../a3_h100m/data/lee-wave_test_a3_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_3 = "../a3_h100m/data/DensBZ_a3_zcoordinate.nc"   ! 入力ファイル名
!!$  character(128), parameter :: InFileName_1 = "../output_data/a3/lee-wave_test_a3_VelX_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_2 = "../output_data/a3/lee-wave_test_a3_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_3 = "../output_data/a3/DensBZ_a3_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_1 = "../output_data/a3-2/lee-wave_test_a3_VelX_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_2 = "../output_data/a3-2/lee-wave_test_a3_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_3 = "../output_data/a3-2/DensBZ_a3_zcoordinate.nc"
!!$  real(8) :: xrange = 40d3
!!$  real(8) :: factor = 1.0d0
!!$  integer :: out_kindex = 26
  ! A3' (10 m)
!!$  character(128), parameter :: InFileName_1 = "../a3_h010m/data/lee-wave_test_a3_VelX_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_2 = "../a3_h010m/data/lee-wave_test_a3_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_3 = "../a3_h010m/data/DensBZ_a3_zcoordinate.nc"   ! 入力ファイル名
!!$  real(8) :: xrange = 40d3
!!$  real(8) :: factor = 1.0d0
!!$  integer :: out_kindex = 26
  ! A4
!!$  character(128), parameter :: InFileName_1 = "/home/shiraham/flux_ver2/output_data/lee-wave_test_a4_VelX_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_2 = "/home/shiraham/flux_ver2/output_data/lee-wave_test_a4_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_3 = "/home/shiraham/flux_ver2/output_data/DensBZ_a4_zcoordinate.nc"   ! 入力ファイル名
!!$  real(8) :: xrange = 10d3
!!$  real(8) :: factor = 1.0d0
!!$  integer :: out_kindex = 101


  character(128), parameter :: InFileName_1 = &
    & "../output_data_joho03_sponge/a3_L080km_H020km_SpongeVelX10Zt10kmTau1800s/lee-wave_test_a3_VelX_zcoordinate.nc"
  character(128), parameter :: InFileName_2 = &
    & "../output_data_joho03_sponge/a3_L080km_H020km_SpongeVelX10Zt10kmTau1800s/lee-wave_test_a3_VelZ_zcoordinate.nc"
  character(128), parameter :: InFileName_3 = &
    & "../output_data_joho03_sponge/a3_L080km_H020km_SpongeVelX10Zt10kmTau1800s/DensBZ_a3_zcoordinate.nc"
  real(8), parameter :: xrange = 80d3
  real(8), parameter :: factor = xrange / 40d3
  integer :: out_kindex = 26

!!$  character(128), parameter :: InFileName_1 = "../output_data_joho08/a3_L160km/lee-wave_test_a3_VelX_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_2 = "../output_data_joho08/a3_L160km/lee-wave_test_a3_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: InFileName_3 = "../output_data_joho08/a3_L160km/DensBZ_a3_zcoordinate.nc"
!!$  real(8) :: xrange = 160d3
!!$  real(8) :: factor = 160d3 / 40d3
!!$  integer :: out_kindex = 26


!!$  character(128), parameter :: InFileName_4 = "input_data/flux_a3.nc"                  ! 入力ファイル名
  character(128), parameter :: OutFileName = "momfluxspe_deepconv.nc"

  character(128), parameter :: InValName_1 = "VelX"          ! 入力変数名
  character(128), parameter :: InValName_2 = "VelZ"          ! 入力変数名
  character(128), parameter :: InValName_3 = "DensBZ"          ! 入力変数名
  
  character(128), parameter :: OutputTitle = "Momentum Flux" ! 出力ファイルのタイトル
  character(128), parameter :: OutValName = "Mflux"          ! 出力変数名
  character(128), parameter :: ValUnit = ""                  ! 出力変数の単位

  integer :: Nx, Ny, Nz
  integer :: NTime

  !%%%%%%%%%%%%%%%%%%%%%%%%%%%

  real(8), allocatable   :: x(:)
  real(8), allocatable   :: y(:)
  real(8), allocatable   :: z(:)
  real(8), allocatable   :: t_Time(:)
  real(8), allocatable   :: Bval_1  (:, :, :)  ! before
  real(8), allocatable   :: Bval_2  (:, :, :)  ! before
  real(8), allocatable   :: DensBZ  (:, :, :)  ! before
  real(8), allocatable   :: Aval    (:, :)       ! after
  real(8), allocatable   :: Aval_2  (:, :)       ! after
  real(8), allocatable   :: AFlux(:)
  
  real(8)            :: sum
  real(8)            :: avs
  real(8)            :: missing
  
  real(8)            :: a, b, ans
  real(8)            :: mounttop

  integer            :: it
  integer            :: i, j, k, l

  character(64)      :: range
  real(8), pointer   :: a_Axis(:)

  real(8) :: ave1, ave2

  real(8), allocatable   :: xz_ur(:,:)
  real(8), allocatable   :: xz_ui(:,:)
  real(8), allocatable   :: xz_wr(:,:)
  real(8), allocatable   :: xz_wi(:,:)
  integer                :: ita( 5, 2 )
  real(8), allocatable   :: dta(:,:)
  real(8), allocatable   :: w_fft_UWbar(:)
  real(8), allocatable   :: fft_UWbar
  real(8), allocatable   :: wz_fft_UWbar(:,:)
  real(8), allocatable   :: z_fft_UWbar(:)
  real(8), allocatable   :: w_wn(:)
  integer :: lneg
  real(8) :: pi = acos(-1.0d0)


  write( 6, * ) "mounttop = "
  read( 5, * ) mounttop

  call HistoryGetPointer( InFileName_1, 'x', a_Axis )
  Nx = size( a_Axis )
  Nullify( a_Axis )

  call HistoryGetPointer( InFileName_1, 'y', a_Axis )
  Ny = size( a_Axis )
  Nullify( a_Axis )

  call HistoryGetPointer( InFileName_1, 'z', a_Axis )
  Nz = size( a_Axis )
  Nullify( a_Axis )

  call HistoryGetPointer( InFileName_1, 't', a_Axis )
  NTime = size( a_Axis )
  Nullify( a_Axis )

  allocate( x(Nx) )
  allocate( y(Ny) )
  allocate( z(Nz) )
  allocate( t_Time(NTime) )
  allocate( Bval_1(Nx, Ny, Nz) )
  allocate( Bval_2(Nx, Ny, Nz) )
  allocate( DensBZ(Nx, Ny, Nz) )
  allocate( Aval  (Ny, Nz) )
  allocate( Aval_2(Ny, Nz) )
  allocate( AFlux(500) )



  allocate( xz_ur(nx,nz) )
  allocate( xz_ui(nx,nz) )
  allocate( xz_wr(nx,nz) )
  allocate( xz_wi(nx,nz) )
  allocate( dta( nx, 2 ) )
  allocate( wz_fft_UWbar(nx,nz) )
  allocate( z_fft_UWbar(nz) )
  allocate( w_wn(nx) )

  do i = 1, nx
    w_wn(i) = 2.0d0 * pi / xrange * dble(i-1)
  end do

  call fft8_fft_init( nx, ita, dta )



  call HistoryGet(InFileName_1, 'x'  , x)      ! [入力] x
  call HistoryGet(InFileName_1, 'y'  , y)      ! [入力] y
  call HistoryGet(InFileName_1, 'z'  , z)      ! [入力] zstar
  call HistoryGet(InFileName_1, 't'  , t_Time) ! [入力] t


  ! [出力] 出力ファイル作成 (開く)
  call HistoryCreate(                                                 &
    & file=OutFileName, title=OutputTitle,                            &
    & source='NCEP/NCAR Reanalysis data',                             &
    & institution='Department of Planetology, Kobe University',       &
    & dims=(/'x ','y ','z ','t ','wn'/),                                   &
    & dimsizes=(/Nx,Ny,Nz,NTime,Nx/),                                &
    & longnames=(/'x   ','y   ','z   ','time', 'wn  '/),              &
    & units=(/'m  ',                     &
    &         'm  ',                     &
    &         'm  ',                     &
    &         's  ',                     &
    &         'm-1'/)                  &
    & )

  call HistoryPut('x' ,x       )            ! [出力] x
  call HistoryPut('y' ,y       )            ! [出力] y
  call HistoryPut('z' ,z       )            ! [出力] z
  call HistoryPut('t' ,t_Time  )            ! [出力] t
  call HistoryPut('wn',w_wn    )            ! [出力] wn

  call HistoryAddVariable(                              &
       & varname=OutValName, dims=(/'y ','z ', 't '/),   &
       & longname=OutputTitle, units=ValUnit, xtype='real' &
       & )

!!$  call HistoryAddVariable(                              &
!!$       & varname="Normflux", dims=(/'y ','z ', 't '/),   &
!!$       & longname=OutputTitle, units=ValUnit, xtype='real' &
!!$       & )

  call HistoryAddVariable(                              &
       & varname='MfluxSpe', dims=(/'wn','z ', 't '/),  &
       & longname=OutputTitle, units=ValUnit, xtype='real' &
       & )


  call HistoryAddattr( &  !属性を与える, 欠損値
    & varname=OutValName, attrname="missing_value", &
    & value=-999.0 )
  call HistoryAddattr( &  !ncdump で欠損値にするため
    & varname=OutValName, attrname="_FillValue", &
    & value=-999.0 )

  call HistoryAddattr( &  !ncdump で欠損値にするため
    & varname='MfluxSpe', attrname="missing_value", &
    & value=-999.0 )
  call HistoryAddattr( &  !ncdump で欠損値にするため
    & varname='MfluxSpe', attrname="_FillValue", &
    & value=-999.0 )

  call HistoryGet(InFileName_3, InValName_3, DensBZ)
!!$  call HistoryGet(InFileName_4, "flux", AFlux)

  ! 時間ループの開始
  do it = 1, NTime
!!$  do it = NTime, NTime
    write( 6, * ) it, t_Time(it)

    ! 出力用の時刻情報文字列の作成
    write( range, '(a,i5.5)' ) 't=^', it
    write( 6, * ) range

    ! [入力] it 番目のデータを読む
    call HistoryGet(InFileName_1, InValName_1, Bval_1, range = range)
    call HistoryGet(InFileName_2, InValName_2, Bval_2, range = range)

    do k = 1, nz
      do j = 1, ny
        ave1 = 0.0d0
        ave2 = 0.0d0
        do i = 1, nx
          ave1 = ave1 + Bval_1(i,j,k)
          ave2 = ave2 + Bval_2(i,j,k)
        end do
        ave1 = ave1 / dble( nx )
        ave2 = ave2 / dble( nx )
        Bval_1(:,j,k) = Bval_1(:,j,k) - ave1
        Bval_2(:,j,k) = Bval_2(:,j,k) - ave2
      end do
    end do

    j = 1
    do k = 1, nz
      xz_ur(:,k) = Bval_1(:,j,k)
      xz_wr(:,k) = Bval_2(:,j,k)
    end do
    xz_ui = 0.0d0
    xz_wi = 0.0d0
    call fft8_fft( nx, nz, ita, dta, xz_ur, xz_ui, 1 )
    call fft8_fft( nx, nz, ita, dta, xz_wr, xz_wi, 1 )

    do k = 1, nz
      l = 1
      wz_fft_UWbar(l,k) = 0.0d0
      do l = 1+1, nx
        lneg = nx-(l-2)
        wz_fft_UWbar(l,k) = xz_ur(l,k) * xz_wr(lneg,k) - xz_ui(l,k) * xz_wi(lneg,k)
      end do
    end do
    do k = 1, nz
      wz_fft_UWbar(:,k) = wz_fft_UWbar(:,k) * DensBZ(:,j,k)
    end do
    wz_fft_UWbar = wz_fft_UWbar * factor
!!$    wz_fft_UWbar = wz_fft_UWbar &
!!$      & / ( (2.0d0*pi/xrange)*(2.0d0*pi/xrange) * 2.0d0*pi/(2.0d0*pi/xrange) / xrange )
    z_fft_UWbar = 0.0d0
    do k = 1, nz
      do l = 1, nx
!!$        z_fft_UWbar(k) = z_fft_UWbar(k) &
!!$          & + wz_fft_UWbar(l,k) &
!!$          &   * ( (2.0d0*pi/xrange)*(2.0d0*pi/xrange) * 2.0d0*pi/(2.0d0*pi/xrange) / xrange )
        z_fft_UWbar(k) = z_fft_UWbar(k) + wz_fft_UWbar(l,k)
      end do
    end do
    do k = 1, nz
      if ( z(k) <= mounttop ) then
        z_fft_UWbar (k)   = -999.0d0
        wz_fft_UWbar(:,k) = -999.0d0
      end if
    end do
    do k = 1, out_kindex
!!$    do k = 1, nz
!!$      if ( ( z(k) > mounttop ) .and. ( z(k) <= mounttop + 8*(z(2)-z(1)) ) ) then
        write( 6, * ) 'SUM (in spe. space ) : ', k, z_fft_UWbar(k)
!!$      end if
    end do
    if ( it == NTime ) then
      k = 6
      k = 8
      k = out_kindex
!!$      do k = 1, nz
!!$        if ( z(k) > mounttop ) then
!!$          exit
!!$        end if
!!$      end do
!!$      k = k + 4
      open( 81, file = 'momfluxspe_deepconv.txt', status = 'unknown' )
      write( 81, * ) '# mounttop = ', mounttop
      write( 81, * ) '# k = ', k
      do j = 1, k+1
        write( 81, * ) '# z_fft_UWbar : ', j, z_fft_UWbar(j)
      end do
      write( 81, * ) '# <index>   <wn (m-1)>   <uw spe. dens. (kg s-1)>'
      do l = 1, nx
        write( 81, * ) l, (2.0d0*pi/xrange)*(l-1), &
          & wz_fft_UWbar(l,k) / (2.0d0*pi/xrange)
!!$        write( 81, * ) l, (2.0d0*pi/xrange)*(l-1), &
!!$          & wz_fft_UWbar(l,k)
      end do
      close( 81 )
      !
      open( 81, file = 'momfluxspe_deepconv_rm.txt', status = 'unknown' )
      write( 81, * ) '# mounttop = ', mounttop
      write( 81, * ) '# k = ', k
      do j = 1, k+1
        write( 81, * ) '# z_fft_UWbar : ', j, z_fft_UWbar(j)
      end do
      write( 81, * ) '# <index>   <wn (m-1)>   <uw spe. dens. (kg s-1)>'
      do l = 1+1, nx-2, 3
        write( 81, * ) l, (2.0d0*pi/xrange)*((l-1)+l+(l+1))/3.0, &
          & (wz_fft_UWbar(l,k)+wz_fft_UWbar(l+1,k)+wz_fft_UWbar(l+2,k))/3 / (2.0d0*pi/xrange)
      end do
!!$      do l = 1+1, nx-1, 2
!!$        write( 81, * ) l, (2.0d0*pi/xrange)*((l-1)+l)/2.0, &
!!$          & (wz_fft_UWbar(l,k)+wz_fft_UWbar(l+1,k))/2 / (2.0d0*pi/xrange)
!!$      end do
      close( 81 )


      z_fft_UWbar = 0.0d0
      do k = 1, nz
        do l = 1, 15
          z_fft_UWbar(k) = z_fft_UWbar(k) + wz_fft_UWbar(l,k)
        end do
      end do
      k = 6
      l = 15
      write( 6, * )  (2.0d0*pi/xrange)*(l-1)
      write( 6, * )  z_fft_UWbar(k)


    end if


    j = 1
    do k = 1, Nz-1
      sum = 0.0d0
      avs = 0.0d0
      missing = 0.0d0

      do i = 1, Nx
         sum = sum + DensBZ(i,j,k)*Bval_1(i,j,k)*Bval_2(i,j,k)
      end do

      avs = sum/(Nx - missing)

      if ( z(k) <= mounttop ) then
         Aval(j,k) = -999.0d0
      else
         do i = 1, Nx
!!$            Aval(j,k) = avs
!!$!            Aval_2(j,k) = avs/AFlux(1)
            Aval(j,k) = avs * factor
         end do
      end if

    end do

    do k = 1, out_kindex
!!$    do k = 1, nz
!!$      if ( ( z(k) > mounttop ) .and. ( z(k) <= mounttop + 8*(z(2)-z(1)) ) ) then
        write( 6, * ) 'SUM (in phys. space) : ', k, Aval(j,k)
!!$      end if
    end do

    call HistoryPut(OutValName, Aval, range=range )
!!$   call HistoryPut("Normflux", Aval_2, range=range )

    do k = 1, nz
      if ( z(k) <= mounttop ) then
        wz_fft_UWbar(:,k) = -999.0d0
      else
        wz_fft_UWbar(:,k) = wz_fft_UWbar(:,k) / (2.0d0*pi/xrange)
      end if
    end do
    call HistoryPut('MfluxSpe', wz_fft_UWbar, range=range )

  end do

  ! [出力] 出力ファイルを閉じる
  call HistoryClose

  stop

end program flux
