  ! deepconv で出力した z* 座標系でのデータを
  ! z 座標系のデータに書き換えるためのプログラム

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

program zstar_to_z

  use gtool_history ! モジュール指定

  implicit none

  !%%%%%%%%%%%%%%%%%%%%%%%%%%
  !- 実行前に確認する変数
  !--- 入出力に関連すること
!!$  character(128), parameter :: InFileName = "input_data/lee-wave_test_a3_VelZ.nc"  ! 入力ファイル名
!!$  character(128), parameter :: OutFileName = "output_data/lee-wave_test_a3_VelZ_zcoordinate.nc"      ! 出力ファイル名
!!$  character(128), parameter :: OutputTitle = "VelZ"          ! 出力ファイルのタイトル
!!$  character(128), parameter :: ValName = "VelZ"              ! 出力変数名


  character(128), parameter :: InFileName = "data_tmp/a3_L080km/lee-wave_test_a3_VelX.nc"
  character(128), parameter :: OutFileName = "output_data/a3_L080km/lee-wave_test_a3_VelX_zcoordinate.nc"
  character(128), parameter :: OutputTitle = "VelX"
  character(128), parameter :: ValName = "VelX"
  !
!!$  character(128), parameter :: InFileName = "data_tmp/a3_L080km/lee-wave_test_a3_VelZ.nc"
!!$  character(128), parameter :: OutFileName = "output_data/a3_L080km/lee-wave_test_a3_VelZ_zcoordinate.nc"
!!$  character(128), parameter :: OutputTitle = "VelZ"
!!$  character(128), parameter :: ValName = "VelZ"

!!$  character(128), parameter :: InFileName = "data_tmp/a3_L080km/lee-wave_test_a3_PTemp.nc"
!!$  character(128), parameter :: OutFileName = "output_data/a3_L080km/lee-wave_test_a3_PTemp_zcoordinate.nc"
!!$  character(128), parameter :: OutputTitle = "PTemp"
!!$  character(128), parameter :: ValName = "PTemp"

  !---地形に関連すること
  real(8), parameter :: H   = 10.0d+3  ! 領域の上端 (check)
  real(8), parameter :: a_0 = 1.0d+2  ! 山の半値幅
  real(8), parameter :: h_0 = 1.0d+2  ! 山の高さ
  real(8), parameter :: x_0 = 40.0d+3 ! 山の中央の x 座標
!!$  real(8), parameter :: x_0 = 100.0d+3 ! 山の中央の x 座標




  character(128), parameter :: ValUnit = ""                  ! 出力変数の単位

  integer :: Nx                ! x 方向の格子点数 (check)
  integer :: Ny                ! y 方向の格子点数 (check)
  integer :: Nzstar            ! z 方向の格子点数
  integer :: NTime             ! 時間ステップ数
!!$  integer, parameter :: NTimeStart = 1       ! はじめの時間ステップ
!!$  integer, parameter :: NTimeEnd =  121      ! 終わりの時間ステップ
  integer :: NTimeStart
  integer :: NTimeEnd
  integer, parameter :: NTimeStep =  1
  
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%

  real(8)            :: z_s        !地形の指定のために用意した変数

  real(8), allocatable  :: x          (:)
  real(8), allocatable  :: y          (:)
  real(8), allocatable  :: z          (:)
  real(8), allocatable  :: t_Time     (:)
  real(8), allocatable  :: Bval       (:, :, :)  ! z* 座標系の変数
  real(8), allocatable  :: Aval       (:, :, :)  ! z  座標系の変数
  real(8), allocatable  :: Bz         (:, :, :)  ! x-z* 座標系での z の格子
  real(8), allocatable  :: Az         (:)          ! x-z  座標系での z の格子

  real(8)            :: dz

  real(8)            :: a, b, ans

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

  character(64)      :: range
  
  real(8), pointer       :: a_Axis(:) !一次的な変数

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

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

  call HistoryGetPointer( InFileName, 'z', a_Axis )
  Nzstar = size( a_Axis )
  Nullify( a_Axis )

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

  NTimeStart = 1
