c--------( PLAMAG : Plot AMAG map with surrounding mask)--------
      parameter (MAXSZ=1201*1201, LXO=1201, LYO=1201)
c
      dimension f(LXO,LYO), g(MAXSZ)
      character*80 wdr, flog, fnin, fobs, fnam
      character area*8, area1*8, s4*4, yn*1, spl*3/'   '/
      character out*72, buf*80, sdtm*20
      external cviken
c
      call prompt('PLAMAG : plot AMAG map with surrounding mask')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnin = wdr
      fobs = wdr
      fnam = wdr
      call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
      call gparma('Enter Input AMAG data filename  ==> ',
     *            81-ldr, fnin(ldr:80))
      open(10,file=fnin,status='old')
    1 read(10,'(a)') buf
      if(buf(1:1).eq.'#') goto 1
      read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb
      if(ncc.ge.200) then
        mc = 1
        nc = ncc - 200
      else
        mc = 0
        nc = ncc
      endif
      if(nc.ge.1.and.nc.le.62) then
        iorg = 0
        ispa = 0
        ispb = 0
        korg = (nc-30)*360 - 180
        if(nc.gt.60) korg = 0
        if(nc.eq.61) iorg = +5400
        if(nc.eq.62) iorg = -5400
      endif
      read(10,*) ixs, iys, mszx, mszy, nmx, nmy, vnul
      xs = float(ixs)/1000.
      ys = float(iys)/1000.
c
      call gparma('Enter obs.anomaly data filename ==> ',
     *            81-ldr, fobs(ldr:80))
      open(11,file=fobs,status='old')
    2 read(11,'(a)') buf
      if(buf(1:1).eq.'#') goto 2
      read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
      if(ncc1.ge.200) then
        nc1 = ncc1 - 200
      else
        nc1 = ncc1
      endif
      if(nc1.ge.1.and.nc1.le.62) then
        iorg1 = 0
        ispa1 = 0
        ispb1 = 0
        korg1 = (nc1-30)*360 - 180
        if(nc1.gt.60) korg1 = 0
        if(nc1.eq.61) iorg1 = +5400
        if(nc1.eq.62) iorg1 = -5400
      endif
      if (area1.ne.area .or. ncc1.ne.ncc .or.
     *    iorg1.ne.iorg. or. korg1.ne.korg .or.
     *    ispa1.ne.ispa .or. ispb1.ne.ispb) then
        call abendm('Header inconsistent')
      endif
      read(11,*) ixs1, iys1, mszx1, mszy1, nmx1, nmy1, vnul1
      if(nmx1.gt.LXO.or.nmy1.gt.LYO) call abendm('too large size')
      if(nmx1.le.0.or.nmy1.le.0) call abendm('illegal array size')
      xs1 = float(ixs1)/1000.
      ys1 = float(iys1)/1000.
      read(11,*) ((f(i,j),i=1,nmx1),j=1,nmy1)
      close(11)
c
      call premsg(' ')
      write(out,'(a,a8,a,i4,4x,4i8)')
     *     '--- Area : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      call premsg(out)
      write(out,'(14x,2a10,2x,2a8,2a6)')
     *     'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my'
      call premsg(out)
      write(out,'(6x,a4,6x,2f10.2,2i8,2i6)')
     *     'AMAG', xs, ys, mszx, mszy, nmx, nmy
      call premsg(out)
      write(out,'(6x,a4,6x,2f10.2,2i8,2i6)')
     *     'obs.', xs1, ys1, mszx1, mszy1, nmx1, nmy1
      call premsg(out)
      call premsg('  Specify widths to limit drawing range')
      call gparmi2('    No. of grids on S-bound, N-bound ==> ',
     *             nks,nkn)
      if (nks.lt.0.or.nkn.lt.0) call abendm('invalid value')
      npx = nmx - nks - nkn
      if (npx.lt.2) call abendm('invalid value')
      call gparmi2('    No. of grids on W-bound, E-bound ==> ',
     *             nkw,nke)
      if (nkw.lt.0.or.nke.lt.0) call abendm('invalid value')
      npy = nmy - nkw - nke
      if (npy.lt.2) call abendm('invalid value')
      if((npx*npy).gt.MAXSZ) then
        call abendm('too many grids')
      endif
      read(10,*)
     *     ((dmy,i=1,nmx),j=1,nkw),
     *     ((dmy,i=1,nks),
     *      (g(i+j),i=0,(npx-1)*npy,npy),
     *      (dmy,i=1,nkn),j=1,npy),
     *     ((dmy,i=1,nmx),j=1,nke)
      close(10)
