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

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

program flux

  use gtool_history ! モジュール指定

  implicit none

  !%%%%%%%%%%%%%%%%%%%%%%%%%%
  !- 実行前に確認する変数
  !--- 入出力に関連すること
  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_4 = "input_data/flux_a3.nc"                  ! 入力ファイル名
  character(128), parameter :: OutFileName = "out.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


  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) )

  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 '/),                                   &
    & dimsizes=(/Nx,Ny,Nz,NTime/),                                &
    & longnames=(/'x   ','y   ','z   ','time'/),                      &
    & units=(/'m',                     &
    &         'm',                     &
    &         'm',                     &
    &         's'/)                    &
    & )

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

  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 HistoryAddattr( &  !属性を与える, 欠損値
    & varname=OutValName, attrname="missing_value", &
    & value=-999.0 )

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

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

  ! 時間ループの開始
  do it = 1, 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-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)
         end do
      end if

    end do

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

  end do

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

  stop

end program flux
