!
!  バイナリファイルから netCDF データを生成するプログラム
!
!  バイナリ : /global/pr_preci/2000/monthly/gpr00100.zip
!  
!  方針 : 格子点, 次元の取り方等, 構造はそのまま
!         netCDF コンベンションは Gtool4 - NetCDF Convention に
!         従うものとする.
!  2002/10/19  塚原 大輔 作成
!

program test_netcdf

!    モジュールの呼び出し
use netcdf          !netcdfモジュール呼び出し

  implicit none

!   元データ(TRMM/GANAL cdrom)から変数を取り出すための
!                                           容れ物変数の宣言
!

  character(len = 80) :: fname          !ファイル名(ex.gpr00100)

!ファイルに格納される各変数

  character(len = 80) :: header(7)      !ファイルの簡単な情報が書かれている変数.

  integer(2)          :: pr(720,153)      !地上における層状性, 対流性降雨量の合計
  integer(2)          :: pr_2(720,153)    !高度2kmにおける層状性, 対流性降雨量の合計
  integer(2)          :: pr_4(720,153)    !高度4kmにおける層状性. 対流性降雨量の合計
  integer(2)          :: pr_6(720,153)    !高度6kmにおける層状性, 対流性降雨量の合計

  integer(2)          :: cpr_2(720,153)  !高度2kmにおける対流性降雨量
  integer(2)          :: cpr_4(720,153)  !高度4kmにおける対流性降雨量
  integer(2)          :: cpr_6(720,153)  !高度6kmにおける対流性降雨量

  integer(2)          :: spr_2(720,153)  !高度2kmにおける層状性降雨量
  integer(2)          :: spr_4(720,153)  !高度4kmにおける層状性降雨量
  integer(2)          :: spr_6(720,153)  !高度6kmにおける層状性降雨量

!   netcdf 出力用変数の宣言

  integer   ::  ncid, i               
  integer   ::  dim_lonid, dim_latid, dim_altid  !次元ID. 左から緯度, 経度, 高度.
  integer   ::  var_lonid, var_latid, var_altid  !変数ID. 同上.
  integer   ::  var_prid, var_cprid, var_sprid   !降雨量. 左から合計, 対流性, 層状性.
  integer   ::  stat                             !エラー判断
  integer   ::  prmid(3)                         !降雨量の次元IDベクトル.
  real, dimension(720,153) :: &
                rhValues = 0                     !全ての格子において値が0である配列.
                                                 !対流性, 層状性降雨量が共に地上での値がないので
                                                 !高度の次元をそろえるために用意.
!データ読み取り
  fname='gpr00100'
  open(80,file=fname,form='unformatted',status='old')
  read(80) header
  read(80) pr
  read(80) pr_2
  read(80) pr_4
  read(80) pr_6
    
  read(80) cpr_2
  read(80) cpr_4
  read(80) cpr_6

  read(80) spr_2
  read(80) spr_4
  read(80) spr_6
  close(80)
!    netCDFファイルの生成 命令文
  stat = nf90_create('gpr00100.nc',nf90_noclobber,ncid)
  if(stat /= nf90_noerr)  call handle_err(stat)

!    次元の定義

  stat = nf90_def_dim(ncid, 'lon', 720,dim_lonid)     !経度の定義
    if(stat /= nf90_noerr) call handle_err(stat)

  stat = nf90_def_dim(ncid, 'lat', 153, dim_latid)     !緯度の定義
    if(stat /= nf90_noerr) call handle_err(stat)

  stat = nf90_def_dim(ncid,'alt',4,dim_altid)     !鉛直方向ベクトルの定義
    if(stat /= nf90_noerr) call handle_err(stat)

!    変数の定義

  stat = nf90_def_var(ncid,'lon',nf90_double,&
                                   &dim_lonid,var_lonid)   !緯度座標変数の定義
    if(stat /= nf90_noerr) call handle_err(stat)


  stat = nf90_def_var(ncid,'lat',nf90_double,&
                                   &dim_latid,var_latid)   !緯度座標変数の定義
    if(stat /= nf90_noerr) call handle_err(stat)


  stat = nf90_def_var(ncid,'dim',nf90_double,&
                                   &dim_altid,var_altid)   !鉛直座標変数の定義
    if(stat /= nf90_noerr) call handle_err(stat)

  prmid = (/dim_lonid,dim_latid,dim_altid/)              ! 降雨量の次元配列

  stat = nf90_def_var(ncid,'pr',nf90_double,&
                                     &prmid,var_prid)  ! 降雨強度の定義
    if(stat /= nf90_noerr) call handle_err(stat)

  stat = nf90_def_var(ncid,'pr_c',nf90_double,&
                                     &prmid,var_cprid)  ! 降雨強度の定義
    if(stat /= nf90_noerr) call handle_err(stat)

  stat = nf90_def_var(ncid,'pr_s',nf90_double,&
                                     &prmid,var_sprid)  ! 降雨強度の定義
    if(stat /= nf90_noerr) call handle_err(stat)

