fortran/NetCDFファイルを読む

NetCDF ファイルを Fortran で読む †

他の言語では NetCDF 用のパッケージを使えば簡単に変数を読み込めるし， Fortran も余裕…と思っていたけど甘くなかった．
例として，生存圏研究所のデータベースにある気象庁MSMデータの NetCDF を使う．
netcdf-c と netcdf-fortran がインストールされている前提．Ubuntu でのインストールについては NetCDFのインストールに書いた．
ここで扱うのは NetCDF4 の話で，旧バージョンの NetCDF3 については全くわからない．

↑

コンパイル †

まずは nf-config や nc-config が実行できるように，環境変数 PATH に NetCDF のインストールパスを追加する．
コンパイルはコマンドをうつのが面倒くさいので Makefile に書いておく．

FC = $(shell nf-config --fc) 
FFLAGS = $(shell nf-config --fflags)
FLIBS = $(shell nf-config --flibs)
EXE = xread_msmnc
FILE = test_read.f90

.PHONY: exe clean

exe:
	$(FC) $(FILE) $(FFLAGS) $(FLIBS) -o $(EXE)

clean:
	-rm -f $(EXE)

これを test_read.f90 と同じディレクトリに入れておけば make をうつだけで実行ファイルが生成される．

↑

Fortran ソース †

MSM の .nc ファイルから lon, lat や海面更正気圧 psea などを読み取る．
引数で .nc ファイルのファイル名を指定する設定にしている (iargc のところ)．
ほとんど下記参考欄にあるリンク先のコードを参考にして書いた．

program read_msmnc
  use netcdf
  implicit none

  ! dimension
  integer :: nx, ny, nt
  character(len=1024) :: f_in
  real(kind=4), allocatable :: lon(:), lat(:), timelap(:)
  real(kind=4), allocatable :: psea(:,:,:)
  real(kind=4), allocatable :: u10(:,:,:), v10(:,:,:)
  ! parameters
  integer :: ncid, varid, dimid
  integer :: start_nc(3), count_nc(3)
  real(kind=4) :: scale_factor, add_offset
  ! loop counter
  integer :: k

  ! netcdf filename
  if(iargc()>0)then
    call getarg(1,f_in)
  else    
    stop 'No ncfile specified'
  end if
 
  ! open
  call check_ncstatus( nf90_open( trim( f_in ), nf90_nowrite, ncid) )

  ! number of array
  ! -- lon
  call check_ncstatus( nf90_inq_dimid(ncid, 'lon', dimid) )
  call check_ncstatus( nf90_inquire_dimension(ncid, dimid, len=nx) )
  allocate(lon(nx))
  call check_ncstatus( nf90_get_var(ncid, dimid, lon) )
  ! -- lat
  call check_ncstatus( nf90_inq_dimid(ncid, 'lat', dimid) )
  call check_ncstatus( nf90_inquire_dimension(ncid, dimid, len=ny) )
  allocate(lat(ny))
  call check_ncstatus( nf90_get_var(ncid, dimid, lat) )
  ! -- timelap
  call check_ncstatus( nf90_inq_dimid(ncid, 'time', dimid) )
  call check_ncstatus( nf90_inquire_dimension(ncid, dimid, len=nt) )
  allocate(timelap(nt))
  call check_ncstatus( nf90_get_var(ncid, dimid, timelap) )

  ! allocate
  allocate(psea(nx,ny,nt))
  allocate(u10(nx,ny,nt), v10(nx,ny,nt))

  ! indices
  start_nc = [1,1,1]
  count_nc = [nx,ny,nt]

  ! read variables
  ! -- psea
  call check_ncstatus( nf90_inq_varid(ncid, "psea", varid) )
  call check_ncstatus( nf90_get_var(ncid, varid, psea, start=start_nc, count=count_nc) )
  call check_ncstatus( nf90_get_att(ncid, varid, "scale_factor", scale_factor) )
  call check_ncstatus( nf90_get_att(ncid, varid, "add_offset", add_offset) )
  psea(:,:,:) = psea(:,:,:)*scale_factor + add_offset
  ! -- u10
  call check_ncstatus( nf90_inq_varid(ncid, "u", varid) )
  call check_ncstatus( nf90_get_var(ncid, varid, u10, start=start_nc, count=count_nc) )
  call check_ncstatus( nf90_get_att(ncid, varid, "scale_factor", scale_factor) )
  call check_ncstatus( nf90_get_att(ncid, varid, "add_offset", add_offset) )
  u10(:,:,:) = u10(:,:,:)*scale_factor + add_offset
  ! -- v10
  call check_ncstatus( nf90_inq_varid(ncid, "v", varid) )
  call check_ncstatus( nf90_get_var(ncid, varid, v10, start=start_nc, count=count_nc) )
  call check_ncstatus( nf90_get_att(ncid, varid, "scale_factor", scale_factor) )
  call check_ncstatus( nf90_get_att(ncid, varid, "add_offset", add_offset) )
  v10(:,:,:) = v10(:,:,:)*scale_factor + add_offset

  ! close nc file
  call check_ncstatus( nf90_close(ncid) )

  ! --- print and check
  write(*,*) nx, ny, nt
  write(*,*) timelap(lbound(timelap)), timelap(ubound(timelap))
  write(*,*) minval(lon), maxval(lon), minval(lat), maxval(lat)
  do k = 1, nt
    write(*,*) k, maxval(psea(:,:,k)), minval(psea(:,:,k)), &
  &               maxval(u10(:,:,k)), minval(u10(:,:,k)), maxval(v10(:,:,k)), minval(v10(:,:,k))
  end do

  contains
