c  FFT_DECOMP.FOR
c
c  Decompensative correction to isostatic gravity anomaly data, to extend
c  isostatic correction to geologic, as well as topographic, loads.  See
c  Cordell, Zorin and Keller, 1991, The decompensative gravity anomaly and
c  deep structure of the region of the Rio Grande rift: Journal of
c  Geophysical Research v. 96, p. 6557-6568.
c  
c  Calculate the decompensative gravity correction in descrete fft
c  domain, given an existing .cof file from program FFTFIL containing
c  the DFT of the isostatic gravity anomaly, and, as requested parameters: 
c  the assumed compensation depth (in km) and low- and high-
c  wavelength cutoffs (in 1/km). Note that if isostatic anomaly data are 
c  not handily available, a servicable isostatic anomaly can be calculated
c  by modified FFT of the topography, as described in the cited paper.
c
c  Note also that this procedure applies to very large regions, where, for
c  reasons discussed in the paper, the transform of the isostatic anomaly
c  approaches zero in the lowestmost wave numbers.  Normally the DC term is
c  set to zero and, if necessary, the spectrum is dampened for wavelengths
c  larger than a specified wavelength (in km, i.e.) filter threshold.  The
c  analysis might be extended to smaller sized areas,  if special treatment
c  were made of the DC term.  Some prelimimary efforts in this direction are
c  incorporated in the program, but these would not normally be used.
c
c  Note that program F_DECOMP calculates the transform of the decompensative
c  correction (i.e., the negative of the gravity effect of the root), which,
c  when inverse-transformed back to the space domain, is to be added to the
c  isostatic anomaly to give the decompensative anomaly.
c
c  Program modified from FFT_COF.FOR:  
c  In complex frequency domain, modifies fft coefficient file "fftfil.cof"
c  (obtained  from FFTFIL using the icoef option) according to user-supplied
c  code, folds coefficients into lower quadrants  to preserve appropriate
c  symmetries for real-valued functions in time (space) domain, and outputs
c  new "fftfil.cof" file, which can then be inverse-transformed
c  by means of FFTFIL and use of the icoef=-1 option.   This program follows
c  the blocking of 'NEWFFTFIL'.
c
c  Decompensative anomaly calculations: Cordell 26 Feb 89.
c  PC version by Cordell Jan, 1992.
c
c  Fold into lower two quadrants.

c  To low-pass with radial exponential taper.
c
c  In complex frequency domain, modifies fft coefficient file "fftfil.cof"
c  (obtained  from FFTFIL using the icoef option) according to user-supplied
c  code, folds coefficients into lower quadrants  to preserve appropriate
c  symmetries for real-valued functions in time (space) domain, and outputs
c  new "fftfil.ncf" file, which can be renamed "fftfil.cof" and then
c  inverse-transformed by means of FFTFIL and use of the icoef=-1 option.
c  This program follows the blocking of 'NEWFFTFIL'.
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, 9 Dec '87
c  Repairs 1 Feb 91.
c  Revised for PC 21 Aug '90.
c
      character title*56,pgm*8
      character*50 ingrd
      dimension g(2,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. nx+2*nadd.
c  m=no. of rows, l=no. from 9-16, k=interger.
      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
c  Blocking parameters for fft are set; call subroutine.
      n22=2*n2
      open(14,access='direct',status='old',
     1 form='unformatted',file='fftfil.cof',
     2 recl=n22*4)
      open(12,access='direct',status='new',
     1 form='unformatted',file='fftfil.ncf',
     2 recl=n22*4)
      call modify(g,n2,n1,n22,dy,dx)
      print*,' Output in file FFTFIL.NCF'
      print*,' Dont forget to rename this FFTFIL.COF before rerunning'
      print*,' fftfil.'
      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,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)
      ncp1=n2+1
      cn=1./(n2*dy)
      rn=1./(n1*dx)
      scale=float(n1)*float(n2)
      nrnq=float(n1)/2.+1.0000001
      ncnq=float(n2)/2.+1.0000001
      pi=3.1415927e0
