!--------( PLXIM : Arbitrary-section contour map of 3D image )--------
  implicit none
  integer,parameter :: KVER = 2018, LVER = 4
  integer :: lwkdir, lrtrim
  real,allocatable :: f(:,:), v(:,:), g(:,:), z(:), ds(:), dc(:)
  real,allocatable :: spx(:), spy(:)
  character(80) :: wdr, flog, fnin, fnam, buf
  character(72) :: out;  character(40) :: str;  character(20) :: sdtm
  character(8) :: area, area1;  character(3) :: spl = '   '
  integer :: ldr, i, j, l, n, nhx, nhy, nly, kly, nver, nds
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: ixs, iys, mszx, mszy, mx, my
  integer :: kds, kpl, nlp, np, nsl, lcap, kx, ky
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs1, iys1, mszx1, mszy1, mx1, my1, nly1, kly1, nlp1
  real :: alt, dl, dll, dz, f1, f2, high, hlen, hscl
  real :: csz, step, thick, tmin, vmax, vmin, astp
  real :: vnul, vscl, vsum, wide, xs, ys, zmax, zmin, zt
  real :: dx, dy, px, py, s1x, s1y, s2x, s2y, sf1, sf2, szmx, szmy, xt, yt
  real :: vnul1, alt1

  call premsg('PLXIM : Arbitrary-section contour map of 3D image')
  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 AIMG data filename  ==> ', 81-ldr, fnin(ldr:80))
  open(10,file=fnin,status='old')
  n = 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));  allocate(dc(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, nlp
  if (alt /= -1.) call abendm('invalid Alt. value')

  allocate(f(mx,my))

  xs = float(ixs) / 1000.;  xt = xs + float(mszx*(mx-1))/1000.
  ys = float(iys) / 1000.;  yt = ys + float(mszy*(my-1))/1000.
  read(10,*) ((f(i,j),i=1,mx),j=1,my)

  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,'(16x,2f10.2,2i8,3i6)') xs,ys, mszx,mszy, mx,my, nly
  call premsg(out)
  call premsg('Specify Both ends of section (in km)')
  call gparmf2('  Northing(X),Easting(Y) of Start position ==> ', s1x, s1y)
  if ((s1x <= xs) .or. (s1x >= xt)) call abendm('out of range')
  if ((s1y <= ys) .or. (s1y >= yt)) call abendm('out of range')
  call gparmf2('  Northing(X),Easting(Y) of End position   ==> ', s2x, s2y)
  if ((s2x <= xs) .or. (s2x >= xt)) call abendm('out of range')
  if ((s2y <= ys) .or. (s2y >= yt)) call abendm('out of range')
  hlen = sqrt((s2x-s1x)*(s2x-s1x) + (s2y-s1y)*(s2y-s1y)) * 1000.
  call gparmi(' How many points to conform the section ? ==> ', np)
  if (np < 2) call abendm('invalid value')
  write(str,'(2(a,f8.3,1h,,f8.3),a)') ' (', s1x,s1y, ')-(', s2x,s2y, ')'

  allocate(z(np));  allocate(spx(np));  allocate(spy(np))
  allocate(v(np,nly))

  dx = (s2x-s1x) / float(np-1);  dy = (s2y-s1y) / float(np-1)
  spx(1) = s1x;  spy(1) = s1y
  do i=2,np
    spx(i) = spx(i-1) + dx;  spy(i) = spy(i-1) + dy
  enddo
  szmx = float(mszx) / 1000.;  szmy = float(mszy) / 1000.
  do i=1,np
    px = (spx(i)-xs) / szmx;  py = (spy(i)-ys) / szmy
    kx = ifix(px);  dx = px - float(kx);  ky = ifix(py);  dy = py - float(ky)
    if ((f(kx+1,ky+1) == vnul) .or. (f(kx+2,ky+1) == vnul) .or.   &
        (f(kx+1,ky+2) == vnul) .or. (f(kx+2,ky+2) == vnul)) then
      z(i) = vnul
    else
      sf1 = f(kx+1,ky+1)*(1.-dx) + f(kx+2,ky+1)*dx
      sf2 = f(kx+1,ky+2)*(1.-dx) + f(kx+2,ky+2)*dx
      z(i) = -(sf1*(1.-dy) + sf2*dy)            !---  z():depth <== f():alt
    endif
  enddo
  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')

  do l=1,nly
    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)')   &
          area1, ncc1, iorg1, korg1, ispa1, ispb1, nly1, kly1
      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,8x,2i4)')   &
          area1, ncc1, iorg1, korg1, ispa1, ispb1, nly1, kly1
      nc1 = ncc1
      if ((ncc1 >= 800) .and. (ncc1 < 1000)) nc1 = ncc1 - 800
    endif
    if ((nly1 <= 0) .or. (kly1 /= l)) call abendm('AIMG file broken')
    read(10,'(a)') buf
    read(buf,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1, alt1, nlp1
    read(10,*) ((f(i,j),i=1,mx),j=1,my)
    do i=1,np
      px = (spx(i)-xs) / szmx;  py = (spy(i)-ys) / szmy
      kx = ifix(px);  ky = ifix(py)
      dx = px - float(kx);  dy = py - float(ky)
      if ((f(kx+1,ky+1) == vnul) .or. (f(kx+2,ky+1) == vnul) .or.   &
          (f(kx+1,ky+2) == vnul) .or. (f(kx+2,ky+2) == vnul)) then
        v(i,l) = vnul
      else
        sf1 = f(kx+1,ky+1)*(1.-dx) + f(kx+2,ky+1)*dx
        sf2 = f(kx+1,ky+2)*(1.-dx) + f(kx+2,ky+2)*dx
        v(i,l) = (sf1*(1.-dy) + sf2*dy) / 100.
      endif
    enddo
  enddo
  close(10)

  n = -1
  do l=1,nly
    do i=1,np
      if (v(i,l) /= vnul) then
        n = 0;  vsum = 0.;  vmin = v(i,l);  vmax = v(i,l);  exit
      endif
    enddo
    if (n == 0) exit
  enddo
  if (n < 0) call abendm('no valid data')
  do l=1,nly; do i=1,np
    if (v(i,l) == vnul) cycle
    n = n + 1;  vsum = vsum + v(i,l)
    if (v(i,l) < vmin) vmin = v(i,l)
    if (v(i,l) > vmax) vmax = v(i,l)
  enddo; enddo
  vsum = vsum / float(n)
  zmin = z(1);  zmax = z(1)
  do i=2,np
    if (z(i) < zmin) zmin = z(i)
    if (z(i) > zmax) zmax = z(i)
  enddo
!---  zmin: depth of highest point of Terr Surf
!---  zmax: depth of lowest point of Terr Surf
  if (kds == 0) then
    zmax = ds(nly)
  else
    zmax = zmax + ds(nly)
  endif
!---  zmax: depth of deepest point of Bottom Layer bottom
  call premsg(' ')
  write(out,'(4x,a,f7.2,a,f7.2,a,f7.2)')   &
    'Data Range :  [', vmin, ' ,', vmax, ' ]    Average :', vsum
  call premsg(out)
  write(out,'(4x,a,f7.1,a,f7.1,a)')   &
    'Depth Range : [', zmin, ' ,', zmax, ' ] (m)'
  call premsg(out)
  call premsg('')
!---  kds   Layer config          bottom ds()  center dc()  ds(0)
!     0   Var.Depth Horizon         depth below Sea Level  -99999.
!     1   Const.Thick below Terr    depth below Terr Surf      0.
!---  2   Var.Thick below Terr      depth below Terr Surf      0,
  if (kds == 0) then
    tmin = ds(2) - ds(1)                !---  dc(1) undefined
  else
    tmin = ds(1);  dc(1) = ds(1) / 2.
  endif
  do l=2,nly
    thick = ds(l) - ds(l-1)
    if (thick < tmin) tmin = thick
    dc(l) = (ds(l-1)+ds(l)) / 2.
  enddo
  nsl = ifix((zmax-zmin)*3. / tmin) + 2
  dz = (zmax-zmin) / float(nsl-1)
!---  Define evenly spaced depth nodes for section mapping,
!      so that even the thinnest layer would be devided into 3.
!---    (total <nsl> nodes (including both ends) with the spacing <dz>)

  allocate(g(nsl,np))

  do i=1,np
    j = 0;  zt = zmin;  dl = z(i)
    do while (zt < dl)
      j = j + 1;  g(j,i) = vnul;  zt = zt + dz
    enddo
    if (kds == 0) then
      dl = (z(i)+ds(1)) / 2.
    else
      dl = z(i) + dc(1)
    endif
    do while (zt < dl)
      j = j + 1;  g(j,i) = v(i,1);  zt = zt + dz
    enddo
    do l=2,nly
      dll = dl
      if (kds == 0) then
        dl = dc(l)
      else
        dl = z(i) + dc(l)
      endif
      do while (zt < dl)
        f1 = zt - dll;  f2 = dl - zt;  j = j + 1
        if ((v(i,l-1) == vnul) .or. (v(i,l) == vnul)) then
          g(j,i) = vnul
          if (f1 == 0.) g(j,i) = v(i,l-1)
        else
          g(j,i) = (v(i,l-1)*f2 + v(i,l)*f1) / (f1 + f2)
        endif
        zt = zt + dz
      enddo
    enddo
    if (kds == 0) then
      dl = ds(nly)
    else
      dl = z(i) + ds(nly)
    endif
    do while (zt < dl)
      j = j + 1;  g(j,i) = v(i,nly);  zt = zt + dz
    enddo
    do while (j < nsl)
      j = j + 1;  g(j,i) = vnul
    enddo
  enddo

  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 gparmf('  Width of Drawing (in cm) ?  ==> ', wide)
  if ((wide <= 0.) .or. (wide > 25.)) call abendm('invalid width')
  if ((kpl == 1) .and. (wide > 18.)) call abendm('invalid width')
  hscl = hlen * 100. / wide
  call gparmf('  Height of Drawing (in cm) ?  ==> ', high)
  if ((high <= 0.) .or. (high > 25.)) call abendm('invalid height')
  if ((kpl == 2) .and. (high > 18.)) call abendm('invalid height')
  vscl = (zmax-zmin) * 100. / high
  call premsg('  Specify contour interval (A/m)')
  call gparmf('   (minus to draw even crowded contour) ==> ', step)
  if (step == 0.) call abendm('invalid parameter')
  astp = abs(step);  lcap = 5
  if ((astp == 5.) .or. (astp == 0.5) .or. (astp == 0.05)) lcap = 4
  if ((astp == 2.5) .or. (astp == 0.25) .or. (astp == 0.025)) lcap = 4
  if (step < 0.) lcap = -lcap
  call premsg('  Char-size of contour value (in cm)')
  call gparmf('   (If zero (0.), no values drawn.)  ==> ', csz)
  if (csz < 0.) call abendm('invalid character size')
  call clspin()
  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'PLXIM : Arbitrary-section contour map of 3D image'
  write(6,'(a,a)') ' Input AIMG 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,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,a,4x,a,i4)')   &
    'Section :  ', str, 'np=', np
  write(6,'(8x,a,f7.1,a,f7.1,a,4x,a,i4)')   &
    'Depth range : [', zmin, ' ,', zmax, ' ] (m)', 'nsl=', nsl
  write(6,'(6x,2a,4x,a,f8.0,4x,a,f8.0)')   &
    'Sheet: ', spl, 'H-scale  1 /', hscl, 'V-scale  1 /', vscl
  write(6,'(6x,a,f6.3,1x,a,4x,a,i3,4x,f8.1,a1,f8.1,a1)')   &
    'step=', step, 'A/m', 'lcapt=', lcap, wide, 'W', high, 'H'
  write(6,'(8x,a,f6.2,a)') 'char.size of cont.value : ', csz, ' cm'

  call psopn(fnam, spl)
  call plots(1.5, 1.5)
  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 conts(0., high, wide, -high, -np, -nsl, csz)
  call contr(g, vnul, astp, lcap, -99.999, 0., 11)
  call contr(g, vnul, astp, lcap, -astp/2., 99.999, 10)
  call newpen(2)
  call wrect(0., 0., wide, high)
  call newpen(1)
  call ptext(str, 40, 0., high+1., 0)
  call plote()
  call pscls()

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'PLXIM : Arbitrary-section contour map of 3D image'
    write(99,'(a,a)') ' Input AIMG filename : ', fnin(1:lrtrim(fnin))
    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,a,4x,a,i4)') 'Section :  ', str, 'np=', np
    write(99,'(8x,a,f7.1,a,f7.1,a,4x,a,i4)')   &
      'Depth range : [', zmin, ' ,', zmax, ' ] (m)', 'nsl=', nsl
    write(99,'(6x,2a,4x,a,f8.0,4x,a,f8.0)')   &
      'Sheet: ', spl, 'H-scale  1 /', hscl, 'V-scale  1 /', vscl
    write(99,'(6x,a,f6.3,1x,a,4x,a,i3,4x,f8.1,a1,f8.1,a1)')   &
      'step=', step, 'A/m', 'lcapt=', lcap, wide, 'W', high, 'H'
    write(99,'(8x,a,f6.2,a)') 'char.size of cont.value : ', csz, ' cm'
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end