! Sample program for gtool_history/gtool5
!
program nc_io

  use gtool_history ! モジュール指定

  implicit none

  character(128)     :: InFileName               ! 入力ファイル名
  character(128)     :: OutFileName              ! 出力ファイル名

  integer, parameter :: NLon   = 144             ! 経度格子点数
  integer, parameter :: NLat   =  73             ! 緯度格子点数
  integer, parameter :: NLevel =  17             ! 高度(圧力)格子点数
  integer, parameter :: NTime  = 365             ! 時間格子点数

  real(4)            :: x_Lon      (NLon  )             ! 経度配列
  real(4)            :: y_Lat      (NLat  )             ! 緯度配列
  real(4)            :: z_Level    (NLevel)             ! 高度(圧力)配列
  real(4)            :: t_Time     (NTime )             ! 時間配列
  real(4)            :: xyz_Temp   (NLon, NLat, NLevel) ! 温度配列
  real(4)            :: xyz_PotTemp(NLon, NLat, NLevel) ! 温位配列

  real(4)            :: Kappa = 0.29     ! R/Cp
  real(4)            :: P0    = 1000.0   ! 基準圧力

  integer            :: it
  integer            :: i
  integer            :: j
  integer            :: k

  character(64)      :: range


  InFileName  = "air.2019.nc"   ! 入力ファイル名の設定
  OutFileName = "output.nc"     ! 出力ファイル名の設定


  call HistoryGet(InFileName, 'lon'  , x_Lon  ) ! [入力] 経度を読む
  call HistoryGet(InFileName, 'lat'  , y_Lat  ) ! [入力] 緯度を読む
  call HistoryGet(InFileName, 'level', z_Level) ! [入力] 圧力を読む
  call HistoryGet(InFileName, 'time' , t_Time ) ! [入力] 時間を読む


  ! [出力] 出力ファイル作成 (開く)
  call HistoryCreate(                                                 &
    & file=OutFileName, title='atmospheric temperature',              &
    & source='NCEP/NCAR Reanalysis data',                             &
    & institution='Department of Planetology, Kobe University',       &
    & dims=(/'lon  ','lat  ','level','time '/),                       &
    & dimsizes=(/NLon,NLat,NLevel,NTime/),                            &
    & longnames=(/'longitude','latitude ','level    ', 'time     '/), &
    & units=(/'degrees_east                    ',                     &
    &         'degrees_north                   ',                     &
    &         'hPa                             ',                     &
    &         'hours since 1800-01-01 00:00:0.0'/)                    &
    & )

  call HistoryPut('lon'  ,x_Lon  )            ! [出力] 次元変数出力, 経度
  call HistoryPut('lat'  ,y_Lat  )            ! [出力] 次元変数出力, 緯度
  call HistoryPut('level',z_Level)            ! [出力] 次元変数出力, 高度(圧力)
  call HistoryPut('time' ,t_Time )            ! [出力] 次元変数出力, 時間


  ! [出力] 変数定義, 温位
  call HistoryAddVariable(                                         &
    & varname='pottemp', dims=(/'lon  ','lat  ','level','time '/), &
    & longname='potential temperature', units='K', xtype='real'    &
    & )

  do it = 1, NTime
    write( 6, * ) it, t_Time(it)

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

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

    ! 温位の計算
    do k = 1, NLevel
      do j = 1, NLat
        do i = 1, NLon
          xyz_PotTemp(i,j,k) = xyz_Temp(i,j,k) * ( P0 / z_Level(k) )**Kappa
        end do
      end do
    end do

    ! [出力] it 番目のデータを出力, 温位
    call HistoryPut('pottemp', xyz_PotTemp, range=range)
  end do

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

  stop
end program nc_io