c  If necessary, SUPPLY parameter read here.
c
10      write(6,11)
11      format(1x,'Enter compensation depth:'/' *'$)
      read(5,*)z
      if(z.le.0.0)go to 900
      c=2.0*pi*z
      write(6,12)
12    format(' Enter wavelength filter threshold'
     $' (wl2 >= 0, in km):'/' *'$)
      read(5,*)wl2
      if(wl2.lt.0.0)go to 900
      w1=0.0
      w2=1.e-10
      if(wl2.ne.0.0) w2=1./wl2
      if(w2.gt.ncnq*.5)go to 900
c  Read fftfil.coef file and modify.
c   First row, DC term treated separately.
      jr=1
      read(14,rec=jr)g
c   DC term.
      if(wl2.eq.0)then
        print*, 'Enter G(0) >= 0 (normally 0):'
        read*, g(1)
      else
        g(1)=0.0
      endif
      print*,'The first few terms of |G(w)|:'
      print*,'     w           |G(w)|        log|G(w)|'
      w=0
      ag=g(1)
      if(ag.ne.0)then
        agl=log(abs(ag))
        write(5,50)w,ag,agl
      else
        write(5,51)w,ag
      endif
50    format(1x,f8.5,2e16.5)
51    format(1x,f8.5,e16.5)
      wlim=2.0*(max(cn,rn))
c
c   Rest of row 1.
      do 1 i=2,n2
      ic=2*i
      ii=i-1
      if(i.gt.ncnq) ii=-(ncp1-i)
      v2=float(ii)*cn
c    (If necessary, here consider effect of sign of v2.)
      Gr=g(ic-1)
      Gi=g(ic)
      w=abs(v2)
      tf=1.0/(exp(c*w)-1.0)
      if(w.lt.w1)f=0.0
      if(w.gt.w1)f=(w-w1)/(w2-w1)
      if(w.gt.w2)f=1.0
      g(ic-1)=Gr*f*tf
      g(ic)=Gi*f*tf
      if(w.le.wlim)then
        ag=sqrt(g(ic-1)**2+g(ic)**2)
        agl=log(ag)
        write(5,50)w,ag,agl
      endif
1     continue
      write(12,rec=jr)g
c
c  Read rest of fftfil.cof file, row-wise, and modify.
      do 3 jr=2,nrnq
      jj=jr-1
      v1=float(jj)*rn
      v1sq=v1*v1
      read(14,rec=jr)g
      do 2 i=1,n2
      ic=2*i
      Gr=g(ic-1)
      Gi=g(ic)
      ii=i-1
      if(i.gt.ncnq) ii=-(ncp1-i)
      v2=float(ii)*cn
      v2sq=v2*v2
      w=sqrt(v1sq+v2sq)
c  (If necessary, consider effect of sgn of v1, v2.)
      tf=1.0/(exp(c*w)-1.0)
      if(w.lt.w1)f=0.0
      if(w.gt.w1)f=(w-w1)/(w2-w1)
      if(w.gt.w2)f=1.0
      g(ic-1)=Gr*f*tf
      g(ic)=Gi*f*tf
      if(w.le.wlim)then
        ag=sqrt(g(ic-1)**2+g(ic)**2)
        agl=log(ag)
        write(5,50)w,ag,agl
      endif
2     continue
      write(12,rec=jr)g
3     continue
c
c  Fold into lower two quadrants.
      nrp2=n1+2
      ncp2=n2+2
      do 5 i=2,nrnq-1
      ir=i
      read(12,rec=ir)g
      ik=nrp2-i
      g(2)=-g(2)
      do 4 j=2,ncnq
      jc=2*j
      jk=ncp2-j
      jkc=2*jk
      gr=g(jc-1)
      gi=-g(jc)
      g(jc-1)=g(jkc-1)
      g(jc)=-g(jkc)
      g(jkc)=gi
      g(jkc-1)=gr
4     continue
      write(12,rec=ik)g
5     continue
      go to 999
900   print*, 'Cant handle it.'
      stop
999   continue
      return
      end
