c--------( SHADE : Plot shaded relief map for Grid data )--------
      parameter (MAXSZ=3001*3001)
      parameter (rad=(180./3.14159))
c
      dimension g(MAXSZ), s(MAXSZ)
      character area*8, s4*4, yn*1, spl*3/'   '/
      character out*72, buf*80, sdtm*20
      character*80 wdr, flog, fnin, fnam
      external cviken
c
      call premsg('SHADE : plot shaded relief map for Grid data')
      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(1,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(1,'(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(1,*) ixs, iys, mszx, mszy, mx, my, 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, mx, my
      call premsg(out)
      if((mx*my).gt.MAXSZ) call abendm('too large array size')
      read(1,*) ((g(i+j),i=0,(mx-1)*my,my),j=1,my)
      do 4 i=1,mx*my
        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,mx*my
        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)
      n = 0
      gdev = 0.
      do 7 i=1,mx*my
        if (g(i).eq.vnul) goto 7
        dg = g(i) - gsum
        n = n + 1
        gdev = gdev + dg*dg
    7 continue
      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('')
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
   10 call prompt('Process this set ?')
      call gparma('  "y" (Yes) or "s" (Skip) ==> ', 1, yn)
      if(yn.ne.'y'.and.yn.ne.'s') goto 10
      if(yn.eq.'s') goto 3
      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((my-1)*mszy) / wide
        high = float((mx-1)*mszx) / fscl
      else
        call gparmf('  Height of Drawing (in cm) ?  ==> ', high)
        fscl = float((mx-1)*mszx) / high
        wide = float((my-1)*mszy) / fscl
      endif
      fscl = fscl * 100.
      ccm = fscl/100000.
      kscl = 0
      sclr = 0.
      sclu = 0.
      klll = 0
      kshl = 0
c     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
c     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, mx, my
      write(9,*) az, el, vs, fscl, high, wide
      write(9,*) kscl, sclr, sclu, klll, kshl
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'SHADE : plot shaded relief 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, mx, my
      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 100 i=1,mx
      do 100 j=1,my
        n = (i-1)*my + j
        if(g(n).eq.vnul) then
          s(n) = 1.
        else
          if(i.eq.1.or.g(n-my).eq.vnul) then
            if(i.eq.mx.or.g(n+my).eq.vnul) then
              dx = 0.
            else
              dx = g(n+my) - g(n)
            endif
          else
            if(i.eq.mx.or.g(n+my).eq.vnul) then
              dx = g(n) - g(n-my)
            else
              dx = (g(n+my)-g(n-my)) / 2.
            endif
          endif
          if(j.eq.1.or.g(n-1).eq.vnul) then
            if(j.eq.my.or.g(n+1).eq.vnul) then
              dy = 0.
            else
              dy = g(n+1) - g(n)
            endif
          else
            if(j.eq.my.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
  100 continue
c
      call plots(1.5, 1.5)
      call dftone(-1)
      call dframe(0., 0., wide, high, my, mx)
      call paintm(s, -1., +1., 9.9999)
      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)
      call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb))
      if(klll.ne.0.or.kshl.ne.0) then
        xt = xs + float((mx-1)*mszx)/1000.
        yt = ys + float((my-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 300 kdf=ksf+1,ktf
            call cviken(float(idf), float(kdf), y, x)
            call plot((y-ys)/ccm, (x-xs)/ccm, 2)
  300     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 500 idf=isf+1,itf
            call cviken(float(idf), float(kdf), y, x)
            call plot((y-ys)/ccm, (x-xs)/ccm, 2)
  500     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., 0)
      call ptext(fnin, 50, 2., -1., 0)
      call plote()
    3 read(1,'(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(1)
      endfile(9)
      rewind(9)
      call pscls()
      call clspin()
c
      if(flog(ldr:ldr).ne.' ') then
        open(99,file=flog,access='append')
        write(99,'(//a,$)') 'SHADE :'
        write(99,'(a)')   ' plot shaded relief map for Grid data'
        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, mx, my
        read(9,*) 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, mx, my
        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   call strdtm(sdtm)
        write(99,'(a,a,a)') '=============== ', sdtm, '==============='
        close(99)
      endif
      close(9)
      stop
      end
