c   XIAM_SAV
c
c   attempts to recover a magnetization file resulting from premature
c   termination of XIAM.
c
c   Jeff Phillips Oct 1994
c
      dimension work(4096)
      character*80 ifile1,ifile2,ofile
      character*80 filet,fileo,filem
c      character*80 ofilet,ofileo
c      character*1 ans,drive
      character*1 ans
      character*56 id
c
      common/mane/id,pgm(2),dx,dy,nc,nr,nz,xo,yo
      common/other/sincl,sdecl,fincl,fdecl,nstop,errn,mstop,errm
c      common/files/drive,filet,fileo,ans,filem
      common/files/filet,fileo,ans,filem
      common/surf/smed,smin,smax,elevel
c
c      data eps/1.e-7/
      data dval/1.e37/
c
      open(12,file='xia.rec',status='old',form='formatted')
c get the magnetic anomaly
c      print'(a\)',' Standard file containing the observed anomaly: '
      read(12,1) ifile1
    1 format(a)
      open(10,file=ifile1,form='unformatted',status='old')
        read(10)id,pgm,nc,nr,nz,xo,dx,yo,dy
c        print*,' title = ',id
        nctest=nc
        nrtest=nr
c test for dvals, get mean
        tm=0.
        do 20 j=1,nr
        call rowio(work,yr,nz,nc,10,1)
        do 20 i=1,nc
        if(work(i).ge.dval) stop 'dvals found'
        tm=tm+work(i)
   20   continue
        tm=tm/float(nc*nr)
c
c      print'(a\)',' Standard file containing the observation surface: '
      read(12,1) ifile2
        open(11,file=ifile2,status='old',form='unformatted')
        read(11)id,pgm,nc,nr,nz,xo,dx,yo,dy
c        print*,' title = ',id
c test for registration with observed anomaly
          if(nc.ne.nctest.or.nr.ne.nrtest) stop 'grids do not overlay'
c get smin,smax and test for dvals
        smin=1.e37
        smax=-1.e37
        do 21 j=1,nr
        call rowio(work,yr,nz,nc,11,1)
        do 21 i=1,nc
        if(work(i).ge.dval) stop 'dvals found'
        if(work(i).lt.smin) smin=work(i)
        if(work(i).gt.smax) smax=work(i)
   21   continue
c get scale factor for observation surface
c        print*,'grid dx/dy = ',dx,dy
c        print*,'surface min/max = ',smin,smax
c        print*,'What single number should these values be multiplied by'
c        print*,'  so that they increase vertically down and have the'
c        print'(a\)','   same units as dx?:'
        read(12,*) sscale
        if(sscale.gt.0.0) then
          smin=smin*sscale
          smax=smax*sscale
        else
          temp=smin
          smin=smax*sscale
          smax=temp*sscale
        endif
        sscale=sscale/dx
c        print*,'the bounds of the observation surface are:'
c        print*,'        upper = ',smin
c        print*,'        lower = ',smax
c        smed=(smin+smax)/2.
        smed=(smin+smax)/(2.*dx)
c        print'(a\)',' Elevation of equivalent source layer: '
        read(12,*) elevel
        elevel=elevel/dx
c get dimensions of expanded magnetization grid
        read(12,*) ncol,nrow
c        call check_dims(nc,nr,ncol,nrow)
c
      call get_dims(nrow,ncol,nx,ny,nri,id2,nxa)
c
c get output file and set up calculations
c      print'(a\)',' Standard file to receive magnetization:'
      read(12,1) ofile
      open(13,file=ofile,form='unformatted',status='unknown')
      pgm='xiam_sav'
c      print'(a\)',' Title: '
      read(12,1) id
      write(13) id,pgm,ncol,nrow,nz,xo,dx,yo,dy,tm,smed,elevel,nc,nr
c
c      print'(a\)',' Enter the (RAM) drive letter for temporary files: '
c      read(5,1) drive
       close(12)
c
      call fmag3dsav(work,nrow,ncol,id2)
c
      stop
      end
c **************************************************************
      subroutine rowio(a,y,nz,ncol,ld,key)
      dimension a(nz,ncol)
      go to (1,2),key
1     read(ld) y,a
      go to 90
2     write(ld)y,a
90    return
      end
c **************************************************************
      subroutine check_dims(nc,nr,ncol,nrow)
        print*,nc,nr
10      call fftdims(nc,nr)
        print*, 'Enter (new) ncol and nrow:'
        read*, ncol,nrow
c  Check to make sure SFFTMG will accept these.
c  Added 19 June 91.
c        kr=11
c        iw=6
c        ir=5
c        nadd=0
c  note that x & y have been switched from normal (usgs) grid specficati
c  ...convention used is x-north y-east and z-down
        n1=nrow
        n2=ncol
        call get_dims(n1,n2,nx,ny,nri,id2,nxa)
        itest=n1-nx
        if(itest.ne.0) then
          print*,'Select another nrow; SFFTMG cannot handle this one.'
          go to 10
        endif
        if(nc.gt.ncol.or.nr.gt.nrow)then
          print*,'New ncol, nrow are too small.'
          go to 10
        endif
        return
        end