c
      do ii=1,npx
        i = (ii-1)*npy
        ix = ixs + mszx*(nks+ii-1)
        mx = nint(float(ix-ixs1)/float(mszx1))
        if (mx.le.0 .or. mx.gt.nmx1) then
          do j=1,npy
            g(i+j) = vnul
          enddo
        else
          do j=1,npy
            iy = iys + mszy*(nkw+j-1)
            my = nint(float(iy-iys1)/float(mszy1))
            if (my.le.0 .or. my.gt.nmy1) then
              g(i+j) = vnul
            else
              if (f(mx,my).eq.vnul1) then
                g(i+j) = vnul
              else
                g(i+j) = g(i+j) / 100.
              endif
            endif
          enddo
        endif
      enddo
c
      call gparma('Enter Output PS filename  ==> ',
     *            81-ldr, fnam(ldr:80))
      if (fnam(ldr:ldr).eq.' ') fnam(1:1) = char(0)
      if (fnam(80:80).ne.' ') call abendm('too long filename')
      do i=1,npx*npy
        if(g(i).ne.vnul) goto 3
      enddo
      call abendm('no valid data')
    3 n = 1
      gsum = g(i)
      gmin = g(i)
      gmax = g(i)
      do j=i+1,npx*npy
        if (g(j).ne.vnul) then
          n = n + 1
          gsum = gsum + g(j)
          if(g(j).lt.gmin) gmin = g(j)
          if(g(j).gt.gmax) gmax = g(j)
        endif
      enddo
      gsum = gsum / float(n)
      n = 0
      gdev = 0.
      do i=1,npx*npy
        if (g(i).ne.vnul) then
          dg = g(i) - gsum
          n = n + 1
          gdev = gdev + dg*dg
        endif
      enddo
      gdev = sqrt(gdev / float(n))
      write(out,'(4x,3(a,f7.2),a,f6.2,a)')
     *     'Range (A/m): [', gmin, ' ,', gmax,
     *     ' ]   Mean :', gsum, '   St.Dev.:', gdev
      call premsg(out)
      call premsg('')
c
      call gparmi('  Select  1(portrait) or 2(landscape) ==> ', kpl)
      if(kpl.eq.1) spl = 'A4P'
      if(kpl.eq.2) spl = 'A4L'
      if(spl.eq.'   ') call abendm('invalid value')
      call premsg('  Specify contour interval (A/m)')
      call gparmf('   (minus to draw even crowded contour) ==> ', step)
      if(step.eq.0.) call abendm('illegal contour interval')
      astp = abs(step)
      lcap = 5
      if (astp.eq.5. .or. astp.eq.0.5 .or. astp.eq.0.05) lcap = 4
      if (astp.eq.2.5 .or. astp.eq.0.25 .or. astp.eq.0.025) lcap = 4
      if (step.lt.0.) lcap = -lcap
      if(kpl.eq.1) then
        call gparmf('  Width of Drawing (in cm) ?  ==> ', wide)
        fscl = float((npy-1)*mszy) / wide
        high = float((npx-1)*mszx) / fscl
      else
        call gparmf('  Height of Drawing (in cm) ?  ==> ', high)
        fscl = float((npx-1)*mszx) / high
        wide = float((npy-1)*mszy) / fscl
      endif
      fscl = fscl * 100.
      ccm = fscl/100000.
      kscl = 0
      sclr = 0.
      sclu = 0.
      klll = 0
      kshl = 0
c     if(fscl.lt.3333334..and.fscl.ge.20000.) then
        km = nint(fscl)/25000
        if(km.eq.0)  km = 1
        if(km.gt.50) km = 50
        if(km.gt.2.and.km.lt.5)   km = 2
        if(km.gt.5.and.km.lt.10)  km = 5
        if(km.gt.10.and.km.lt.20) km = 10
        if(km.gt.20.and.km.lt.50) km = 20
        iinc = km
        kinc = km + km/2
        if(km.eq.50) iinc = 60
        if(km.eq.1.or.km.eq.5.or.km.eq.50) kinc = iinc
        dscl = float(km)/ccm
   11   call prompt('  Draw  Scale Bar  ?  ')
        call gparma('  "y" (Yes) or "n" (No) ==> ', 1, yn)
        if(yn.ne.'y'.and.yn.ne.'n') goto 11
        if(yn.eq.'y') then
          kscl = 1
          write(out,'(a,f4.1,a)')
     *          '      size :', dscl, 'cm Wide, 0.8cm High'
          call premsg(out)
          call prompt('    Lower-Left Pos. ? (right and up)')
          call gparmf2(' in cm ==> ', sclr, sclu)
        endif
   12   call prompt('  Draw Lat.Lon.Lines ?')
        call gparma('  "y" (Yes) or "n" (No) ==> ', 1, yn)
        if(yn.ne.'y'.and.yn.ne.'n') goto 12
        if(yn.eq.'y') klll = 1
