!--------( PLMAPS : Plot contour map with shading )--------
  implicit none
  integer :: lwkdir, lrtrim
  real,allocatable :: g(:), s(:)
  real,parameter :: RAD = (180./3.141593)
  character(80) :: wdr, flog, fnin, 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, lcap, 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
  real :: vnul, xs, ys, gsum, gmin, gmax, gdev, dg
  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
  real :: az, cx, cy, cz, dx, dy, el, fx, fy, p, q, r, sn, vs
  external cviken

  call premsg('PLMAPS : Plot contour map with shading')
  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) == ' ') 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));  allocate(s(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 (minus to draw even crowded)')
    call gparmi('    ( enter 0 to skip countouring ) ==> ', istep)

    if (istep /= 0) then
      istp = iabs(istep)
      lcap = 5
      if ((istp == 25) .or. (istp == 50) .or. (istp == 250)) lcap = 4
      if (istep < 0) lcap = -lcap
      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 == 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
      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, lcap, az, el, vs, fscl, high, wide
      write(9,*) kscl, sclr, sclu, klll, kshl
      write(6,'(a)') '----------------------------------------'
      write(6,'(a)') 'PLMAPS : Plot contour map with shading'
      write(6,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin))
      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, 'lcapt=', lcap
      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 /= 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)
      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 plots(1.5, 1.5)
      call dftone(-1)
      call dframe(0., 0., wide, high, nmy, nmx)
      call paintm(s, -1., +1., 9.9999)
      call wrect(0., 0., wide, high)
      call conts(0., 0., wide, high, nmy, nmx, 0.)
      call contx(g, vnul, istp, lcap, -99999, -1, 11)
      call contx(g, vnul, istp, lcap,  0,  99999, 10)
      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 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); deallocate(s)

    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)') 'PLMAPS : Plot contour map with shading'
    write(99,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin))
    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, lcap, az, el, vs, fscl, high, wide
      read(9,*) kscl, sclr, sclu, klll, kshl
      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,'(6x,a,i4,4x,a,i4)') 'step=', istp, 'lcapt=', lcap
      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 /= 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