!--------( GGRIDS : gen. surf.grid from DPAM line data or equiv.)-------- !-------- This program must be linked with C subroutine implicit none integer :: lwkdir, lrtrim real :: hgeoid real(8) :: v(30) character(80) :: wdr, flog, fnam1, fnam2 character(150) :: buf; character(20) :: sdtm, sgmdl, sgmdln character(24) :: ca, can; character(9) :: strim character(8) :: area; character(6) :: cu; character(3) :: cg integer :: ldr, ktyp, kg, ku, k1, k2, k3, k4, kmx integer :: i, lc, ka, kan, l, IER, i1, i2 integer :: nc, ncc, iorg, korg, ispa, ispb, ixl, iyl, mszx, mszy, mx, my real :: sr, xl, yl, smx, smy, amag, calt, valtc real(8) :: dxn, dye, vlat, vlon, valt, tlat, tlon, talt, wlat, wlon, walt call premsg('GGRIDS : gen. surf.grid from DPAM line data or equiv.') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr; fnam1 = wdr; fnam2 = wdr call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80)) call premsg('Select AM Line data type :') call premsg(' DPAM Line data (1), AMDB-GSJ Line data (2),') call premsg(' AMDB-NEDO Line data (3), StdLIN data (4),') call gparmi(' or general Random data (0) ==> ', ktyp) if ((ktyp < 0) .or. (ktyp > 4)) call abendm('invalid value') call gparma('Enter Input line/random data filename ==> ', & 81-ldr, fnam1(ldr:80)) open(1,file=fnam1,status='old') call gparmi(' Geodetic ? in WGS(1) or TokyoDatum(2) ==> ', kg) if (kg == 1) then cg = 'WGS' call gparmi(' Alt. ? above Ellipsoid(1) or Geoid(2) ==> ', ka) if (ka == 1) then ca = 'above ITRF-GRS' else if (ka == 2) then sgmdl = ' ' call gparma(' Geoid model ? 0:NGA, or others [Enter NAME] ==> ', & 20, sgmdl) if (sgmdl == '0') sgmdl = 'NGA' call sgeoid(sgmdl) l = lrtrim(sgmdl) if (l <= 12) then ca = 'above Geoid ' // sgmdl(1:l) else ca = 'above Geoid ' // sgmdl(1:9) // '...' endif else call abendm('Alt. reference unknown') endif else if (kg == 2) then cg = 'TYO' call premsg(' Input Alt. data is treated as "above Geoid NGA"') ca = 'above Geoid NGA'; sgmdl = 'NGA'; call sgeoid(sgmdl) else call abendm('Geodetic reference unknown') endif ku = 1; cu = ' Min. ' if (ktyp == 1) then k1 = 5; k2 = 6; k3 = 7; k4 = 9; kmx = 9; ku = 2; cu = ' Deg. ' else if (ktyp == 2) then k1 = 2; k2 = 3; k3 = 0; k4 = 4; kmx = 4 read(1,'(a)') buf if ((buf(1:2) /= '##') .or. (buf(27:28) /= 'ft')) & call abendm('seems non AMDB-GSJ data') read(buf,'(20x,f6.0)') calt valtc = calt * 0.3048 else if (ktyp == 3) then k1 = 3; k2 = 4; k3 = 10; k4 = 8; kmx = 10 read(1,'(a)') buf if (buf(1:8) /= '## NEDO') call abendm('seems non AMDB-NEDO data') else if (ktyp == 4) then k1 = 1; k2 = 2; k3 = 3; k4 = 4; kmx = 4 else !------ 1st data-line display (to answer the following query) do read(1,'(a)') buf call premsg(buf(1:lrtrim(buf))) if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) cycle if (buf(1:1) /= '#') exit enddo rewind(1) !--- call gparmi('LatLon unit ? in Min.(1), Deg.(2) or D:Min.(3) ==> ', ku) if (ku == 1) then cu = ' Min. ' else if (ku == 2) then cu = ' Deg. ' else if (ku == 3) then cu = 'D:Min.' else call abendm('invalid value') endif call gparmi2(' data positions of Lat., Lon. : ', k1,k2) if ((k1 <= 0) .or. (k2 <= 0) .or. (k2==k1)) call abendm('invalid value') call gparmi2(' data positions of Alt., Mag. : ', k3,k4) if ((k3 <= 0) .or. (k3==k1) .or. (k3==k2) .or. & (k4 <= 0) .or. (k4==k1) .or. (k4==k2) .or. (k4==k3)) & call abendm('invalid value') kmx = max(k1, k2, k3, k4) if (ku == 3) then if (((k1+1)==k2) .or. ((k1+1)==k3) .or. ((k1+1)==k4) .or. & ((k2+1)==k1) .or. ((k2+1)==k3) .or. ((k2+1)==k4)) & call abendm('invalid value') kmx = max(k1+1, k2+1, k3, k4) endif if (kmx > 30) call abendm('too large data position') endif call premsg('Specify trim radius in km') call gparmf(' ( 0 for no-trimming ) ==> ', sr) if (sr < 0.) call abendm('invalid value') write(strim,'(f6.2,1x,a2)') sr, 'km' if (sr == 0.) strim = ' no-trim' call gparma('Enter Name-label for grid-data ==> ', 8, area) call premsg('Specify gridding parameters') call gparmi(' Coordinate Number (+800 for Old-TOKYO) : ', ncc) nc = ncc if (ncc >= 800) nc = ncc - 800 ispa = 0; ispb = 0 if ((nc >= 1) .and. (nc <= 62)) then iorg = 0 korg = (nc-30)*360 - 180 if (nc > 60) korg = 0 if (nc == 61) iorg = +5400 if (nc == 62) iorg = -5400 else if (nc == 65) then iorg = 0 call gparmi2(' Central Longitude ? (Deg. Min.) : ', i1,i2) korg = i1*60 + i2 else if (nc == 72) call abendm('nc=72 not implemented') call gparmi2(' Origin Latitude ? (Deg. Min.) : ', i1,i2) iorg = i1*60 + i2 call gparmi2(' Origin Longitude ? (Deg. Min.) : ', i1,i2) korg = i1*60 + i2 endif if (ncc < 800) then call gparmi(' Output Alt. ? above Ellipsoid(1) or Geoid(2) ==> ', kan) if (kan == 1) then can = 'above ITRF-GRS ' else if (kan == 2) then sgmdln = ' ' call gparma(' Geoid model ? 0:NGA, or others [Enter NAME] ==> ', & 20, sgmdln) if (sgmdln == '0') sgmdln = 'NGA' call sgeoid(sgmdln) l = lrtrim(sgmdln) if (l <= 12) then can = 'above Geoid ' // sgmdln(1:l) else can = 'above Geoid ' // sgmdln(1:9) // '...' endif else call abendm('Alt. reference unknown') endif else call premsg(' Output Alt. data is "above Geoid NGA"') can = 'above Geoid NGA'; sgmdln = 'NGA'; call sgeoid(sgmdln) endif 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 : ', mszx) if (mszx <= 0) call abendm('invalid value') smx = float(mszx) / 1000. call gparmi(' Number of Mesh to N : ', mx) if (mx <= 0) call abendm('invalid value') call premsg(' Mesh size and numbers to E ?') call gparmi(' Unit Mesh size in m : ', mszy) if (mszy <= 0) call abendm('invalid value') smy = float(mszy) / 1000. call gparmi(' Number of Mesh to E : ', my) if (my <= 0) call abendm('invalid value') call gparma('Enter grid-data output filename ==> ', 81-ldr, fnam2(ldr:80)) call clspin() call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb)) write(6,'(a)') '----------------------------------------' write(6,'(a)') 'GGRIDS : gen. surface grid from DPAM line data or equiv.' write(6,'(a,a)') ' Input line/random data file : ', fnam1(1:lrtrim(fnam1)) write(6,'(8x,i2,5a,i2,2a)') & k1, ': ', cg, '-Lat.(', cu, ') ', k3, ': Alt.(m) ', ca, & k2, ': ', cg, '-Lon.(', cu, ') ', k4, ': Mag.(nT)', '' write(6,'(a,a)') ' trim radius :', strim write(6,'(a,a)') ' Surface grid output filename : ', fnam2(1:lrtrim(fnam2)) write(6,'(8x,a,a)') 'Mag. anomaly, and Surface elevation ', can write(6,'(a12,a8,6x,a6,i4,4x,4i8)') & 'Areaname : ', area, 'Coord.', ncc, iorg, korg, ispa, ispb write(6,'(4x,a)') & ' (mesh) xs ys msz-x msz-y mx my alt ' write(6,'(4x,a12,2f10.2,2i8,2i6,a8)') & 'surface grid', xl, yl, mszx, mszy, mx, my, ' var.' open(2,file=fnam2,status='new') write(2,'(a,a,a)') '# output from GGRIDS ( Alt.: ', can, ' )' write(2,'(a,a)') '# input line/random data file : ', fnam1(1:lrtrim(fnam1)) write(2,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixl,iyl, mszx,mszy, mx,my, 9999.9, 0. close(2) call gsurf0(xl, yl, smx, smy, mx, my) do read(1,'(a)',iostat=IER) buf if (IER < 0) exit if ((buf(1:1)=='&') .or. (buf(1:1)=='%') .or. (buf(1:1)=='#')) then call premsg(buf(1:lrtrim(buf))); cycle endif if (ktyp == 4) then i = index(buf,'N') if (i /= 0) buf(i:i) = ' ' i = index(buf,'E') if (i /= 0) buf(i:i) = ' ' i = index(buf,'m') if (i /= 0) buf(i:i) = ' ' i = index(buf,'nT') if (i /= 0) buf(i:i+1) = ' ' read(buf,*) (v(i),i=1,4) else if (ktyp == 0) then do lc = index(buf,':') if (lc == 0) exit buf(lc:lc) = ' ' enddo endif read(buf,*) (v(i),i=1,kmx) endif amag = real(v(k4)) if (ku == 1) then vlat = v(k1); vlon = v(k2) else if (ku == 2) then vlat = v(k1) * 60.; vlon = v(k2) * 60. else vlat = v(k1)*60. + v(k1+1); vlon = v(k2)*60. + v(k2+1) endif ! [Let dxn/dye be in Output side spec (WGS/TYO)] if ((kg == 1) .and. (ncc >= 800)) then call xw84td(vlat, vlon, 0.d0, tlat, tlon, talt) call cvdiken(tlat, tlon, dye, dxn) ! [now dxn/dye are in TYO] else if ((kg == 2) .and. (ncc < 800)) then call xtw84d(vlat, vlon, 0.d0, wlat, wlon, walt) call cvdiken(wlat, wlon, dye, dxn) ! [now dxn/dye are in WGS] else call cvdiken(vlat, vlon, dye, dxn) ! [now dxn/dye are in common spec of IN/OUT (WGS or TYO)] endif call gsurf2(1, sngl(dxn), sngl(dye), amag) enddo rewind(1) call gsurf2(0, 0., 0., 0.) call gsurf3(sr, 9999.9) call gsurf4(fnam2) call gsurf5() open(2,file=fnam2,access='append') write(2,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixl,iyl, mszx,mszy, mx,my, 9999.9, -1. close(2) call gsurf0(xl, yl, smx, smy, mx, my) do read(1,'(a)',iostat=IER) buf if (IER < 0) exit if ((buf(1:1)=='&') .or. (buf(1:1)=='%') .or. (buf(1:1)=='#')) then call premsg(buf(1:lrtrim(buf))); cycle endif if (ktyp == 4) then i = index(buf,'N') if (i /= 0) buf(i:i) = ' ' i = index(buf,'E') if (i /= 0) buf(i:i) = ' ' i = index(buf,'m') if (i /= 0) buf(i:i) = ' ' i = index(buf,'nT') if (i /= 0) buf(i:i+1) = ' ' read(buf,*) (v(i),i=1,4) else if (ktyp == 0) then do lc = index(buf,':') if (lc == 0) exit buf(lc:lc) = ' ' enddo endif read(buf,*) (v(i),i=1,kmx) endif if (ku == 1) then vlat = v(k1); vlon = v(k2) else if (ku == 2) then vlat = v(k1) * 60.; vlon = v(k2) * 60. else vlat = v(k1)*60. + v(k1+1); vlon = v(k2)*60. + v(k2+1) endif if (ktyp == 2) valt = valtc if (ktyp /= 2) valt = v(k3) ! [Let vlat/vlon/valt be in Output side spec] if (kg == 1) then if (ncc < 800) then ! (vlat/vlon are in WGS, common IN/OUT) if (ka == 1) then if (kan == 1) then ! [valt is above WGS Ellips, common IN/OUT] else valt = valt - hgeoid(sngl(vlat), sngl(vlon)) ! [now valt is above WGS Geoid (conv. from WGS Ellips)] endif else if (kan == 1) then valt = valt + hgeoid(sngl(vlat), sngl(vlon)) ! [now valt is above WGS Ellips (conv. from WGS geoid)] else if (sgmdl.ne.sgmdln) then call sgeoid(sgmdl) valt = valt + hgeoid(sngl(vlat), sngl(vlon)) call sgeoid(sgmdln) valt = valt - hgeoid(sngl(vlat), sngl(vlon)) ! [now valt is above WGS Geoid (conv. from another model)] endif endif endif else ! (Convert vlat/vlon from WGS to TYO) if (ka == 1) then valt = valt - hgeoid(sngl(vlat), sngl(vlon)) ! [now valt is above Geoid [TYO] (conv. from WGS Ellips)] else if (sgmdl.ne.sgmdln) then call sgeoid(sgmdl) valt = valt + hgeoid(sngl(vlat), sngl(vlon)) call sgeoid(sgmdln) valt = valt - hgeoid(sngl(vlat), sngl(vlon)) ! [now valt is above Geoid [TYO] (conv. from another model)] endif endif call xw84td(vlat, vlon, valt, tlat, tlon, talt) vlat = tlat; vlon = tlon ! [now vlat/vlon are in TYO (conv. from WGS)] endif else if (ncc < 800) then ! (Convert vlat/vlon from TYO to WGS) call xtw84d(vlat, vlon, valt, wlat, wlon, walt) vlat = wlat; vlon = wlon ! [now vlat/vlon are in WGS (conv. from TYO)] if (kan == 1) then valt = valt + hgeoid(sngl(vlat), sngl(vlon)) ! [now valt is above WGS Ellips (conv. from Geoid [TYO]) else if (sgmdl.ne.sgmdln) then call sgeoid(sgmdl) valt = valt + hgeoid(sngl(vlat), sngl(vlon)) call sgeoid(sgmdln) valt = valt - hgeoid(sngl(vlat), sngl(vlon)) ! [now valt is above WGS Geoid (conv. from Geoid [TYO])] endif endif else ! (vlat/vlon are in TYO, common IN/OUT) ! (valt is above Geoid [TYO}, common IN/OUT) endif endif call cvdiken(vlat, vlon, dye, dxn) call gsurf2(1, sngl(dxn), sngl(dye), sngl(valt)) enddo close(1) call gsurf2(0, 0., 0., 0.) call gsurf3(sr, 9999.9) call gsurf4(fnam2) call gsurf5() if (flog(ldr:ldr) /= ' ') then open(99,file=flog,access='append') write(99,'(//a)') & 'GGRIDS : gen. surface grid from DPAM line data or equiv.' write(99,'(a,a)') & ' Input line/random data file : ', fnam1(1:lrtrim(fnam1)) write(99,'(8x,i2,5a,i2,2a)') & k1, ': ', cg, '-Lat.(', cu, ') ', k3, ': Alt.(m) ', ca, & k2, ': ', cg, '-Lon.(', cu, ') ', k4, ': Mag.(nT)', '' write(99,'(a,a)') ' trim radius :', strim write(99,'(a,a)') & ' Surface grid output filename : ', fnam2(1:lrtrim(fnam2)) write(99,'(a12,a8,6x,a6,i4,4x,4i8)') & 'Areaname : ', area, 'Coord.', ncc, iorg, korg, ispa, ispb write(99,'(4x,a)') & ' (mesh) xs ys msz-x msz-y mx my alt ' write(99,'(4x,a12,2f10.2,2i8,2i6,a8)') & 'surface grid', xl, yl, mszx, mszy, mx, my, ' var.' call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif stop end