c     endif
   13 call prompt('  Draw Coastlines ?   ')
      call gparma('  "y" (Yes) or "n" (No) ==> ', 1, yn)
      if(yn.ne.'y'.and.yn.ne.'n') goto 13
      if(yn.eq.'y') then
        kshl = 2
   14   call prompt('       Rivers  ?      ')
        call gparma('  "y" (Yes) or "n" (No) ==> ', 1, yn)
        if(yn.ne.'y'.and.yn.ne.'n') goto 14
        if(yn.eq.'y') kshl = 12
   15   call prompt('       Pref.Boundary ?')
        call gparma('  "y" (Yes) or "n" (No) ==> ', 1, yn)
        if(yn.ne.'y'.and.yn.ne.'n') goto 15
        if(yn.eq.'y') then
          kshl = kshl + 200
          if(kshl.eq.202) kshl = 102
        endif
      endif
      call clspin()
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'PLAMAG : plot AMAG map with surrounding mask'
      write(6,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin))
      write(6,'(a,a)') ' Obs. data filename  : ', fobs(1:lrtrim(fobs))
      write(6,'(a,a)') ' Output PS filename  : ', fnam(1:lrtrim(fnam))
      write(6,'(a,a8,a,i4,4x,4i8)')
     *     'Areaname : ', area, '   Coord.', ncc, iorg,korg, ispa,ispb
      write(6,'(14x,2a10,2x,2a8,2a6)')
     *     'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my'
      write(6,'(16x,2f10.2,2i8,2i6)') xs, ys, mszx, mszy, nmx, nmy
      write(6,'(6x,a,2x,2i4,2x,2i4)')
     *     'range limit (in grids) [S,N, W,E] :', nks,nkn, nkw, nke
      write(6,'(6x,a,f6.2,4x,a,i4)') 'step=', astp, 'lcapt=', lcap
      write(6,'(6x,a,f8.0,4x,a3,f8.1,a1,f8.1,a1)')
     *     'scale  1 /', fscl, spl, high, 'H', wide, 'W'
      if(kscl.ne.0.or.klll.ne.0.or.kshl.ne.0) then
        write(6,'(a,$)') '      with '
        if(kshl.ne.0) write(6,'(a,$)') '  Coast'
        if(mod(kshl,100)/10.ne.0) write(6,'(a,$)') '+River'
        if(kshl/100.ne.0) write(6,'(a,$)') '+PrefB'
        if(klll.ne.0) write(6,'(a,$)') '  LatLonLines'
        if(kscl.ne.0) then
          write(6,'(a,f5.1,$)') '  ScaleBar [', sclr
          write(6,'(a,f5.1,a,$)') ',', sclu, ' ]'
        endif
        write(6,'(a)')
      endif
