grib2のかきかた

まだGrADSですが、grib2データの書き方です。データはいつものように京大生存研からいただきます。

まず、wgrib2とg2ctlのインストールが必要です。

$ wget ftp://ftp.cpc.ncep.noaa.gov/wd51we/wgrib2/wgrib2.tgz
$ tar xzvf wgrib2
$ cd grib2/
$ make
# cp wgrib2/wgrib2 /usr/local/bin/

$ wget ftp://ftp.cpc.ncep.noaa.gov/wd51we/g2ctl/g2ctl
$ chmod a+x g2ctl
(最初の方にある$wgrib2を$wgrib2='/usr/local/bin/wgrib2';にかきかえる)
# cp g2ctl /usr/local/bin

データを取得した後、以下のシェルスクリプト（コマンドラインでもよい）

#!/bin/sh
g2ctl -0 $1 > tmp.ctl
gribmap -0 -i tmp.ctl

を実行します。$1なのでシェル実行時にGrADSで読みたいgrib2形式のファイルを引数に指定します。

あとは、Gradsでこのtmp.ctlをopenすれば描画できます。

GSMの描画

下準備として、

wget http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2014/10/20/Z__C_RJTD_20141020000000_GSM_GPV_Rgl_FD0312_grib2.bin
g2ctl -0 $1 > tmp.ctl
gribmap -0 -i tmp.ctl

としておいて、GrADSで

'reinit'
'open tmp.ctl'
'set mpdset hires'
'set lat 23 44'
'set lon 120 150'
'set gxout grfill'
'set clevs  5 10 20 40 80 160 320'
'set ccols 0 4 11  10 7 12 8 2 6'
'd apcpsfc'
'cbarn'
'set gxout contour'
'set cint 2'
'set ccolor 1'
'd prmslmsl*0.01'
'print -R'
'!convert grads.eps grads.png'

1週間予報[EPSW]対応GrADS

#!/bin/sh
day=20131125120000
hgt=500hgt
temp=850temp
slp=slp
wgrib2 'Z__C_RJTD_'$day'_EPSW_GPV_Rgl_FD00-08_grib2.bin' -match ":HGT:500" -grib $hgt.grib2
wgrib2 'Z__C_RJTD_'$day'_EPSW_GPV_Rgl_FD00-08_grib2.bin' -match ":TMP:850" -grib $temp.grib2
wgrib2 'Z__C_RJTD_'$day'_EPSW_GPV_Rgl_FD00-08_grib2.bin' -match ":PRMSL:" -grib $slp.grib2

for i in $hgt $temp $slp
do
g2ctl -0 $i.grib2 > $i.ctl
gribmap -0 -i $i.ctl
done

とりあえず500高度と850気温とslpがみれればよいか。

異常天候早期警戒情報[EPS1]

wgrib2 Z__C_RJTD_20131125120000_GSM_GPV_Rgl_FD1100_grib2.bin -match ":HGT:" -no_header -ieee hgt.grd

とすれば変数HGTの情報が単精度実数となってhgt.grdに書き出されます。これはfortranで書式なし直接形式big_endianで読めます。
fortranも含めるとこんなかんじです

echo "***WORKING ON "$1" HGT500 ***"
wgrib2 $1 -match ":HGT:" -no_header -ieee hgt500.grd > hgt500.log
echo "***WORKING ON "$2" TMP850 ***"
wgrib2 $2 -match ":TMP:" -no_header -ieee tmp850.grd > tmp850.log

program readgribgrid
implicit none
integer :: i,j,k,l,nrec
integer,parameter :: nx=144,ny=73,nz=850
real(4),dimension(nx,ny,nz) :: h,hmean

open(10,file="hgt500.grd",form="unformatted",access="direct",convert='big_endian',recl=4*nx*ny)
open(11,file="tmp2")
open(12,file="hgt500member.grd",form="unformatted",access="direct",convert="little_endian",recl=4*nx*ny)
open(13,file="hgt500member.ctl")
open(14,file="hgt500mean.grd",form="unformatted",access="direct",convert="little_endian",recl=4*nx*ny)
open(15,file="hgt500mean.ctl")
do k=1,nz
read(10,rec=k) ((h(i,j,k),i=1,nx),j=1,ny)
end do

do l=1,17
 do j=1,ny
  do i=1,nx
   do k=l,l+833,17
    hmean(i,j,l)=hmean(i,j,l)+h(i,j,k)/50.0
   end do
  end do
 end do
end do


nrec=1
do k=1,50
 do l=1,17
 write(12,rec=nrec) ((h(i,j,nrec),i=1,nx),j=1,ny)
 nrec=nrec+1
 end do
end do

do k=1,17
 write(14,rec=k) ((hmean(i,j,k),i=1,nx),j=1,ny)
 do j=1,ny
 do i=1,nx
! write(*,*) hmean(i,j,k)
 end do
 end do
end do

write(13,*) "dset ^hgt500member.grd"
write(13,*) "undef 9.999E+20"
write(13,*) "title 500hgt.1monthfrcst.ensemble"
write(13,*) "xdef 144 linear 0.0 2.5"
write(13,*) "ydef 73 linear -90.0 2.5"
write(13,*) "zdef 17 linear 1 1"
write(13,*) "tdef 50 linear 00Z01JAN0001 1DY"
write(13,*) "VARS 1"
write(13,*) "hgt 17 0 height 500hpa"
write(13,*) "ENDVARS"
write(13,*) ""
write(15,*) "dset ^hgt500mean.grd"
write(15,*) "undef 9.999E+20"
write(15,*) "title 500hgt.1monthfrcst.ensemble.mean"
write(15,*) "xdef 144 linear 0.0 2.5"
write(15,*) "ydef 73 linear -90.0 2.5"
write(15,*) "zdef 1 levels 1"
write(15,*) "tdef 17 linear 00Z01JAN0001 1DY"
write(15,*) "VARS 1"
write(15,*) "hgt 1 0 height 500hpa"
write(15,*) "ENDVARS"
write(15,*) ""
stop
end program

back