c **************************************************************
      subroutine fftdims(ncol, nrow)
c
c     suggests new dimensions for fourier transform
c
c     Jeff Phillips 1992
c
      dimension nlist(54), ngood(54)
      data nlist / 9, 10, 11, 12, 13, 16
     &, 18, 20, 22, 24, 26, 32
     &, 36, 40, 44, 48, 52, 56, 64
     &, 72, 80, 88, 96, 104, 112, 128
     &, 144, 160, 176, 192, 208, 224, 256
     &, 288, 320, 352, 384, 416, 448, 512
     &, 576, 640, 704, 768, 832, 896, 1024
     &, 1152, 1280, 1408, 1536, 1664, 1792, 2048 /
c
      write(unit=6, fmt=901) ncol, nrow
  901 format(' ncol = ',i5,' nrow = ',i5)
      n = 0
      do 10 i = 1, 54
      if (nlist(i) .lt. ncol) goto 10
      if (nlist(i) .gt. ncol+ncol/2) goto 20
      n = n + 1
      ngood(n) = nlist(i)
   10 continue
   20 write(unit=6, fmt=902)
  902 format(' Suggested new values:')
      write(unit=6, fmt=101) (ngood(i),i = 1, n)
  101 format('  ncol = ',14i5)
      n = 0
      do 30 i = 1, 54
      if (nlist(i) .lt. nrow) goto 30
      if (nlist(i) .gt. nrow+nrow/2) goto 40
      n = n + 1
      ngood(n) = nlist(i)
   30 continue
   40   write(unit=6, fmt=102) (ngood(i),i = 1, n)
  102 format('  nrow = ',14i5)
      return 
      end
c***********************************************************************
      subroutine get_dims(n1,n2,nx,ny,nri,id2,nxa)
      data lnri/16/,iw/6/
c
c  nri = blocking factor for rows (l)
c  nxa = number of rows added (n1-nx)
c  id2 = 2*max(n1,n2)
c
      nx=n1
      ny=n2
      l=lnri+1
      lnri21=lnri/2+1
c  set no. of rows for fft: need m=l*2**k, m.gt. or .eq. nx.
c  m=no. of rows, l=no. from 9-16, k=integer.
  190 l=l-1
      if(l.lt.lnri21) go to 260
      mr=n1/l+0.0000001
      k=1
      idiv=2
  200 iquot=mr/idiv+0.0000001
      if(iquot.lt.idiv) go to 210
      k=k+1
      mr=iquot
      go to 200
  210 k=k+1
      m=l*2**k
  220 mtest=l*2**(k-1)
      if(mtest.lt.n1) go to 230
      k=k-1
      m=mtest
      if(k.eq.0) go to 230
      go to 220
  230 lnxa=m-n1
      if(l.ne.lnri) go to 250
      nxa16=lnxa
  240 nri=l
      nxa=lnxa
      go to 190
  250 if(lnxa.ge.nxa) go to 190
      go to 240
  260 n1=n1+nxa
c  check to see if row block size of 16 will be more efficient
      n116=n1-nxa+nxa16
      ntest=0.9*n116
      if(ntest.gt.n1.or.n116.gt.2048) go to 270
      n1=n1-nxa+nxa16
      nxa=nxa16
      nri=16
  270 continue
      if(n1.gt.2048.or.n2.gt.2048) go to 320
      id2=n1
      if(n2.gt.n1)id2=n2
      id2=2*id2
      return
  320 write(iw,330)nx,ny,n1,n2
  330 format(' #no. of extended rows or columns exceeds 2048:'/
     1' input no. of rows and columns='2i4/
     2' no. of rows and columns required for fft =',2i4,/)
      stop
      end
c **************************************************************
      subroutine fmag3dsav(work[huge],nrow,ncol,id2)
      dimension work(id2)
c      dimension fact(31)
c      real mx,my,mz,maxd
      character*80 fileo,filet,filem
c      character*1 drive,ans
      character*1 ans
c
      common/mane/id(14),pgm(2),dx,dy,nc,nr,nz,xo,yo
      common/other/sincl,sdecl,fincl,fdecl,nstop,errn,mstop,errm
c      common/files/drive,filet,fileo,ans,filem
      common/files/filet,fileo,ans,filem
      common/surf/smed,smin,smax,elevel
      data iunitm/26/,dum/0.0/
cc
        open(iunitm,access='direct',status='old',
c     1  form='unformatted',file=drive//':\slave5.tmp',recl=ncol*4)
     1  form='unformatted',file='slave5.tmp',recl=ncol*4)
      do 2220 jr=1,nrow
      call srdda(iunitm,jr,work,ncol)
      write(13) dum,(work(i),i=1,ncol)
 2220 continue
      return
      end
c***********************************************************************
       subroutine srdda(no,ipos,dat[huge],n)
c       subroutine 'srrda' reads keyed sequential files.
      dimension dat(n)
      read(no,rec=ipos)dat
      return
      end