!!$  NTimeStart = NTime
  NTimeEnd   = NTime
  if( NTime < NTimeEnd ) then
    write( 6, * ) "NTime < NTimeEnd"
    stop
  end if
  
  allocate( x(Nx) )
  allocate( y(Ny) )
  allocate( z(Nzstar) )
  allocate( t_Time(NTime) )
  allocate( Bval(Nx, Ny, Nzstar) )  ! z* 座標系の変数
  allocate( Aval(Nx, Ny, Nzstar) )  ! z  座標系の変数
  allocate( Bz(Nx, Ny, Nzstar) )    ! x-z* 座標系での z の格子
  allocate( Az(Nzstar) )            ! x-z  座標系での z の格子

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

  ! 鉛直格子間隔を定義
  dz = H/Nzstar

  ! 絵を描くための格子点を作る
  do k = 1, Nzstar
    Az(k) = dz*k
  end do

  ! [出力] 出力ファイル作成 (開く)
  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,Nzstar,0/),                                &
    & 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 HistoryAddVariable(                              &
    & varname=ValName, dims=(/'x ','y ','z ', 't '/),   &
    & longname=OutputTitle, units=ValUnit, xtype='real' &
    & )

  call HistoryAddattr( &  !属性を与える, 欠損値
    & varname=ValName, attrname="missing_value", &
    & value=-999.0 )

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

  ! 時間ループの開始
  do it = NTimeStart, NTimeEnd, NTimeStep

    if ( NTimeEnd < it ) then
      write( 6, *) "NTimeEnd < it"
      stop
    end if
    
    write( 6, * ) it, t_Time(it)

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

    ! [入力] it 番目のデータを読む
    call HistoryGet(InFileName, ValName, Bval, range = range)

    j = 1
    iloop: do i = 1, Nx
      kloop: do k = 1, Nzstar
        ! 地形を指定
        z_s       = h_0 * (a_0**2)/(a_0**2 + (x(i) - x_0)**2 )
        Bz(i,j,k) = z_s + (z(k) * (H-z_s)/H)
      end do kloop
    end do iloop

    j = 1
    ticks_iloop: do i = 1, Nx
      !格子点を作る
      ticks_kloop: do k = 1, Nzstar

        !格子点が今あるデータ点に対してどこにあるか確認する
        !データの上端より上
        if ( Az(k) >= Bz(i,j,Nzstar) ) then
          ans = -999
          Aval(i,j,k)=ans

          !データの下端より下
        else if ( Az(k) <= Bz(i,j,1) ) then
          ans = -999
          Aval(i,j,k)=ans

          !データがあるところの間
        else
          dummy_k_lloop: do l = Nzstar, 1, -1
            if ( Az(k) >= Bz(i,j,l) ) then
              !前後のデータが同じだった場合の処理(傾き a が定義できなくなることの対策)
              if ( Bval(i,j,l) == Bval(i,j,l+1) )then
                ans = Bval(i,j,l)
                Aval(i,j,k)=ans
                exit

              else
                a =   ( Bz(i,j,l+1) - Bz(i,j,l) )/( Bval(i,j,l+1) - Bval(i,j,l) )
                b = - ( Bz(i,j,l+1) - Bz(i,j,l) ) * Bval(i,j,l)/ &
                  & ( Bval(i,j,l+1) - Bval(i,j,l) ) + Bz(i,j,l)
                ans         = ( Az(k) - b )/a
                Aval(i,j,k) = ans
                exit
              end if

            else
            end if

            ! if ( Az(k) >= Bz(i,j,l) ) then
            !   ans  =    ( Bval(i,j,l+1) - Bval(i,j,l) ) &
            !     &     / ( Bz  (i,j,l+1) - Bz  (i,j,l) ) &
            !     &     * ( Az  (k)       - Bz  (i,j,l) ) &
            !     &     + Bval(i,j,l) 
            !   Aval(i,j,k) = ans
            ! end if
            
          end do dummy_k_lloop

        end if

      end do ticks_kloop
    end do ticks_iloop
    
    call HistoryPut('t' ,t_Time(it) )
    call HistoryPut(ValName, Aval )
    
  end do

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

  stop

end program zstar_to_z
