風速の平面図

出来上がった図の例 （500hPa面の風向風速、ジオポテンシャル高度）

データフォーマットの変換

grib2形式のGSMデータをGrADS形式に変換する

プログラムは、防災科研の清水さん作成のread_gsmを使用する。詳細は以下に記載。もしくは、防災科研の清水さんのwikiサイトにある、モデル関連メモを参照のこと。

プログラムの入手

防災科研の清水さんのwikiサイトにある、モデル関連メモというページより、read_gsm.tar.gzをダウンロードする。

プログラムの書庫ファイルの解凍

$ tar xzvf read_gsm.tar.gz

で解凍して、できたファイルをread_gsmというディレクトリに入れる。

プログラムのコンパイル

read_gsmディレクトリに移り、プログラムをコンパイルする。

$ cd read_gsm

$ make

(注：付属のMakefileではコンパイラにgccを指定している→大抵のOSで動きそう)

GSMデータの入手

inputというディレクトリを作成する。

$ mkdir input

inputディレクトリ内に、 下記のget.gsm.global.shという名称のファイルを作成する。

cat > get.gsm.global.sh

#!/bin/bash

# Z__C_RJTD_yyyymmddhh0000_GSM_GPV_Rgl_FDtttt_grib2.bin

# yyyymmddhh ... 初期時刻(UTC), hh=00,06,12,18

# tttt ... 0000: 予報時刻 000h (=初期値)

yyyy=2012

mm=06

dd=15

tttt=0000

url=http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/$yyyy/$mm/$dd

hhlist="00 06 12 18"

for hh in $hhlist; do

file=${url}/Z__C_RJTD_${yyyy}${mm}${dd}${hh}0000_GSM_GPV_Rgl_FD${tttt}_grib2.bin

wget $file

done

Ctl+D

get.gsm.global.shを実行

このサンプルでは、2012年6月15日午前00時（UTC)のデータがダウンロードされる。

プログラム実行用シェルスクリプトの作成

inputディレクトリから親ディレクトリに移動する

$ cd ..

親ディレクトリ内に、 下記のgsm2bin.run.shという名称のファイルを作成する。

cat > gsm2bin.run.sh

#!/bin/bash

exe=read_gsm/gsm2bin

if [ ! -f $exe ]; then

echo Error in $0: No such file, $exe

exit 1

fi

indir=input

outdir=output

mkdir -vp $outdir

yyyy=2012

mm=06

dd=15

hh=00

tttt=0000

input=Z__C_RJTD_${yyyy}${mm}${dd}${hh}0000_GSM_GPV_Rgl_FD${tttt}_grib2.bin

in=${indir}/${input}

if [ ! -f $in ]; then

echo Error in $0: No such file, $in

exit 1

fi

yyyymmdd=${yyyy}${mm}${dd}

output=GSM${yyyymmdd}${hh}Z.bin

out=${outdir}/${output}

$exe $in $out

ctl=${outdir}/$(basename $output .bin).ctl

imm=$mm

if [ ${imm} -eq 1 ]; then

mmm=Jan

elif [ ${imm} -eq 2 ]; then

mmm=Feb

elif [ ${imm} -eq 3 ]; then

mmm=Mar

elif [ ${imm} -eq 4 ]; then

mmm=Apr

elif [ ${imm} -eq 5 ]; then

mmm=May

elif [ ${imm} -eq 6 ]; then

mmm=Jun

elif [ ${imm} -eq 7 ]; then

mmm=Jul

elif [ ${imm} -eq 8 ]; then

mmm=Aug

elif [ ${imm} -eq 9 ]; then

mmm=Sep

elif [ ${imm} -eq 10 ]; then

mmm=Nov

elif [ ${imm} -eq 11 ]; then

mmm=Dec

fi

cat <<EOF > $ctl

dset ^GSM%y4%m2%d2%h2Z.bin

options template

undef -999

xdef 720 LINEAR 0.0 0.5

ydef 361 LINEAR -90.0 0.5

zdef 17 LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 70 50

30 20 10

tdef 2 LINEAR ${hh}:00Z${dd}Jun${yyyy} 6hr

vars 17

slp 0 0 sea level pressure

sp 0 0 surface pressure

su 0 0 surface westerly wind

sv 0 0 surface southerly wind

st 0 0 surface temperature

srh 0 0 surface relative humidity

lca 0 0 LowCloudAmount

mca 0 0 MidCloudAmount

hca 0 0 HighCloudAmount

tca 0 0 TotalCloudAmount

z 17 0 height

u 17 0 westerly wind

v 17 0 southerly wind

t 17 0 temperature

rh 17 0 relative humidity

w 17 0 p-velocity

pre 0 0 precipitation

endvars

EOF

echo

echo $in

echo $ctl

echo $out

echo

CTL+D

プログラムの実行

gsm2bin.run.shを実行

$ bash gsm2bin.run.sh

結果の作図

この例では、作図にはGrADSが必要

親ディレクトリ内に、次の2つのファイル

GSM.snapshot.run.sh

GSM.snapshot.gs

を作成する。

$ cat > GSM.snapshot.run.sh

#!/bin/sh

exe=GSM.snapshot.gs

if [ ! -f $exe ]; then

echo Error in $0 : No such file, $exe

exit 1

fi

yyyy=2012

