付録
1.1 標高データ作成プログラム(1600×1600)
1.2 標高データ作成プログラム(4800×4800)
2 WGRIBを用いた外部気象データの作成
1.1 標高データ作成プログラム(1600×1600)
使用方法(Windows,UNIX)
1.プログラム文中のファイルパスを使用に応じてかえる。
2.CD『数値地図50mメッシュ(標高)』を入れる。
3.標高地域メッシュコード一覧を参考に地域番号を入力する(4桁)。
4.確認データを作るか作らないかを選択する。
5.直接使うならファイル名を変更、4800×4800に使うならファイル名はそのままで。
***** 標高データをCReSSの使用形式に変換するプログラム(1600×1600) ******
character mn1*4, mn2*2
integer i,j,m,n
real a(1600,1600), e(200,200),s
write(*,*)'地域番号(4桁入れてください。)'
read(*,*)mn1
************************************************************************
***** 番号指定 *********************************************************
do i=0,77
j=i
if(j.eq. 8.or.j.eq. 9)go to 2
if(j.eq.18.or.j.eq.19)go to 2
if(j.eq.28.or.j.eq.29)go to 2
if(j.eq.38.or.j.eq.39)go to 2
if(j.eq.48.or.j.eq.49)go to 2
if(j.eq.58.or.j.eq.59)go to 2
if(j.eq.68.or.j.eq.69)go to 2
write(mn2,'(I2.2)') j
call input(mn1,mn2,e)
call import(a,e,i)
2 end do
do n=1,1600
do m=1,1600
if(a(n,m).lt.0)go to 3
go to 4
3 a(n,m)=-1
4 end do
end do
open(unit=7,file='c:\Fortran\'//mn1//'.bin',access='direct'
+ ,form='unformatted',recl=6400,err=952)
do n=1,1600
write(7,rec=n,err=951)(a(n,m)/10,m=1,1600)
end do
close(unit=7)
go to 600
************************************************************************
951 write(*,*)'出力エラー'
go to 999
952 write(*,*)'最後のエラー'
go to 999
************************************************************************
600 write(*,*)'どちらかを番号により選択してください。'
write(*,*)'1:確認ファイルを作成してから終了'
write(*,*)'2:すぐに終了'
read(5,*)s
if(s.eq.1)call kakunin(mn1)
if(s.ne.1)go to 999
go to 999
999 end
***** 標高データを読み込むプログラム(input) ****************************
subroutine input(mn1,mn2,e)
integer c(200),d(200)
integer xdeg, xmin, ydeg, ymin
character line*1011,fname*1024
character mn1*4 , mn2*2
real e(200,200)
real xsec, ysec
fname = 'd:\data\'//mn1//'\'//mn1//mn2//'.mem'
open(1, file=fname, access='direct',err=12,
+ form='unformatted', recl=1011)
read(1,rec=1) line
read(line(766:773), '(I3,I2,F3.1)') ydeg, ymin, ysec
read(line(774:781), '(I3,I2,F3.1)') xdeg, xmin, xsec
do 11 m=1,200
read(1,rec=m+1) line(1:1009)
read(line, '(i6,i3,200f5.0)')c(m),d(m),(e(m,n),n=1,200)
11 continue
go to 13
12 do m=1,200
do n=1,200
e(m,n)=-1
end do
end do
13 close(1)
return
end
************************************************************************
***** 標高データを配列に代入するプログラム(import) *********************
subroutine import(a,e,n)
real a(1600,1600), e(200,200)
integer i,j,n,nn,x,y
if(n.ge. 0.and.n.le. 7)nn=0
if(n.ge.10.and.n.le.17)nn=1
if(n.ge.20.and.n.le.27)nn=2
if(n.ge.30.and.n.le.37)nn=3
if(n.ge.40.and.n.le.47)nn=4
if(n.ge.50.and.n.le.57)nn=5
if(n.ge.60.and.n.le.67)nn=6
if(n.ge.70.and.n.le.77)nn=7
x=n-10*nn
y=nn
do i=1,200
do j=1,200
a(j+200*y,i+200*x)=e(201-j,i)
end do
end do
return
end
************************************************************************
***** 確認プログラム ***************************************************
subroutine kakunin(mn)
integer i, j
real k(1600,1600)
character mn*4
open(unit=77,file='c:\Fortran\'//mn//'.bin',access='direct',
+ form='unformatted', recl=6400,err=996)
do 53 i=1,1600
read(77,rec=i)(k(i,j),j=1,1600)
53 continue
close(77)
open(unit=78,file='c:\gmt-ground\check\check.dat')
do 51 i=1,1600
do 52 j=1,1600
write(78,100)j, i, k(i,j)
100 format(' ',i4,' ',i4,' ',f10.3)
52 continue
51 continue
close(78)
go to 999
***** プログラム作成エラー *********************************************
* *
996 write(*,*) 'You are bad!'
go to 999
* *
************************************************************************
999 write(*,*)'確認終了'
return
end
************************************************************************
1.2 標高データ作成プログラム(4800×4800)
使用方法(UNIX)
作成したい領域に入る1600×1600のデータを9個作る。
9個のデータを左下から順番に4桁の地域番号を入力する。
************************************************************************
***** 標高データを配列に代入するプログラム(import)
*********************
c
subroutine(a,e,n)
real a(4800,4800),k(1600,1600)
integer i, j,
n,x,y
character mn*4
************************************************************************
do n=1,9
write(*,*) '地域番号を入れてください。'
read(*,*)mn
open(unit=77,file=mn//'.bin',access='direct',
+
form='unformatted',
recl=6400)
do
i=1,1600
read(77,rec=i)(k(i,j),j=1,1600)
end do
close(77)
************************************************************************
if(n.ge.1.and.n.le.3)y=0
if(n.ge.4.and.n.le.6)y=1
if(n.ge.7.and.n.le.9)y=2
if(n.eq.1
.or. n.eq.4 .or. n.eq.7)x=0
if(n.eq.2
.or. n.eq.5 .or. n.eq.8)x=1
if(n.eq.3
.or. n.eq.6 .or. n.eq.9)x=2
do
i=1,1600
do
j=1,1600
a(i+1600*y,j+1600*x)=k(i,j)
end
do
end do
end do
open(unit=7,file='data.terrain.bin',access='direct'
+
,form='unformatted',recl=19200)
do
n=1,4800
write(7,rec=n)(a(n,m),m=1,4800)
end do
close(unit=7)
c write(*,*)a(3000,4000)
************************************************************************
c close(78)
************************************************************************
end
2.WGRIBを用いた外部気象データの作成
注):太字は任意のファイル名および値を入力する。
1.プログラムをコンパイルする。
コマンド・・・cc
wgrib.c > wgrib
2.GRIBファイルの内容物を確認する。
コマンド・・・wgrib GRIBファイル or wgrib –s GRIBファイル
3.地表面以外のデータを取り出す
コマンド・・・wgrib
GRIBファイル | tail –n 91 | wgrib –i GRIBファイル –o 出力ファイル名
4.各解析値をヘッダー無しのバイナリ形式で出力する。
コマンド・・・wgrib GRIBファイル | grep “:HGT:” | wgrib –i –bin –nh GRIBファイル –o出力ファイル名
5.足りない解析値を用意する。
6.CReSSの各解析値を読み込む順番にcatコマンド(UNIX)でつなげる。