! ------------------------------------------------------------------------------
  subroutine check_ncstatus( status )
    integer, intent (in) :: status
    if(status /= nf90_noerr) then 
       print *, trim(nf90_strerror(status))
       stop "Something went wrong while reading ncfile."
    end if
  end subroutine check_ncstatus
! ------------------------------------------------------------------------------
end program read_msmnc

基本的には nf90_inq_varid で変数名に対応する id を取得して， nf90_get_var で読むという流れ．
nf90_* 系の関数は全て実行ステータスを返してくれるので，失敗していないか subtoutine check_ncstatus を使って都度チェックを行っている.

Python や Julia のパーサでは勝手に scale_factor や offset を処理した上で読んでくれるが， Fortran では自力で処理するコードを書く必要があるみたい．
また，Fortran ソース内で宣言した変数の方が real(kind=4) や real(kind=8) で，NetCDF の中に格納されている変数の型が integer(kind=4) などであっても，そこはうまいことやってくれる(エラーが出たり変な値にはならない)らしい．

↑

実行と確認 †

Makefile に書いた通り xread_msmnc という実行ファイルができたので，ファイル名を引数に指定して，読めたかどうかテスト．

$ ./xread_msmnc  MSM2021111100S.nc
         481         505          34
   0.00000000       33.0000000    
   120.000000       150.000000       22.3999996       47.5999985    
           1   102380.734       98720.1797       18.3058109      -7.87767601       15.5412846      -13.8960247    
           2   102347.703       98708.7188       18.3119278      -7.69419003       16.2874622      -14.1896029    
           3   102314.219       98696.7891       18.5932732      -7.63302755       15.7798166      -14.3425083    
           4   102247.703       98677.5234       19.1314983      -7.31498480       16.0856266      -13.6819572    
           5   102155.047       98647.2500       19.0336399      -7.00917435       15.1681967      -12.8256884    
           6   102072.016       98661.0078       18.7217140      -7.03975582       14.6605511      -11.7614679    
           7   102021.562       98694.9531       18.5382271      -7.19877720       14.2629976      -11.4617739    
           8   102016.516       98707.3359       18.2752304      -7.60244656       14.1284409      -12.2935781    
           9   102048.164       98722.9375       18.5382271      -8.39143753       14.3792057      -12.7461777    
          10   102117.891       98858.7188       17.4984722      -9.30886841       14.2935781      -13.0214071    
          11   102174.773       98965.1406       17.3577995      -10.3302755       15.0764532      -12.9602451    
          12   102246.789       99013.7578       17.0825691      -11.0275230       15.4311934      -13.3516827    
          13   102270.188       99098.1641       16.7645264      -11.2660551       14.5871563      -14.3302755    
          14   102309.633       99187.6172       16.8195724      -11.3577986       14.2446489      -14.9785938    
          15   102355.508       99236.6953       19.7431202      -11.4128447       13.5351686      -14.6360865    
          16   102393.117       99319.2656       20.0795116      -11.5657492       16.6360855      -14.2263002    
          17   102394.953       99373.8516       16.1162090      -11.6513767       16.0550461      -13.8042822    
          18   102404.586       99407.7969       16.3302765      -11.7492361       13.7003059      -13.1743126    
          19   102414.219       99425.6875       17.4311924      -11.9082575       14.5688076      -12.8134565    
          20   102423.391       99457.7969       17.5229359      -12.0550461       14.7155972      -12.7339458    
          21   102444.039       99516.0547       17.8899097      -12.0611629       15.9633036      -12.7461777    
          22   102467.891       99557.7969       18.5871563      -12.0000000       16.6911316      -12.7828751    
          23   102502.750       99610.5469       18.9663620      -12.2446489       17.2905197      -12.8195724    
          24   102567.891       99638.0703       18.8807354      -12.4036703       16.1773701      -12.4648323    
          25   102638.070       99628.8984       17.3516827      -12.5565758       16.7706432      -13.1070337    
          26   102672.477       99633.0312       16.6972485      -12.6238537       16.1467896      -12.9663610    
          27   102647.703       99627.5234       17.0030594      -12.3425083       15.8715601      -12.6666670    
          28   102578.898       99597.2500       17.4556580      -12.0428143       15.2354746      -12.0305815    
          29   102489.453       99599.0859       16.8868504      -11.6207952       13.1314993      -12.0550461    
          30   102404.125       99616.0547       17.1987782      -11.3088684       13.2966366      -11.9388380    
          31   102348.164       99651.8359       17.1070347      -11.2599392       12.4097862      -11.7370033    
          32   102325.688       99692.6641       17.2966366      -11.4495420       12.4892969      -11.7798166    
          33   102337.156       99728.4375       18.2813454      -11.5535173       14.3119268      -11.6330280    
          34   102387.156       99764.6797       17.7798176      -11.7737007       14.3058109      -11.5474014

できた．気圧や風速の値も妥当な範囲内に収まっている．

↑

参考 †

FortranでNetCDFを読み書きする - Qiita

Takuya Miyashita

MENU

recent(10)

NetCDF ファイルを Fortran で読む †

コンパイル †

Fortran ソース †

実行と確認 †

参考 †