ASCATのデータをFortranで読む

メインプログラム: 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