Fortran　NetCDFファイル　読み込み

本記事では、FortranでNetCDFファイルを読み込む方法を紹介します。

NetCDFファイルとは

NetCDFとは、バイナリファイルフォーマットの一つであり、その拡張子は”.nc”である。 NetCDFは気象、海洋、気候変動などの分野で国際的に広く使われている。 Wikipedia参照

また、NetCDFファイルは、ncdumpコマンドを使うことで、どのようなデータが格納されているのか調べることができます。

2021年10月24日 ncdumpコマンド　使い方

使用するデータ

さっそくコードをNetCDFファイルを読み込むコードを見ていきます。

例では、このような4次元データを読み込んでいます。

netcdf PRES_0020.pe000000 {
dimensions:
        z = 75 ;
        x = 512 ;
        y = 512 ;
        time = UNLIMITED ; // (21 currently)
variables:
        float PRES(time, z, y, x) ;
                PRES:long_name = "pressure" ;
                PRES:units = "Pa" ;
                PRES:standard_name = "air_pressure" ;
                PRES:_FillValue = -9.9999e+30f ;
                PRES:missing_value = -9.9999e+30f ;
                PRES:coordinates = "lon lat height time" ;
                PRES:cell_measures = "volume: cell_volume" ;
                PRES:grid = "grid" ;
                PRES:location = "face" ;

上記の結果は、ncdumpコマンドを用いています。

コード

program netcdf_read
use netcdf
implicit none
!!!!read
character(10),parameter::var_pres4D="PRES"
character(len=*),parameter::FILE_pres4D = "./"//trim(var_pres4D)//"_0020.pe000000.nc"

integer::ncid
integer,parameter::ndims = 4
integer::xx_dim,yy_dim,zz_dim,tt_dim
real,allocatable::xx(:),yy(:),zz(:),tt(:)
integer::var_pres
real,allocatable::press(:,:,:)
integer::start4D(ndims),count4D(ndims)
integer::x,y,z,t
real::var_mid
integer::err

call check( nf90_open(FILE_pres4D,nf90_nowrite,ncid) )
call get_dimension(ncid,"x",xx_dim,xx)
call get_dimension(ncid,"y",yy_dim,yy)
call get_dimension(ncid,"z",zz_dim,zz)
call get_dimension(ncid,"time",tt_dim,tt)

allocate(press(xx_dim,yy_dim,zz_dim))
call check( nf90_inq_varid( ncid, var_pres4D, pres ) )
count3D = (/ xx_dim, yy_dim, zz_dim, 1 /)
start3D = (/ 1, 1, 1, 1 /)
do t = 1, tt_dim
   start4D(4) = t
   call check( nf90_get_var( ncid, var_pres, press, start = start4D, count = count4D ))
   do z = 1, zz_dim
   do y = 1, yy_dim
   do x = 1, xx_dim
!!!処理
   enddo
   enddo
   enddo
enddo
call check( nf90_close(ncid) )
contains
subroutine get_dimension( ncid, name, dim, dims )
  integer,         intent( in)::ncid
  character(len=*),intent( in)::name
  integer,         intent(out)::dim
  real,allocatable, intent(out)::dims(:)

  integer::varid,dimid

  call check( nf90_inq_dimid( ncid, name, dimid ) )
  call check( nf90_inquire_dimension( ncid, dimid, len=dim ) )
  allocate(dims(dim), stat=error)
  if ( error /= 0 ) print *,name,":Allocation request denied"
  call check( nf90_inq_varid( ncid, name, varid ))
  call check( nf90_get_var( ncid, varid, dims ))
end subroutine get_dimension

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

FortranでNetCDFファイルを読み込むと、書くコード量が多くなります。

解説

program netcdf_read
use netcdf
implicit none

use netcdf でnetcdfライブラリを使います。

call check( nf90_open(FILE_pres4D,nf90_nowrite,ncid) )
call get_dimension(ncid_pres4D,"x",xx_dim,xx)
call get_dimension(ncid_pres4D,"y",yy_dim,yy)
call get_dimension(ncid_pres4D,"z",zz_dim,yy)
call get_dimension(ncid_pres4D,"time",tt_dim,tt)

nf90_open関数で NetCDF ファイルを開きます。

nf90_nowrite としているの読み込み専用を意味します。

これでファイルが開かれ、ncid にファイルと対応した識別子が代入されます。

get_dimension関数でNetCDFファイルの各次元の長さを抽出しています。

”x”は、今回のNetCDFファイルでは、512ですので、xx_dimに512が代入されます。

”y”は、今回のNetCDFファイルでは、512ですので、yy_dimに512が代入されます。

"z"は、今回のNetCDFファイルでは、75ですので、zz_dimに75が代入されます。

”time”は、今回のNetCDFファイルでは、21ですので、tt_dimに21が代入されます。

allocate(press(xx_dim,yy_dim,zz_dim))

pressの配列に先ほど取得したx方向y方向z方向の長さを代入しています。

call check( nf90_inq_varid( ncid, var_pres4D, var_pres ) )

nf90_inq_varid関数で変数の識別子を取得します。

count4D = (/ xx_dim, yy_dim, zz_dim, 1 /)
start4D = (/ 1, 1, 1, 1 /)

start4Dとcount4Dで変数の次元毎に読み込む開始と一回に読み込む量を取得します。

今回は、x=1,y=1,z=1,t=1から開始し、一回のループでx方向に512,y方向に512,z方向に75のデータを読み込みます。

do t = 1, tt_dim
   start4D(4) = t
   call check( nf90_get_var( ncid, var_pres, sfc_press, start = start4D, count = count4D ))
   do z = 1, zz_dim
   do y = 1, yy_dim
   do x = 1, xx_dim
!!!処理
   enddo
   enddo
   enddo
enddo

一気にデータを読み込むとメモリを多量に使うことから、次元timeで繰り返し処理をして、値を取得しています。

start4Dの4つ目の要素をtにすることで、その時刻のデータを取り出すことができます。

実行

このまま、gfortranでコンパイルすることはできません。

use netcdf
      1
Fatal Error: Can't open module file 'netcdf.mod' for reading at (1): No such file or directory

このようなエラーメッセージが出てしまいます。

netcdf.modとlibnetcdff.aをインストールする必要があります。

この方法は、こちらのサイトにて詳しく書かれています。

$ gfortran ファイル名.f90 -I INCLUDE -L LIBROOT -lnetcdff

netcdf.modがある場所をINCLUDE、libnetcdff.aがある場所をLIBROOTとします。

これで実行できます。

また、毎回、このコマンドを打つのはしんどいので、シェルスクリプトを書いて、そのファイルをコマンドとして使いましょう。

#!/bin/bash
LIBROOT=-L/usr/local/lib64
INCLUDE=-I/usr/local/include
LIBS=-lnetcdff
echo gfortran $@ $INCLUDE $LIBROOT $LIBS
gfortran $@ $INCLUDE $LIBROOT $LIBS

LIBROOTとINCLUDEは、環境に合わせて変えてください。

このファイルをgfortncとすると、