GDMP プログラム群中の gtopo プログラムは，DEM40 データから
地形高度グリッドデータを作成するが，その DEM40 データとは，国土地理院発行の
数値地図50mメッシュ (標高) の CD-ROM 3枚（日本-Ⅰ，日本-Ⅱ，日本-Ⅲ）
[訂正情報を含む] を用いて，日本の地表高データを再構成したものである．
その再構成の内容は，
(1) 旧東京測地系から世界測地系へ，旧東京測地系の歪み補正を含めて 変換，
(2) Mesh size を，緯度 1/40 分，経度 3/80 分 から 緯経度とも 1/40分 に，
(3) 緯度範囲・経度範囲とも 1度 の区画ごとのデータファイルに，
(4) 細分Mesh中心位置の値 から 細分Meshの格子点での値へ，
(5) 負の陸域標高値を 0. に置換，海域を示すデータを -1. に置換，
である．
この処理を行う，FORTRANプログラムのソースコードを 下記に示す．
旧東京測地系の歪み補正については，
国土地理院が提供している「拡張パラメータファイル」
http://www.gsi.go.jp/MAP/CD-ROM/sekaitaiou/2500Conv/Enhance_Par.exe
の情報を参考に，必要な精度を考慮して，数多くない地域に区分して 測地系変換
サブプログラム 'xtw84' の結果に必要な補正を加えるようにしてある．
このプログラムの入力データは，国土地理院による元データにならって，
1次メッシュコードが 'xxxx' で 2次メッシュコードが 'yy'
で表される範囲のデータが ファイルのpath名
'GSI50dem/xxxx/xxxxyy.MEM' で表現されるものとし，
出力データは，緯度 'ii' 度・経度 'kkk' 度 を南西角とする
緯経度とも 1度 の範囲のデータが，path名で 'wm40/iikkk.alt'
のファイルとなる．（すなわち，GSJ50dem サブディレクトリからデータ入力し，
wm40 サブディレクトリへデータ出力する．）
gtopo プログラムでの利用のためには， 上記の wm40 サブディレクトリへ出力される 120個のファイル （1ファイル約45MBの固定サイズ）を，絶対path が '/home/SHARE/data/DEM/wm40/' なるディレクトリ内に配置すればよい．
c---- Source program of generating DEM40 data files
c
dimension mh(200,200), h(2800,2000), a(2401,2401)
character fnam1*24, fnam2*14, cmesh*6, code*6, area*8
logical exist
c
do 100 lli=20,45
do 100 llk=122,153
llc = lli*1000 + llk
write(fnam2,'(a5,i5,a4)') 'wm40/', llc, '.alt'
if (lli.lt.24 .and. llc.ne.20136) goto 100
if (llk.gt.145 .and. llc.ne.24153) goto 100
c
kc = 0
if (llc.eq.25131 .or. llc.eq.32139) kc = 1
if (llc.ge.24122 .and. llc.le.24125) kc = 1
do 200 i=1,2800
do 200 k=1,2000
h(i,k) = -1.
200 continue
key = 0
ins = lli*12 - 1
kns = (llk-100)*8 - 1
do 300 in=ins,ins+13
do 300 kn=kns,kns+9
ima = in / 8
imb = in - ima*8
kma = kn / 8
kmb = kn - kma*8
write(fnam1,'(a9,2i2,a1,2i2,2i1,a4)')
* 'GSI50dem/', ima,kma, '/', ima,kma,imb,kmb, '.MEM'
inquire(file=fnam1,exist=exist)
if (exist) then
open(1,file=fnam1,status='old')
read(1,'(a6)') cmesh
if (cmesh.ne.fnam1(15:20)) call abendm('mesh-code error')
do 400 n=1,200
read(1,'(a6,i3,200i5)') code, nn, (mh(k,201-n),k=1,200)
if (code.ne.cmesh) call abendm('mesh-code conflict')
if (nn.ne.n) call abendm('format error')
400 continue
close(1)
is = (in-ins) * 200
ks = (kn-kns) * 200
do 500 i=1,200
do 500 k=1,200
if (mh(k,i).eq.-9999) goto 500
if (mh(k,i).le.0) then
h(is+i,ks+k) = 0.
else
h(is+i,ks+k) = float(mh(k,i)) / 10.
endif
500 continue
key = 1
endif
300 continue
if (key.eq.0) goto 100
c
key = 0
wlat0 = float(lli*60)
wlon0 = float(llk*60)
tlat0 = wlat0 - 5.0125
tlon0 = wlon0 - 7.51875
do 600 kw=1,2401
do 600 iw=1,2401
wlat = wlat0 + float(iw-1)/40.
wlon = wlon0 + float(kw-1)/40.
call xw84t(wlat, wlon, 0., tlat, tlon, talt)
if (kc.eq.1) then
if (tlat.ge.1920. .and. tlat.le.1960. .and.
* tlon.ge.8340. .and. tlon.le.8400.) then
tlat = tlat - 0.0404
tlon = tlon - 0.0023
else if (tlat.ge.1520. .and. tlat.le.1560. .and.
* tlon.ge.7860. .and. tlon.le.7920.) then
tlat = tlat + 0.2016
tlon = tlon - 0.3138
else if (tlat.ge.1440. .and. tlat.le.1500. .and.
* tlon.ge.7320. .and. tlon.le.7560.) then
if (tlon.ge.7500.) then
tlat = tlat + 0.0310
tlon = tlon - 0.0401
else if (tlon.ge.7470.) then
tlat = tlat + 0.1527
tlon = tlon - 0.2867
else if (tlon.ge.7443.) then
tlat = tlat - 0.0784
tlon = tlon - 0.1221
else if (tlon.le.7410.) then
tlat = tlat - 0.0782
tlon = tlon - 0.1196
endif
endif
endif
dlat = tlat - tlat0
ip = ifix(dlat / 0.025)
flat = dlat - float(ip)*0.025
dlon = tlon - tlon0
kp = ifix(dlon / 0.0375)
if (ip.le.0.or.ip.ge.2800.or.kp.le.0.or.kp.ge.2000) then
write(6,*) llc
write(6,*) ' ip=',ip, ' kp=',kp, ' iw=',iw, ' kw=',kw
write(6,*) ' tlat=',tlat, ' tlon=',tlon
write(6,*) ' wlat=',wlat, ' wlon=',wlon
write(6,*) ' lat0=',tlat0, ' lon0=',tlon0
call abendm('out of range')
endif
c
flon = dlon - float(kp)*0.0375
if (h(ip,kp).lt.0..or.h(ip+1,kp+1).lt.0..or.
* h(ip,kp+1).lt.0..or.h(ip+1,kp).lt.0.) then
a(iw,kw) = -1.
else
h0 = h(ip,kp) + (h(ip+1,kp)-h(ip,kp))*flat
h1 = h(ip,kp+1) + (h(ip+1,kp+1)-h(ip,kp+1))*flat
a(iw,kw) = h0 + (h1-h0)*flon
key = 1
endif
600 continue
if (key.eq.0) goto 100
c
write(6,*) fnam2
area = 'w' // fnam2(6:10) // ' '
open(2,file=fnam2,status='new')
write(2,'(a8,i4,4x,4i8)') area, 399, lli*60, llk*60, 0, 0
write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')
* 0,0, 25,25, 2401,2401, 9999.9, -1.
do 700 kw=1,2401
write(2,'((f7.1,9(1x,f7.1)))') (a(iw,kw),iw=1,2401)
700 continue
close(2)
c
100 continue
stop
end