FortranでNetCDFを読み書きする

FortranでNetCDFを読み書きする。

# 基本

### netCDFとは

自己記述的なデータフォーマットの一種。気象・気候分野の配布データとしてよく使われている。
Fortran, C, C++, pythonなどのAPIが用意されている
ここではFortran90での読み書きを対象とする。
[netcdf-fortranの本家サイト](https://www.unidata.ucar.edu/software/netcdf/docs-fortran/)

### インストール
ソースからコンパイルしても良いし、apt-getやanacondaでインストールしてもよいと思う。

### コンパイル
インストールするとnc-configコマンドも同時にインストールされる。
パスの指定等はこのコマンド用いると便利。

```bash
　 $ `nc-config --fc` ncread.f90 `nc-config --fflags` `nc-config --flibs`
```

### netCDFモジュールを読む。

```fortran
use netcdf
```

### netCDFの関数の引数を受け取るsubroutineを用意しておく。

```fortran:test.f90
program main 
use netcdf
implicit none
!! main code
contains
subroutine check(status)
 integer, intent (in) :: status
 if(status /= nf90_noerr) then 
   print *, trim(nf90_strerror(status))
   stop "Stopped"
 end if
end subroutine check
end program main
```

# 読む

基本的にncid：netcdfのファイル番号。
varid:netcdf内の変数の番号
を指定して読めば良い。

## ncidを得る 。

```fortran
character(len=1024) :: filename
integer :: ncid
call check( nf90_open(filename, nf90_nowrite, ncid) )
```
## varid を得る（要ncid)

```fortran
integer :: ncid, varid
call check ( nf90_inq_varid( ncid, "varname", varid) )
```
## 変数を読む。

```fortran
integer :: start_nc(4), count_nc(4)
integer :: nx,ny,nz
start_nc=[1,1,1,1]
count_nc=[nx,ny,nz,1]
! data に読み取ったデータを入れる
! dataの型に注意
call check ( nf90_get_var(ncid, varid, data, start=start_nc, count=count_nc) )
```

ここでは読み取るデータの型に注意する。
ncdump -h などで事前に必要なデータの型を確認すること
| ncdump での変数 | Fortranの変数 |
|:-:|:-:|
|double   | real(8)  |
|float   | real(4)  |
|int   | integer(4)  |
|short   |integer(2)   |

となる。

また、多次元配列の場合ncdump等のコマンドで表示されるものと、逆に読まないといけない。これはFortranの配列アクセスが列優先で、c言語などと異なるからである。
例えばncdump で表示されるのが、

``` 
float air(time, lat, lon) ;
```

である場合には

``` Fortran
integer,parameter :: ntime=12, nlat=73, nlon=144
real(4) :: air(nlon,nlat,ntime)
```

のようにする。ntime,nlat,nlonはそれぞれtime,lat,lonの要素数


### さらに、次元の情報を得るやり方。

```fortran 
  ! 次元の大きさを得る
  call check( nf90_inq_dimid(ncid, 'lon', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nx) )
```

## サンプルコード
対象は経度、緯度、高度（気圧）、時間の4次元あるデータを読み、各気圧層の最大最小を出力するプログラム。
サンプルデータは以下から得た。
[https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html]

```bash:サンプルデータの中身を確認する
$ ncdumo -h air.mon.nc
```

### サンプルコード本体

```fortran:read.f90 
program read
  use netcdf
  implicit none
  ! dimension
  integer :: nx, ny, nz
  character( len=1024 ) :: f_in
  real(4), allocatable :: data(:,:,:)
  ! for netcdf
  integer :: ncid, varid, dimid
  integer :: start_nc(4), count_nc(4)
  ! doloop
  integer :: k, n
  !
  !
  f_in = "air.mon.mean.nc"
  call check( nf90_open( trim( f_in ), nf90_nowrite, ncid) )
  ! 次元の大きさを得る
  call check( nf90_inq_dimid(ncid, 'lon', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nx) )
  call check( nf90_inq_dimid(ncid, 'lat', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=ny) )
  call check( nf90_inq_dimid(ncid, 'level', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nz) )
  allocate( data(nx,ny,nz) )
  ! varid を得る
  call check( nf90_inq_varid(ncid, "air", varid))
  n = 1
  ! n-step目の変数を読む
  start_nc = [1,1,1,n] 
  count_nc = [nx,ny,nz,1]
  call check( nf90_get_var(ncid, varid, data, start=start_nc, count=count_nc) )
  do k = 1, nz
     write(*,*) k, maxval(data(:,:,k)), minval(data(:,:,k))
  end do
  ! file を閉じる
  call check( nf90_close(ncid) )
contains
  subroutine check( status )
    integer, intent (in) :: status
    if(status /= nf90_noerr) then 
       print *, trim(nf90_strerror(status))
       stop "Stopped"
    end if
  end subroutine check
end program read
```


# 書き込む

## 書き込むファイルを開く(ncidを得る)。
特にこだわりがなければnetcdf4形式で良いと思う。

```fortran
  call check( nf90_create( trim(f_out), NF90_HDF5, ncido) ) 
```

## 変数を定義する(define mode)

### 次元の定義
次元を定義しその次元のdimidを得る。次の変数の定義で用いる。

```fortran
call check( nf90_def_dim(ncido, 'lat', ny, lat_dim_id) )
! 1次元のみNF90_UNLIMITEDとして次元長を指定しなくても良い。RECORD次元という。
call check( nf90_def_dim(ncido, 'time', NF90_UNLIMITED, time_dim_id) )
```

### 変数の定義

```fortran
 ! 1次元量の場合
 call check( nf90_def_var(ncido, 'lat', NF90_REAL, lat_dim_id, lat_id) )
 !　多次元量は配列でdimidを指定する。
! deflate_levelを指定すると圧縮がかかる。netcdf4のみ
 call check( nf90_def_var(ncido, "air", NF90_REAL, (/ lat_dim_id, level_dim_id, time_dim_id/), &
& varido, deflate_level = 1 ) )
```

### 変数にattributionを付ける
attribution(属性)を付ける。変数の説明(単位、名前 etc.）を付加できる。
ncid、 varidを指定し書き込む。下の例は"units" として 'degrees_north'を書き込む例。

```fortran
call check( nf90_put_att(ncido, lat_id, 'units', 'degrees_north') )
```

### 変数定義モードを終える。

```fortran
  call check( nf90_enddef(ncido) )
```

## 変数を書き込む。
ncidとvaridを指定して、変数を入れれば良い。
時間ステップ毎に入れる際などはstartとcountを指定する。

```fortran
 call check( nf90_put_var(ncido, lat_id, lat ) )
 start_nco = [1,1,n]
 count_nco = [ny,nz,1]
 call check( nf90_put_var(ncido, varido, out, start=start_nco, count=count_nco) )
```
## 書き込むファイルを閉じる。

```fortran
 call check( nf90_close(ncido) )
```
### サンプルコード
対象は経度、緯度、高度（気圧）、時間の4次元あるデータを読み、経度方向に平均をとって、緯度、高度（気圧）、時間の3次元のデータに変換する。

```fortran:write.f90
program read
  use netcdf
  implicit none
  ! dimension
  integer :: nx, ny, nz, nt
  ! file
  character( len=1024 ) :: f_in, f_out
  real(4), allocatable :: data(:,:,:)
  real(4), allocatable :: out(:,:)
  real(4), allocatable :: lat(:), level(:)
  real(8) :: time(1)
  ! for netcdf input
  integer :: ncid, varid, dimid
  integer :: start_nc(4), count_nc(4)
  ! for netcdf output
  integer :: ncido, varido
  integer :: lat_dim_id, level_dim_id, time_dim_id
  integer :: lat_id, level_id, time_id
  integer :: start_nco(3), count_nco(3)
  ! doloop
  integer :: k, n
  !
  !
  f_in = "air.mon.mean.nc"
  f_out = "air.mon.zonmean.nc"
  !
  call check( nf90_open( trim( f_in ), nf90_nowrite, ncid) )
  ! 次元の大きさを得る
  call check( nf90_inq_dimid(ncid, 'lon', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nx) )
  call check( nf90_inq_dimid(ncid, 'lat', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=ny) )
  allocate(lat(ny))
  call check( nf90_get_var(ncid, dimid, lat) )
  call check( nf90_inq_dimid(ncid, 'level', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nz) )
  allocate(level(nz))
  call check( nf90_get_var(ncid, dimid, level) )
  call check( nf90_inq_dimid(ncid, 'time', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nt) )
  ! 時間次元だけ読み込む
  !
  allocate( data(nx,ny,nz) )
  ! varid を得る
  call check( nf90_inq_varid(ncid, "air", varid))
  ! 読む設定はここまで。
  ! 書く設定
  ! 
  call check( nf90_create( trim(f_out), NF90_HDF5, ncido) )
  ! 次元を定義する
  call check( nf90_def_dim(ncido, 'lat', ny, lat_dim_id) )
  Call check( nf90_def_dim(ncido, 'level', nz, level_dim_id) )
  call check( nf90_def_dim(ncido, 'time', NF90_UNLIMITED, time_dim_id) )
  ! 変数を定義する
  call check( nf90_def_var(ncido, 'lat', NF90_REAL, lat_dim_id, lat_id) )
  call check( nf90_def_var(ncido, 'level', NF90_REAL, level_dim_id, level_id) )
  call check( nf90_def_var(ncido, 'time', NF90_DOUBLE, time_dim_id, time_id) )
  call check( nf90_def_var(ncido, "air", NF90_REAL, (/ lat_dim_id, level_dim_id, time_dim_id/),&
       & varido, deflate_level = 1 ) )
  ! 変数に attributionをつける　
  call check( nf90_put_att(ncido, lat_id, 'units','degrees_north') )
  call check( nf90_put_att(ncido, level_id, 'units','millibar') )
  call check( nf90_put_att(ncido, time_id, 'units', 'hours since 1800-01-01 00:00:0.0' ) )
  ! define モードを終了する。
  call check( nf90_enddef(ncido) )
  call check( nf90_put_var(ncido, lat_id, lat ) )
  call check( nf90_put_var(ncido, level_id, level ) )
  allocate( out(ny,nz) )
  !
  do n = 1,nt 
  ! n-step目の変数を読む
     start_nc = [1,1,1,n]
     count_nc = [nx,ny,nz,1]
     call check( nf90_get_var(ncid, varid, data, start=start_nc, count=count_nc) )
     call check( nf90_get_var(ncid, dimid, time, start=[n], count=[1]) )
     out = sum(data,dim=1)/nx
     start_nco = [1,1,n]
     count_nco = [ny,nz,1]
     call check( nf90_put_var(ncido, varido, out, start=start_nco, count=count_nco) )
     call check( nf90_put_var(ncido, time_id, time, start=[n], count=[1]) )
  end do
  ! file を閉じる
  call check( nf90_close(ncid) )
  call check( nf90_close(ncido) )
  ! deallocate
  deallocate( data, out, lat, level )

contains

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

end program read
```

#おまけ

上記はcdo(Climate Data Operators)を使うと一瞬である。
cdoのページ[https://code.mpimet.mpg.de/projects/cdo/]
気候の計算で使うような処理の多くはこのコマンドで扱える。

```bash
$ cdo zonmean air.mon.mean.nc air.mon.zonmean.nc
```