c------( GGRIDS : gen. surface grid from DPAM line data or equiv.)------
c------ This program must be linked with C language subroutine <gsurf.o>
c
      dimension v(30)
      character*80 wdr, flog, fnam1, fnam2
      character cu*4, cg*3, ca*14, can*14, area*8, strim*9, buf*150
c
      call premsg('GGRIDS : '
     *            // 'gen. surface 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)
      ku = 1
      cu = 'Min.'
      if (ktyp.eq.1) then
        k1 = 5
        k2 = 6
        k3 = 7
        k4 = 9
        kmx = 9
        ku = 2
        cu = 'Deg.'
      else if (ktyp.eq.2) then
        k1 = 2
        k2 = 3
        k3 = 0
        k4 = 4
        kmx = 4
      else if (ktyp.eq.3) then
        k1 = 3
        k2 = 4
        k3 = 10
        k4 = 8
        kmx = 10
      else if (ktyp.eq.4) then
        k1 = 1
        k2 = 2
        k3 = 3
        k4 = 4
        kmx = 4
      else if (ktyp.eq.0) then
      else
        call abendm('invalid Alt. value')
      endif
      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.eq.1) then
        cg = 'WGS'
        call premsg(' Alt. ?   above ITRF/GRS Ellipsoid (1)')
        call gparmi('            or Geoid / Tokyo-Datum (2) ==> ', ka)
        if (ka.eq.1) then
          ca = 'above ITRF-GRS'
        else if (ka.eq.2) then
          ca = 'above Geoid   '
        else
          call abendm('invalid Alt. value')
        endif
      else if (kg.eq.2) then
        cg = 'TYO'
        ka = 2
        ca = 'above Geoid   '
      else
        call abendm('invalid Alt. value')
      endif
      if (ktyp.eq.0) then
   10   read(1,'(a)',end=2) buf
        call premsg(buf(1:lrtrim(buf)))
        if (buf(1:1).eq.'&' .or. buf(1:1).eq.'%'
     *                      .or. buf(1:1).eq.'#') goto 10
        rewind(1)
        call gparmi('LatLon unit ?   in Min.(1) or Deg.(2) ==> ', ku)
        if (ku.eq.1) then
          cu = 'Min.'
        else if (ku.eq.2) then
          cu = 'Deg.'
        else
          call abendm('invalid Alt. value')
        endif
        call gparmi2('  data positions of  Lat., Lon. : ', k1,k2)
        if (k1.le.0.or.k2.le.0.or.k1.eq.k2)
     *       call abendm('invalid value')
        if (k1.gt.30.or.k2.gt.30)
     *       call premsg('too large data position')
        call gparmi2('  data positions of  Alt., Mag. : ', k3,k4)
        if (k3.le.0.or.k3.eq.k1.or.k3.eq.k2)
     *       call abendm('invalid value')
        if (k4.le.0.or.k4.eq.k1.or.k4.eq.k2.or.k4.eq.k3)
     *       call abendm('invalid value')
        kmx = max(k1, k2, k3, k4)
        if (kmx.gt.30) call abendm('too large data position')
      else if (ktyp.eq.2) then
        read(1,'(a)') buf
        if (buf(1:2).ne.'##' .or. buf(27:28).ne.'ft')
     *       call abendm('seems non AMDB-GSJ data')
        read(1,'(20x,f6.0,a2)') calt, buf
        valt = calt * 0.3048
      else if (ktyp.eq.3) then
        read(1,'(a)') buf
        if (buf(1:8).ne.'##  NEDO')
     *       call abendm('seems non AMDB-NEDO data')
      endif
c
      call premsg('Specify trim radius in km')
      call gparmf('     ( 0 for no-trimming ) ==> ', sr)
      if(sr.lt.0.) call abendm('invalid value')
      write(strim,'(f6.2,1x,a2)') sr, 'km'
      if(sr.eq.0.) strim = '  no-trim'
      call gparma('Enter Name-label for grid-data   ==> ', 8, area)
      call premsg('Specify gridding parameters')
      call gparmi('  Coordinate Number (+200 for WGS) : ', ncc)
      nc = ncc
      if(ncc.ge.200) nc = ncc - 200
      ispa = 0
      ispb = 0
      if(nc.ge.1.and.nc.le.62) then
        iorg = 0
        korg = (nc-30)*360 - 180
        if(nc.gt.60) korg = 0
        if(nc.eq.61) iorg = +5400
        if(nc.eq.62) iorg = -5400
      else
        call gparmi2('    Origin  Lat. & Long. in minutes  : ',
     *               iorg, korg)
        if(nc.eq.72) then
          call gparmi2('    two Standard Latitudes (in min.) : ',
     *                 ispa, ispb)
        endif
      endif
      if (ncc.ge.200) then
        call premsg(' Alt. ?   above ITRF/GRS Ellipsoid (1)')
        call gparmi('            or Geoid / Tokyo-Datum (2) ==> ', kan)
      else
        kan = 2
      endif
      if (kan.eq.1) then
        can = 'above ITRF-GRS'
      else if (kan.eq.2) then
        can = 'above Geoid   '
      else
        call abendm('invalid value')
      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.le.0) call abendm('invalid value')
      smx = float(mszx) / 1000.
      call gparmi('    Number of Mesh to N :  ', mx)
      if(mx.le.0) call abendm('invalid value')
      call premsg('  Mesh size and numbers to E ?')
      call gparmi('    Unit Mesh size in m :  ', mszy)
      if(mszy.le.0) call abendm('invalid value')
      smy = float(mszy) / 1000.
      call gparmi('    Number of Mesh to E :  ', my)
      if(my.le.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))
