c  F_DX
c
c  Horizontal derivative, d f(x,y)/dx, of a standard-file grid.
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, 25 September 1992
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)
      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
      tpi=2.*pi
c
c  Read fftfil.coef file and modify.
c   First row, DC term treated separately.
      jr=1
      read(14,rec=jr)g
c   DC term.
      g(1)=0.
      g(2)=0.
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
      t=tpi*v2
c    (If necessary, here consider effect of sign of v2.)
      Gr=g(ic-1)
      Gi=g(ic)
      g(ic-1)=-t*Gi
      g(ic)=t*Gr
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
      read(14,rec=jr)g
      do 2 i=1,n2
      ic=2*i
      ii=i-1
      if(i.gt.ncnq) ii=-(ncp1-i)
      v2=float(ii)*cn
      t=tpi*v2
c  (If necessary, consider effect of sgn of v1, v2.)
      Gr=g(ic-1)
      Gi=g(ic)
      g(ic-1)=-t*Gi
      g(ic)=t*Gr
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
      return
      end
