#author("2021-12-05T12:29:03+09:00","default:Miyashita","Miyashita") #author("2021-12-07T23:37:15+09:00","default:Miyashita","Miyashita") * NetCDF ファイルを Fortran で読む [#i6b5c4de] 他の言語では NetCDF 用のパッケージを使えば簡単に変数を読み込めるし， Fortran も 余裕…と思っていたけど甘くなかった．~ 例として，[[生存圏研究所のデータベース>http://database.rish.kyoto-u.ac.jp/index-e.html]]にある気象庁MSMデータの NetCDF を使う．~ netcdf-c と netcdf-fortran がインストールされている前提．Ubuntu でのインストールについては [[NetCDFのインストール>Ubuntu/Install_NetCDF]] に書いた．~ ここで扱うのは NetCDF4 の話で， 旧バージョンの NetCDF3 については全くわからない． ~ ** コンパイル [#y4e3ac03] まずは nf-config や nc-config が実行できるように，環境変数 PATH に NetCDF のインストールパスを追加する．~ コンパイルはコマンドをうつのが面倒くさいので Makefile に書いておく．~ #codeprettify{{ 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 ソース [#icd32bae] MSM の .nc ファイルから lon, lat や 海面更正気圧 psea などを読み取る．~ 引数で .nc ファイルのファイル名を指定する設定にしている (iargc のところ)．~ ほとんど下記参考欄にあるリンク先のコードを参考にして書いた．~ #codeprettify(lang-fortran){{ 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 for check ! --- 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) などであっても，そこはうまいことやってくれる(エラーが出たり変な値にはならない)らしい．~ **実行と確認 [#b129e789] Makefile に書いた通り xread_msmnc という実行ファイルができたので，ファイル名を引数に指定して，読めたかどうかテスト． #codeprettify{{ $ ./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 }} できた．気圧や風速の値も妥当な範囲内に収まっている． ** 参考 [#z37ffca4] -[[FortranでNetCDFを読み書きする - Qiita>https://qiita.com/convection/items/14eea7b4638780fa8d55]]