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

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

program zstar_to_z

  use gtool_history ! モジュール指定

  implicit none

  !%%%%%%%%%%%%%%%%%%%%%%%%%
  !- 実行前に確認する変数
  !--- 入出力に関連すること
!!$  character(128), parameter :: InFileName = "../test_a3_x200km_basic/BasicZ.nc"       ! 入力ファイル名
!!$  character(128), parameter :: OutFileName = "output_data/DensBZ_a3_zcoordinate.nc"          ! 出力ファイル名
!!$  character(128), parameter :: OutputTitle = "density"         ! 出力ファイル名
  !
  character(128), parameter :: InFileName = "../test_a3_L080km_basic/BasicZ.nc"
  character(128), parameter :: OutFileName = "output_data/a3_L080km/DensBZ_a3_zcoordinate.nc"
  !
  character(128), parameter :: OutputTitle = "density"         ! 出力ファイル名
  character(128), parameter :: ValName = "DensBZ"              ! 変換する変数名
  character(128), parameter :: ValUnit = ""                    ! 変数の単位
  

  !--- 地形に関連すること
  real(8), parameter :: H   = 10.0d+3  ! 領域の上端 deepconv の Zmaxと同じにしておく
  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 座標

  integer            :: Nx, Nzstar
  integer, parameter :: Ny     =  1
  integer, parameter :: dz     =  20           ! z方向の格子間隔

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

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

  real(8), allocatable :: x(:), y(:), z(:)
  real(8), allocatable :: Bval(:,:,:), Aval(:,:,:)
  real(8), allocatable :: Bz(:,:,:), Az(:)

  real(8)            :: a, b, ans

  integer            :: i, j, k, l

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

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

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

  allocate( x(Nx) )
  allocate( y(Ny) )
  allocate( z(Nzstar) )
  allocate( Bval(Nx, Ny, Nzstar) )
  allocate( Aval(Nx, Ny, Nzstar) )
  allocate( Bz(Nx, Ny, Nzstar) )
  allocate( Az(Nzstar) )
  
  call HistoryGet(InFileName, 'x'  , x)      ! [入力] x
  call HistoryGet(InFileName, 'y'  , y)      ! [入力] y
  call HistoryGet(InFileName, 'z'  , z)      ! [入力] zstar

  ! 絵を描くための格子点を作る
  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'/),                                       &
    & dimsizes=(/Nx,Ny,Nzstar/),                                  &
    & longnames=(/'x','y','z'/),                                  &
    & units=(/'m',                     &
    &         'm',                     &
    &         'm'/)                    &
    & )

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

  call HistoryAddVariable(                              &
    & varname=ValName, dims=(/'x','y','z'/),            &
    & 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 )

  ! 変換後の変数の出力
  call HistoryGet(InFileName, ValName, Bval)

  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
        end do dummy_k_lloop

      end if

    end do ticks_kloop
  end do ticks_iloop

  call HistoryPut(ValName, Aval )
  

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

  stop

end program zstar_to_z
