MANAL を読んでテキスト出力する tcsh スクリプトのサンプル

研究室への説明文 (2009.1)

研究室の皆さん、特に山田君

玉川です。気象庁客観解析データの MANAL を
読み出して、経度、緯度、データの形のテキストファイルをつくる
スクリプトのサンプルを送ります。

このまま、 192.168.1.12 の PC 上へ置いて
chmod +x Textdump.tcsh
./Textdump.tcsh
とすると、最初だけ wgrib.c を取りに行きますが、
2008 8 27 06Z の 1000mb の気温データを読み出して
経度 緯度 気温 のデータを並べたデータセットを作成します。
このデータの中から、例えば欲しい場所の周囲のデータを
下の行のように、経度、緯度を使って (上限 下限をそれぞれ x1, x2, y1, y2 と書いていますが
実際には、数値をいれてください )

awk '{if( x1 <= $1 && $1 <= x2 && y1 <= $2 && $2 <= y1 ) print $1, $2, $3}' < 2005082706_TMP_1000mb.txt

などとして、とりだして使ってください。

最初の方の読み出すデータ名や、WGRIB の行の
TMP 1000 mb の付近を適当に変更すると、他のデータ、他の項目も
読み出せます。

他のPCで使う時も、最初の方の FC や GRBFILE を適切に設定しなおせば
OKです。

--------------------------------------------------------------------------------
〒 501-1193 岐阜市柳戸1-1
岐阜大学 流域圏科学研究センター
玉川一郎 (e-mail: tama@green.gifu-u.ac.jp)
phone/fax 兼用: 058-293-2430
WWW: 	http://tama.green.gifu-u.ac.jp/~tama
--------------------------------------------------------------------------------

#!/bin/tcsh
set WGRIB = ./wgrib
set GRBFILE = /work/JMA/manal/2005/08/MANAL_2005082706.grb
set OUTFILE = 2005082706_TMP_1000mb
set FC     = gfortran


if( -e wgrib.c)  then
echo ok
else
 echo no
   wget ftp://ftp.cpc.ncep.noaa.gov/wd51we/wgrib/wgrib.c
endif


cc $WGRIB.c -o $WGRIB

$WGRIB $GRBFILE |  grep TMP | grep '1000 mb' | $WGRIB $GRBFILE -i -bin -nh -o ${OUTFILE}.bin

# ./wgrib ~/JMA_DATA/MANAL_200508/MANAL_2005082706.grb \
#   | grep TMP | grep '1000 mb' |\
#   ./wgrib ~/JMA_DATA/MANAL_200508/MANAL_2005082706.grb -i -bin -nh -o 2005082706_TMP_1000mb.bin2


cat << ENDPROGRAM > bin2txt_manal.f90

program bin2txt_manal
implicit none
  integer, parameter :: m=361, n=289
  real :: x(m), lat, lon
  character(len=1024) :: infile, outfile
  integer :: argc, i, j

  argc = iargc()
  if( argc < 2 ) then
   write(0,*) 'Please give me filename to read and write'
   write(0,*) './bin2txt_manal  infile.bin outfile.txt'
   stop
  end if
  call getarg(1, infile)
  call getarg(2, outfile)

  open(unit=10, file=infile, err=900, &
    access='direct', form='unformatted', recl=4*m )
  open(unit=11, file=outfile, err=910)

  do j=1, n
    read(unit=10, rec=j) (x(i), i=1,m)
    do i=1, m
     call ijtolatlon( i, j, lon, lat, m, n )
     write(11,*) lon, lat, x(i)
    end do
  end do

  close(10); close(11)


 stop

900  write(0,*) 'Open error ' ,infile
 stop

910  write(0,*) 'Open error ' ,outfile
 stop
 end

subroutine ijtolatlon( i, j, lon, lat, m, n )
!   for MANAL from JMA (10km grid)

  integer :: i, j
  real    :: lat, lon

  integer :: m, n
  integer :: j1
  integer :: xs=245, ys=205
  real*8 :: xp, yp
  real*8, parameter    :: fact=180.0/3.14159265358979323846
  real, parameter    :: dx=1.0e4, dy=1.0e4

 j1 = n+1 -j

   xp = (i-xs)*dx
   yp = 7.71061e6 + (j1-ys)*dy

!   write(*,*) 'xp, yp:', xp, yp

   lat = 90.0 - 2.0*fact*atan((1.37003e-10)*((xp**2 + yp**2 )**0.69875))
   lon = 1.39750*fact*atan2(xp, yp)+140.0
end subroutine



ENDPROGRAM

$FC bin2txt_manal.f90 -o bin2txt_manal

./bin2txt_manal ${OUTFILE}.bin ${OUTFILE}.txt

echo ${OUTFILE}.txt is produced
echo each line show 'longitude lattude value'