!--------( CMAGF : calc. anomaly of <amag> fine model )--------
  implicit none
  integer :: lwkdir, lrtrim;  real :: calma
  real,allocatable :: f(:,:), h(:,:), s(:,:), t(:,:)
  character(80) :: wdr, flog, fnam1, fnam2, fnam3, fnam4, buf
  character(72) :: out;  character(20) :: strs, sdtm
  character(8) :: area, area1, area2
  integer :: ldr, i, j, ir, jr, ki, kj, nmsx, nmsy, is, js, it, jt
  integer :: nc, ncc, iorg, korg, ispa, ispb
  integer :: nc1, ncc1, iorg1, korg1, ispa1, ispb1
  integer :: ixs, iys, mszx, mszy, mx, my, ixt, iyt
  integer :: ixl, iyl, mszxm, mszym,  mxm, mym, ixh, iyh
  integer :: ixst, iyst, mszxt, mszyt, mxt, myt, ixtt, iytt
  integer :: klp, nhx, nhy, ks, ixc, iyc, ia, ja, kic, kjc
  real :: hnul, xs, ys, smx, smy, anul, dmy, xl, yl, smxm, smym, tnul
  real :: xst, yst, smxt, smyt, hw, hwm, botm, dip, dec, the, phi, am
  real :: ht, x, y, hc, sz, za, hnmsx, hnmsy

  call premsg('CMAGF : calc. anomaly of <amag> fine model')
  call opnpin()
  ldr = lwkdir(50,wdr) + 1
  flog = wdr;  fnam1 = wdr;  fnam2 = wdr;  fnam3 = wdr;  fnam4 = wdr
  call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80))
  call gparma('Enter Calc.Alt. data 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,*) ixs, iys, mszx, mszy, mx, my, hnul
  if ((mx <= 0) .or. (my <= 0)) call abendm('illegal array size')

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

  xs = float(ixs)/1000.;  ys = float(iys)/1000.
  ixt = ixs + mszx*(mx-1);  iyt = iys + mszy*(my-1)
  read(1,*) ((h(i,j),i=1,mx),j=1,my)
  close(1)
  smx = float(mszx);  smy = float(mszy)

  call gparma('Enter AMAG model input filename ==> ', 81-ldr, fnam2(ldr:80))
  open(2,file=fnam2,status='old')
  do
    read(2,'(a)') buf
    if (buf(1:1) /= '#') exit
  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.')
    nc1 = ncc1
    if (ncc1 < 200) then
      nc1 = ncc1; ncc1 = nc1 + 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 = (nc1-30)*360 - 180
    if (nc1 > 60) korg1 = 0
    if (nc1 == 61) iorg1 = +5400
    if (nc1 == 62) iorg1 = -5400
  endif
  read(2,'(a)') buf
  read(buf,*) ixl, iyl, mszxm, mszym, mxm, mym, anul, dmy, klp
  if ((ncc1 /= ncc) .or. (iorg1 /= iorg) .or. (korg1 /= korg)   &
      .or. (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
    call abendm('Header inconsistent')
  if ((mxm < 10) .or. (mym < 10)) call abendm('illegal array size')

  allocate(s(mxm,mym))

  xl = float(ixl)/1000.;  yl = float(iyl)/1000.
  ixh = ixl + mszxm*(mxm-1);  iyh = iyl + mszym*(mym-1)
  if ((ixl > ixs) .or. (iyl > iys) .or.   &
      (ixh < ixt) .or. (iyh < iyt)) call abendm('ranges conflict')
  read(2,*) ((s(i,j),i=1,mxm),j=1,mym)
  close(2)
  smxm = float(mszxm);  smym = float(mszym)

  call gparma('Enter Source Alt. data filename  ==> ', 81-ldr, fnam3(ldr:80))
  open(10,file=fnam3,status='old')
  do
    read(10,'(a)') buf
    if (buf(1:1) /= '#') exit
  enddo
  if (buf(13:16) == '    ') then        !-- v2005
    read(buf,'(a8,i4,4x,4i8)') area2, ncc1, iorg1,korg1, ispa1,ispb1
    if ((ncc1<0) .or. (ncc1>=400)) call abendm('invalid coord.no.')
    nc1 = ncc1
    if (ncc1 < 200) then
      nc1 = ncc1; ncc1 = nc1 + 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 = (nc1-30)*360 - 180
    if (nc1 > 60) korg1 = 0
    if (nc1 == 61) iorg1 = +5400
    if (nc1 == 62) iorg1 = -5400
  endif
  read(10,'(a)') buf
  read(buf,*) ixst, iyst, mszxt, mszyt, mxt, myt, tnul
  if ((ncc1 /= ncc) .or. (iorg1 /= iorg) .or. (korg1 /= korg)   &
      .or. (ispa1 /= ispa) .or. (ispb1 /= ispb))   &
    call abendm('Header inconsistent')
  if ((mxt < 10) .or. (myt < 10)) call abendm('illegal array size')

  allocate(t(mxt,myt))

  nmsx = mszxm / mszxt
  if ((nmsx*mszxt) /= mszxm) call abendm('mesh size illegal')
  nmsy = mszym / mszyt
  if ((nmsy*mszyt) /= mszym) call abendm('mesh size illegal')
  hnmsx = float(nmsx-1)/2.;  hnmsy = float(nmsy-1)/2.
  xst = float(ixst)/1000.;  yst = float(iyst)/1000.
  ixtt = ixst + mszxt*(mxt-1);  iytt = iyst + mszyt*(myt-1)
  if ((ixst > ixl) .or. (iyst > iyl) .or.   &
      (ixtt < ixh) .or. (iytt < iyh)) call abendm('ranges conflict')
  read(10,*) ((t(i,j),i=1,mxt),j=1,myt)
  close(10)
  smxt = float(mszxt);  smyt = float(mszyt)

  write(out,'(a12,a8,6x,a6,i4,4x,4i8)')   &
    '--- Area :  ', area, 'Coord.', ncc, iorg,korg, ispa,ispb
  call premsg(out)
  write(out,'(a14,2a10,2x,2a8,2a6)')   &
    '---   (mesh)  ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my'
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6)')   &
    '--- Calc. Alt.  ', xs, ys, mszx,mszy, mx,my
  call premsg(out)
  write(out,'(a16,2f10.2,2i8,2i6)')   &
    '--- AMAG model  ', xl, yl, mszxm,mszym, mxm,mym
  write(out,'(a16,2f10.2,2i8,2i6)')   &
    '--- Source Alt. ', xst, yst, mszxt,mszyt, mxt,myt
  call premsg(out)
  call premsg('  Truncation of Source effect ?')
  call gparmf('    Half width of calc.window (in km) : ', hw)
  hwm = hw * 1000.;  nhx = nint(hwm/smx);  nhy = nint(hwm/smy)

  call premsg('')
  call premsg('Specify Source-Bottom configuration')
  call gparmi('    Constant Thickness (1) or Flat Bottom (0) : ', ks)
  if (ks == 0) then;  strs = ' below sea level    '
  else if (ks == 1) then; strs = ' below terr.surface '
  else;  call abendm('invalid value')
  endif
  call gparmf('  Bottom-depth (in m)' // strs // ': ', botm)
  call premsg('Specify Direction Parameters (in degrees)')
  call gparmf2('  Inclin., Declin. of Ambient field : ', dip, dec)
  if ((dip < -90.) .or. (dip > 90.) .or.   &
      (dec < -180.) .or. (dec >= 360.)) call abendm('invalid value')
  call gparmf2('  Inclin., Declin. of Magnetization : ', the, phi)
  if ((the < -90.) .or. (the > 90.) .or.   &
      (phi < -180.) .or. (phi >= 360.)) call abendm('invalid value')

  call gparma('Enter Calc.Anom. output filename ==> ', 81-ldr, fnam4(ldr:80))
  call clspin()
  call strdtm(sdtm)
  open(3,file=fnam4,status='new')
  write(6,'(a)')   '----------------------------------------'
  write(6,'(a)')   'CMAGF : calc. anomaly of <amag> fine model'
  write(6,'(a,a)') ' Calc.Alt. data infile : ', fnam1(1:lrtrim(fnam1))
  write(6,'(a,a)') ' AMAG model infile     : ', fnam2(1:lrtrim(fnam2))
  write(6,'(a,a)') ' Source Alt. infile    : ', fnam3(1:lrtrim(fnam3))
  write(6,'(a,a)') ' Calc.Anom. outfile    : ', fnam4(1:lrtrim(fnam4))
  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'
  write(6,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
    'Calc. Alt.  ', xs, ys, mszx, mszy, mx, my
  write(6,'(4x,a12,2f10.2,2i8,2i6)')   &
    'AMAG model  ', xl, yl, mszxm, mszym, mxm, mym
  write(6,'(4x,a12,2f10.2,2i8,2i6)')   &
    'Source Alt  ', xst, yst, mszxt, mszyt, mxt, myt
  write(6,'(a,f6.2,a,i4,a,i4,a)')   &
    ' Half width of calc.window :', hw, ' km (', nhx, ' *', nhy, '  calc.mesh)'
  write(6,'(3x,a,f8.1,a,a20)') 'bottom depth  : ', botm, ' m', strs
  write(6,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
    'ambient field : ', 'dip=', dip, 'dec=', dec
  write(6,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
    'magnetization : ', 'dip=', the, 'dec=', phi

  do j=1,my; do i=1,mx
    f(i,j) = 0.;  if (h(i,j) == hnul) f(i,j) = 9999.9
  enddo; enddo
  call dpcini('... calc. model anomaly   : ')
  call magafd(the, phi, dip, dec)
  am = smxt * smyt
  do j=1,mym
    iyc = iyl + mszym*(j-1)
    ja = nint(float(iyc-iyst)/smyt - hnmsy)
    kjc = nint(float(iyc-iys)/smy)
    do i=1,mxm
      ixc = ixl + mszxm*(i-1)
      ia = nint(float(ixc-ixst)/smxt - hnmsx)
      kic = nint(float(ixc-ixs)/smx)
      do js=1,nmsy
        jt = ja + js
        if ((jt <= 0) .or. (jt > myt)) cycle
        iyt = iyst + (jt-1)*mszyt
        do is=1,nmsx
          it = ia + is
          if ((it <= 0) .or. (it > mxt)) cycle
          ixt = ixst + (it-1)*mszxt
          ht = t(ia+1,ja+1)
          if (ht == tnul) call abendm('source alt. undefined')
          do kj=-nhy,nhy
            jr = kjc + kj
            if ((jr < 0) .or. (jr >= my)) cycle
            y = float(iyt - iys - jr*mszy)
            do ki=-nhx,nhx
              ir = kic + ki
              x = float(ixt - ixs - ir*mszx)
              if ((ir < 0) .or. (ir >= mx)) cycle
              hc = h(ir+1,jr+1)
              if (hc == hnul) cycle
              za = hc - ht;  sz = botm
              if (ks == 0) sz = ht + botm
              if ((iabs(ki) <= 10) .and. (iabs(kj) <= 10)) then
                call mprism(s(i,j)/100., x, smxt, y, smyt, za, sz)
              else
                call mvline(s(i,j)/100., x, y, za, sz, am)
              endif
              f(ir+1,jr+1) = f(ir+1,jr+1) + calma(0.,0.,0.)
            enddo
          enddo
        enddo
      enddo
      call dpcent((j-1)*mxm+i, mxm*mym)
    enddo
  enddo

  write(3,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
     ixs, iys, mszx, mszy, mx, my, 9999.9, 0.
  do j=1,my
    write(3,'((f7.1,9(1x,f7.1)))') (f(i,j),i=1,mx)
  enddo
  write(3,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb
  write(3,'(2i12,4i6,1x,f7.1,1x,f7.0)')   &
     ixs, iys, mszx, mszy, mx, my, hnul, -1.
  do j=1,my
    write(3,'((f7.1,9(1x,f7.1)))') (h(i,j),i=1,mx)
  enddo
  close(3)

  if (flog(ldr:ldr) /= ' ') then
    open(99,file=flog,access='append')
    write(99,'(//a)') 'CMAGF : calc. anomaly of <amag> fine model'
    write(99,'(a,a)') ' Calc.Alt. data infile : ', fnam1(1:lrtrim(fnam1))
    write(99,'(a,a)') ' AMAG model infile     : ', fnam2(1:lrtrim(fnam2))
    write(99,'(a,a)') ' Source Alt. infile    : ', fnam3(1:lrtrim(fnam3))
    write(99,'(a,a)') ' Calc.Anom. outfile    : ', fnam4(1:lrtrim(fnam4))
    write(99,'(a,a,a)') '--------------- ', sdtm, '---------------'
    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'
    write(99,'(4x,a12,2f10.2,2i8,2i6,a8)')   &
      'Calc. Alt.  ', xs, ys, mszx, mszy, mx, my
    write(99,'(4x,a12,2f10.2,2i8,2i6)')   &
      'AMAG model  ', xl, yl, mszxm, mszym, mxm, mym
    write(99,'(4x,a12,2f10.2,2i8,2i6)')   &
      'Source Alt  ', xst, yst, mszxt, mszyt, mxt, myt
    write(99,'(a,f6.2,a,i4,a,i4,a)') ' Half width of calc.window :',   &
      hw, ' km (', nhx, ' *', nhy, '  calc.mesh)'
    write(99,'(3x,a,f8.1,a,a20)') 'bottom depth  : ', botm, ' m', strs
    write(99,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
      'ambient field : ', 'dip=', dip, 'dec=', dec
    write(99,'(3x,a,4x,a,f6.1,4x,a,f7.1)')   &
      'magnetization : ', 'dip=', the, 'dec=', phi
    call strdtm(sdtm)
    write(99,'(a,a,a)') '=============== ', sdtm, '==============='
    close(99)
  endif
  stop
  end