GrADS station data
〜Fortran90での作成から作図まで〜

　　GrADSでは，格子点（グリッド）データだけではなく，AMeDASのような不規則な配置であっても
　Station Dataファイルを作成することにより使用可能となります．Station Dataの作成方法は，
　公式サイト，また日本語でも幾つかインターネット上に掲載されています．
　
　　しかし，Fortranでの作成の場合，出力ファイルをopenする際の，form，recordtype，access等が
　各解説サイトで異なることが多く，また，作成（Fortran90）から作図（ctlファイル）の手順・完成図の
　流れで解説されているサイトは少ないです．公式サイトのサンプルデータ，Fortranプログラムで
　作成すると，tdef（タイムステップ数）が2のはずなのに，tdef 1でしかctlファイルを開くことが
　できません．
　
　　自分は普段あまりGrADSを解析に使用しないので，Station Dataの作成・描図の際にエラーが
　数多く発生したのですが，作成方法に問題があるのか，ctlファイルに問題があるのか，原因が特定
　できずかなり苦戦しました．

　　そこで，複数要素のGrADSのStation DataファイルをFortran90で作成し，このファイルを開くための
　ctlファイルの作成，完成図・・・の流れで下に解説します．

　※使用環境：Intel(R) Fortran Compiler Version 8.0 (Windows版)，GrADS Version 1.8 (Windows版)
＜解説＞
1. Fortran90でのStation Dataの作成
1-1. データの準備
　　まずは，データの準備です．ここでは，記録的な猛暑だった2008年8月の平均気温，平均
　最高/最低気温を使用します．データの並びは以下の通りです．左から順に地点番号，経度，
　緯度，高度，年，月，平均気温（以下Tmean），平均最高気温（以下Tmax），平均最低気温
　（以下Tmin）です．

　40041 140.6367 36.8683 370 2008 8 999.00 999.00 999.00
　40046 140.7700 36.8417 45 2008 8 23.48 26.38 21.09
　40061 140.3450 36.7783 120 2008 8 24.17 29.18 20.93

　データのダウンロード→amedas_kanto_2008_aug_dat


1-2　Station Dataの構造
　　Station Dataはバイナリーで書き込みます．
　ファイルの構成は原則以下の通りです．
　地点番号のようなID（character 8byte），緯度（real 4byte），経度（real 4byte），時（real 4byte），
　鉛直層（zdef）とデータの終了フラグ（integer 4byte），時間層（tdef）フラグ（integer 4byte），
　データ1（型とbyteは任意），データ2（型とbyteは任意），データ3（型とbyteは任意）
　下図のようなイメージ．




　id: 地点番号や地点名をcharacter 8byteで．地点番号のような整数でもcharacterで入力．
　lat: 緯度をreal 4byteで入力．
　lon: 経度をreal 4byteで入力．
　t: 時間．普通の時間とは異なるので注意！データが毎正時のものであれば0.0，毎正15分の
　　ものであれば4分の1時間なので0.25・・・のように入力．したがって，tは-0.5〜0.5の範囲に．
　nlev: データ入力時は1を入力．ある時間のデータが終了，もしくはデータそのものの終了
　　（＝ファイルの終わり）のとき，次の基本情報部id，lat，lonは適当で良いのでに0を
　　入力する（下図のイメージ参照）．
　flag:: データ入力時は1を入力．データそのものの終了（＝ファイルの終わり）のとき，0を入力する．
　
　以上をまとめると今回のサンプルで作成するファイルは下図のようなイメージとなる．
　ファイルの最後にnlev=0，flag=0の基本情報部が必要…ということが重要です．



1-3　Station Dataの作成
　　1-2で解説した構造を，Fortran90で作成します．サンプルプログラムは以下のようになります．
　サンプルプログラムのダウンロード→make_grads_station_data_sample.f90

! programここから---------------------------------------------------------------------
!!
!! 2008年8月の平均気温，平均最高/最低気温データをGrADS Station Dataに変換する
!! サンプルプログラム．
!!
!!-----------------------------------------------------------------------------

implicit none

　integer,parameter :: iwcmax=2000
　integer :: i,k,io,iwc

! file(10)
　character*8 :: stnum(iwcmax)
　integer(4) :: alt(iwcmax),year(iwcmax),month(iwcmax)
　real(4) :: lon(iwcmax),lat(iwcmax)
　real(4) :: tmean(iwcmax),tmax(iwcmax),tmin(iwcmax)
　integer(4) :: nlev,nflag
　real(4) :: t

! 入力データのopen & read
　open(10,file='./amedas_kanto_2008_aug.dat', status='old')
　read(10,*)

　iwc=0

　do i=1,iwcmax

