c--------( GTOPO : generate topography grid-data from DEM )-------- parameter (LXO=2401, LYO=2401, LXN=4001, LYN=4001) character*26 DEMDIR parameter (DEMDIR = '/home/SHARE/data/DEM/wm40/') c dimension f1(LXO,LYO), f2(LXO,LYO), f3(LXO,LYO), g(LXN,LYN) character*80 wdr, flog, fnam1, fnam2 character area*8, narea*8, buf*80, sdtm*20 logical exist c call premsg('GTOPO : generate topography grid-data from DEM') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr fnam2 = wdr call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80)) call premsg('') call gparma('Enter Name-label for topo grid ==> ', 8, narea) call gparmi('Select Projection-Number ==> ', nccn) ncn = nccn if(nccn.ge.200) ncn = nccn - 200 ispan = 0 ispbn = 0 if(ncn.ge.1.and.ncn.le.62) then iorgn = 0 korgn = (ncn-30)*360 - 180 if(ncn.gt.60) korgn = 0 if(ncn.eq.61) iorgn = +5400 if(ncn.eq.62) iorgn = -5400 else if(ncn.eq.65) then iorgn = 0 call gparmi(' Central Longitude ? (Deg.) : ', i1) korgn = i1*60 else call gparmi2(' Origin Latitude ? (Deg. Min.) : ', i1,i2) iorgn = i1*60 + i2 call gparmi2(' Origin Longitude ? (Deg. Min.) : ', i1,i2) korgn = i1*60 + i2 if(ncn.eq.72) then call gparmi2(' St.Parallel-1 Lat. ? (Deg. Min.) : ', i1,i2) ispan = i1*60 + i2 call gparmi2(' St.Parallel-2 Lat. ? (Deg. Min.) : ', i1,i2) ispbn = i1*60 + i2 endif endif call premsg('Specify gridding parameters') call premsg(' Southwest corner Coord. in km ?') call gparmf(' Northing : ', xl) ixl = nint(xl*1000.) call gparmf(' Easting : ', yl) iyl = nint(yl*1000.) call premsg(' Mesh size and numbers to N ?') call gparmi(' Unit Mesh size in m : ', mszxn) if(mszxn.le.0) call abendm('invalid value') call gparmi(' Number of Mesh to N : ', mxn) if(mxn.gt.LXN) call abendm('too large array size') call premsg(' Mesh size and numbers to E ?') call gparmi(' Unit Mesh size in m : ', mszyn) if(mszyn.le.0) call abendm('invalid value') call gparmi(' Number of Mesh to E : ', myn) if(myn.gt.LYN) call abendm('too large array size') call gparma('Enter new grid output filename ==> ', * 81-ldr, fnam2(ldr:80)) call clspin() open(2,file=fnam2,status='new') write(6,'(a)') '----------------------------------------' write(6,'(a)') 'GTOPO : generate topography grid-data from DEM' write(6,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2)) write(6,1000) narea, nccn, iorgn,korgn, ispan,ispbn, * xl, yl, mszxn, mszyn, mxn, myn c call cvinit(nccn, float(iorgn), float(korgn), * float(ispan), float(ispbn)) nik1 = 0 nik2 = 0 nik3 = 0 ksel = 3 do 10 j=1,myn iy = iyl + (j-1)*mszyn ye = float(iy) / 1000. do 10 i=1,mxn ix = ixl + (i-1)*mszxn xn = float(ix) / 1000. if (nccn .lt. 200) then call cvenik(ye, xn, fit, fkt) call xtw84(fit, fkt, 0., fi, fk, alt) else call cvenik(ye, xn, fi, fk) endif ido = ifix(fi/60.) fif = fi - float(ido*60) kdo = ifix(fk/60.) fkf = fk - float(kdo*60) nik = ido*1000 + kdo key = 0 if (nik .eq. nik1) then ksel = 1 else if (nik .eq. nik2) then ksel = 2 else if (nik .eq. nik3) then ksel = 3 else write(fnam1,'(a,i5,a)') DEMDIR, nik, '.alt' inquire(file=fnam1,exist=exist) if (exist) then open(1,file=fnam1,status='old') 1 read(1,'(a)') buf if (buf(1:1).eq.'#') goto 1 read(buf,'(a8,i4,4x,4i8)') area,ncc, iorg,korg,ispa,ispb if (ncc.ne.399) call abendm('Illegal DEM file spec.') if (iorg.ne.(ido*60) .or. korg.ne.(kdo*60)) * call abendm('DEM file broken') read(1,*) ixs,iys, mszx,mszy, mx,my, vnul, alt if (ixs.ne.0.or.iys.ne.0) call abendm('DEM file broken') if (mx.gt.LXO.or.my.gt.LYO) * call abendm('too large array size') if (ksel.eq.3) then ksel = 1 nik1 = nik uif1 = float(mszx) / 1000. ukf1 = float(mszy) / 1000. mmi1 = mx - 1 mmk1 = my - 1 read(1,*) ((f1(ii,jj),ii=1,mx),jj=1,my) else if (ksel.eq.1) then ksel = 2 nik2 = nik uif2 = float(mszx) / 1000. ukf2 = float(mszy) / 1000. mmi2 = mx - 1 mmk2 = my - 1 read(1,*) ((f2(ii,jj),ii=1,mx),jj=1,my) else ksel = 3 nik3 = nik uif3 = float(mszx) / 1000. ukf3 = float(mszy) / 1000. mmi3 = mx - 1 mmk3 = my - 1 read(1,*) ((f3(ii,jj),ii=1,mx),jj=1,my) endif close(1) else key = 1 g(i,j) = -1. endif endif if (key.eq.0) then if (ksel .eq. 1) then si = fif / uif1 mi = ifix(si) if (mi.eq.mmi1) mi = mmi1 - 1 dif = si - float(mi) sk = fkf / ukf1 mk = ifix(sk) if (mk.eq.mmk1) mk = mmk1 - 1 dkf = sk - float(mk) a0 = f1(mi+1,mk+1) a1 = f1(mi+1,mk+2) a2 = f1(mi+2,mk+1) a3 = f1(mi+2,mk+2) else if (ksel.eq.2) then si = fif / uif2 mi = ifix(si) if (mi.eq.mmi2) mi = mmi2 - 1 dif = si - float(mi) sk = fkf / ukf2 mk = ifix(sk) if (mk.eq.mmk2) mk = mmk2 - 1 dkf = sk - float(mk) a0 = f2(mi+1,mk+1) a1 = f2(mi+1,mk+2) a2 = f2(mi+2,mk+1) a3 = f2(mi+2,mk+2) else si = fif / uif3 mi = ifix(si) if (mi.eq.mmi3) mi = mmi3 - 1 dif = si - float(mi) sk = fkf / ukf3 mk = ifix(sk) if (mk.eq.mmk3) mk = mmk3 - 1 dkf = sk - float(mk) a0 = f3(mi+1,mk+1) a1 = f3(mi+1,mk+2) a2 = f3(mi+2,mk+1) a3 = f3(mi+2,mk+2) endif if (a0.lt.0..or.a1.lt.0..or.a2.lt.0..or.a3.lt.0.) then g(i,j) = -1. else g(i,j) = (a0*(1.-dif) + a2*dif) * (1.-dkf) * + (a1*(1.-dif) + a3*dif) * dkf endif endif 10 continue write(2,'(a8,i4,4x,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixl, iyl, mszxn,mszyn, mxn,myn, 9999.9, -1. do 100 j=1,myn write(2,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mxn) 100 continue close(2) c if(flog(ldr:ldr).ne.' ') then open(99,file=flog,access='append') write(99,'(//a,$)') 'GTOPO :' write(99,'(a)') ' generate topography grid-data from DEM' write(99,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2)) write(99,1000) narea, nccn, iorgn,korgn, ispan,ispbn, * xl, yl, mszxn, mszyn, mxn, myn call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif stop c 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/ * 3x, ' (mesh) ', ' xs ', ' ys ', * ' msz-x', ' msz-y', ' mx', ' my'/ * 3x, 'gen.topo grid', 2f10.2, 2i8, 2i6) end