Linux-Note

http://www.gfd-dennou.org/library/gtool4/gt4f90io/...

http://www.sci.hokudai.ac.jp/grp/poc/top/old/softw...

GrADSのバージョン2.0では、sdfwriteというコマンドでNetCDFファイルの作成が可能。
ただ、2.0は64bit PCにしか入らない。

open cmap_precip.ctl
set x 1 144
set y 1 73
set lev 1000
set t 1 12
set e 1 
a = precip
set sdfwrite precip.nc
sdfwrite a

sdfwrite

set sdfwrite

set sdfattr

http://hydro.iis.u-tokyo.ac.jp/~agata/archive/ALMA...

http://www.rain.hyarc.nagoya-u.ac.jp/laboratory/OB...

基本的に

status=nf90_***(******)

という文を使って変数や次元の設定を行っていく。

 nf90_create:ファイル作成宣言とファイルのオープン
 nf90_def_dim:要素の次元（緯度、経度、時間等）を宣言
 nf90_def_var:データを入力する要素を宣言（実数、整数なども）
 (次元のデータを入れる場合、次元も要素として再宣言）
 nf90_put_att:要素の説明を入力
 nf90_enddef:要素等の宣言、説明の終了
 nf90_put_var:要素ごとにデータの入力
 nf90_close:ファイルを閉じる

コンパイル時に同じディレクトリにNETCDF.modというファイルが必要。
NETCDFのディレクトリのinclude/にあると思われる。
コンパイル時にnetcdfのライブラリを参照する。パスは適宜書き換える。

ifort test.f90 -lnetcdf -L/usr/local/netcdf-3.5.1/lib

./a.out　で実行

入力：../med-res/20c3m/y1980/z{700,500}.grd
出力：TEST.nc

     program main
     use netcdf
 
     implicit none
     integer,parameter::nx=128,ny=33,nts=366,nz=2
     integer::varid(5),lonid,latid,levid,timid
     real(4)::lon(nx),lat(ny),lev(nz),vdat(nx,ny,nz,nts),lat_val(ny)
     integer(4)::tim(nts)
 
     integer::ncid,status,irec,ix,iy,iz,it
 !　緯度軸の設定（等間隔でない場合）
     data lat_val /-87.8638,-85.0965,-82.3129,-79.5256,-76.7369, &
        -73.9475,-71.1578,-68.3678,-65.5776,-62.7874,-59.9970, &
        -57.2066,-54.4162,-51.6257,-48.8352,-46.0447,-43.2542, &
        -40.4636,-37.6731,-34.8825,-32.0919,-29.3014,-26.5108, &
        -23.7202,-20.9296,-18.1390,-15.3484,-12.5578,-9.7671, &
        -6.9765,-4.1859,-1.3953,1.3953/
 
     character(len=50)::FLOUT
 
 !　出力ファイル
     FLOUT='TEST.nc'
 
     tim(1)=0
     do it=2,nts
        tim(it)=it-1.0
     enddo
 !　経度軸の設定（等間隔の場合）
     do ix=1,nx
         lon(ix)=(ix-1)*2.8125+1.40625
     enddo
     do iy=1,ny
         lat(iy)=lat_val(iy)
     enddo
     lev(1)=700
     lev(2)=500
 
 !　入力ファイル
     open(100,file='../med-res/20c3m/y1980/z700.grd',&
           access='direct', form='unformatted', recl=4*nx*ny)
     open(200,file='../med-res/20c3m/y1980/z500.grd',&
           access='direct', form='unformatted', recl=4*nx*ny)
     do it=1,nts
     read(100,rec=it) ((vdat(ix,iy,1,it),&
           ix=1,nx),iy=1,ny)
     read(200,rec=it) ((vdat(ix,iy,2,it),&
           ix=1,nx),iy=1,ny)
     enddo
 
 !------------------------------------------------------------------
      status=nf90_create(trim(FLOUT),nf90_clobber,ncid)
      call err(status)
 !------------------------------------------------------------------
 !     definition of dimensions
 
      status=nf90_def_dim(ncid,"lon",nx,lonid)
      call err(status)
 
      status=nf90_def_dim(ncid,"lat",ny,latid)
      call err(status)
 
      status=nf90_def_dim(ncid,"lev",nz,levid)
      call err(status)
 
      status=nf90_def_dim(ncid,"time",nts,timid)
      call err(status)
 
 !--------------------------------------------------------------------
 !     definition of variables
 
      status=nf90_def_var(ncid,"lon",nf90_float,lonid,varid(1))
      call err(status)
 
      status=nf90_def_var(ncid,"lat",nf90_float,latid,varid(2))
      call err(status)
 
      status=nf90_def_var(ncid,"lev",nf90_float,levid,varid(3))
      call err(status)
 
 !     status=nf90_def_var(ncid,"time",nf90_int,timid,varid(4))
      status=nf90_def_var(ncid,"time",nf90_double,timid,varid(4))
      call err(status)
 
      status=nf90_def_var(ncid,"Z",nf90_float,&
         (/lonid,latid,levid,timid/),varid(5))
      call err(status)
 
 !------------------------------------------------------------------------
 !     putting attributions
 
      status=nf90_put_att(ncid,varid(1),"long_name","longitutde")
      call err(status)
 
      status=nf90_put_att(ncid,varid(1),"units","degree_east")
      call err(status)
 
      status=nf90_put_att(ncid,varid(1),"missing_value",-999.0)
      call err(status)
 
      status=nf90_put_att(ncid,varid(2),"long_name","latitutde")
      call err(status)
 
      status=nf90_put_att(ncid,varid(2),"units","degree_north")
      call err(status)
 
      status=nf90_put_att(ncid,varid(2),"missing_value",-999.0)
      call err(status)
 
      status=nf90_put_att(ncid,varid(3),"long_name","pressure level")
      call err(status)
 
      status=nf90_put_att(ncid,varid(3),"units","hPa")
      call err(status)
 
      status=nf90_put_att(ncid,varid(3),"missing_value",-999.0)
      call err(status)
 
      status=nf90_put_att(ncid,varid(4),"long_name","time")
      call err(status)
 
      status=nf90_put_att(ncid,varid(4),"units",&
         "days since 1980-1-1 00:00:0.0")
      call err(status)
      
      status=nf90_put_att(ncid,varid(4),"delta_t",&
         "0000-00-01 00:00:00")
      call err(status)
 
      status=nf90_put_att(ncid,varid(4),"missing_value",-999.0)
      call err(status)
 
      status=nf90_put_att(ncid,varid(5),"long_name",&
         "geopotential height")
      call err(status)
 
      status=nf90_put_att(ncid,varid(5),"missing_value",-999.0)
      call err(status)
 
      status=nf90_put_att(ncid,nf90_global,"model","MIROC-midres")
      call err(status)
 
 !--------------------------------------------------------------------------
 
      status=nf90_enddef(ncid)
      call err(status)
 
 !-------------------------------------------------------------------------- 
 !     putting data
 
      status=nf90_put_var(ncid,varid(1),lon,(/1/),(/nx/))
      call err(status)
 
      status=nf90_put_var(ncid,varid(2),lat,(/1/),(/ny/))
      call err(status)
 
      status=nf90_put_var(ncid,varid(3),lev,(/1/),(/nz/))
      call err(status)
 
      status=nf90_put_var(ncid,varid(4),tim,(/1/),(/nts/))
      call err(status)
 
      status=nf90_put_var(ncid,varid(5),vdat,(/1,1,1,1/),&
         (/nx,ny,nz,nts/))
      call err(status)
 
      status=nf90_close(ncid)
      call err(status)
 
      END
 
      subroutine err(status) 
      use netcdf 
      integer::status
 
      if (status /= nf90_noerr ) call handle_err (status)
 
      end subroutine err
 
      subroutine handle_err (status)
      use netcdf
      integer,intent( in):: status
 
      if(status /= nf90_noerr) then
      print *, trim(nf90_strerror(status))
      stop "stopped"
      endif
 
      end subroutine handle_err