!     属性定義
!!    大域属性

  stat = nf90_put_att(ncid,nf90_global,'Conventions',&
      &'http://www.gfd-dennou.org/arch/gtool4/conventions/')! お決まり
    if(stat /= nf90_noerr) call handle_err(stat)


  stat = nf90_put_att(ncid, nf90_global,'title',&
      & 'TRMM/PR Precipitation ( Monthly, Jan. 2000 ), Global ') ! タイトル
    if(stat /= nf90_noerr) call handle_err(stat)

  stat = nf90_put_att(ncid,nf90_global,'source',& 
      &'For study of global climate change::TRMM Rain Data Set')!出典
    if(stat /= nf90_noerr) call handle_err(stat)

  stat = nf90_put_att(ncid,nf90_global,'history',& 
      &'2002-11-27T04:28:20+0900 tsukahara')!履歴
    if(stat /= nf90_noerr) call handle_err(stat)

!!    各変数の属性@単位
  stat = nf90_put_att(ncid,var_lonid,'units','degrees_east')!緯度
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_latid,'units','degrees_north')!経度
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_altid,'units','km')
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_prid,'units','mm')
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_cprid,'units','mm')
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_sprid,'units','mm')
    if(stat /= nf90_noerr) call handle_err(stat)
  

!!    各変数の属性@longname

  stat = nf90_put_att(ncid,var_lonid,'long_name','longitude')
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_lonid,'modulo','360.0')
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_latid,'long_name','latitude')
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_altid,'long_name','altitude')
    if(stat /= nf90_noerr) call handle_err(stat)

  stat = nf90_put_att(ncid,var_prid,'long_name',&
                              &  ' precipitaion ')
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_cprid,'long_name',&
                              &'convective precipitation ')
    if(stat /= nf90_noerr) call handle_err(stat)
  stat = nf90_put_att(ncid,var_sprid,'long_name',&
                              & 'stratiform precipitation')
    if(stat /= nf90_noerr) call handle_err(stat)


!    定義モードから抜ける

  stat = nf90_enddef(ncid)

    if(stat /= nf90_noerr) call handle_err(stat) 


!    変数にデータを代入する.

  stat = nf90_put_var(ncid,var_lonid,(/((i-1)*0.5,i=1,720)/) )!緯度方向の座標値
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_latid,(/(((i-1)*0.5-38),i=1,153)/) )!経度方向の座標値
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_altid,(/(i,i=0,6,2)/) )        !高度方向の座標値
    if(stat /= nf90_noerr) call handle_err(stat) 

!   通常降雨

  stat = nf90_put_var(ncid,var_prid,pr,&
                      & start=(/1,1,1/)) !通常降雨 地表面
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_prid,pr_2,&
                      & start=(/1,1,2/)) !通常降雨 地上から2km
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_prid,pr_4,&
                      & start=(/1,1,3/)) !通常降雨 地上から4km   
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_prid,pr_6,&
                      & start=(/1,1,4/)) !通常降雨 地上から6km
    if(stat /= nf90_noerr) call handle_err(stat) 

!  対流性降雨
!
  stat = nf90_put_var(ncid,var_cprid,rhValues,&
                      & start=(/1,1,1/)) !対流性降雨 地表面はなし
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_cprid,cpr_2,&
                      & start=(/1,1,2/)) !対流性降雨 地上から2km
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_cprid,cpr_4,&
                      & start=(/1,1,3/)) !対流性降雨 地上から4km   
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_cprid,cpr_6,&
                      & start=(/1,1,4/)) !対流性降雨 地上から6km
    if(stat /= nf90_noerr) call handle_err(stat) 

!  層状性降雨
!
  stat = nf90_put_var(ncid,var_sprid,rhValues,&
                      & start=(/1,1,1/)) !対流性降雨 地表面はなし
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_sprid,spr_2,&
                      & start=(/1,1,2/)) !層状降雨 地上から2km
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_sprid,spr_4,&
                      & start=(/1,1,3/)) !層状降雨 地上から4km   
    if(stat /= nf90_noerr) call handle_err(stat) 
  stat = nf90_put_var(ncid,var_sprid,spr_6,&
                      & start=(/1,1,4/)) !層状降雨 地上から6km
    if(stat /= nf90_noerr) call handle_err(stat) 



!    作成したnetCDFファイルを閉じる

  stat = nf90_close(ncid)
  if(stat /= nf90_noerr)    call handle_err(stat) 


!    エラー表示用のサブルーチン部分. 杉山氏作成.

contains

  subroutine handle_err(status)
    integer, intent (in) :: status
    if(status /= nf90_noerr) then
      print *, trim(nf90_strerror(status))
      stop "Stopped"
    end if
  end subroutine handle_err

end program test_netcdf
