netcdfをfortranで読む

これまでのやり方でGrADSやGMTで読めるようにはなりましたが、自分でデータを好きなように加工する場合はやはりfortranで読みたくなります。

準備として、netcdfライブラリのインストール

まず、netcdfのライブラリをインストールする必要があります。

netcdfのパッケージをダウンロードし、（例netcdf-4.3.0.tar.gz）、インストールします。

やり方は、概ね以下のとおりです。

gunzip netcdf-4.3.0.gz
tar xvf netcdf-4.3.0.tar
cd netcdf-4.3.0
su
./configure --prefix=/usr/local/netcdf-4.3.0/           (ひょっとしたら--disable-netcdf4などいるかも)
make 
mkdir /usr/local/netcdf-4.3.0
make install

とし、/etc/profileに

PATH="$PATH:/usr/local/netcdf-4.3.0/bin/"
LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/netcdf-4.3.0/lib"

を追加。

ようやく本編です。GrADSやGMTのときと同じように、

ncdump -h MSM2013102421P.nc

としてヘッダー情報を読み取ります。これをもとにctlファイルを作ります。
たとえば、MSMのPressure level（x,y,t,z方向に4次元持っています）のデータをncdump -hすると、長いファイルが出てきますが、最初の方を見ると、

netcdf MSM2013102421P {
dimensions:
	lon = 241 ;
	lat = 253 ;
	p = 16 ;
	time = 12 ;
variables:

と書いてあります。これはすなわちx軸[lon]に241個、y軸[lat]に253個、z軸[p](鉛直)に16個、t軸[time]に12個の長さを持っていることを示しています。
さらに下の方を読んでいくと、

variables:
(中略)
	float z(time, p, lat, lon) ;
		z:long_name = "geopotential height" ;
		z:units = "m" ;
		z:standard_name = "geopotential_height" ;
	float w(time, p, lat, lon) ;
		w:long_name = "vertical velocity in p" ;
		w:units = "Pa/s" ;
		w:standard_name = "lagrangian_tendency_of_air_pressure" ;
（略）

などと書いてあります。これは、time,p,lat,lonの4つの軸で定義される変数[z](ジオポテンシャル高度)が、float、すなわち単精度実数で格納されていることを示しています。この様な変数がたくさん並んでいます。

サンプルプログラムは、

!
! to compile this fortran program, 
! gfortran readfort92.f90 -L/usr/local/netcdf-4.3.0/lib -lnetcdf -L/usr/lib -lnetcdff

program netcdf_example
  implicit none
!  include '/usr/local/netcdf-3.6.0/include/netcdf.inc'
  include '/usr/include/netcdf.inc'

  integer, parameter :: &
nx=241, &
ny=253, &
nz=16 , &
nt=12 ! 次元長
  integer :: ncid, status, rhid
  real(4) :: hz(nx,ny,nz,nt),w(nx,ny,nz,nt)  ! データ
  character(30) :: filename, var,var2
 integer :: i,j,k,l,m,n

  var='z' ! 変数設定名。''にはnetcdfに入っている変数名を書く。
!var2='w' !二つ以上読む場合は並べます
  filename='MSM2013102421P.nc'

  status = nf_open(filename, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得
  status = nf_inq_varid(ncid, var, rhid) ! 変数のID(rhid)取得
  status = nf_get_var_real(ncid, rhid, hz)! 変数読み込み。3個目の引数には自分で変数名を設定できる。
!  status = nf_inq_varid(ncid, var2, rhid) ! 
!  status = nf_get_var_real(ncid, rhid, w)! 複数読む場合は並べます
  status = nf_close(ncid) ! ファイルのclose
! do l=1,nt
! do k=1,nz
! do j=1,ny
! do i=1,nx
! write(*,*) i,j,k,l,hz(i,j,k,l)
! end do
! end do
! end do
! end do
  stop
end program

です。コンパイルする時には、gfortran環境なら、

gfortran readnetcdf.f90 -L/usr/local/netcdf-4.3.0/lib -lnetcdf -L/usr/lib -lnetcdff

でうまく行きました。インストールした場所に応じてオプションをうまく指定してください。
やりやすいミスとして、fortranで定義する配列hz(xsz,ysz,zsz,tsz)の配列x,y,z,tの順序は、ncdumpしたときにでてくる

float h(time, p, lat, lon)

の並び順t,z,y,xの逆にしなければなりません。GrADSとちょうど逆ですね。

もう1つの例として、world ocean atlasによる9月のSST。

#!/bin/sh
file=temperature_monthly_1deg.nc
region=225/-40/45/-40r
proj=S0/-90/7i
#region=0/360/-90/90
#region=120/150/20/50
#proj=q210/35/1:125000000
#proj=q135/0/1:20000000
out=test.eps
fortran=readotemp.f90
time=2
level=1
echo $times
echo program netcdf_example >$fortran
echo  implicit none >>$fortran
echo  include "'"/usr/include/netcdf.inc"'" >>$fortran
echo  integer, parameter :: nx=360, ny=180, nz=24 ,nt=12  >>$fortran
echo  integer :: ncid, status, rhid >>$fortran
echo " real(4) :: otemp(nx,ny,nz,nt)"  >>$fortran
echo " real(4) :: tmp1 " >>$fortran
echo character :: filename*60 >>$fortran
echo "character(20) :: var1,var2,var3,var4" >>$fortran
echo integer :: i,j,k,l,m,n >>$fortran
echo var1="'"t_an"'" >>$fortran
echo  filename="'"$file"'" >>$fortran
echo " status = nf_open(filename, nf_nowrite, ncid)">>$fortran
echo " status = nf_inq_varid(ncid, var1, rhid)" >>$fortran
echo " status = nf_get_var_real(ncid, rhid, otemp)">>$fortran
echo " status = nf_close(ncid)" >>$fortran
echo "open(91,file="'"otemp.dat"'")" >>$fortran
echo "open(92,file="'"land.dat"'")" >>$fortran
echo do j=1,ny >>$fortran
echo do i=1,nx >>$fortran
echo "tmp1=otemp(i,j,"$level","$time")" >>$fortran
echo "if(tmp1>100.0) then" >>$fortran
echo "write(92,'(a,f10.2)') "'"> -Z "'",-999.0" >>$fortran
echo "write(92,*) (i-1)*1.0,(j-91)*1.0" >>$fortran
echo "write(92,*)     i*1.0,(j-91)*1.0" >>$fortran
echo "write(92,*)     i*1.0,(j-90)*1.0" >>$fortran
echo "write(92,*) (i-1)*1.0,(j-90)*1.0" >>$fortran
echo "write(92,*) (i-1)*1.0,(j-91)*1.0" >>$fortran
echo else >>$fortran
echo "write(91,'(a,f10.2)') "'"> -Z "'",tmp1" >>$fortran
echo "write(91,*) (i-1)*1.0,(j-91)*1.0" >>$fortran
echo "write(91,*)     i*1.0,(j-91)*1.0" >>$fortran
echo "write(91,*)     i*1.0,(j-90)*1.0" >>$fortran
echo "write(91,*) (i-1)*1.0,(j-90)*1.0" >>$fortran
echo "write(91,*) (i-1)*1.0,(j-91)*1.0" >>$fortran
echo end if >>$fortran
echo end do >>$fortran
echo end do >>$fortran
echo  stop>>$fortran
echo end program>>$fortran
gfortran $fortran -L/usr/local/netcdf-4.3.0/lib -lnetcdf -L/usr/lib -lnetcdff
./a.out
cpt=tmp.cpt
gmtset ANNOT_FONT_SIZE_PRIMARY 20
makecpt -T-2/6/0.25 -Chaxby > $cpt
gray=179
scalebar=18/8.9/17/0.5
echo 0 -90 | psxy -J$proj -R$region -G219/219/219 -Ss32 -K -X2 -Y2>$out
psxy -J$proj -R$region otemp.dat -C$cpt -M -L -K -O >> $out
psxy -J$proj -R$region land.dat -G$gray/$gray/$gray -M -L -K -O -BNWSeg60a60/g10a10>> $out
psscale -D$scalebar -C$cpt -K -O -N -E -B1>>$out

back