メインプログラム: read_ascat
program read_ascat
! Description:
!
! Author: satsuki
!
! Host: aofd30
! Directory: /work3/satsuki/ASCAT/AVE
!
! Revision history:
! This file is created by /usr/local/mybin/nff.sh at 17:24 on 08-03-2012.
! use
! implicit none
! write(*,'(a)')'Program read_ascat starts.'
! write(*,*)''
! 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, ofle2*500, ofle2_base*500, ofle3*500
character var*50
real,allocatable :: lon(:),lat(:)
integer(2),allocatable :: eastward_wind(:,:,:,:), northward_wind(:,:,:,:)
! real,allocatable :: sdata(:,:,:)
integer,allocatable :: dimids(:),nshape(:)
character*70 err_message
character*100 dimname !dimnameは各次元の名前。
real,parameter::RMISS=-999.9 !欠損値
real,allocatable :: u(:,:),v(:,:)
write(*,'(a)')'Program read_ncep starts.'
! write(*,*)''
im=1440 !192
jm=641 ! !94
nm=1 !366
allocate(lon(im),lat(jm),eastward_wind(im,jm,1,1)&
& ,northward_wind(im,jm,1,1))
infle="../2007063000_2007070100_daily-ifremer-L3-MWF-GLO-20110529093951-01.0.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='eastward_wind'
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_u)
print *,'scale_u= ',scale_u
print *,'sta= ',sta
print *
! オフセットの 値を得る(*2)
sta = nf_get_att_real(ncid,varid,'add_offset',offset_u)
print *,'offset_u=', offset_u
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, eastward_wind)
print *,'stat= ',stat
var='northward_wind'
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_v)
print *,'scale= ',scale_v
print *,'sta= ',sta
print *
! オフセットの 値を得る(*2)
sta = nf_get_att_real(ncid,varid,'add_offset',offset_v)
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
deallocate(dimids,nshape)
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, northward_wind)
print *,'stat= ',stat
var='longitude'
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='latitude'
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_ascat.asc"
ofle2="output/test_ascat.bin"
ofle3="output/test_ascat.ctl"
! print '(A,A)','Output: ',ofle(1:lnblnk(ofle))
! print *
allocate(u(im,jm),v(im,jm))
! open(20,file=ofle)
! Note: time index (n) is fixed and set to 1 here.
n=1 !366
do j=1,jm
do i=1,im
u(i,j)=RMISS
v(i,j)=RMISS
if(eastward_wind(i,j,1,1) > -32760 )then
u(i,j)=eastward_wind( i,j,1,1 )*scale_u + offset_u
endif
if(northward_wind(i,j,1,1) > -32760 )then
v(i,j)=northward_wind(i,j,1,1 )*scale_v + offset_v
endif
! write(20,'(2f12.6,2f12.4,1x,2i5)') &
! lon(i),lat(j),u(i,j),v(i,j),i,j
enddo !i
enddo !j
close(20)
print '(A,A)','Output: ',ofle2(1:lnblnk(ofle2))
print *
open(20,file=ofle2,form="unformatted")
write(20)((u(i,j),i=1,im),j=1,jm)
write(20)((v(i,j),i=1,im),j=1,jm)
close(20)
print '(A,A)','Output: ',ofle3(1:lnblnk(ofle3))
print *
open(20,file=ofle3)
call basename(ofle2, ofle2_base, idxbi)
write(20,'(A,A)')'DSET ',ofle2_base(1:lnblnk(ofle2_base))
write(20,'(A)')'OPTIONS sequential'
write(20,'(A)')'UNDEF -999.9'
write(20,'(A,i5,A,500f12.3)')'XDEF ',im,' LINEAR ' ,-179.875, 0.25
write(20,'(A,i5,A,500f12.6)')'YDEF ',jm,' LINEAR ', -79.875, 0.25
write(20,'(A,i5,A,500f10.3)')'ZDEF ',1,' LEVELS ',10
write(20,'(A,i7,A,A)')'TDEF ',1,' LINEAR 00z30Jun2007 ','1dy'
write(20,'(A)')'VARS 2'
write(20,'(A,i5,A)')'u ',1,' 99 eastward_wind [m/s]'
write(20,'(A,i5,A)')'v ',1,' 99 northward_wind [m/s]'
write(20,'(A)')'ENDVARS'
write(*,'(a)')'Done program read_ascat.'
write(*,*)
end program read_ascat
サブルーチン: basename.f
!***********************************************************************
! .f
! NOTE
! Reference
!***********************************************************************
subroutine basename(infle, bifle, idxbi)
!------------------------- DECLARATION ---------------------------------
! implicit none
! implicit double precision(a-h,o-z)
!------------------------- INCLUDE FILE --------------------------------
! include 'comblk.h'
!-------------------------- PARAMETERS ---------------------------------
! parameter()
!------------------------- COMMON VARIABLES ----------------------------
! common //
!----------------------------ARGUMENTS ---------------------------------
! double precision
! real
! integer
! dimension
! Input
character infle*(*)
! Output
character bifle*(*)
integer idxbi
! Input/output
!------------------------- LOCAL VARIABLES -----------------------------
character testchar*1
!-----------------------------------------------------------------------
!---+$|--1----+----2----+----3----+----4----+----5----+----6----+----7--
last=index(infle,' ')-1
do i=last, 1, -1
testchar=infle(i:i)
id=index(testchar,'/')
if(id.ne.0)then
idx=i+1
goto 1000
end if
end do !i
! If there is no slash in infle
bifle=infle(1:last)
idxbi=last
return
! If there is a slash in infle
1000 bifle=infle(idx:last)
idxbi=last-idx+1
return
end
コンパイル
注意
この例でにおけるnetCDFライブラリのインストール先: /usr/local/lib
コンパイルオプションの説明
-CB -traceback -fpe0: ifortのデバッグ用
-module ../obj: ifortでmoduleファイルの置き場所を指定している
[satsuki@aofd30 AVE]$ make
if [ ! -d ../obj ]; then \
mkdir -p ../obj ; \
fi
ifort -c -CB -traceback -fpe0 -module ../obj -c -o ../obj/read_ascat.o read_ascat.f90
ifort -o read_ascat ../obj/read_ascat.o -module ../obj -L/usr/local/lib -lnetcdf
[satsuki@aofd30 AVE]$ ./read_ascat
Program read_ncep starts.
Input: ../2007063000_2007070100_daily-ifremer-L3-MWF-GLO-20110529093951-01.0.nc
ncid= 3
status= 0.0000000E+00
eastward_wind
varid= 6
sta= 0.0000000E+00
scale= 9.9999998E-03
sta= 0.0000000E+00
offset= 0.0000000E+00
sta= 0.0000000E+00
ndims= 4
dimids( 1)= 4
dimids( 2)= 3
dimids( 3)= 2
dimids( 4)= 1
nshape( 1)= 1440
nshape( 2)= 641
nshape( 3)= 1
nshape( 4)= 1
dimname( 1)= longitude
dimname( 2)= latitude
dimname( 3)= depth
dimname( 4)= time
Read eastward_wind
stat= 0.0000000E+00
longitude
varid= 4
sta= 0.0000000E+00
Read longitude
stat= 0.0000000E+00
latitude
varid= 3
sta= 0.0000000E+00
Read latitude
stat= 0.0000000E+00
Output: output/test.asc
[satsuki@aofd30 ASCAT]$ ncdump -h 2007063000_2007070100_daily-ifremer-L3-MWF-GLO-20110529093951-01.0.nc
-h: ヘッダのみ出力する
netcdf 2007063000_2007070100_daily-ifremer-L3-MWF-GLO-20110529093951-01.0 {
dimensions:
time = 1 ;
depth = 1 ;
latitude = 641 ;
longitude = 1440 ;
variables:
int time(time) ;
time:long_name = "time" ;
time:units = "hours since 1900-01-01 00:00:00.0 00:00" ;
time:valid_min = 942276. ;
time:valid_max = 942276. ;
time:axis = "T" ;
time:standard_name = "time" ;
float depth(depth) ;
depth:long_name = "depth" ;
depth:units = "m" ;
depth:valid_min = 10. ;
depth:valid_max = 10. ;
depth:axis = "Z" ;
depth:positive = "up" ;
depth:standard_name = "depth" ;
float latitude(latitude) ;
latitude:long_name = "latitude" ;
latitude:units = "degrees_north" ;
latitude:valid_min = -79.875 ;
latitude:valid_max = 80.125 ;
latitude:axis = "Y" ;
latitude:standard_name = "latitude" ;
float longitude(longitude) ;
longitude:long_name = "longitude" ;
longitude:units = "degrees_east" ;
longitude:valid_min = -179.875 ;
longitude:valid_max = 179.875 ;
longitude:axis = "X" ;
longitude:standard_name = "longitude" ;
short wind_speed(time, depth, latitude, longitude) ;
wind_speed:long_name = "wind speed" ;
wind_speed:units = "m/s" ;
wind_speed:scale_factor = 0.01 ;
wind_speed:add_offset = 0. ;
wind_speed:valid_min = 0. ;
wind_speed:valid_max = 6000. ;
wind_speed:_FillValue = -32768s ;
wind_speed:standard_name = "wind_speed" ;
short eastward_wind(time, depth, latitude, longitude) ;
eastward_wind:long_name = "eastward wind speed" ;
eastward_wind:units = "m/s" ;
eastward_wind:scale_factor = 0.01 ;
eastward_wind:add_offset = 0. ;
eastward_wind:valid_min = -6000. ;
eastward_wind:valid_max = 6000. ;
eastward_wind:_FillValue = -32768s ;Fillvalueとは欠損値のこと
eastward_wind:standard_name = "eastward_wind" ;
eastward_wind:positive = "meteorological convention" ;
short northward_wind(time, depth, latitude, longitude) ;
northward_wind:long_name = "northward wind speed" ;
northward_wind:units = "m/s" ;
northward_wind:scale_factor = 0.01 ;
northward_wind:add_offset = 0. ;
northward_wind:valid_min = -6000. ;
northward_wind:valid_max = 6000. ;
northward_wind:_FillValue = -32768s ;
northward_wind:standard_name = "northward_wind" ;
northward_wind:positive = "meteorological convention" ;
ncdumpの下記の出力結果と
short eastward_wind(time, depth, latitude, longitude) ;
下記のFortranプログラムの配列eastward_windの添え字の順番が異なるのに注意すること。
Fortranプログラムの配列eastward_windでは(経度, 緯度, 深さ, 時間)の順番になっている。
おそらくdimidsの値が重要だと思われる
dimids( 1)= 4
dimids( 2)= 3
dimids( 3)= 2
dimids( 4)= 1