c--------( GGRID : gen. grid from DPAM line data or equiv.)--------
c-------- This program must be linked with C subroutine <gsurf.o>
c
      dimension v(30)
      character*80 wdr, flog, fnam1, fnam2
      character cu*4, cg*3, area*8, strim*9, buf*150, sdtm*20
c
      call premsg('GGRID : gen. 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 = 9
        kmx = 9
        ku = 2
        cu = 'Deg.'
      else if (ktyp.eq.2) then
        k1 = 2
        k2 = 3
        k3 = 4
        kmx = 4
      else if (ktyp.eq.3) then
        k1 = 3
        k2 = 4
        k3 = 8
        kmx = 8
      else if (ktyp.eq.4) then
        k1 = 1
        k2 = 2
        k3 = 4
        kmx = 4
      else if (ktyp.eq.0) then
      else
        call abendm('invalid 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'
      else if (kg.eq.2) then
        cg = 'TYO'
      else
        call abendm('invalid value')
      endif
      if (ktyp.eq.0) then
   10   read(1,'(a)') 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 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 abendm('too large data position')
        call gparmi('  data position of Mag. data : ', k3)
        if (k3.le.0.or.k3.eq.k1.or.k3.eq.k2)
     *       call abendm('invalid value')
        if (k3.gt.30) call abendm(' too large data position')
        kmx = max(k1, k2, k3)
      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
      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)') 'GGRID : gen. 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,5a,i2,a)')
     *     k1, ': ', cg, '-Lat.(', cu, ')   ',
     *     k2, ': ', cg, '-Lon.(', cu, ')   ',
     *     k3, ': Mag.(nT)'
      write(6,'(a,a)') '  trim radius :', strim
      write(6,'(a,a)')
     *     ' Gen. grid output filename : ', fnam2(1:lrtrim(fnam2))
      write(6,1000) area, ncc, iorg,korg, ispa,ispb,
     *              xl, yl, mszx, mszy, mx, my, '   undef'
      write(2,'(a)') '# output from GGRID'
      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, -1.
      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
        if (ktyp.eq.0) then
   21     lc = index(buf,':')
          if (lc.eq.0) goto 22
          buf(lc:lc) = ' '
          goto 21
        endif
   22   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, talt)
        call cviken(tlat, tlon, ye, xn)
      else if (kg.eq.2 .and. ncc.ge.200) then
        call xtw84(vlat, vlon, 0., wlat, wlon, walt)
        call cviken(wlat, wlon, ye, xn)
      else
        call cviken(vlat, vlon, ye, xn)
      endif
      call gsurf2(1, xn, ye, v(k3))
      goto 1
    2 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,access='append')
        write(99,'(//a)')
     *       'GGRID : gen. 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,5a,i2,a)')
     *       k1, ': ', cg, '-Lat.(', cu, ')   ',
     *       k2, ': ', cg, '-Lon.(', cu, ')   ',
     *       k3, ': Mag.(nT)'
        write(99,'(a,a)') '  trim radius :', strim
        write(99,'(a,a)')
     *       ' Surface grid output filename : ', fnam2(1:lrtrim(fnam2))
        write(99,1000) area, ncc, iorg,korg, ispa,ispb,
     *                 xl, yl, mszx, mszy, mx, my, '   undef'
        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', '    alt '/
     *       3x, 'gen. grid   ', 2f10.2, 2i8, 2i6, a8)
      end