c
      open(2,file=fnam2,status='new')
      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,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xl, yl, mszx, mszy, mx, my, '    var.'
      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,i4,4x,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)
c
      call gsurf0(xl, yl, smx, smy, mx, my)
    1 read(1,'(a)',end=2) buf
      if (buf(1:1).eq.'&' .or. buf(1:1).eq.'%'
     *                    .or. buf(1:1).eq.'#') then
        call premsg(buf(1:lrtrim(buf)))
        goto 1
      endif
      if (ktyp.eq.4) then
        read(buf,'(f10.4,1x,f11.4,1x,f8.1,1x,f8.1,2x)') (v(i),i=1,4)
      else
        read(buf,*) (v(i),i=1,kmx)
      endif
      vlat = v(k1)
      vlon = v(k2)
      if (ku.eq.2) then
        vlat = vlat * 60.
        vlon = vlon * 60.
      endif
      if (kg.eq.1 .and. ncc.lt.200) then
        call xw84t(vlat, vlon, 0., tlat, tlon, alt)
        call cviken(tlat, tlon, ye, xn)
      else if (kg.eq.2 .and. ncc.ge.200) then
        call xtw84(vlat, vlon, 0., wlat, wlon, alt)
        call cviken(wlat, wlon, ye, xn)
      else
        call cviken(vlat, vlon, ye, xn)
      endif
      call gsurf2(1, xn, ye, v(k4))
      goto 1
    2 rewind(1)
      call gsurf2(0, 0., 0., 0.)
      call gsurf3(sr, 9999.9)
      call gsurf4(fnam2)
      call gsurf5()
c
      open(2,file=fnam2,status='old')
    3 read(2,*,end=4)
      goto 3
    4 write(2,'(a8,i4,4x,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)
c
      call gsurf0(xl, yl, smx, smy, mx, my)
    5 read(1,'(a)',end=6) buf
      if (buf(1:1).eq.'&' .or. buf(1:1).eq.'%'
     *                    .or. buf(1:1).eq.'#') then
        call premsg(buf(1:lrtrim(buf)))
        goto 5
      endif
      if (ktyp.eq.4) then
        read(buf,'(f10.4,1x,f11.4,1x,f8.1,1x,f8.1,2x)') (v(i),i=1,4)
      else
        read(buf,*) (v(i),i=1,kmx)
      endif
      vlat = v(k1)
      vlon = v(k2)
      if (ktyp.ne.2) valt = v(k3)
      if (ku.eq.2) then
        vlat = vlat * 60.
        vlon = vlon * 60.
      endif
      if (kg.eq.1) then
        if (ncc.ge.200) then
          if (ka.ne.kan) then
            call xw84t(vlat, vlon, valt, tlat, tlon, talt)
            if (ka.eq.1) then
              valt = talt
            else
              valt = valt+valt - talt
            endif
          endif
        else
          call xw84t(vlat, vlon, valt, tlat, tlon, talt)
          vlat = tlat
          vlon = tlon
          if (ka.eq.1) valt = talt
        endif
      else
        if (ncc.ge.200) then
          call xtw84(vlat, vlon, valt, wlat, wlon, walt)
          vlat = wlat
          vlon = wlon
          if (kan.eq.1) valt = walt
        endif
      endif
      call cviken(vlat, vlon, ye, xn)
      call gsurf2(1, xn, ye, valt)
      goto 5
    6 close(1)
      call gsurf2(0, 0., 0., 0.)
      call gsurf3(sr, 9999.9)
      call gsurf4(fnam2)
      call gsurf5()
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,status='old',iostat=ios)
        if(ios.eq.0) then
   90     read(99,*,end=91)
          goto 90
   91     backspace(99)
        else
          open(99,file=flog,status='new')
        endif
        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,'(8x,2a)') 'Mag. anomaly, and Surface elevation ', can        
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xl, yl, mszx, mszy, mx, my, '    var.'
        close(99)
      endif
      stop
c
 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/
     *       3x, '   (mesh)   ', '      xs  ', '      ys  ',
     *       '   msz-x', '   msz-y', '    mx', '    my', '    alt '/
     *       3x, 'surface grid', 2f10.2, 2i8, 2i6, a8)
      end
