!--------( XPLMAPCS : Extended plotting of color-shading map )--------
!-------- This program must be linked with C-subroutine "xplgobj.o"
  implicit none
  integer :: lwkdir, lrtrim
  integer,parameter :: LBUF = 256, MXPT = 150
  real,parameter :: RAD = (180./3.14159)
  real,allocatable :: g(:), s(:)
  real :: xp(MXPT), yp(MXPT)
  real :: robj(8);  integer :: kobj(6)
  character(LBUF) :: sbf
  character(80) :: wdr, flog, fnin, fnobj, fncap, fnam, buf
  character(72) :: out;  character(40) :: slbl;  character(20) :: sdtm
  character(10) :: cunit;  character(8) :: area;  character(6) :: s6
  character(4) :: sht;  character(3) :: fnt;  character(1) :: yn
  integer :: ldr, nseq, ns, n, i, j, k, kpolyg, m, nct, IER
  integer :: nc, ncc, mc, iorg, korg, ispa, ispb
  integer :: ixs, iys, mszx, mszy, nmx, nmy
  integer :: istep, istp, km, iinc, kinc, kscl, klll, kshl, klbl
  integer :: ibcol, itcol, idf, isf, itf, kdf, ksf, ktf, ksp, mark
  integer :: kcgl, mval, nras
  real :: vnul, dmy, xs, ys, gsum, gmin, gmax, szx, szy, swd, shi
  real :: fscl, am, bm, high, wide, rs, ts, ccm, dscl, sclr, sclu
  real :: ca, fi0, fi1, fi2, fi3, fi4, fi5, fk0, fk1, fk2, fk3, fk4, fk5
  real :: wd, hi, angle, sa, size, x, xo, xt, y, yo, yt
  real :: cglf, cglr, cglu, vhi, vlo
  real :: az, el, cx, cy, cz, dx, dy, fx, fy, vs, sn, p, q, r
  external cviken

  call premsg('XPLMAPCS : Extended plotting of color-shading map')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnin = wdr;  fnobj = wdr;  fncap = 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 gparmi(' dataset sequence number  ==> ', nseq)
  if (nseq <= 0) call abendm('invalid value')
  ns = 0
  do
    read(10,'(a)') buf
    if (buf(1:1) == '#') cycle
    read(10,*) ixs, iys, mszx, mszy, nmx, nmy, vnul
    ns = ns + 1
    if (nseq == ns) exit
    read(10,*) ((dmy,i=1,nmx),j=1,nmy)
  enddo
  call gparma('Specify Overlay Obj. file ==> ', 81-ldr, fnobj(ldr:80))
  if (fnobj(ldr:ldr) /= ' ') open(11,file=fnobj,status='old')
  call gparma('Specify Caption data file ==> ', 81-ldr, fncap(ldr:80))
  if (fncap(ldr:ldr) /= ' ') open(12,file=fncap,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')

  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
  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));  allocate(s(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)
  write(out,'(4x,a,f7.1,a,f7.1,a,f7.1)')   &
    'Data Range : [', gmin, ' ,', gmax, ' ]    Average :', gsum
  call premsg(out)

  call gparma(' Sheet Size ?  (A4P, B4L, A1L, B0P, etc.) : ', 4, sht)
  if (sht(1:1) == 'a') sht(1:1) = 'A'
  if (sht(1:1) == 'b') sht(1:1) = 'B'
  if (sht(1:2) == 'A4') then
    szx = 19.7;  szy = 27.9
  else if (sht(1:2) == 'B4') then
    szx = 24.1;  szy = 34.3
  else if (sht(1:2) == 'A3') then
    szx = 27.9;  szy = 39.4
  else if (sht(1:2) == 'B3') then
    szx = 34.3;  szy = 48.3
  else if (sht(1:2) == 'A2') then
    szx = 39.4;  szy = 55.9
  else if (sht(1:2) == 'B2') then
    szx = 48.3;  szy = 68.6
  else if (sht(1:2) == 'A1') then
    szx = 55.9;  szy = 78.7
  else if (sht(1:2) == 'B1') then
    szx = 68.6;  szy = 96.5
  else if (sht(1:2) == 'A0') then
    szx = 78.7;  szy =111.8
  else if (sht(1:2) == 'B0') then
    szx = 96.5;  szy =137.2
  else
    call abendm('unknown sheet size')
  endif
  if (sht(3:3) == 'p') sht(3:3) = 'P'
  if (sht(3:3) == 'l') sht(3:3) = 'L'
  if (sht(3:3) == 'P') then
    swd = szx;  shi = szy
  else if (sht(3:3) == 'L') then
    swd = szy;  shi = szx
  else
    call abendm('unknown sheet orientation')
  endif
  if (sht(4:4) == 'e') sht(4:4) = 'E'
  if (sht(4:4) /= 'E') sht(4:4) = ' '
  call gparmf(' Scale ?  (in reciprocal) : ', fscl)
  call gparmf(' Left Margin ?  (in cm)   : ', am)
  call gparmf(' Bottom Margin ?  (in cm) : ', bm)
  call gparmi(' Color grading interval (1/20 of full-range) : ', istep)
  call gparmi(' Median value ?  (for color-grading) : ', mval)
  istp = iabs(istep)
  vlo = float(mval - istp*10)
  vhi = float(mval + istp*10)
  call premsg(' Raster resolution of Color-shading')
  call gparmi('  <default=2> (in 600dpi) [1/2/4/8-] : ', nras)
  if (nras < 1) nras = 2
  call premsg('Specify illuminant direction (in Degrees)')
  call gparmf('    Azimuth angle   : ', az)
  call gparmf('    Elevation angle : ', el)
  call premsg('Assign vertical magnification')
  call gparmf('    Value to be scaled to 1km : ', vs)
  fx = vs * float(mszx)/1000.
  fy = vs * float(mszy)/1000.
  cx = cos(el/rad)*cos(az/rad)
  cy = cos(el/rad)*sin(az/rad)
  cz = sin(el/rad)

  high = float((nmx-1)*mszx)*100. / fscl
  wide = float((nmy-1)*mszy)*100. / fscl
  rs = swd - wide - am
  ts = shi - high - bm
  write(out,'(a,i4,4x,a,f8.0,4x,a4)')   &
    '---   step=', istep, 'scale  1 /', fscl, sht
  call premsg(out)
  write(out,'(8x,a,f5.1,a,f5.1)')   &
    'Drawing   Width:', wide, '   Height:', high
  call premsg(out)
  write(out,'(8x,a,f5.1,a,f5.1)')   &
    'Margin    Left :', am, '   Bottom:', bm
  call premsg(out)
  write(out,'(8x,a,f5.1,a,f5.1)')   &
    'Rspace    Right:', rs, '   Top   :', ts
  call premsg(out)
  write(out,'(8x,a,i2)') 'Raster resolution :  (600dpi) *', nras
  call premsg(out)

  ccm = fscl/100000.
  km = nint(fscl)/25000
  if (km < 1) km = 1
  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
  if ((km >= 50) .and. (km < 120)) km = 60
  if ((km > 120) .and. (km < 300)) km = 120
  if ((km > 300) .and. (km < 600)) km = 300
  if (km > 600) km = 600
  iinc = km
  kinc = km + km/2
  if (km > 50) km = 50
  if ((km == 5) .or. (km == 50)) kinc = iinc
  dscl = float(km)/ccm

  kcgl = 0;  cglf = 1.;  cglr = 0.;  cglu = 0.
  call prompt('  Draw Color-Shading Legend ?')
  call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
  if ((yn == 'y') .or. (yn == 'Y')) then
    kcgl = 1
    call premsg('      st.size :  11.5cm Wide, 22.5cm High')
    call gparmf('    Size factor ?  (ratio to st.size) : ', cglf)
    call prompt('    Lower-Left Pos. ? (right and up)')
    call gparmf2(' in cm : ', cglr, cglu)
    call gparma('    Unit notation ?  (max. 10 chars) : ', 10, cunit)
  endif
  kscl = 0;  sclr = 0.;  sclu = 0.;  klll = 0;  kshl = 0
  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 gparma('Enter "y" to continue O.K. [otherwise Break] ==> ', 1, yn)
  if ((yn /= 'y') .and. (yn /= 'Y')) call abendm('Process terminated')
  call prompt('  Write Identifier-Label ?')
  call gparma('  "y" for Yes, othewise No ==> ', 1, yn)
  if ((yn == 'y') .or. (yn == 'Y')) then
    klbl = 1
    call gparma('    Label ?  (max. 40 chars) : ', 40, slbl)
  else
    klbl = 0
  endif
  call clspin()

  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'XPLMAPCS: plot color-shading map for Grid data'
  write(6,'(a,a)') ' Input data filename   : ', fnin(1:lrtrim(fnin))
  write(6,'(a,i2)') '   dataset seq. number : ', nseq
  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,'(a,a)') ' Output PS filename  : ', fnam(1:lrtrim(fnam))
  write(6,'(6x,a,i4,4x,a,i4,4x,a,f8.0,4x,a4)')   &
    'step=', istep, 'mval=', mval, 'scale  1 /', fscl, sht
  write(6,'(6x,a,f7.1,a,f7.1)')   &
    'illuminant :   Azimuth', az, ',   Elevation', el
  write(6,'(8x,a,f8.2)') 'value to be scaled to 1km :', vs
  write(6,'(6x,a,i2)') 'Raster resolution :  DPI = 600 /', nras
  write(6,'(6x,a,f5.1,a,f5.1)') 'Margin      Left :', am, '   Bottom:', bm
  write(6,'(8x,a,f5.1,a,f5.1)') 'Drawing   Width:', wide, '   Height:', high
  write(6,'(8x,a,f5.1,a,f5.1)') 'Rspace    Right:', rs, '   Top   :', ts
  if (kcgl /= 0) then
    write(6,'(a,$)') '      with Color-Shading Legend  '
    write(6,'(a,f6.3,a)') '[ Size factor =', cglf, ' ]'
    write(6,'(a,f5.1,$)') '           at [ ', cglr
    write(6,'(a,f5.1,a,$)') ', ', cglu, ' ]'
    write(6,'(5x,a,a10,a)') '[ Unit notation = ', cunit, ']'
  endif
  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) write(6,'(a,f5.1,a,f5.1,a,$)')   &
                     '  ScaleBar [', sclr, ',', sclu, ' ]'
    write(6,'(a)')
  endif
  if (klbl /= 0) write(6,'(a,a)') '      Identifier : ', slbl
  if (fnobj(ldr:ldr) /= ' ')   &
    write(6,'(a,a)') ' Overlay obj. filename : ', fnobj(1:lrtrim(fnobj))
  if (fncap(ldr:ldr) /= ' ')   &
    write(6,'(a,a)') ' Caption data filename : ', fncap(1:lrtrim(fncap))

  do i=1,nmx;  do j=1,nmy
    n = (i-1)*nmy + j
    if (g(n) == vnul) then
      s(n) = 1.
    else
      if (i == 1.or.g(n-nmy) == vnul) then
        if (i == nmx.or.g(n+nmy) == vnul) then
          dx = 0.
        else
          dx = g(n+nmy) - g(n)
        endif
      else
        if (i == nmx.or.g(n+nmy) == vnul) then
          dx = g(n) - g(n-nmy)
        else
          dx = (g(n+nmy)-g(n-nmy)) / 2.
        endif
      endif
      if (j == 1.or.g(n-1) == vnul) then
        if (j == nmy.or.g(n+1) == vnul) then
          dy = 0.
        else
          dy = g(n+1) - g(n)
        endif
      else
        if (j == nmy.or.g(n+1) == vnul) then
          dy = g(n) - g(n-1)
        else
          dy = (g(n+1)-g(n-1)) / 2.
        endif
      endif
      p = dx/fx;  q = dy/fy;  r = sqrt(p*p + q*q + 1.)
      sn = (p*cx + q*cy - cz) / r
