c--------( PLMAPCS : Plot color-shading map without contours )--------
      parameter (MAXSZ=3001*3001)
      parameter (rad=(180./3.14159))
c
      dimension g(MAXSZ), s(MAXSZ)
      character*80 wdr, flog, fnin, fnam
      character out*72, buf*80, area*8, s4*4, yn*1, spl*3/'   '/
      external cviken
c
      call premsg('PLMAPCS : plot color-shading without contours')
      call opnpin()
      ldr = lwkdir(50,wdr) + 1
      flog = wdr
      fnin = 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 Output PS filename  ==> ',
     *            81-ldr, fnam(ldr:80))
      if(fnam(ldr:ldr).eq.' ') fnam(1:1) = char(0)
      if (fnam(80:80).ne.' ') then
        call abendm('too long filename')
      endif
      open(9,status='scratch')
c
    1 read(10,'(a)',end=8) buf
      if(buf(1:1).eq.'#') goto 1
    2 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.
      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)
      if((nmx*nmy).gt.MAXSZ) then
        call abendm('too many grids')
      endif
      read(10,*) ((g(i+j),i=0,(nmx-1)*nmy,nmy),j=1,nmy)
      do 4 i=1,nmx*nmy
        if(g(i).ne.vnul) goto 5
    4 continue
      call abendm('no valid data')
    5 n = 1
      gsum = g(i)
      gmin = g(i)
      gmax = g(i)
      do 6 j=i+1,nmx*nmy
        if(g(j).eq.vnul) goto 6
        n = n + 1
        gsum = gsum + g(j)
        if(g(j).lt.gmin) gmin = g(j)
        if(g(j).gt.gmax) gmax = g(j)
    6 continue
      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 premsg('')
c
      if(spl.eq.'   ') then
        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 psopn(fnam, spl)
      endif
      call premsg('  Specify color grading interval '
     *            // '(1/20 of full-range)')
      call gparmi('    ( enter 0 to skip mapping ) ==> ', istep)
      if(istep.eq.0) goto 3
      istp = iabs(istep)
      call gparmi('  Median value of color-grading ?  ==> ', mval)
      mval = mval/istp*istp
      vlo = float(mval - istp*10)
      vhi = float(mval + istp*10)
      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)
      if(kpl.eq.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
      fscl = fscl * 100.
      ccm = fscl/100000.
      kscl = 0
      sclr = 0.
      sclu = 0.
      klll = 0
      kshl = 0
      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
      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
      write(9,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb
      write(9,*) xs, ys, mszx, mszy, nmx, nmy
      write(9,*) istp, mval, az, el, vs, fscl, high, wide
      write(9,*) kscl, sclr, sclu, klll, kshl
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'PLMAPCS : plot color-shading map for Grid data'
      write(6,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin))
      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=', istp, 'mval=', mval
      write(6,'(6x,a,f7.1,a,f7.1)')
     *      'illuminant :   Azimuth', az, ',   Elevation', el
      write(6,'(6x,a,f8.2)') ' value to be scaled to 1km :', vs
      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
      do 50 i=1,nmx
      do 50 j=1,nmy
        n = (i-1)*nmy + j
        if (g(n).eq.vnul) then
          s(n) = 1.
        else
          if (i.eq.1.or.g(n-nmy).eq.vnul) then
            if (i.eq.nmx.or.g(n+nmy).eq.vnul) then
              dx = 0.
            else
              dx = g(n+nmy) - g(n)
            endif
          else
            if (i.eq.nmx.or.g(n+nmy).eq.vnul) then
              dx = g(n) - g(n-nmy)
            else
              dx = (g(n+nmy)-g(n-nmy)) / 2.
            endif
          endif
          if (j.eq.1.or.g(n-1).eq.vnul) then
            if (j.eq.nmy.or.g(n+1).eq.vnul) then
              dy = 0.
            else
              dy = g(n+1) - g(n)
            endif
          else
            if (j.eq.nmy.or.g(n+1).eq.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
c          s(n) = (1.-sn) / 2.
          if (sn.lt.0.) then
            s(n) = 0.75 - sn/4.
          else
            s(n) = (1.-sn) * 0.75
          endif
        endif
   50 continue
c
      call plots(2., 2.)
      call dfcols(20, 0.85, 1., -1)
      call dresol(2)
      call dframe(0., 0., wide, high, nmy, nmx)
      call paintw(g, vlo, vhi, vnul, s)
      call wrect(0., 0., wide, high)
      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)
      if(klll.ne.0.or.kshl.ne.0) then
        call cvinit(ncc, float(iorg), float(korg),
     *                   float(ispa), float(ispb))
        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.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-ys)/ccm, (x-xs)/ccm, 3)
          do 100 kdf=ksf+1,ktf
            call cviken(float(idf), float(kdf), y, x)
            call plot((y-ys)/ccm, (x-xs)/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-ys)/ccm, (x-xs)/ccm, 3)
          do 300 idf=isf+1,itf
            call cviken(float(idf), float(kdf), y, x)
            call plot((y-ys)/ccm, (x-xs)/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(-ys/ccm, -xs/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.5, 0)
      call ptext(fnin, 50, 2., -1.5, 0)
      call plote()
    3 read(10,'(a)',end=8) buf
      if(buf(1:1).eq.'#') goto 3
   20 call prompt('  Process Next set ?')
      call gparma('  "y" (Yes) or "n" (No) ==> ', 1, yn)
      if(yn.ne.'y'.and.yn.ne.'n') goto 20
      if(yn.eq.'y') goto 2
c
    8 close(10)
      endfile(9)
      rewind(9)
      call pscls()
      call clspin()
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,status='old',iostat=ios)
        if(ios.eq.0) then
   90     read(99,*,end=91)
          goto 90
   91     backspace(99)
        else
          open(99,file=flog,status='new')
        endif
        write(99,'(//a,$)') 'PLMAPCS :'
        write(99,'(a)')   ' plot color-shading without contours'
        write(99,'(a,a)')
     *       ' Input data filename : ', fnin(1:lrtrim(fnin))
        write(99,'(a,a)')
     *       ' Output PS filename  : ', fnam(1:lrtrim(fnam))
   92   read(9,'(a8,i4,4x,4i8)',end=93) area, ncc, iorg,korg, ispa,ispb
        read(9,*) xs, ys, mszx, mszy, nmx, nmy
        read(9,*) istp, mval, az, el, vs, fscl, high, wide
        read(9,*) kscl, sclr, sclu, klll, kshl
        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=', istp, 'mval=', mval
        write(99,'(6x,a,f7.1,a,f7.1)')
     *      'illuminant :   Azimuth', az, ',   Elevation', el
        write(99,'(6x,a,f8.2)') ' value to be scaled to 1km :', vs
        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
        goto 92
   93   close(99)
      endif
      close(9)
      stop
      end
