!--------( PLMAPCL : Plot color-contour map with Line loc. )--------
  implicit none
  integer :: lwkdir, lrtrim
  real,allocatable :: g(:)
  character(80) :: wdr, flog, fnin, fnln, fnam, buf
  character(72) :: out;  character(20) :: sdtm;  character(8) :: area;
  character(4) :: s4;  character(3) :: spl = '   ';  character(1) :: yn
  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, kl, mval, np
  real :: vnul, xs, ys, gsum, gmin, gmax, gdev, dg, alat, alon, vhi, vlo
  real :: fscl, high, wide, ccm, sclr, sclu, dscl, xo, yo, xt, yt, x, y
  real :: fi0, fi1, fi2, fi3, fi4, fi5, fk0, fk1, fk2, fk3, fk4, fk5
  external cviken

  call premsg('PLMAPCL : plot color-contour map with Line loc.')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnin = wdr;  fnln = wdr;  fnam = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Input data filename ==> ', 81-ldr, fnin(ldr:80))
  open(10,file=fnin,status='old')
  call gparma('Enter StdLIN filename     ==> ', 81-ldr, fnln(ldr:80))
  open(11,file=fnln,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) /= ' ') call abendm('too long filename')

  open(9,status='scratch')
  read(10,'(a)') buf
  do while (buf(1:1) == '#')
    read(10,'(a)') buf
  enddo
  do
    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 premsg(' ')
    write(out,'(a12,a8,6x,a6,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)
    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 premsg('  Specify contour interval (1/2 of color-grading interval)')
    call gparmi('    ( enter 0 to skip contouring ) ==> ', istep)

    if (istep /= 0) then
      istp = iabs(istep)
      call gparmi('  Median value of color-grading ?  ==> ', mval)
      vlo = float(mval - istp*20)
      vhi = float(mval + istp*20)
      if (kpl == 1) then
        call gparmf('  Width of Drawing (in cm) ?  ==> ', wide)
        fscl = float((nmy-1)*mszy) / wide
        high = float((nmx-1)*mszx) / fscl
      else
        call gparmf('  Height of Drawing (in cm) ?  ==> ', high)
        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.
      kscl = 0;  sclr = 0.;  sclu = 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
      write(9,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
      write(9,*) xs, ys, mszx, mszy, nmx, nmy
      write(9,*) istp, mval, fscl, high, wide, np
      write(9,*) kscl, sclr, sclu, klll, kshl
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'PLMAPCL : plot color-contour map with Line location'
      write(6,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin))
      write(6,'(a,a)') ' StdLIN filename     : ', fnln(1:lrtrim(fnln))
      write(6,'(a,a)') ' Output PS filename  : ', fnam(1:lrtrim(fnam))
      write(6,'(a12,a8,6x,a6,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=', istp, 'mval=', mval
      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 ((kscl /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
        write(6,'(a,$)') '      with '
        if (kshl /= 0) write(6,'(a,$)') '  Coast'
        if (mod(kshl/10,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 flush(6)
      call plots(1.5, 1.5)
      call dfrgbt(0)
      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 (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 cviken(float(iorg), float(korg), yo, xo)
        call cvenik(yo, xs, fi0, fk0)
        call cvenik(ys, xs, fi1, fk1)
        call cvenik(yt, xs, fi2, fk2)
        call cvenik(ys, xt, fi3, fk3)
        call cvenik(yt, xt, fi4, fk4)
        call cvenik(yo, xt, fi5, fk5)
        isf = nint(amin1(fi0,fi1,fi2) - 0.5)
        itf = nint(amax1(fi3,fi4,fi5) + 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)
      kl = 0
      do
        if (kl == 0) then
          read(11,'(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 plot((y-ys)/ccm, (x-xs)/ccm, 3)
          kl = 1
        endif
        read(11,'(a)',iostat=IER) buf
        if (IER < 0) exit
        if ((buf(1:1)=='&') .or. (buf(1:1)=='%') .or. (buf(1:1)=='#')) then
          kl = 0; 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, 2)
      enddo
      rewind(11)
      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()
    endif
    deallocate(g)

    read(10,'(a)',iostat=IER) buf
    do while (buf(1:1) == '#')
      if (IER < 0) exit
      read(10,'(a)',iostat=IER) buf
    enddo
    if (IER < 0) exit
    call prompt('  Process Next set ?')
    call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
    if ((yn /= 'y') .and. (yn /= 'Y')) exit
  enddo
  close(10)
  endfile(9)
  rewind(9)
  call pscls()
  call clspin()

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'PLMAPCL : plot color-contour map with Line location'
    write(99,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin))
    write(99,'(a,a)') ' StdLIN filename     : ', fnln(1:lrtrim(fnln))
    write(99,'(a,a)') ' Output PS filename  : ', fnam(1:lrtrim(fnam))
    do
      read(9,'(a8,4x,i4,4i8)',iostat=IER) area, ncc, iorg, korg, ispa, ispb
      if (IER < 0) exit
      read(9,*) xs, ys, mszx, mszy, nmx, nmy
      read(9,*) istp, mval, fscl, high, wide, np
      read(9,*) kscl, sclr, sclu, klll, kshl
      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,i4,4x,a,i4)') 'step=', istp, 'mval=', mval
      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 ((kscl /= 0) .or. (klll /= 0) .or. (kshl /= 0)) then
        write(99,'(a,$)') '      with '
        if (kshl /= 0) write(99,'(a,$)') '  Coast'
        if (mod(kshl/10,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
    enddo
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  close(9)
  stop
  end