!--------( CALMAS : calc. model anomaly on a given surface )--------
  implicit none
  integer :: lwkdir, lrtrim; real :: calma
  real,allocatable :: f(:,:), h(:,:)
  real :: pv(7)
  character(5) :: nmdl(4) = (/'Block', 'hRect', 'vLine', 'Point'/)
  character(80) :: wdr, flog, fnam1, fnam2, buf
  character(20) :: sdtm;  character(8) :: area
  integer :: ldr, n, i, j, kl, ks, IER
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: ixl, ixh, iyl, iyh, mszx, mszy, mx, my
  real :: am, cx, cy, dip, dec, hnul, phi, the, vol, rx, ry, rz
  real :: xc, xl, xp, yc, yl, yp, zc, za, zp

  call premsg('CALMAS : calc. model anomaly on a given surface')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = wdr;  fnam2 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Calc.Surface Alt. filename ==> ', 81-ldr, fnam1(ldr:80))
  open(1,file=fnam1,status='old')
  do
    read(1,'(a)') buf
    if (buf(1:1) /= '#') exit
  enddo
  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
    else
      ncc = ncc - 200; nc = ncc
    endif
  else                                  !-- v2018
    read(buf,'(a8,4x,i4,4i8)') area, ncc, iorg, korg, ispa, ispb
    nc = ncc
    if ((ncc >= 800) .and. (ncc < 1000)) nc = ncc - 800
  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(1,'(a)') buf
  read(buf,*) ixl, iyl, mszx, mszy, mx, my, hnul

  if ((mx <= 0) .or. (my <= 0)) call abendm('illegal array size')
  allocate(f(mx,my));  allocate(h(mx,my))

  read(1,*) ((h(i,j),i=1,mx),j=1,my)
  close(1)
  xl = float(ixl)/1000.;  yl = float(iyl)/1000.
  do j=1,my; do i=1,mx;  f(i,j) = 0.;  enddo; enddo
  ixh = ixl + mszx*(mx-1);  iyh = iyl + mszy*(my-1)
  call gparma('Enter Calc. MA output filename ==> ', 81-ldr, fnam2(ldr:80))
  open(2,file=fnam2,status='new')

  write(6,'(a)') '----------------------------------------'
  write(6,'(a)') 'CALMAS : calc. model anomaly on a given surface'
  write(6,'(a,a)')   &
    ' Calc. alt. data infile : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)')   &
    ' Model anomaly outfile  : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    'Areaname :  ', area, 'Coord.', ncc, iorg, korg, ispa, ispb
  write(6,'(4x,a)')   &
    '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'calc. surf. ', xl, yl, mszx, mszy, mx, my, '    --- '
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'model anom. ', xl, yl, mszx, mszy, mx, my, '    var.'
  call premsg('')
  call premsg('How to specify the source location')
  call premsg('  0: Coordinates values   or')
  call gparmi('  1: Distance from SW corner ==> ', kl)
  if ((kl < 0) .or. (kl > 1)) call abendm('invalid selection')
  call premsg('Specify Ambient field Direction')
  call gparmf2('  Inclin., Declin. (in degrees) : ', dip, dec)
  if ((dip < -90.) .or. (dip > 90.) .or.   &
      (dec < -180.) .or. (dec >= 360.)) call abendm('invalid value')

  open(9,status='scratch')
  n = 0
  do
    call premsg('Select source model type:')
    call gparmi('  1:Block 2:hRect 3:vLine 4:Point or 0:End ==> ', ks)
    if ((ks < 0) .or. (ks > 4)) call abendm('invalid selection')
    if (ks == 0) exit
    n = n + 1

    call premsg('  Horizontal location:')
    if ((ks == 1) .or. (ks == 2)) then
      call prompt('   Northing Center (cx), N-S size (rx) ')
      call gparmf2('(in km) ==> ', cx, rx)
      if (kl == 1) cx = cx + xl
      pv(1) = cx;  pv(2) = rx
      call prompt('   Easting Center (cy), E-W size (ry)  ')
      call gparmf2('(in km) ==> ', cy, ry)
      if (kl == 1) cy = cy + yl
      pv(3) = cy;  pv(4) = ry
    else
      call gparmf('   Northing (xc) (in km) ==> ', xc)
      if (kl == 1) xc = xc + xl
      pv(1) = xc;  pv(2) = 0.
      call gparmf('   Easting (yc) (in km)  ==> ', yc)
      if (kl == 1) yc = yc + yl
      pv(3) = yc;  pv(4) = 0.
    endif
    call premsg('  Depth below Sea-level:')
    if ((ks == 1) .or. (ks == 3)) then
      call prompt('   Top depth (za), and Thickness (rz)  ')
      call gparmf2('(in km) ==> ', za, rz)
      pv(5) = za;  pv(6) = rz
    else
      call gparmf('    (in km) ==> ', zc)
      pv(5) = zc;  pv(6) = 0.
    endif
    if (ks == 4) then
      call gparmf('  Equivalent source volume (in cubic-km) ==> ', vol)
    else if (ks == 3) then
      call gparmf('  Equivalent source area (in square-km) ==> ', vol)
    else if (ks == 2) then
      call gparmf('  Equivalent source thickness (in km) ==> ', vol)
    else
      vol = 0.
    endif
    pv(7) = vol
    call gparmf('  Magnetization Intensity (A/m) : ', am)
    call premsg('  Magnetization Direction ?')
    call gparmf2('    Inclin., Declin. (in degrees) : ', the, phi)
    if ((the < -90.) .or. (the > 90.) .or.   &
        (phi < -180.) .or. (phi >= 360.)) call abendm('invalid value')
    write(6,'(i4,2x,a5,f8.2,2f8.1)') n, nmdl(ks), am, the, phi
    write(6,'(6x,3(f9.2,f7.2),f9.4)') (pv(i),i=1,7)
    write(9,*) ks, am, the, phi, pv

    call magafd(the, phi, dip, dec)
    if (ks == 4) then
      call mpoint(am, xc, yc, zc, vol)
    else if (ks == 3) then
      call mvline(am, xc, yc, za, rz, vol)
    else if (ks == 2) then
      call mhrect(am, cx, rx, cy, ry, zc, vol)
    else
      call mprism(am, cx, rx, cy, ry, za, rz)
    endif
    do j=1,my; do i=1,mx
      if (h(i,j) == hnul) then
        f(i,j) = 9999.9
      else
        xp = float(ixl + mszx*(i-1)) / 1000.
        yp = float(iyl + mszy*(j-1)) / 1000.
        zp = - h(i,j) / 1000.
        f(i,j) = f(i,j) + calma(xp, yp, zp)
      endif
    enddo; enddo
  enddo
  call clspin()
  write(2,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
    ixl, iyl, mszx, mszy, mx, my, 9999.9, 0.
  do j=1,my
    write(2,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  enddo
  write(2,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
    ixl, iyl, mszx, mszy, mx, my, hnul, -1.
  do j=1,my
    write(2,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  enddo
  close(2)
  endfile(9)
  rewind(9)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'CALMAS : calc. model anomaly on a given surface'
    write(99,'(a,a)') ' Calc. alt. data infile : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' Model anomaly outfile  : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a12,a8,6x,a6,i4,4x,4i8)')   &
      'Areaname :  ', area, 'Coord.', ncc, iorg, korg, ispa, ispb
    write(99,'(4x,a)')   &
      '  (mesh)          xs        ys     msz-x   msz-y    mx    my    alt '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'calc. surf. ', xl, yl, mszx, mszy, mx, my, '    --- '
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'model anom. ', xl, yl, mszx, mszy, mx, my, '    var.'
    write(99,'(4x,a,f6.1,a,f6.1)')   &
      'Ambient field direction:    dip=', dip, '    dec=', dec
    n = 0
    do
      read(9,*,iostat=IER) ks, am, the, phi, pv
      if (IER < 0) exit
      n = n + 1
      write(99,'(i4,2x,a5,f8.2,2f8.1)') n, nmdl(ks), am, the, phi
      write(99,'(6x,3(f9.2,f7.2),f9.4)') (pv(i),i=1,7)
    enddo
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  close(9)
  stop
  end