!--------( PLIMVC : Plot 3D image view color map )--------
  implicit none
  integer,parameter :: KVER = 2018, LVER = 4
  integer :: lwkdir, lrtrim
  real,allocatable :: f(:,:), g(:,:), h(:,:), gg(:), ds(:)
  character(80) :: wdr, flog, fnin, fobs, fnam, buf
  character(72) :: out;  character(20) :: sdtm
  character(8) :: area, area1, area2
  integer :: ldr, i, j, n, l, ix, iy
  integer :: kds, nds, nhx, nhy, nly, kly, nly2, nver
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ncc2, iorg2, korg2, ispa2, ispb2
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1
  integer :: ixs2, iys2, mszx2, mszy2, mx2, my2
  integer :: klp, klp2, kls, klt, kly2, kx, ky
  integer :: nkn, nke, nks, nkw, nlp, npx, npy
  real :: alt, alt2, cint, dh, fx, fy, gmax, gmin, gsum
  real :: step, vhi, vlo, vmed, vn, ve, vnul, vnul1, vnul2
  real :: x0, x1, x2, xs, xs1, y0, y1, y2, ys, ys1

  call premsg('PLIMVC : plot 3D image view color map')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnin = wdr;  fobs = wdr;  fnam = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Input AIMG data filename  ==> ', 81-ldr, fnin(ldr:80))
  open(10,file=fnin,status='old')
  n = 0;  l = 0;  nver = 0
  read(10,'(a)') buf
  do while (buf(1:1) == '#')
    if (buf(1:6) == '#AIMG:') then
      if (n == 0) then
        read(buf(7:80),*) kds, nds, nver, nhx, nhy
        if ((kds<0) .or. (kds>2) .or. (nds<1)) then
          call abendm('invalid AIMG file')
        endif
        l = (nds+7) / 8;  n = 1
        allocate(ds(nds))

      else if (n < l) then
        read(buf(7:80),*) (ds(i),i=n*8-7,n*8)
        n = n + 1
      else if (n == l) then
        read(buf(7:80),*) (ds(i),i=n*8-7,nds)
        n = l + 1
      endif
    endif
    read(10,'(a)') buf
  enddo
  if ((nver/=KVER) .and. (nver/=LVER)) call abendm('AIMG Version mismatch')
  if (n /= (l+1)) call abendm('AIMG file header broken')
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8,8x,2i4)')   &
        area, ncc, iorg, korg, ispa, ispb, nly, kly
    if ((ncc<0) .or. (ncc>=400)) call abendm('invalid coord.no.')
    if (ncc < 200) then
      nc = ncc; ncc = ncc + 800
    else
      ncc = ncc - 200; nc = ncc
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8,8x,2i4)')   &
        area, ncc, iorg, korg, ispa, ispb, nly, kly
    nc = ncc
    if ((ncc >= 800) .and. (ncc < 1000)) nc = ncc - 800
  endif
  if ((nly /= nds) .or. (kly /= 0)) call abendm('AIMG file broken')
  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,'(a)') buf
  read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt, klp
  if (alt /= -1.) call abendm('invalid Alt. value')

  allocate(g(mx,my));  allocate(h(mx,my))

  xs = float(ixs) / 1000.;  ys = float(iys) / 1000.
  read(10,*) ((h(i,j),i=1,mx),j=1,my)

  call gparma('Enter obs.anomaly data filename ==> ', 81-ldr, fobs(ldr:80))
  open(11,file=fobs,status='old')
  read(11,'(a)') buf
  do while (buf(1:1) == '#')
    read(11,'(a)') buf
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
    if ((ncc1<0) .or. (ncc1>=400)) call abendm('invalid coord.no.')
    if (ncc1 < 200) then
      nc1 = ncc1; ncc1 = ncc1 + 800
    else
      ncc1 = ncc1 - 200; nc1 = ncc1
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1
    nc1 = ncc1
    if ((ncc1 >= 800) .and. (ncc1 < 1000)) nc1 = ncc1 - 800
  endif
  if ((nc1 >= 1) .and. (nc1 <= 62)) then
    iorg1 = 0;  ispa1 = 0;  ispb1 = 0
    korg1 = (nc-30)*360 - 180
    if (nc1 > 60) korg1 = 0
    if (nc1 == 61) iorg1 = +5400
    if (nc1 == 62) iorg1 = -5400
  endif
  if ((area1 /= area) .or. (ncc1 /= ncc) .or.   &
      (iorg1 /= iorg) .or. (korg1 /= korg) .or.   &
      (ispa1 /= ispa) .or. (ispb1 /= ispb)) then
    call abendm('Header inconsistent')
  endif
  read(11,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1
  if ((mx1 <= 0) .or. (my1 <= 0)) call abendm('illegal array size')

  allocate(f(mx1,my1))

  xs1 = float(ixs1)/1000.;  ys1 = float(iys1)/1000.
  read(11,*) ((f(i,j),i=1,mx1),j=1,my1)
  close(11)

  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,a8)')   &
    'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'Layers'
  call premsg(out)
  write(out,'(6x,a4,6x,2f10.2,2i8,3i6)')   &
    'AIMG', xs, ys, mszx, mszy, mx, my, nly
  call premsg(out)
  write(out,'(6x,a4,6x,2f10.2,2i8,2i6)')   &
    'obs.', xs1, ys1, mszx1, mszy1, mx1, my1
  call premsg(out)
  call premsg('  Specify widths to limit drawing range')
  call gparmi2('    No. of grids on S-bound, N-bound ==> ', nks,nkn)
  if ((nks < 0) .or. (nkn < 0)) call abendm('invalid value')
  npx = mx - nks - nkn
  if (npx < 2) call abendm('invalid value')
  call gparmi2('    No. of grids on W-bound, E-bound ==> ', nkw,nke)
  if ((nkw < 0) .or. (nke < 0)) call abendm('invalid value')
  npy = my - nkw - nke
  if (npy < 2) call abendm('invalid value')

  allocate(gg(npx*npy))

  call gparmi('  Topmost layer no. to draw ==> ', kls)
  if ((kls <= 0) .or. (kls >= (nly-1)))   &
    call abendm('illegal layer no.')
  call gparmi('  Bottommost layer no. to draw  ==> ', klt)
  if ((klt < (kls+2)) .or. (klt > nly))   &
    call abendm('illegal layer no.')
  nlp = klt-kls+1
  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')

  n = 0
  do
    read(10,'(a)') buf
    do while (buf(1:1) == '#')
      read(10,'(a)') buf
    enddo
    kly = kly + 1
    if (buf(13:16) == '    ') then        !-- v2005
      read(buf,'(a8,i4,4x,4i8,8x,2i4)')   &
          area2, ncc2, iorg2, korg2, ispa2, ispb2, nly2, kly2
      if ((ncc2<0) .or. (ncc2>=400)) call abendm('invalid coord.no.')
      if (ncc2 < 200) then
        ncc2 = ncc2 + 800
      else
        ncc2 = ncc2 - 200
      endif
    else                                  !-- v2018
      read(buf,'(a8,4x,i4,4i8,8x,2i4)')   &
          area2, ncc2, iorg2, korg2, ispa2, ispb2, nly2, kly2
    endif
    if ((area2 /= area) .or. (ncc2 /= ncc))   &
      call abendm('Header1 inconsistent')
    if ((nly2 /= nly) .or. (kly2 /= kly)) call abendm('AIMG file broken')
    read(10,'(a)') buf
    read(buf,*) ixs2,iys2, mszx2,mszy2, mx2,my2, vnul2, alt2, klp2
    if ((ixs2 /= ixs) .or. (iys2 /= iys) .or. (mszx2 /= mszx) .or.   &
        (mszy2 /= mszy) .or. (mx2 /= mx) .or. (my2 /= my) .or.   &
        (vnul2 /= vnul)) call abendm('Header2 inconsistent')
    read(10,*) ((g(i,j),i=1,mx),j=1,my)
    if (kly < kls) cycle
    do i=nks+1,mx-nkn
      ix = ixs + mszx*(i-1)
      kx = nint(float(ix-ixs1)/float(mszx1))
      if ((kx <= 0) .or. (kx > mx1)) cycle
      do j=nkw+1,my-nke
        iy = iys + mszy*(j-1)
        ky = nint(float(iy-iys1)/float(mszy1))
        if ((ky <= 0) .or. (ky > my1)) cycle
        if (f(kx,ky) == vnul1) cycle
        if ((kds == 0) .and. ((h(i,j)+ds(kly)) <= 0.)) cycle
        if (n == 0) then
          n = 1;  gsum = g(i,j);  gmin = g(i,j);  gmax = g(i,j)
        else
          n = n + 1;  gsum = gsum + g(i,j)
          if (g(i,j) < gmin) gmin = g(i,j)
          if (g(i,j) > gmax) gmax = g(i,j)
        endif
      enddo
    enddo
    if (kly >= klt) exit
  enddo
  gsum = gsum / float(n)
  rewind(10)

  call premsg(' ')
  write(out,'(12x,a,i3,a,i3,a)')   &
    'plot layers   from', kls, '  ( to', klt, ' )'
  call premsg(out)
  write(out,'(4x,a,f7.1,a,f7.1,a,f7.1)')   &
    'Data Range : [', gmin/100., ' ,', gmax/100., ' ]    Average :', gsum/100.
  call premsg(out)
  call premsg('')
  call gparmf('  Specify color-grading interval (A/m) ==> ', step)
  if (step <= 0.) call abendm('invalid parameter')
  call gparmf('  Median value of color-grading (A/m) ? ==> ', vmed)
  vlo = vmed - step*10.;  vhi = vmed + step*10.
  call premsg('  Specify contour-line interval (A/m)')
  call gparmf('    ( 0 for no line contours)  ==> ', cint)
  if (cint < 0.) call abendm('invalid parameter')
  call clspin()
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'PLIMVC : plot 3D image view color map'
  write(6,'(a,a)') ' Input AIMG filename : ', fnin(1:lrtrim(fnin))
  write(6,'(a,a)') ' Obs. anom. filename : ', fobs(1:lrtrim(fobs))
  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,a8)')   &
    'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', 'Layers'
  write(6,'(16x,2f10.2,2i8,3i6)') xs,ys, mszx,mszy, mx,my, nly
  write(6,'(6x,a,i3,a,i3,4x,a,i4,4x,a,i4)')   &
    'plot layers : ', kls, ' to ', klt, 'npx=', npx, 'npy=', npy
  write(6,'(34x,a,i3,3x,a,i3,3x,a,i3,3x,a,i3)')   &
    'LmS=', nks, 'LmN=', nkn, 'LmW=', nkw, 'LmE=', nke
  write(6,'(6x,a,f6.3,1x,a,4x,a,f6.3)')   &
    'step=', step, 'A/m', 'mval=', vmed
  write(6,'(6x,a,f6.3,1x,a)')   &
    'line-contour interval :', cint, 'A/m'

  read(10,'(a)') buf
  do while (buf(1:1) == '#')
    read(10,'(a)') buf
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8,8x,2i4)')   &
        area2, ncc2, iorg2, korg2, ispa2, ispb2, nly2, kly
    if ((ncc2<0) .or. (ncc2>=400)) call abendm('invalid coord.no.')
    if (ncc2 < 200) then
      ncc2 = ncc2 + 800
    else
      ncc2 = ncc2 - 200
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8,8x,2i4)')   &
        area2, ncc2, iorg2, korg2, ispa2, ispb2, nly2, kly
  endif
  if ((area2 /= area) .or. (ncc2 /= ncc))   &
    call abendm('Header1 inconsistent')
  if ((nly2 /= nly) .or. (kly /= 0)) call abendm('AIMG file broken')
  read(10,'(a)') buf
  read(buf,*) ixs2,iys2, mszx2,mszy2, mx2,my2, vnul2, alt2, klp2
  if ((ixs2 /= ixs) .or. (iys2 /= iys) .or. (mszx2 /= mszx) .or.   &
      (mszy2 /= mszy) .or. (mx2 /= mx) .or. (my2 /= my) .or.   &
      (vnul2 /= vnul)) call abendm('Header2 inconsistent')
  if (alt2 /= -1.) call abendm('invalid Alt. value')
  read(10,*) ((h(i,j),i=1,mx),j=1,my)

  vn = float(mszx*(npx-1))/1000.;  ve = float(mszy*(npy-1))/1000.
  call psopn(fnam, 'A4P')
  call plots(0., 0.)
  call lstyle('T', 0.4, 0., 0, -255)
  call ptext(area, 8,  1., 0.5, 0)
  call ptext(fnin, 50, 4., 0.5, 0)
  dh = 27.72 / float(nlp+1)
  call plot(10., dh*float(nlp), -3)
  if (ve >= vn) then
    fx = (vn*0.6 + ve*0.8) / 11.2;  fy = fx * 10.5 / dh
    x0 = -5.6;  y0 = (ve*0.3 - vn*0.4) / fy
    x1 = x0 + ve*0.8 / fx;  y1 = y0 - ve*0.6 / fy
    x2 = x0 + vn*0.6 / fx;  y2 = y0 + vn*0.8 / fy
  else
    fx = (vn*0.8 + ve*0.6) / 11.2;  fy = fx * 10.5 / dh
    x0 = -5.6;  y0 = (ve*0.4 - vn*0.3) / fy
    x1 = x0 + ve*0.6 / fx;  y1 = y0 - ve*0.8 / fy
    x2 = x0 + vn*0.8 / fx;  y2 = y0 + vn*0.6 / fy
  endif

  do
    read(10,'(a)') buf
    do while (buf(1:1) == '#')
      read(10,'(a)') buf
    enddo
    kly = kly + 1
    if (buf(13:16) == '    ') then        !-- v2005
      read(buf,'(a8,i4,4x,4i8,8x,2i4)')   &
          area2, ncc2, iorg2, korg2, ispa2, ispb2, nly2, kly2
      if ((ncc2<0) .or. (ncc2>=400)) call abendm('invalid coord.no.')
      if (ncc2 < 200) then
        ncc2 = ncc2 + 800
      else
        ncc2 = ncc2 - 200
      endif
    else                                  !-- v2018
      read(buf,'(a8,4x,i4,4i8,8x,2i4)')   &
          area2, ncc2, iorg2, korg2, ispa2, ispb2, nly2, kly2
    endif
    if ((area2 /= area) .or. (ncc2 /= ncc))   &
      call abendm('Header1 inconsistent')
    if ((nly2 /= nly) .or. (kly2 /= kly)) call abendm('AIMG file broken')
    read(10,'(a)') buf
    read(buf,*) ixs2,iys2, mszx2,mszy2, mx2,my2, vnul2, alt2, klp2
    if ((ixs2 /= ixs) .or. (iys2 /= iys) .or. (mszx2 /= mszx) .or.   &
        (mszy2 /= mszy) .or. (mx2 /= mx) .or. (my2 /= my) .or.   &
        (vnul2 /= vnul)) call abendm('Header2 inconsistent')
    if (alt2 /= -1.) call abendm('invalid Alt. value')
    read(10,*) ((g(i,j),i=1,mx),j=1,my)
    if (kly < kls) cycle
    n = 0
    do i=nks+1,mx-nkn
      ix = ixs + mszx*(i-1)
      kx = nint(float(ix-ixs1)/float(mszx1)) + 1
      if ((kx <= 0) .or. (kx > mx1)) then
        do j=nkw+1,my-nke
          n = n + 1;  gg(n) = vnul
        enddo
      else
        do j=nkw+1,my-nke
          iy = iys + mszy*(j-1)
          ky = nint(float(iy-iys1)/float(mszy1)) + 1
          n = n + 1
          if ((ky <= 0) .or. (ky > my1)) then
            gg(n) = vnul
          else if (f(kx,ky) == vnul1) then
            gg(n) = vnul
          else if ((kds == 0) .and. ((h(i,j)+ds(kly)) <= 0.)) then
            gg(n) = vnul
          else
            gg(n) = g(i,j) / 100.
          endif
        enddo
      endif
    enddo
    call dfrgbt(0)
    if (cint == 0.) call dresol(2)
    call dframo(x0, y0, x1, y1, x2, y2, npy, npx)
    call paintm(gg, vlo, vhi, vnul)
    if (cint /= 0.) then
      call contso(x0, y0, x1, y1, x2, y2, npy, npx)
      call contr(gg, vnul, cint, 1, -99.999, 99.999, 0)
    endif
    call newpen(2)
    call plot(x0, y0, 3)
    call plot(x1, y1, 2)
    call plot(x1+x2-x0, y1+y2-y0, 2)
    call plot(x2, y2, 2)
    call plot(x0, y0, 2)
    if (kly == klt) exit
    call plot(0., -dh, -3)
  enddo
  close(10)
  call plote()
  call pscls()

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'PLIMVC : plot 3D image view color map'
    write(99,'(a,a)') ' Input AIMG filename : ', fnin(1:lrtrim(fnin))
    write(99,'(a,a)') ' Obs. anom. filename : ', fobs(1:lrtrim(fobs))
    write(99,'(a,a)') ' Output PS filename  : ', fnam(1:lrtrim(fnam))
    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,i3,a,i3,4x,a,i4,4x,a,i4)')   &
      'plot layers : ',kls, ' to ',klt, 'npx=',npx, 'npy=',npy
    write(99,'(34x,a,i3,3x,a,i3,3x,a,i3,3x,a,i3)')   &
      'LmS=', nks, 'LmN=', nkn, 'LmW=', nkw, 'LmE=', nke
    write(99,'(6x,a,f6.3,1x,a,4x,a,f6.3)')   &
      'step=', step, 'A/m', 'mval=', vmed
    write(99,'(6x,a,f6.3,1x,a)')   &
      'line-contour interval :', cint, 'A/m'
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end