　　read(10,'(a5,f9.4,f8.4,2i5,i3,3f7.2)',iostat=io) &
　　stnum(i),lat(i),lon(i),alt(i),year(i),month(i),&
　　tmean(i),tmax(i),tmin(i)

　　if(io < 0) exit ! end
　　iwc=iwc+1

　end do ! i


! 出力ファイルのopen
　open(55,file='./amedas_kanto_2008_aug_grads_station.dat', form='unformatted',recordtype='stream')

　do i=1,iwc

　　nlev = 1
　　nflag = 1
　　t =0.0

　　write(55) stnum(i),lon(i),lat(i),t,nlev,nflag
　　write(55) tmean(i),tmax(i),tmin(i)

　end do


! 最後にnlev=0，flag=0の基本情報部をwrite
　nlev = 0
　nflag = 0

　write(55) stnum(i),lon(i),lat(i),t,nlev,nflag

　close(55)


stop
end program
 ! programここまで--------------------------------------------------------------------


　これをコンパイルして，実行．

　プログラムのポイント（上のプログラム中の赤文字部分）としては，
　●変数定義の段階で，基本情報部の型（character, integer, real）とバイト数を指定すること．
　●formはunformatted，recordtypeはstream．
　●データの最後にnlev=0，flag=0の基本情報部のみ（データ部は不要）をwriteすること．
　　その際，stnum, lon, latはなんでも良い．特に記述しなくてもバイト数を定義してあれば
　　適当に埋め込まれるはず．



2. ctl（コントロール）ファイルの作成
　　1の手順により，station dataファイルが完成（のはず・・・）．次はGrADSで描図できるように
　ctl（コントロールファイル）を作成する．サンプルctlファイルは以下のようになります．
　サンプルctlファイルのダウンロード→amedas_kanto_2008_aug_sample.ctl

! ctlファイルここから--------------------------------------------------------------------
dset ^amedas_kanto_2008_aug_grads_station.dat
dtype station
stnmap ^amedas_kanto_2008_aug_grads_station.map
undef 999.00
title AMeDAS 2008 Aug
tdef 1 linear 00Z01Aug2008 01yr
vars 3
tmean 0 99 Tmean
tmax 0 99 Tmax
tmin 0 99 Tmin
endvars
 ! ctlファイルここまで--------------------------------------------------------------------

　ctlファイルのポイント（上のプログラム中の赤文字部分）としては，
　●dtypeをstationに設定．
　●station mapを指定（3-1で作成します）．
　●変数はwriteした順（今回のサンプルではTmean→Tmax→Tmin）に記述すること．



3. GrADSで作図
3-1. station mapの作成
　station dataを使用する場合，ctlファイルを開く前にstation mapを作成する必要があります．

　作成方法は
　1. GrADSを起動．
　2. 画面に先程作成したctlファイルをstnmapコマンドで実行．
　　※windowsの場合はUNIXコマンド実行時はコマンドの前に!が必要．cygwin上やlinuxでは!は不要．
　　!stnmap -i （ctlファイルのあるディレクトリ）/amedas_kanto_2008_aug_sample.ctl



　station mapの作成が成功すると下のような画面が表示される．作業フォルダにstation mapが作成
　されているかを確認する．




3-2. 作図（ctlファイルのopen）→お好みの図に
　station mapが作成できれば，あとは通常と同様にctlファイルをopenして，図を好きなように調整する．
　open （ctlファイルのあるディレクトリ）/amedas_kanto_2008_aug_sample.ctl



　鉛直方向1層でstation dataを作成すると，ctlファイルをopenしたときにlevel set to 500 500に
　なるが，作図は問題なくできるのでここではスルー．

　関東領域の緯度，経度，表示する数値の大きさ・小数点の桁数,，地図の解像度，背景色を白に
　指定するなどお好みの設定で作図．

set lat 34.5 37.2
set lon 138 141
set digsiz 0.15
set dignum 1
set mpdset hires
set display color white
d tmean



　以上で完成．あとは内挿したり，色付けしたり作図者の好みで．

　「無事使えた！」「ここ違うのでは？」等，BBSにお書き込みいただけると，今後のプログラム公開の
モチベーションアップになりますので，是非書き込んでいただけると有難いです．
＜参考文献＞
　1. http://www.iges.org/grads/gadoc/usingstationdata.html
　2. http://www.mets.dyndns.org/misc/gadoc/aboutstationdata.html
　3. http://www28.atwiki.jp/fmemo/pages/20.html
　4. http://nattoumaruo.com/blog/archives/32
　5. http://hydro.iis.u-tokyo.ac.jp/~kei/?IT%20memo%2FGrADS%20memo%2FStation%20data
　6. http://mc39.genin.jp/GrADS_Station_Data.html
＜検索用タグ＞
　GrADS, station data, ステーションデータ, Fortran90, ctl, コントロールファイル