c  F_A_AMP.FOR
c
c  To make grid of top two quads of area-weighted Fourier amplitude
c  spectrum, i.e., log(w|F(w)|).
c
c  See comments in subroutine MODIFY regarding switch of coordinate
c  axes from x=east, y=north; to x=north, y=east.
c
c  Lin Cordell, 17 Jan 90.
c  Modified for pc version of fftfil 20 Aug 90
c
      character title*56
      character*50 ingrd,ofile
      dimension g(2,2048),pgm(2),s(2048)
      data lnri/16/
      ierror=0
      write(6,10)
10    format(1x,'Enter grid from which Fourier coeff.s created:'/1x'*'$)
      read(5,100)ingrd
100   format(a50)
      open(10,file=ingrd,status='old',form='unformatted')
      read(10)title,pgm,n2,n1,nz,yo,dy,xo,dx
c  Note: n2,n1 and dx,dy rotate sense of rows and cols. Now x=north,
c   y=east.
      close(10)
      write(6,7)
7     format(' Enter no. of rows & cols added to grid(nadd):'/1x'*'$)
      read(5,*)nadd
      icmd=10
      l=lnri+1
      lnri21=lnri/2+1.
c  set no. of rows for fft: need m=l*2**k, m.gt. or .eq. nrows+2*nadd.
c  m=no. of rows, l=no. from 9-16, k=integer.
      n1=n1+2*nadd
  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 n2=n2+2*nadd
      write(6,275)nri
  275 format(' blocking rows =',i3)
      if(n1.gt.2048.or.n2.gt.2048) go to 320
      nxa=nxa+2*nadd
      id2=n1
      if(n2.gt.n1)id2=n2
      id2=2*id2
      n22=2*n2
      open(14,access='direct',status='old',
     1 form='unformatted',file='fftfil.cof',
     2 recl=n22*4)
      write(6,11)
11    format(1x,'Enter output |FFT| grid file name:'/1x'*'$)
      read(5,100)ofile
      nrnq=float(n1)/2.+1.0000001
      open(11,file=ofile,status='new',form='unformatted')
      dyy=1./(n2*dy)
      dxx=1./(n1*dx)
      yo=-1./(2*dy)
      xo=0.
      write(11)title,pgm,n2,nrnq,nz,yo,dyy,xo,dxx
      call modify(g,s,n2,n1,n22,dy,dx)
      go to 99
320   write(6,321)
321   format(1x,'Extended rows/columns exceeds 2048; run aborted.')
99    continue
      close(14)
      close(12)
      stop
      end
      subroutine modify(g,s,n2,n1,n22,dy,dx)
c  Note:  In contrast with standard-grid notation, the x,dx,v1
c  direction is north; y,dy,v2 direction is east, and angles are
c  measured clockwise from the +x (north) direction.
      dimension g(n22),s(n2)
      data dval/0.170412e+39/
      ncp1=n2+1
      cn=1./(n2*dy)
      rn=1./(n1*dx)
      rave=0.5*(rn+cn)
      scale=float(n1)*float(n2)
      scalesq=scale*scale
      nrnq=float(n1)/2.+1.0000001
      ncnq=float(n2)/2.+1.0000001
      pi=3.1415927e0
      tpi=2.*pi
c  Note that in FFTFIL nrow and ncol are now even. 
c
c  Read fftfil.cof file row-wise.
      do 3 jr=1,nrnq
      v1=jr*rave
      read(14,rec=jr)g
      do 2 i=1,ncnq
      ic=2*i
      l=i+ncnq-3
      lc=2*l
      arg=v1*(g(lc-1)**2+g(lc)**2)/scalesq
      amp=dval
      if(arg.gt.0.)amp=0.5*log(arg)
      s(i)=amp
      arg=v1*(g(ic-1)**2+g(ic)**2)/scalesq
      amp=dval
      if(arg.gt.0.)amp=0.5*log(arg)
      s(i+ncnq-2)=amp
2     continue
      write(11)v1,(s(i),i=1,n2)
3     continue
      return
      end