!!!!      s(n) = (1.-sn) / 2.
      if (sn < 0.) then
        s(n) = 0.75 - sn/4.
      else
        s(n) = (1.-sn) * 0.75
      endif
!!!!
    endif
  enddo;  enddo

  call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
  call psopn(fnam, sht)
  call plots(am, bm)
  call premsg('--- Drawing Color Contour')
  call dfcols(20, 0.85, 1., -1)
  call dresol(nras)
  call dframe(0., 0., wide, high, nmy, nmx)
  call paintw(g, vlo, vhi, vnul, s)
  call wrect(0., 0., wide, high)
  if (klbl /= 0) then
    call lstyle('CFN', 0.3, 0., 0, -999)
    call ptext(slbl, 40, -am, -bm, 0)
  endif

  if (fnobj(ldr:ldr) /= ' ') then
    call premsg('--- Drawing Overlay Obj.')
    ns = 0;  kpolyg = 0
    do
      ns = ns + 1
      read(11,'(a)',iostat=IER) sbf(1:LBUF-1)
      if (IER < 0) exit
      if (sbf(1:1) == '#') cycle
      j = LBUF
      do while (j > 0)
        j = j-1
        if (sbf(j:j) /= ' ') exit
      enddo
      if (j == 0) cycle
      sbf(j+1:j+1) = char(0)

      if (kpolyg == 0) then
        call getobj(sbf, kobj, robj)
        k = kobj(1);  nct = kobj(2)
        if (((nct /= -1) .and. (nct /= ncc)) .or. (k <= 0) .or. (k > 12)) then
          write(buf,'(a17,i5,a8)') '  (Obj.) Error at', ns, '-th Line'
          call premsg(buf(1:30)); cycle
        endif
        if ((k <= 7) .or. (k == 9)) then
          if (nct == -1) then
            call cviken(robj(1), robj(2), y, x)
            robj(1) = (x-xs)/ccm;  robj(2) = (y-ys)/ccm
            if ((k >= 3) .and. (k <= 6)) then
              call cviken(robj(3), robj(4), y, x)
              robj(3) = (x-xs)/ccm;  robj(4) = (y-ys)/ccm
            endif
          else
            robj(1) = (robj(1) - xs) / ccm
            robj(2) = (robj(2) - ys) / ccm
            if ((k >= 3) .and. (k <= 6)) then
              robj(3) = (robj(3) - xs) / ccm
              robj(4) = (robj(4) - ys) / ccm
            endif
          endif
        else if ((k == 10) .or. (k == 11)) then
          if ((kobj(3) < 3) .or. (kobj(3) > MXPT))   &
            call abendm('too few/many corners (polygon)')
          kpolyg = 1;  i = 1;  m = 1;  cycle
        else if ((k == 12) .and. (nct /= -1)) then
          call cvenik(robj(2), robj(1), robj(1), robj(2))
          call cvenik(robj(4), robj(3), robj(3), robj(4))
        endif
      else
        j = 1
        do
          if (m == 1) then
            if (nct == -1) then
                call scandm(j, sbf, xp(i))
            else
                call scanfv(j, sbf, xp(i))
            endif
          else if (m == 2) then
            if (nct == -1) then
                call scandm(j, sbf, yp(i))
            else
                call scanfv(j, sbf, yp(i))
            endif
          else
            if (k == 11) call scanfv(j, sbf, dmy)
          endif
          if (j == 0) exit
          if (j < 0) then
            write(buf,'(a17,i5,a8)') '  (Obj.) Error at', ns, '-th Line'
            call premsg(buf(1:30));  k = 0;  exit
          endif
          m = m + 1
          if (m <= 3) cycle
          m = 1;  i = i + 1
          if (i > kobj(3)) then
            kpolyg = 0; exit
          endif
        enddo
        if (kpolyg == 1) cycle
        if (nct == -1) then
          do i=1,kobj(3)
            call cviken(xp(i), yp(i), y, x)
            xp(i) = (x-xs)/ccm
            yp(i) = (y-ys)/ccm
          enddo
        else
          do i=1,kobj(3)
            xp(i) = (xp(i) - xs) / ccm
            yp(i) = (yp(i) - ys) / ccm
          enddo
        endif
      endif

      if (k == 1) then
        i = kobj(5)
        if ((i < -255) .or. (i > 4095)) i = 0
        call dfpcol(1, i)
        call paintc(robj(2), robj(1), robj(7)/2.)
      else if (k == 2) then
        call penatr(0, kobj(5), kobj(6), robj(8))
        call wcirc(robj(2), robj(1), robj(6)/ccm, 0., 360.)
      else if ((k == 3) .or. (k == 4)) then
        call penatr(0, kobj(5), kobj(6), robj(8))
        call plot(robj(2), robj(1), 3)
        call plot(robj(4), robj(3), 2)
      else if ((k == 5) .or. (k == 6)) then
        wd = robj(4)-robj(2)
        hi = robj(3)-robj(1)
        call penatr(0, kobj(5), kobj(6), robj(8))
        call wrect(robj(2), robj(1), wd, hi)
      else if (k == 7) then
        mark = kobj(3)
        call pmark(mark, robj(2), robj(1), robj(7), robj(8), kobj(5))
      else if (k == 8) then
        fnt = sbf(kobj(3)+1:kobj(4))
        call lstyle(fnt, robj(7), robj(6), kobj(5), kobj(6))
      else if (k == 9) then
        i = kobj(3)+1;  j = kobj(4)
        if (((sbf(i:i)=='''') .and. (sbf(j:j)=='''')) .or.   &
            ((sbf(i:i)=='"') .and. (sbf(j:j)=='"'))) then
          i = i + 1; j = j - 1
        endif
        call ptext(sbf(i:j), j-i+1, robj(2), robj(1), kobj(6))
      else if ((k == 10) .or. (k == 11)) then
        call penatr(0, kobj(5), kobj(6), robj(8))
        call wpolyg(yp, xp, kobj(3), 1)
      else if (k == 12) then
        call penatr(0, kobj(5), kobj(6), robj(8))
        call cviken(robj(1), robj(2), y, x)
        call plot((y-ys)/ccm, (x-xs)/ccm, 3)
        call cviken(robj(3), robj(2), y, x)
        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
        call cviken(robj(3), robj(4), y, x)
        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
        call cviken(robj(1), robj(4), y, x)
        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
        call cviken(robj(1), robj(2), y, x)
        call plot((y-ys)/ccm, (x-xs)/ccm, 2)
      endif
    enddo
    close(11)
    call newpen(1)
  endif

  if (kscl /= 0) then
    call premsg('--- Drawing Scale-Bar')
    write(s6,'(i2,a2)') km, 'km'
    call lstyle('TB ', 0.7, 0., 0, -999)
    call ptext(s6, 4, sclr+dscl/2., sclu+0.5, 2)
    call newpen(2)
    call plot(sclr,      sclu+0.1, 3)
    call plot(sclr+dscl, sclu+0.1, 2)
    call newpen(3)
    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)
    call newpen(1)
  endif
  call scisor(0., 0., wide, high)
  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
    call premsg('--- Drawing LatLonLines')
    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 premsg('--- Drawing ShoreLines etc.')
    call rshore(mc, isf/60, (itf+59)/60, ksf/60, (ktf+59)/60)
    call pshore(-ys/ccm, -xs/ccm, kshl, ccm, cviken)
  endif
  call scisor(0., 0., 0., 0.)
  if (kcgl /= 0) then
    call premsg('--- Drawing Color-Shading Legend')
    g(1) = vlo - float(istp)
    g(2) = vhi + float(istp)
    g(3) = g(1)
    g(4) = g(2)
    s(1) = 1.;  s(2) = 1.;  s(3) = 0.;  s(4) = 0.
    x = cglr + cglf*1.5
    y = cglu
    xt = cglf*10.
    yt = cglf*22.
    call dframe(x, y, xt, yt, -2, -2)
    call paintw(g, vlo, vhi, vnul, s)
    call newpen(2)
    call wrect(x, y, xt, yt)
    call lstyle('TBR', cglf*0.6, 0., 0, -999)
    call ptext(cunit, 10, cglr, y+yt+0.08, 0)
    call lstyle('CFN', cglf*0.48, 0., 0, -999)
    do i=-10,10
      y = y + cglf
      write(s6,'(i6)') mval+istp*2*i
      call ptext(s6, 6, x-cglf*0.06, y-cglf*0.12, 1)
    enddo
  endif

  if (fncap(ldr:ldr) /= ' ') then
    call premsg('--- Drawing Captions')
    xs = 0.;  ys = 0.;  ca = 1.;  sa = 0.;  ksp = 0;  ns = 0
    do
      ns = ns + 1
      read(12,'(a)',iostat=IER) sbf(1:LBUF-1)
      if (IER < 0) exit
      if (sbf(1:1) == '#') cycle
      j = LBUF-1
      do while (sbf(j:j) == ' ')
        if (j == 1) exit
        j = j - 1
      enddo
      sbf(j+1:j+1) = char(0)
      if (sbf(1:1) == ' ') then
        i = 2
        do while (sbf(j:j) == ' ')
          i = i + 1
        enddo
        if (((sbf(i:i)=='''') .and. (sbf(j:j)=='''')) .or.   &
            ((sbf(i:i)=='"') .and. (sbf(j:j)=='"'))) then
          i = i + 1; j = j - 1
        endif
        if (ksp == 0) then
          xs = xs + sa*size*1.2
          ys = ys - ca*size*1.2
        endif
        call ptext(sbf(i:j), j-i+1, xs, ys, 0)
        ksp = 0
      else if (sbf(1:1) == '=') then
        j = 2
        call scanfv(j, sbf, xs)
        if (j > 0) then
          call scanfv(j, sbf, ys)
        endif
        if (j < 0) then
          write(buf,'(a17,i5,a8)') '  (Cap.) Error at', ns, '-th Line'
          call premsg(buf(1:30))
        endif
        ksp = 1
      else
        if (sbf(2:2) == ' ') then
          j = 2
        else if (sbf(3:3) == ' ') then
          j = 3
        else if (sbf(4:4) == ' ') then
          j = 4
        else
          j = -1
        endif
        if (j > 0) then
          fnt = sbf(1:j-1)
          call scanfv(j, sbf, size)
          if (j > 0) then
            angle = 0.;  itcol = 0;  ibcol = 9999
            call scanfv(j, sbf, angle)
            if (j > 0) then
              call scaniv(j, sbf, itcol)
              if (j > 0) then
                call scaniv(j, sbf, ibcol)
              endif
            endif
          endif
        endif
        if (j < 0) then
          write(buf,'(a17,i5,a8)') '  (Cap.) Error at', ns, '-th Line'
          call premsg(buf(1:30))
        else
          call lstyle(fnt, size, angle, itcol, ibcol)
          ca = cos(angle/RAD); sa = sin(angle/RAD)
        endif
      endif
    enddo
    close(12)
  endif

  call plote()
  call pscls()

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'XPLMAPCS: plot color-shading map for Grid data'
    write(99,'(a,a)') ' Input data filename   : ', fnin(1:lrtrim(fnin))
    write(99,'(a,i2)') '   dataset seq. number : ', nseq
    write(99,'(a12,a8,6x,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,'(a,a)') ' Output PS filename  : ', fnam(1:lrtrim(fnam))
    write(99,'(6x,a,i4,4x,a,i4,4x,a,f8.0,4x,a4)')   &
      'step=', istep, 'mval=', mval, 'scale  1 /', fscl, sht
    write(99,'(6x,a,f7.1,a,f7.1)')   &
      'illuminant :   Azimuth', az, ',   Elevation', el
    write(99,'(8x,a,f8.2)') 'value to be scaled to 1km :', vs
    write(99,'(6x,a,i2)') 'Raster resolution :  DPI = 600 /', nras
    write(99,'(6x,a,f5.1,a,f5.1)')   &
      'Margin      Left :', am, '   Bottom:', bm
    write(99,'(8x,a,f5.1,a,f5.1)')   &
      'Drawing   Width:', wide, '   Height:', high
    write(99,'(8x,a,f5.1,a,f5.1)')   &
      'Rspace    Right:', rs, '   Top   :', ts
    if (kcgl /= 0) then
      write(99,'(a,$)') '      with Color-Shading Legend  '
      write(99,'(a,f6.3,a)') '[ Size factor =', cglf, ' ]'
      write(99,'(a,f5.1,$)') '           at [ ', cglr
      write(99,'(a,f5.1,a,$)') ', ', cglu, ' ]'
      write(99,'(5x,a,a10,a)') '[ Unit notation = ', cunit, ']'
    endif
    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) write(99,'(a,f5.1,a,f5.1,a,$)')   &
                       '  ScaleBar [', sclr, ',', sclu, ' ]'
      write(99,'(a)')
    endif
    if (klbl /= 0) write(99,'(a,a)') '      Identifier : ', slbl
    if (fnobj(ldr:ldr) /= ' ')   &
      write(99,'(a,a)') ' Overlay obj. filename : ', fnobj(1:lrtrim(fnobj))
    if (fncap(ldr:ldr) /= ' ')   &
      write(99,'(a,a)') ' Caption data filename : ', fncap(1:lrtrim(fncap))
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif

  call premsg('    completed')
  stop
  end