mm=06

dd=15

hh=00

# Parameters

indir=output

input=GSM${yyyy}${mm}${dd}${hh}Z.ctl

prefix= #MandA2012Leg3

lev=700 # Vertical Level

opt="-indir $indir -input $input"

var="umag" #theta

ci=3 #Contour Interval

it=1

nt=1

opt_grads=bcp

while [ $it -le $nt ]; do

out=$(basename $input .ctl).${var}.${lev}.$(printf %02d $it).eps

grads -${opt_grads} "$exe $opt -t1 $it -lev $lev -var ${var} -ci $ci \

-out $out -q"

it=$(expr $it + 1 )

done

exit 0

CTL+D

$ cat > GSM.snapshot.gs

# undef -999

# xdef 720 linear 0 0.5

# ydef 361 linear -90 0.5

# zdef 17 levels 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10

# tdef 2 linear 00Z15JUN2012 360mn

# vars 17

# slp 0 0 sea level pressure

# sp 0 0 surface pressure

# su 0 0 surface westerly wind

# sv 0 0 surface southerly wind

# st 0 0 surface temperature

# srh 0 0 surface relative humidity

# lca 0 0 LowCloudAmount

# mca 0 0 MidCloudAmount

# hca 0 0 HighCloudAmount

# tca 0 0 TotalCloudAmount

# z 17 0 height

# u 17 0 westerly wind

# v 17 0 southerly wind

# t 17 0 temperature

# rh 17 0 relative humidity

# w 17 0 p-velocity

#

function gsm( args )

*

* Default values

*

indir='output/'

input='GSM2012061500Z.ctl'

lonw=110; lone=160

lats=15 ; latn=55

out='GSM2012061500Z.eps'

quitopt='no'

lev=1000

ci=5

var='umag'

*

* Decode options

*

i = 1

while( 1 )

arg = subwrd( args, i )

i = i + 1;

if( arg = '' ); break; endif

while( 1 )

if( arg = '-indir');indir = subwrd(args,i);i=i+1; break;endif

if( arg = '-input');input = subwrd(args,i);i=i+1; break;endif

if( arg = '-var' ) ; var = subwrd(args,i);i=i+1; break;endif

if( arg = '-lonw' ); lonw = subwrd(args,i);i=i+1; break;endif

if( arg = '-lone' ); lone = subwrd(args,i);i=i+1; break;endif

if( arg = '-lats' ); lats = subwrd(args,i);i=i+1; break;endif

if( arg = '-latn' ); latn = subwrd(args,i);i=i+1; break;endif

if( arg = '-lev' ); lev = subwrd(args,i) ;i=i+1; break;endif

if( arg = '-t1' ); t1 = subwrd(args,i) ;i=i+1; break;endif

if( arg = '-ci' ); ci = subwrd(args,i) ;i=i+1; break;endif

if( arg = '-out' ); out = subwrd(args,i) ;i=i+1;break;endif

if( arg = '-q' ); quitopt=yes ;i=i+1; break; endif

say 'Syntax error : arg= 'arg

return

endwhile

endwhile

'open 'indir'/'input

#'q ctlinfo'

#say result

'set lev 'lev

'set t 't1

#'q file'

#say result

'q dims'

say result

line=sublin(result,5)

datetime=subwrd(line,6)

line=sublin(result,4)

pdb=subwrd(line,6)

pha=pdb

say pha

'set lon 'lonw' 'lone

'set lat 'lats' 'latn

# 'set mproj lambert'

if (var = "ept" )

'es=6.1173*exp(((2.501*pow(10,6))/461.50)*(1/273.16 - 1/t))'

'ws=621.97*(es/(' pha '-es))'

'w=(rh*ws)/(100*1000)'

'ept=(t+((2.501*pow(10,6))/1004)*w)*pow((1000/' pha '),287/1004)'

'set clevs 315 318 321 324 327 330 333 336 339 342 345'

endif

if (var = "theta" )

'theta=t*pow((1000/' pha '),287/1004)'

endif

'set parea 1 6 2 9'

'set xlopts 1 5 0.15'

'set ylopts 1 5 0.15'

'set xlint 10'

'set ylint 10'

'set gxout shade2'

### Vmag

'color 15 35 5 -kind white->navajowhite->red'

'd mag(u,v)'

'cbarn 1 1 6.5 5'

'set strsiz 0.12 0.14'

'set string 1 c 5 0'

'draw string 6.7 6.7 m/s'

### V

'set gxout vector'

'set ccolor 0'

'set cthick 20'

'd skip(u,4,4);skip(v,4,4)'

'set ccolor 1'

'set cthick 5'

'd skip(u,4,4);skip(v,4,4)'

'set strsiz 0.10 0.14'

'set string 1 c 5 0'

'draw string 4.65 2.5 m/s'

### Z

'set gxout contour'

'set cint 50'

'set ccolor 0'

'set cthick 20'

'set clab off'

'd z'

'set ccolor 1'

'set cint 50'

'set clab off'

'set cthick 5'

'd z'

'set ccolor 1'

'set clopts 1 5 0.15'

'set clab on'

'set cint 100'

'set cthick 0'

'd z'

'draw title 'datetime' 'pha' hPa'

'print 'out

if ( quitopt = "yes" ); quit; endif

return

CTL+D