c
      call psopn(fnam, spl)
      call plots(1.5, 1.5)
      call wrect(0., 0., wide, high)
      call conts(0., 0., wide, high, npy, npx, 0.12)
      call contr(g, vnul, astp, lcap, -99.999, 0., 11)
      call contr(g, vnul, astp, lcap, -astp/2., 99.999, 10)
      if(kscl.ne.0) then
        call plot(sclr,      sclu+0.1, 3)
        call plot(sclr+dscl, sclu+0.1, 2)
        call newpen(2)
        call plot(sclr+dscl, sclu+0.3, 3)
        call plot(sclr+dscl, sclu,     2)
        call plot(sclr,      sclu,     2)
        call plot(sclr,      sclu+0.3, 2)
        write(s4,'(i2,a2)') km, 'km'
        call lstyle('T', 0.5, 0., 0, -255)
        call ptext(s4, 4, sclr+dscl/2., sclu+0.5, 2)
        call newpen(1)
      endif
      call scisor(0., 0., wide, high)
      call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
      if(klll.ne.0.or.kshl.ne.0) then
        xs2 = float(ixs + mszx*nks)/1000.
        ys2 = float(iys + mszy*nkw)/1000.
        xt2 = xs2 + float((npx-1)*mszx)/1000.
        yt2 = ys2 + float((npy-1)*mszy)/1000.
        call cvenik(ys2, xs2, fi1, fk1)
        call cvenik(yt2, xs2, fi2, fk2)
        call cvenik(ys2, xt2, fi3, fk3)
        call cvenik(yt2, xt2, fi4, fk4)
        isf = nint(amin1(fi1,fi2) - 0.5)
        itf = nint(amax1(fi3,fi4) + 0.5)
        ksf = nint(amin1(fk1,fk3) - 0.5)
        ktf = nint(amax1(fk2,fk4) + 0.5)
      endif
      if(klll.ne.0) then
        do 200 i=isf/iinc,itf/iinc
          idf = i * iinc
          if(idf.lt.isf.or.idf.gt.itf) goto 200
          call cviken(float(idf), float(ksf), y, x)
          call plot((y-ys2)/ccm, (x-xs2)/ccm, 3)
          do 100 kdf=ksf+1,ktf
            call cviken(float(idf), float(kdf), y, x)
            call plot((y-ys2)/ccm, (x-xs2)/ccm, 2)
  100     continue
  200   continue
        do 400 k=ksf/kinc,ktf/kinc
          kdf = k * kinc
          if(kdf.lt.ksf.or.kdf.gt.ktf) goto 400
          call cviken(float(isf), float(kdf), y, x)
          call plot((y-ys2)/ccm, (x-xs2)/ccm, 3)
          do 300 idf=isf+1,itf
            call cviken(float(idf), float(kdf), y, x)
            call plot((y-ys2)/ccm, (x-xs2)/ccm, 2)
  300     continue
  400   continue
      endif
      if(kshl.ne.0) then
        call rshore(mc, isf/60, (itf+59)/60, ksf/60, (ktf+59)/60)
        call pshore(-ys2/ccm, -xs2/ccm, kshl, ccm, cviken)
      endif
      call scisor(0., 0., 0., 0.)
      call lstyle('T', 0.4, 0., 0, -255)
      call ptext(area, 8, -1., -1., 0)
      call ptext(fnin, 50, 2., -1., 0)
      call plote()
      call pscls()
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a,$)') 'PLAMAG :'
        write(99,'(a)') ' plot AMAG map with surrounding mask'
        write(99,'(a,a)')
     *       ' Input data filename : ', fnin(1:lrtrim(fnin))
        write(99,'(a,a)')
     *       ' Obs. data filename  : ', fobs(1:lrtrim(fobs))
        write(99,'(a,a)')
     *       ' Output PS filename  : ', fnam(1:lrtrim(fnam))
        write(99,'(a,a8,a,i4,4x,4i8)')
     *       'Areaname : ', area, '   Coord.', ncc, iorg,korg,ispa,ispb
        write(99,'(14x,2a10,2x,2a8,2a6)')
     *       'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my'
        write(99,'(16x,2f10.2,2i8,2i6)') xs, ys, mszx, mszy, nmx, nmy
        write(99,'(6x,a,2x,2i4,2x,2i4)')
     *     'range limit (in grids) [S,N, W,E] :', nks,nkn, nkw, nke
        write(99,'(6x,a,f6.2,4x,a,i4)') 'step=', astp, 'lcapt=', lcap
        write(99,'(6x,a,f8.0,4x,a3,f8.1,a1,f8.1,a1)')
     *       'scale  1 /', fscl, spl, high, 'H', wide, 'W'
        if(kscl.ne.0.or.klll.ne.0.or.kshl.ne.0) then
          write(99,'(a,$)') '      with '
          if(kshl.ne.0) write(99,'(a,$)') '  Coast'
          if(mod(kshl,100)/10.ne.0) write(99,'(a,$)') '+River'
          if(kshl/100.ne.0) write(99,'(a,$)') '+PrefB'
          if(klll.ne.0) write(99,'(a,$)') '  LatLonLines'
          if(kscl.ne.0) then
            write(99,'(a,f5.1,$)') '  ScaleBar [', sclr
            write(99,'(a,f5.1,a,$)') ',', sclu, ' ]'
          endif
          write(99,'(a)')
        endif
        call strdtm(sdtm)
        write(99,'(a,a,a)') '=============== ', sdtm, '==============='
        close(99)
      endif
      stop
      end
