!--------( PLMVARC : Plot color-contour of Mag.Variation )--------
!--------(             with Lines and Crossover points   )--------
  implicit none
  integer :: lwkdir, lrtrim
  real,allocatable :: g(:)
  character(80) :: wdr, flog, fnin, fnln, fnpt, fnam, buf
  character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
  character(6) :: s6;  character(4) :: s4
  character(3) :: spl = '   ';  character(1) :: yn
  integer :: kclo = Z'005576', kchi = Z'67006a', kcnv = Z'e8e8e8'
  integer :: mrgb(30)
  integer :: ldr, mc, i, j, n, k, IER
  integer :: istep, istp, kpl, kscl, klll, kshl, km, iinc, kinc
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: ixs, iys, mszx, mszy, nmx, nmy
  integer :: isf, itf, ksf, ktf, idf, kdf, mval, np
  integer :: kleg, kpen, nclev, nspc
  real :: vnul, xs, ys, gsum, gmin, gmax, gdev, dg, vlo, vhi
  real :: fscl, high, wide, ccm, sclr, sclu, dscl, xt, yt, x, y
  real :: fi1, fi2, fi3, fi4, fk1, fk2, fk3, fk4
  real :: alat, alon, dleg, slgr, slgu
  data (mrgb(i),i=1,20) /   &
    Z'00676a', Z'007b5e', Z'009351', Z'00b33f', Z'00e422',   &
    Z'1aff07', Z'46fe00', Z'74ff01', Z'a3ff0d', Z'd0ff16',   &
    Z'ffff00', Z'ffd016', Z'ffa30d', Z'ff7401', Z'fe4600',   &
    Z'ff1a07', Z'e40022', Z'b3003f', Z'930051', Z'7b005e' /
  real :: h(2,2)
  external cviken

  call premsg('PLMVARC : Plot color-contour of Mag.Variation')
  call premsg('            with Lines and Crossover points')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnin = wdr;  fnln = wdr;  fnpt = wdr;  fnam = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Input Mag.Var. filename    ==> ', 81-ldr, fnin(ldr:80))
  open(10,file=fnin,status='old')
  call gparma('Enter LinesLoc. StdLIN filename  ==> ', 81-ldr, fnln(ldr:80))
  open(11,file=fnln,status='old')
  call gparma('Enter ControlPoint data filename ==> ', 81-ldr, fnpt(ldr:80))
  open(12,file=fnpt,status='old')
  call gparma('Enter Output PS filename         ==> ', 81-ldr, fnam(ldr:80))
  if (fnam(ldr:ldr) == ' ') fnam(1:1) = char(0)
  if (fnam(80:80) /= ' ') then
    call abendm('too long filename')
  endif

  read(10,'(a)') buf
  do while (buf(1:1) == '#')
    read(10,'(a)') buf
  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.')
    if (ncc < 200) then
      mc = 0;  nc = ncc;  ncc = ncc + 800
    else
      mc = 1;  ncc = ncc - 200;  nc = ncc
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
    nc = ncc
    if ((ncc >= 800) .and. (ncc < 1000)) nc = ncc - 800
  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 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,'(16x,2f10.2,2i8,2i6)') xs, ys, mszx, mszy, nmx, nmy
  call premsg(out)

  allocate(g(nmx*nmy))
  read(10,*) ((g(i+j),i=0,(nmx-1)*nmy,nmy),j=1,nmy)
  close(10)
  n = 0
  do i=1,nmx*nmy
    if (g(i) == vnul) cycle
    gsum = g(i);  gmin = g(i);  gmax = g(i)
    n = 1;  exit
  enddo
  if (n == 0) call abendm('no valid data')
  do j=i+1,nmx*nmy
    if (g(j) == vnul) cycle
    n = n + 1;  gsum = gsum + g(j)
    if (g(j) < gmin) gmin = g(j)
    if (g(j) > gmax) gmax = g(j)
  enddo
  gsum = gsum / float(n)
  n = 0;  gdev = 0.
  do i=1,nmx*nmy
    if (g(i) == vnul) cycle
    dg = g(i) - gsum
    n = n + 1;  gdev = gdev + dg*dg
  enddo
  gdev = sqrt(gdev / float(n))
  write(out,'(4x,3(a,f7.1),a,f6.1)') 'Data Range : [', gmin,   &
     ' ,', gmax, ' ]   Mean :', gsum, '   St.Dev.:', gdev
  call premsg(out)
  call premsg('')

  if (spl == '   ') then
    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 psopn(fnam, spl)
  endif
  call gparmi('  Specify contour (=color-grading) interval ==> ', istep)
  if (istep == 0) call abendm('invalid contour interval')
  istp = iabs(istep)
  call gparmi('  Median value of color-grading ?  ==> ', mval)
  call gparmi('  Blank zone Half-width ? (in Cont.Intv.) ==> ', nspc)
  if ((nspc < 0) .or. (nspc > 5)) call abendm('invalid Blank zone')
  nclev = 20 + nspc+nspc
  if (nspc /= 0) then
    do i=0,9
      mrgb(20-i+nspc+nspc) = mrgb(20-i)
    enddo
    do i=1,(nspc+nspc)
      mrgb(10+i) = Z'ffffff'
    enddo
  endif
  vlo = float(mval - istep*(10+nspc))
  vhi = float(mval + istep*(10+nspc))
  if (kpl == 1) then
    call gparmf('  Width of Drawing (in cm) ?  ==> ', wide)
    if (wide <= 0.) call abendm('illegal size')
    fscl = float((nmy-1)*mszy) / wide
    high = float((nmx-1)*mszx) / fscl
  else
    call gparmf('  Height of Drawing (in cm) ?  ==> ', high)
    if (high <= 0.) call abendm('illegal size')
    fscl = float((nmx-1)*mszx) / high
    wide = float((nmy-1)*mszy) / fscl
  endif
  call premsg('  Select Pen-No. (1-8) to draw LineLoc')
  call premsg('    1-4:solid (0.1/0.3/0.6/1mm thick) or')
  call gparmi('    5-8:broken(0.1/0.3/0.6/1mm thick) ==> ', np)
  if ((np < 1) .or. (np > 8)) call abendm('invalid pen number')
  fscl = fscl * 100.;  ccm = fscl/100000.
  kleg = 0;  kscl = 0;  sclr = 0.;  sclu = 0.;  klll = 0;  kshl = 0
  if ((fscl < 3333334.) .and. (fscl >= 20000.)) then
    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  ColorGrading Legend  ?  ')
    call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
    if ((yn == 'y') .or. (yn == 'Y')) then
      kleg = 1;  dleg = float(11+nspc) * 0.8
      write(out,'(a,f4.1,a)') '      size :  2cm Wide,', dleg, 'cm High'
      call premsg(out)
      call prompt('    Lower-Left Pos. ? (right and up)')
      call gparmf2(' in cm ==> ', slgr, slgu)
    endif
    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
  endif
  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)')   'PLMVARC : Plot color-contour of Mag.Variation'
  write(6,'(a)')   '            with Lines and Control-points'
  write(6,'(a,a)') ' Input Mag.Var. filename    : ', fnin(1:lrtrim(fnin))
  write(6,'(a,a)') ' LinesLoc. StdLIN filename  : ', fnln(1:lrtrim(fnln))
  write(6,'(a,a)') ' ControlPoint data filename : ', fnpt(1:lrtrim(fnpt))
  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,i4,4x,a,i4)') 'step=CGintv=', istep, 'mval=', mval
  write(6,'(6x,a,i2,a)') 'BlankZone Half-width =', nspc, ' (x Cont.Intv.)'
  write(6,'(6x,a,f8.0,4x,a3,f8.1,a1,f8.1,a7,i2)')   &
     'scale  1 /', fscl, spl, high, 'H', wide, 'W   Np=', np
  if ((kleg /= 0) .or. (kscl /= 0)) then
    write(6,'(a,$)') '      with '
    if (kleg /= 0) then
      write(6,'(a,f5.1,$)') '  ColorLegend [', slgr
      write(6,'(a,f5.1,a,$)') ',', slgu, ' ]'
    endif
    if (kscl /= 0) then
      write(6,'(a,f5.1,$)') '  ScaleBar [', sclr
      write(6,'(a,f5.1,a,$)') ',', sclu, ' ]'
    endif
    write(6,'(a)')
  endif
  if ((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'
    write(6,'(a)')
  endif

  call plots(1.5, 1.5)
  call dfrgbt(nclev, mrgb, kclo, kchi, kcnv)
  call dframe(0., 0., wide, high, nmy, nmx)
  call paintm(g, vlo, vhi, vnul)
  call wrect(0., 0., wide, high)
  call conts(0., 0., wide, high, nmy, nmx, 0.)
  call contx(g, vnul, istp, 2, -99999, 99999, 0)
  if (kleg /= 0) then
    h(1,1) = float(mval - istp*(11+nspc))
    h(2,1) = h(1,1)
    h(1,2) = float(mval + istp*(11+nspc))
    h(2,2) = h(1,2)
    call dframe(slgr+1., slgu, 1., dleg, 2, 2)
    call paintm(h, vlo, vhi, vnul)
    call wrect(slgr+1., slgu, 1., dleg)
    call conts(slgr+1., slgu, 1., dleg, 2, 2, 0.)
    call contx(h, vnul, istp, 1, -99999, 99999, 0)
    do i=-10-nspc,+10+nspc
      write(s6,'(f6.0)') float(i*istp+mval)
      y = float(i+11+nspc)*0.4 + slgu - 0.02
      call pcstr(slgr+1.10, y, 0.2, s6, 0., -6)
    enddo
  endif
  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
    xt = xs + float((nmx-1)*mszx)/1000.
    yt = ys + float((nmy-1)*mszy)/1000.
    call cvenik(ys, xs, fi1, fk1)
    call cvenik(yt, xs, fi2, fk2)
    call cvenik(ys, xt, fi3, fk3)
    call cvenik(yt, xt, 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-ys)/ccm, (x-xs)/ccm, 3)
      do kdf=ksf+1,ktf
        call cviken(float(idf), float(kdf), y, x)
        call plot((y-ys)/ccm, (x-xs)/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-ys)/ccm, (x-xs)/ccm, 3)
      do idf=isf+1,itf
        call cviken(float(idf), float(kdf), y, x)
        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
      enddo
    enddo
  endif
  if (kshl /= 0) then
    call rshore(mc, isf/60, (itf+59)/60, ksf/60, (ktf+59)/60)
    call pshore(-ys/ccm, -xs/ccm, kshl, ccm, cviken)
  endif

  call newpen(np);  kpen = 3
  do
    read(11,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%') .or. (buf(1:1) == '#')) then
      kpen = 3; cycle
    endif
    i = index(buf,'N');  if (i /= 0) buf(i:i) = ' '
    i = index(buf,'E');  if (i /= 0) buf(i:i) = ' '
    read(buf,*) alat, alon
    call cviken(alat, alon, y, x)
    call plot((y-ys)/ccm, (x-xs)/ccm, kpen)
    kpen = 2
  enddo
  close(11)

  call newpen(2)
  do
    read(12,'(a)',iostat=IER) buf
    if (IER < 0) exit
    if ((buf(1:1) == '&') .or. (buf(1:1) == '%') .or. (buf(1:1) == '#')) cycle
    i = index(buf,'N');  if (i /= 0) buf(i:i) = ' '
    i = index(buf,'E');  if (i /= 0) buf(i:i) = ' '
    read(buf,*) alat, alon
    call cviken(alat, alon, y, x)
    call wcirc((y-ys)/ccm, (x-xs)/ccm, 0.05, 0., 360.)
  enddo
  close(12)
  call newpen(1)

  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)') 'PLMVARC : Plot color-contour of Mag.Variation'
    write(99,'(a)')   '            with Lines and Control-points'
    write(99,'(a,a)') ' Input Mag.Var. filename    : ', fnin(1:lrtrim(fnin))
    write(99,'(a,a)') ' LinesLoc. StdLIN filename  : ', fnln(1:lrtrim(fnln))
    write(99,'(a,a)') ' ControlPoint data filename : ', fnpt(1:lrtrim(fnpt))
    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,i4,4x,a,i4)')   &
      'step=CGintv=', istep, 'mval=', mval
    write(99,'(6x,a,i2,a)')   &
      'BlankZone Half-width =', nspc, ' (x Cont.Intv.)'
    write(99,'(6x,a,f8.0,4x,a3,f8.1,a1,f8.1,a7,i2)')   &
      'scale  1 /', fscl, spl, high, 'H', wide, 'W   Np=', np
    if ((kleg /= 0) .or. (kscl /= 0)) then
      write(99,'(a,$)') '      with '
      if (kleg /= 0) then
        write(99,'(a,f5.1,$)') '  ColorLegend [', slgr
        write(99,'(a,f5.1,a,$)') ',', slgu, ' ]'
      endif
      if (kscl /= 0) then
        write(99,'(a,f5.1,$)') '  ScaleBar [', sclr
        write(99,'(a,f5.1,a,$)') ',', sclu, ' ]'
      endif
      write(99,'(a)')
    endif
    if ((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'
      write(99,'(a)')
    endif
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end