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. に置換，
(6) 作成したファイルは，UTM図法のゾーン番号に対応したサブディレクトリ
に区分けして保存する．
この処理を行う，FORTRAN90プログラム 'genDEM40m.f90'
のソースコードを 下記に示す．
旧東京測地系の歪み補正については，
国土地理院が提供している「拡張パラメータファイル」
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度の範囲 (UTMゾーン番号jj ) のデータが，
'dem40/Zjj/NiiEkkk.alt'
の path名 のファイルとなる．
（すなわち，GSJ50dem ディレクトリからデータ入力し，dem40/
ディレクトリ下の UTMゾーン番号別サブディレクトリへデータ出力する．）
なお，実行に先立ち，データ出力先のディレクトリ dem40/ と
サブディレクトリ dem40/Z51/, dem40/Z52/, dem40/Z53/,
dem40/Z54/, dem40/Z55/, dem40/Z56/
を作成しておくものとする．その実行時間は，大容量データファイルの処理のため，
何10分かを要するかも知れない．
gtopo プログラムでの利用のためには， 上記の dem40 ディレクトリの下の 6サブディレクトリに分割して出力された 計120個のファイル（1ファイル約45MBの固定サイズ）を，そのツリー構造を維持した 一つのディレクトリ構造 dem40 として， 絶対path が '/home/SHARE/data/DEM/' なるディレクトリ直下に 配置すればよい．
!---- Source program of generating DEM40 data files
!
integer :: mh(200,200)
real :: h(2800,2000), a(2401,2401)
character(24) :: fnam1, fnam2;
character(8) :: area; character(6) :: cmesh, code
logical :: exist
do lli=20,45; do llk=122,153
llj = llk/6 + 31
llc = lli*1000 + llk
if ((lli < 24) .and. (llc /= 20136)) cycle
if ((llk > 145) .and. (llc /= 24153)) cycle
write(fnam2,'(a7,i2,a2,i2,a1,i3,a4)') &
'dem40/Z', llj, '/N', lli, 'E', llk, '.alt'
kc = 0
if ((llc == 25131) .or. (llc == 32139)) kc = 1
if ((llc >= 24122) .and. (llc <= 24125)) kc = 1
do i=1,2800; do k=1,2000; h(i,k) = -1.; enddo; enddo
key = 0
ins = lli*12 - 1
kns = (llk-100)*8 - 1
do in=ins,ins+13; do 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 /= fnam1(15:20)) call abendm('mesh-code error')
do n=1,200
read(1,'(a6,i3,200i5)') code, nn, (mh(k,201-n),k=1,200)
if (code /= cmesh) call abendm('mesh-code conflict')
if (nn /= n) call abendm('format error')
enddo
close(1)
is = (in-ins) * 200
ks = (kn-kns) * 200
do i=1,200; do k=1,200
if (mh(k,i) == -9999) cycle
if (mh(k,i) <= 0) then
h(is+i,ks+k) = 0.
else
h(is+i,ks+k) = float(mh(k,i)) / 10.
endif
enddo; enddo
key = 1
endif
enddo; enddo
if (key == 0) cycle
key = 0
wlat0 = float(lli*60); tlat0 = wlat0 - 5.0125
wlon0 = float(llk*60); tlon0 = wlon0 - 7.51875
do kw=1,2401; do 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 == 1) then
if ((tlat >= 1920.) .and. (tlat <= 1960.) .and. &
(tlon >= 8340.) .and. (tlon <= 8400.)) then
tlat = tlat - 0.0404
tlon = tlon - 0.0023
else if ((tlat >= 1520.) .and. (tlat <= 1560.) .and. &
(tlon >= 7860.) .and. (tlon <= 7920.)) then
tlat = tlat + 0.2016
tlon = tlon - 0.3138
else if ((tlat >= 1440.) .and. (tlat <= 1500.) .and. &
(tlon >= 7320.) .and. (tlon <= 7560.)) then
if (tlon >= 7500.) then
tlat = tlat + 0.0310
tlon = tlon - 0.0401
else if (tlon >= 7470.) then
tlat = tlat + 0.1527
tlon = tlon - 0.2867
else if (tlon >= 7443.) then
tlat = tlat - 0.0784
tlon = tlon - 0.1221
else if (tlon <= 7410.) then
tlat = tlat - 0.0782
tlon = tlon - 0.1196
endif
endif
endif
dlat = tlat - tlat0; ip = ifix(dlat / 0.025)
dlon = tlon - tlon0; kp = ifix(dlon / 0.0375)
if ((ip <= 0) .or. (ip >= 2800) .or. &
(kp <= 0) .or. (kp >= 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
flat = dlat - float(ip)*0.025
flon = dlon - float(kp)*0.0375
if ((h(ip,kp) < 0.) .or. (h(ip+1,kp+1) < 0.) .or. &
(h(ip,kp+1) < 0.) .or. (h(ip+1,kp) < 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
enddo; enddo
if (key == 1) then
write(6,*) fnam2
area = fnam2(11:17) // ' '
open(2,file=fnam2,status='new')
write(2,'(a8,4x,i4,4i8)') area, 199, 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 kw=1,2401
write(2,'((f7.1,9(1x,f7.1)))') (a(iw,kw),iw=1,2401)
enddo
close(2)
endif
enddo; enddo
stop
end