!--------( PLAMAG : Plot AMAG map with surrounding mask)--------
  implicit none
  integer :: lwkdir, lrtrim
  external cviken
  real,allocatable :: f(:,:), g(:)
  character(80) :: wdr, flog, fnin, fobs, fnam, buf
  character(72) :: out;  character(20) :: sdtm
  character(8) :: area, area1; character(4) :: s4
  character(3) :: spl = '   '; character(1) :: yn
  integer :: ldr, i, j, mc, nks, nkn, nkw, nke, npx, npy
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, nmx, nmy
  integer :: ixs1, iys1, mszx1, mszy1, nmx1, nmy1
  integer :: ii, ix, mx, iy, my, n, kpl, lcap, kscl, km, iinc, kinc
  integer :: klll, kshl, isf, itf, ksf, ktf, idf, kdf, k
  real :: vnul, xs, ys, vnul1, xs1, ys1, dmy, gsum, gmin, gmax, gdev, dg
  real :: step, astp, wide, high, fscl, ccm, dscl, sclr, sclu, x, y
  real :: xs2, ys2, xt2, yt2, fi1, fi2, fi3, fi4, fk1, fk2, fk3, fk4

  call premsg('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')
  do
    read(10,'(a)') buf
    if (buf(1:1) /= '#') exit
  enddo
    if (buf(13:16) == '    ') then      !-- v2005
    read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb
    if ((ncc<0) .or. (ncc>=400)) call abendm('invalid coord.no.')
    nc = ncc
    if (ncc < 200) then
      nc = ncc;  ncc = nc + 800;  mc = 0
    else
      ncc = ncc - 200;  nc = ncc;  mc = 1
    endif
  else                                !-- v2018
    read(buf,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
    nc = ncc;  mc = 1
    if ((ncc >= 800) .and. (ncc < 1000)) then
      nc = ncc - 800;  mc = 0
    endif
  endif
  if ((nc >= 1) .and. (nc <= 62)) then
    iorg = 0;  ispa = 0;  ispb = 0
    korg = (nc-30)*360 - 180
    if (nc > 60) korg = 0
    if (nc == 61) iorg = +5400
    if (nc == 62) iorg = -5400
  endif
  read(10,*) ixs, iys, mszx, mszy, nmx, nmy, vnul
  xs = float(ixs)/1000.;  ys = float(iys)/1000.

  call gparma('Enter obs.anomaly data filename ==> ', 81-ldr, fobs(ldr:80))
  open(11,file=fobs,status='old')
  do
    read(11,'(a)') buf
    if (buf(1:1) /= '#') exit
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
    if ((ncc1<0) .or. (ncc1>=400)) call abendm('invalid coord.no.')
    nc1 = ncc1
    if (ncc1 < 200) then
      nc1 = ncc1; ncc1 = nc1 + 800
    else
      ncc1 = ncc1 - 200; nc1 = ncc1
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
    nc1 = ncc1
    if ((ncc1 >= 800) .and. (ncc1 < 1000)) nc1 = ncc1 - 800
  endif
  if ((nc1 >= 1) .and. (nc1 <= 62)) then
    iorg1 = 0;  ispa1 = 0;  ispb1 = 0
    korg1 = (nc1-30)*360 - 180
    if (nc1 > 60) korg1 = 0
    if (nc1 == 61) iorg1 = +5400
    if (nc1 == 62) iorg1 = -5400
  endif
  if ((area1 /= area) .or. (ncc1 /= ncc) .or. (iorg1 /= iorg)   &
      .or. (korg1 /= korg) .or. (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
    call abendm('Header inconsistent')
  read(11,*) ixs1, iys1, mszx1, mszy1, nmx1, nmy1, vnul1
  if ((nmx1 <= 0) .or. (nmy1 <= 0)) call abendm('illegal array size')
  xs1 = float(ixs1)/1000.;  ys1 = float(iys1)/1000.

  allocate(f(nmx1, nmy1))

  read(11,*) ((f(i,j),i=1,nmx1),j=1,nmy1)
  close(11)

  call premsg(' ')
  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  write(out,'(a14,2a10,2x,2a8,2a6)')   &
    '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my'
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6)')   &
    '---   AMAG      ', xs, ys, mszx, mszy, nmx, nmy
  call premsg(out)
  write(out,'(a16,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 < 0) .or. (nkn < 0)) call abendm('invalid value')
  npx = nmx - nks - nkn
  if (npx < 2) call abendm('invalid value')
  call gparmi2('    No. of grids on W-bound, E-bound ==> ', nkw,nke)
  if ((nkw < 0) .or. (nke < 0)) call abendm('invalid value')
  npy = nmy - nkw - nke
  if (npy < 2) call abendm('invalid value')

  allocate(g(npx*npy))

  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)
  do ii=1,npx
    i = (ii-1)*npy
    ix = ixs + mszx*(nks+ii-1)
    mx = nint(float(ix-ixs1)/float(mszx1))
    if ((mx <= 0) .or. (mx > 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 <= 0) .or. (my > nmy1)) then
          g(i+j) = vnul
        else
          if (f(mx,my) == vnul1) then
            g(i+j) = vnul
          else
            g(i+j) = g(i+j) / 100.
          endif
        endif
      enddo
    endif
  enddo

  call gparma('Enter Output PS filename  ==> ', 81-ldr, fnam(ldr:80))
  if (fnam(ldr:ldr) == ' ') fnam(1:1) = char(0)
  if (fnam(80:80) /= ' ') call abendm('too long filename')
  do i=1,npx*npy
    dmy = g(i);  if (g(i) /= vnul) exit
  enddo
  if (dmy == vnul) call abendm('no valid data')
  n = 1;  gsum = g(i);  gmin = g(i);  gmax = g(i)
  do j=i+1,npx*npy
    if (g(j) /= vnul) then
      n = n + 1;  gsum = gsum + g(j)
      if (g(j) < gmin) gmin = g(j)
      if (g(j) > gmax) gmax = g(j)
    endif
  enddo
  gsum = gsum / float(n)
  n = 0;  gdev = 0.
  do i=1,npx*npy
    if (g(i) /= 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('')

  call gparmi('  Select  1(portrait) or 2(landscape) ==> ', kpl)
  if (kpl == 1) spl = 'A4P'
  if (kpl == 2) spl = 'A4L'
  if (spl == '   ') call abendm('invalid value')
  call premsg('  Specify contour interval (A/m)')
  call gparmf('   (minus to draw even crowded contour) ==> ', step)
  if (step == 0.) call abendm('illegal contour interval')
  astp = abs(step)
  lcap = 5
  if ((astp == 5.) .or. (astp == 0.5) .or. (astp == 0.05)) lcap = 4
  if ((astp == 2.5) .or. (astp == 0.25) .or. (astp == 0.025)) lcap = 4
  if (step < 0.) lcap = -lcap
  if (kpl == 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.;  sclr = 0.;  sclu = 0.
  kscl = 0;  klll = 0;  kshl = 0
  km = nint(fscl)/25000
  if (km == 0)  km = 1
  if (km > 50) km = 50
  if ((km > 2) .and. (km < 5))   km = 2
  if ((km > 5) .and. (km < 10))  km = 5
  if ((km > 10) .and. (km < 20)) km = 10
  if ((km > 20) .and. (km < 50)) km = 20
  iinc = km;  kinc = km + km/2
  if (km == 50) iinc = 60
  if ((km == 1) .or. (km == 5) .or. (km == 50)) kinc = iinc
  dscl = float(km)/ccm
  call prompt('  Draw  Scale Bar  ?  ')
  call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
  if ((yn == 'y') .or. (yn == '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
  call prompt('  Draw Lat.Lon.Lines ?')
  call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
  if ((yn == 'y') .or. (yn == 'Y')) klll = 1
  call prompt('  Draw Coastlines ?   ')
  call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
  if ((yn == 'y') .or. (yn == 'Y')) then
    kshl = 2
    call prompt('       Rivers  ?      ')
    call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
    if ((yn == 'y') .or. (yn == 'Y')) kshl = 12
    call prompt('       Pref.Boundary ?')
    call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
    if ((yn == 'y') .or. (yn == 'Y')) then
      kshl = kshl + 200
      if (kshl == 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 /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
    write(6,'(a,$)') '      with '
    if (kshl /= 0) write(6,'(a,$)') '  Coast'
    if (mod(kshl,100)/10 /= 0) write(6,'(a,$)') '+River'
    if (kshl/100 /= 0) write(6,'(a,$)') '+PrefB'
    if (klll /= 0) write(6,'(a,$)') '  LatLonLines'
    if (kscl /= 0) then
      write(6,'(a,f5.1,$)') '  ScaleBar [', sclr
      write(6,'(a,f5.1,a,$)') ',', sclu, ' ]'
    endif
    write(6,'(a)')
  endif

  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 /= 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 /= 0) .or. (kshl /= 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 /= 0) then
    do i=isf/iinc,itf/iinc
      idf = i * iinc
      if ((idf < isf) .or. (idf > itf)) cycle
      call cviken(float(idf), float(ksf), y, x)
      call plot((y-ys2)/ccm, (x-xs2)/ccm, 3)
      do kdf=ksf+1,ktf
        call cviken(float(idf), float(kdf), y, x)
        call plot((y-ys2)/ccm, (x-xs2)/ccm, 2)
      enddo
    enddo
    do k=ksf/kinc,ktf/kinc
      kdf = k * kinc
      if ((kdf < ksf) .or. (kdf > ktf)) cycle
      call cviken(float(isf), float(kdf), y, x)
      call plot((y-ys2)/ccm, (x-xs2)/ccm, 3)
      do idf=isf+1,itf
        call cviken(float(idf), float(kdf), y, x)
        call plot((y-ys2)/ccm, (x-xs2)/ccm, 2)
      enddo
    enddo
  endif
  if (kshl /= 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()

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'PLAMAG : 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,'(a12,a8,6x,a6,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 /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
      write(99,'(a,$)') '      with '
      if (kshl /= 0) write(99,'(a,$)') '  Coast'
      if (mod(kshl,100)/10 /= 0) write(99,'(a,$)') '+River'
      if (kshl/100 /= 0) write(99,'(a,$)') '+PrefB'
      if (klll /= 0) write(99,'(a,$)') '  LatLonLines'
      if (kscl /= 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