NCEP/NCAR再解析データをFortranで読む

データのダウンロードサイト: ftp.cdc.noaa.gov

ダウンロード方法

$ ftp -n < getflux.txt

getflux.txtの内容

open ftp.cdc.noaa.gov

user anonymous YOUR_EMAIL_ADDRESS

cd Datasets/ncep.reanalysis.dailyavgs/surface_gauss/

binary

get air.2m.gauss.2000.nc

by

YOUR_EMAIL_ADDRESSには、あなたのメールアドレスを入れてください。

データの書式

ncdumpコマンドを使うと、データの書式を調べることが出来ます。

$ ncdump -h air.2m.gauss.2000.nc

netcdf air.2m.gauss.2000 {

dimensions:

lon = 192 ;

lat = 94 ;

time = UNLIMITED ; // (366 currently)

variables:

float lat(lat) ;

lat:units = "degrees_north" ;

lat:actual_range = 88.542f, -88.542f ;

lat:long_name = "Latitude" ;

lat:standard_name = "latitude" ;

lat:axis = "Y" ;

float lon(lon) ;

lon:units = "degrees_east" ;

lon:long_name = "Longitude" ;

lon:actual_range = 0.f, 358.125f ;

lon:standard_name = "longitude" ;

lon:axis = "X" ;

double time(time) ;

time:units = "hours since 1-1-1 00:00:0.0" ;

time:long_name = "Time" ;

time:actual_range = 17522904., 17531664. ;

time:delta_t = "0000-00-01 00:00:00" ;

time:avg_period = "0000-00-01 00:00:00" ;

time:standard_name = "time" ;

time:axis = "T" ;

short air(time, lat, lon) ;

air:long_name = "mean Daily Air temperature at 2 m" ;

air:unpacked_valid_range = 150.f, 400.f ;

air:actual_range = 178.5f, 316.07f ;

air:units = "degK" ;

air:add_offset = 477.65f ;

air:scale_factor = 0.01f ;

air:missing_value = 32766s ;

air:precision = 2s ;

air:least_significant_digit = 1s ;

air:GRIB_id = 11s ;

air:GRIB_name = "TMP" ;

air:var_desc = "Air temperature" ;

air:dataset = "NCEP Reanalysis Daily Averages" ;

air:level_desc = "2 m" ;

air:statistic = "Mean" ;

air:parent_stat = "Individual Obs" ;

air:valid_range = -32765s, -7765s ;

// global attributes:

:Conventions = "COARDS" ;

:title = "mean daily NMC reanalysis (2000)" ;

:history = "created 00/01/30 by Hoop (netCDF2.3)" ;

:description = "Data is from NMC initialized reanalysis\n",

"(4x/day). It consists of T62 variables interpolated to\n",

"pressure surfaces from model (sigma) surfaces." ;

:platform = "Model" ;

:references = "http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html" ;

}

読み込みプログラム: read_ncep.f90

program read_ncep

! Description:

!

! Author: aym

!

! Host: aofd30

! Directory: /work1/aym/11.Work10/15.HeatBudjet/51.jcope2/32.NCEP_FLUX/11.test_read_flux/src

!

! Revision history:

! 2011-02-22 16:38

! Initial Version

! References

! http://www.rain.hyarc.nagoya-u.ac.jp/laboratory/OB/shimizu/HTML/netcdf.html

! Tips: NetCDF. http://iprc.soest.hawaii.edu/users/sasakiyo/Japanese/misc/tips_netcdf.html

! NetCDF Tips. http://www.sci.hokudai.ac.jp/grp/poc/top/software/other/netcdf_tips/index.htm

! use

! implicit none

include '/usr/local/include/netcdf.inc'

! ncid:ファイルのID番号、 varid: 変 数のID番号

integer ncid,varid

character infle*500,ofle*500

character var*50

real,allocatable :: lon(:),lat(:)

real(8),allocatable :: time(:)

integer(2),allocatable :: sdata(:,:,:)

! real,allocatable :: sdata(:,:,:)

integer,allocatable :: dimids(:),nshape(:)

character*70 err_message

character*100 dimname !dimnameは各次元の名前。

write(*,'(a)')'Program read_ncep starts.'

! write(*,*)''

im=192

jm=94

nm=366

allocate(lon(im),lat(jm),time(nm),sdata(im,jm,nm))

infle="/work1/aym/0.data/02.NCEP_NCAR/01.Surface_Flux/surface_data/air.2m.gauss.2000.nc"

print '(A,A)','Input: ',infle(1:lnblnk(infle))

sta = nf_open(infle, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得

print *,'ncid= ',ncid

print *,'status= ',sta

print *

!もし、openが上手く行かなかったら、sta /= 0となる。

if(sta.ne.0) then

! その時のエラー を教えてくれる。

err_message = nf_strerror(sta)

write(6,*) err_message

endif

var='air'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *

! スケールファクターを得る。(*1)

sta = nf_get_att_real(ncid,varid,'scale_factor',scale)

print *,'scale= ',scale

print *,'sta= ',sta

print *

! オフセットの 値を得る(*2)

sta = nf_get_att_real(ncid,varid,'add_offset',offset)

print *,'offset=', offset

print *,'sta= ',sta

print *

! varidからndimsを得る。ndimsとは次元の数。二次元データなら2

stat = nf_inq_varndims(ncid, varid, ndims)

print *,'ndims= ', ndims

print *

! varidからdimidsを得る。dimidsとはそれぞれの次元に対するID

allocate(dimids(ndims))

stat = nf_inq_vardimid(ncid, varid, dimids)

do i=1,ndims

write(6,'("dimids(",i2,")=",i9 )') i,dimids(i)

enddo

print *

allocate(nshape(ndims))

do i=1,ndims

! dimidsからnshapeを得る。

! nshapeとはそれぞれの次元に対する最大データ数 (格子数)

stat = nf_inq_dimlen(ncid, dimids(i), nshape(i))

write(6,'("nshape(",i2,")=",i9 )') i,nshape(i)

enddo !i

print *

do i=1,ndims

stat = nf_inq_dimname(ncid, dimids(i), dimname)

write(6,'("dimname(",i2,")=",1x,a23)') i,dimname

enddo

print *

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_int2(ncid, varid, sdata)

print *,'stat= ',stat

var='lon'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_real(ncid, varid, lon)

print *,'stat= ',stat

var='lat'

print '(A)',var(1:lnblnk(var))

! ファイルID(ncid)からvarで設定した変数のID(varid)を得る。

sta = nf_inq_varid(ncid,var,varid)

print *,'varid=', varid

print *,'sta= ', sta

print *,'Read ',var(1:lnblnk(var))

stat = nf_get_var_real(ncid, varid, lat)

print *,'stat= ',stat

print *

ofle="output/test.asc"

print '(A,A)','Output: ',ofle(1:lnblnk(ofle))

print *

open(20,file=ofle)

! Note: time index (n) is fixed and set to 1 here.

n=366

do j=1,jm

do i=1,im

varout=sdata(i,j,n)*scale + offset

write(20,'(2f12.6,f12.4,1x,2i5)') &

lon(i),lat(j),varout,i,j

enddo !i

enddo !j

close(20)

write(*,'(a)')'Done program read_ncep.'

write(*,*)

end program read_ncep

! *1： scaleとは、ある変数の平均的なオーダー

! *2： offsetとは、ある変数の平均的な値

! data(i,j) = sdata(i,j) * scale + offset

! のようにして、本当のデータを得る

makefile