!--------( GTOPO : generate topography grid-data from DEM )-------- implicit none integer :: lwkdir, lrtrim character(*),parameter :: DEMDIR = '/home/SHARE/data/DEM/dem40/' real,allocatable :: g(:,:) real :: f1(2401,2401), f2(2401,2401), f3(2401,2401) character(80) :: wdr, flog, fnam1, fnam2, buf character(20) :: sdtm character(8) :: area, narea character(11) :: cik integer :: ldr, i, j, ii, jj, key, mi, mk integer :: i1, i2, nik, nik1, nik2, nik3, ksel, ix, iy integer :: ido, kdo, kz, mmi1, mmk1, mmi2, mmk2, mmi3, mmk3 integer :: ncc, iorg, korg, ncn, nccn, iorgn, korgn, ispan, ispbn integer :: ixs, iys, mszx, mszy, mx, my integer :: ixl, iyl, mszxn, mszyn, mxn, myn real :: vnul, alt, xl, yl, xn, ye, fi, fk, fit, fkt, fif, fkf real :: uif1, ukf1, uif2, ukf2, uif3, ukf3, si, sk, dif, dkf real :: a0, a1, a2, a3 logical :: cond 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 TOPOgrid ==> ', 8, narea) call gparmi('Select Coordinate Number ==> ', nccn) ncn = nccn if (nccn >= 800) ncn = nccn - 800 ispan = 0 ispbn = 0 if ((ncn >= 1) .and. (ncn <= 62)) then iorgn = 0 korgn = (ncn-30)*360 - 180 if (ncn > 60) korgn = 0 if (ncn == 61) iorgn = +5400 if (ncn == 62) iorgn = -5400 else if (ncn == 65) then iorgn = 0 call gparmi(' Central Longitude ? (Deg.) : ', i1) korgn = i1*60 else if (ncn == 72) call abendm('nc=72 not supported') call gparmi2(' Origin Latitude ? (Deg. Min.) : ', i1,i2) iorgn = i1*60 + i2 call gparmi2(' Origin Longitude ? (Deg. Min.) : ', i1,i2) korgn = i1*60 + i2 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 <= 0) call abendm('invalid value') call gparmi(' Number of Mesh to N : ', mxn) call premsg(' Mesh size and numbers to E ?') call gparmi(' Unit Mesh size in m : ', mszyn) if (mszyn <= 0) call abendm('invalid value') call gparmi(' Number of Mesh to E : ', myn) allocate(g(mxn,myn)); 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,'(a12,a8,6x,a6,i4,4x,4i8)') & 'Areaname : ', narea, 'Coord.', nccn, iorgn, korgn, ispan, ispbn write(6,'(4x,a)') & ' (mesh) xs ys msz-x msz-y mx my' write(6,'(4x,a12,2f10.2,2i8,2i6)') & 'gen.TOPOgrid', xl, yl, mszxn, mszyn, mxn, myn call cvinit(nccn, float(iorgn), float(korgn), float(ispan), float(ispbn)) nik1 = 0; nik2 = 0; nik3 = 0; ksel = 3 do j=1,myn iy = iyl + (j-1)*mszyn ye = float(iy) / 1000. do i=1,mxn ix = ixl + (i-1)*mszxn xn = float(ix) / 1000. if (nccn >= 800) 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.); if (fi < 0.) ido = ido - 1 fif = fi - float(ido*60) kdo = ifix(fk/60.); if (fk < 0.) kdo = kdo - 1 fkf = fk - float(kdo*60) if (fi == 5400.) then; ido = 89; fif = 60.; endif if (fk >= 10800.) then; ido = ido - 360; endif nik = ido*1000 + kdo kz = (kdo+186) / 6 write(cik,'(a1,i2.2,a1,i3.3,i4.4)') 'Z', kz, '/', ido, kdo if (cik(5:5) == '-') then; cik(5:5) = 'S' else; cik(5:5) = 'N'; endif if (cik(8:8) == '-') then; cik(8:8) = 'W' else; cik(8:8) = 'E'; endif key = 0 if (nik == nik1) then ksel = 1 else if (nik == nik2) then ksel = 2 else if (nik == nik3) then ksel = 3 else fnam1 = DEMDIR // cik // '.alt' inquire(file=fnam1,exist=cond) if (cond) then open(1,file=fnam1,status='old') read(1,'(a)') buf do while (buf(1:1) == '#') read(1,'(a)') buf enddo if (buf(13:16) == ' ') then !-- v2005 read(buf,'(a8,i4,4x,2i8)') area, ncc, iorg, korg if (ncc /= 399) & call abendm('DEM coordinate (nc) invalid: ' // cik // '.alt') else !-- v2018 read(buf,'(a8,4x,i4,2i8)') area, ncc, iorg, korg if (ncc /= 199) & call abendm('DEM coordinate (nc) invalid: ' // cik // '.alt') endif read(1,*) ixs,iys, mszx,mszy, mx,my, vnul, alt if (((iorg*1000+ixs) /= ido*60000) .or. & ((korg*1000+iys) /= kdo*60000) .or. & (mszx*(mx-1) /= 60000) .or. (mszy*(my-1) /= 60000)) & call abendm('DEM data range invalid: ' // cik // '.alt') if (ksel == 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 == 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 == 0) then if (ksel == 1) then si = fif / uif1 mi = ifix(si) if (mi == mmi1) mi = mmi1 - 1 dif = si - float(mi) sk = fkf / ukf1 mk = ifix(sk) if (mk == 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 == 2) then si = fif / uif2 mi = ifix(si) if (mi == mmi2) mi = mmi2 - 1 dif = si - float(mi) sk = fkf / ukf2 mk = ifix(sk) if (mk == 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 == mmi3) mi = mmi3 - 1 dif = si - float(mi) sk = fkf / ukf3 mk = ifix(sk) if (mk == 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 < 0.) .or. (a1 < 0.) .or. (a2 < 0.) .or. (a3 < 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 enddo enddo write(2,'(a8,4x,i4,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 j=1,myn write(2,'((f7.1,9(1x,f7.1)))') (g(i,j),i=1,mxn) enddo close(2) if (flog(ldr:ldr) /= ' ') then open(99,file=flog,access='append') write(99,'(//a)') 'GTOPO : generate topography grid-data from DEM' write(99,'(a,a)') ' Output new file : ', fnam2(1:lrtrim(fnam2)) write(99,'(a12,a8,6x,a6,i4,4x,4i8)') & 'Areaname : ', narea, 'Coord.', nccn, iorgn,korgn, ispan,ispbn write(99,'(4x,a)') & ' (mesh) xs ys msz-x msz-y mx my' write(99,'(4x,a12,2f10.2,2i8,2i6)') & 'gen.TOPOgrid', xl, yl, mszxn, mszyn, mxn, myn call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif stop end