program filter_img
c---5----0----5----0----5----0----5----0----5----0----5----0----5----072
c
c Program to filter an img file by reading it in a strip at a time, and
c FFTing each strip.  Strips are overlapped by half their width and 
c cosine blended.  N and S edges are not obscured by tapering because 
c mirrors are used for first and last strip.
c
c W. H. F. Smith, 24 June 1996
c
C*- Modified by W H F Smith 27 June 1996:  More explanatory comments
C*- about consequences of "nyfft" added, and value changed from 576
C*- (=6336/11) to 704 (=6336/9).
C*- Also added stuff to remove/restore local/global means according
C*- to whether it is lowpass or hipass.
C*- Also changed filt.f to newfilt.f to cut some multiplies out of
C*- a loop.  Added storage for xksq.
c
c These parameters must be set at compile time:
c
c nximg - the number of columns in the imgfile;  for 2' img use 10800 
c nyimg - the number of rows in the imgfile;     for 2' img use 6336
c nyfft - the number of rows in the FFT;         for 2' img use 576
c
c nyimg / nyfft must divide evenly with no remainder.  The program will
c use an array of (nximg/2+1)*nyfft*8 bytes for a real/complex array,
c plus nximg*nyfft*2 bytes for a blending array, plus some other small 
c arrays.
C-* If nximg=10800 and nyfft=576, then 37.3 Mb are used for these arrays.
C-* If nximg=10800 and nyfft=704, then 45.6 Mb are used for these arrays.
C-*
C-* The shortest distance is the y-distance, and it is smallest in the
C-* first and last strips, where the scale is set for 72 degrees latitude.
C-* At this latitude, each pixel is 1.145 km across.  Thus the longest
C-* wavelength resolved in the FFT will be 1.145 * nyfft km, or:
C-* if nyfft=576, this is 660 km, or harmonic degree 60;
C-* if nyfft=704, this is 806 km, or harmonic degree 50.
C-* This isn't well-resolved; it is simply the 1 cycle wavelength in
C-* the y direction.  It isn't fully resolved because of the cosine window.
C-*
C-* Making nyfft reduces the number of FFTs which have to be done, and
C-* the number of strips which have to be read in and out, resulting in
C-* some increase in speed; however, there is a trade-off with increase
C-* in possible scale distortion.  The % change in scale from one solution
C-* to the next is largest at the high latitudes, and is about 15% when
C-* nyfft=576; it is 19 % when nyfft=704.
C-*
C-* I initially had set nyfft=576 but I had some trouble at the edges
C-* of solutions due to unkilled long-wavelengths in the solutions.
C-* I noticed that in paper 2 we had used a cosine taper between degrees
C-* 50-70 to kill long wavelengths (and an explicit kill), and so I
C-* decided to try nyfft=704 so that degree 50 would be present.
c
c  03/24/07 - Changed by dts to work with new grid that extends to +/- 81 deg.
c  The dimensions at 2-min are 10800 by 6840.  Note that 8640 is divisible
c  by 15 which provides strips of height 576.  Although this results
c  in a strip of height 551 km at 75 deg north.  We don't need long
c  wavelengths for the prediction so this is OK.
c
      implicit none
      integer*4 nximg, nyimg, nyfft
c     parameter (nximg=10800, nyimg=8640, nyfft=576)
      parameter (nximg=21600, nyimg=17280, nyfft=1152)
c
      complex savek(nximg/2+1, nyfft), work(nximg)
      real*4 savex(nximg, nyfft), blend(nximg, nyfft/2), window(nyfft)
      equivalence (savex(1,1), savek(1,1))
c
c nyfft2 is used to save the half-width of nyfft; nyfft2 = nyfft/2.
c nstrips is the number of ffts to do; nstrips = 2*(nyimg/nyfft) + 1.
c jrowseek is the number of rows into the imgfile to skip before
c loading the next nyfft rows into the savex array.  It works like
c a C-language array counter convention, being 0 when no seek is done.
c jlatscl is likewise the C-index to the row at which the latitude
c scale is currently being set.
c istrip, i, j are loop counters.
c
      integer*4 nstrips, nyfft2, istrip, i, j, jrowseek, jlatscl
c
      real*8 drlat, pi, winfact
      real*4 xwide, ytall
      character*80 ifilnam, ofilnam
      real params(13), outmul
      integer iwant(6), iouttype
c
      logical debug
c
C-*
      integer*4 iaddgm, ikillm
      real*8 glmean, xksq(nximg/2+1)
C-*
c
	debug = .false.
c
c
c Check that nyimg and nyfft parameters are compatible, and if so,
c initialize nstrips, nyfft2, and window array:
c
      istrip = nyimg/nyfft
      if (istrip * nyfft .ne. nyimg) then
        write(*, 101)
 101    format ('filter_img:  ERROR.  Incompatible parameters.')
        stop
      end if
c
      nstrips = 2 * istrip + 1
      nyfft2 = nyfft / 2
c
      pi = 4.0 * datan(1.0d+00)
      winfact = pi/nyfft2
      do 201 j = 1, nyfft
c       window(j)=float(0.5 * (1.0 + dcos((j - nyfft2 - 0.5)*winfact)))
        window(j)=0.5 * (1.0 + dcos((j - nyfft2 - 0.5)*winfact))
 201  continue
c
      if (debug) then
        do 229 j = 1, nyfft
          i = j + nyfft2
          if (j .gt. nyfft2) i = j - nyfft2
          write(*, 228) j, window(j), window(j)+window(i)
 228      format(i3, 1x, f7.5, 1x, f7.5)
 229    continue
      end if
c
c
c
      call getstuf2(ifilnam,ofilnam,params,iwant,iouttype,outmul)
c
C-*
      if (iwant(1).lt.0 .or. iwant(2).lt.0) then
C-* In this case we have a lowpass filter; 
C-* find global mean, remove it from each window
C-* but do not kill local means, add it at end. 
        iaddgm = 1
        ikillm = 0
        call findgm(glmean, nximg, nyimg, ifilnam)
c       write(*,237) float(glmean)
        write(*,237) glmean
 237    format('Grand mean value is ', f10.3)
      else
        iaddgm = 0
        ikillm = 1
      end if
C-*
c
c
c Here is the main loop:
c
      jrowseek = 0
      jlatscl = 0
      do 901 istrip = 1, nstrips
c
c Load the strip:
c
        if (istrip .gt. 2 .and. istrip .lt. nstrips) then
          jrowseek = jrowseek + nyfft2
        end if
        call fimgr(savex, nximg, nyfft, jrowseek, ikillm, ifilnam)
C-*
        if (iaddgm .eq. 1) call detrend(savex, nximg, nyfft, glmean)
C-*
c
c Window the strip.  This reverses the window and mirrors if
c istrip .ew. 1 .or. istrip .eq. nstrips
c
        call winstrp(savex, nximg, nyfft, istrip, nstrips, window)
c
c Now figure out what latitude we are at and what the array size in km is.
c The calculation in the subroutine getscal assumes that nximg is the 
c number of pixels in 360 degrees of longitude, in order to scale jlatscl
c correctly.  This seems a trivial comment as the FFT must be periodic in
c nximg, but I mention it in case we try to go back and use this program
c with old img files that went to 390 degrees of lon; then everything will
c need modification.
c
        if (istrip .gt. 1) then
          jlatscl = jlatscl + nyfft2
        end if
        call getscal(jlatscl, nximg, nyimg, nyfft, xwide, ytall, drlat)
c       write(*, 259) istrip, float(180.0*drlat/pi), xwide, ytall
        write(*, 259) istrip, 180.0*drlat/pi, xwide, ytall
 259  format('FFT ',i2,' at Lat ',f8.4,' is ',f7.1,' km by ',f7.1,' km') 
c
c Now do the filter
c
C-*        call filt(savex,savek,work,nximg,nyfft,xwide,ytall,iwant,params)
        call newfilt(savex,savek,work,xksq,nximg,nyfft,xwide,ytall,
     +    iwant,params)
c
c Now blend.
c
        if (istrip .eq. 1) then
c
c Copy the first half of the array to the blend strip.
c
          do 302 j = 1, nyfft2
            do 301 i = 1, nximg
              blend(i,j) = savex(i,j)
 301        continue
 302      continue
c
        else if (istrip .ne. nstrips) then
c
c Add the first half to the blend strip; write out, then copy 2nd half.
c
          do 352 j = 1, nyfft2
            do 351 i = 1, nximg
              blend(i,j) = blend(i,j) + savex(i,j)
 351        continue
 352      continue
c
         call fimgw(blend, nximg, nyfft2, iaddgm, glmean, ofilnam)
c
          do 362 j = 1, nyfft2
            do 361 i = 1, nximg
              blend(i,j) = savex(i,j+nyfft2)
 361        continue
 362      continue
c
        else
c
c Add the second half to the blend strip and write out
c
          do 392 j = 1, nyfft2
            do 391 i = 1, nximg
              blend(i,j) = blend(i,j) + savex(i,j+nyfft2)
 391        continue
 392      continue
c
          call fimgw(blend, nximg, nyfft2, iaddgm, glmean, ofilnam)
c
        end if
c
c End of main loop
c
 901  continue
      end
c
c
c
      subroutine winstrp(savex, nximg, nyfft, istrip, nstrips, window)
c
c Window the data in savex with the cosine window stored in window
c Except when istrip .eq. 1 .or. istrip .eq. nstrips, then do a
c special windowing with a mirror.
c
      real*4 savex(nximg, nyfft), window(nyfft)
      integer*4 istrip, nstrips
c
      integer*4 i, j, nyfft2, jw, jmirror
c
      nyfft2 = nyfft/2
c
      if (istrip .eq. 1) then
c
        do 702 j = 1, nyfft2
          jmirror = nyfft + 1 - j
          jw = j + nyfft2
          do 701 i = 1, nximg
            savex(i,j) = savex(i,j) * window(jw)
            savex(i,jmirror) = savex(i,j)
 701      continue
 702    continue
c
      else if (istrip .ne. nstrips) then
c
        do 752 j = 1, nyfft
          do 751 i = 1, nximg
            savex(i,j) = savex(i,j) * window(j)
 751      continue
 752    continue
c
      else
c
        do 792 j = 1, nyfft2
          jmirror = nyfft2 + 1 - j
          jw = j + nyfft2
          do 791 i = 1, nximg
            savex(i,jw) = savex(i,jw) * window(j)
            savex(i,jmirror) = savex(i,jw)
 791      continue
 792    continue
c
      end if
      return
      end
c
c
c
      subroutine getscal(jlatscl,nximg,nyimg,nyfft,xwide,ytall,drlat)
c
c Given nximg, the number of pixels in 360 degrees of longitude,
c and also the number of pixels in the horizontal direction of the FFT
c nyimg, such that nyimg/2 is the origin of the counter jlatscl,
c nyfft, the number of pixels in the vertical direction of the FFT,
c return xwide, ytall, and drlat, the width and height of the FFT in
c km and the double precision radian latitude for this scale factor.
c
      integer*4 jlatscl, nximg, nyimg, nyfft
      real*4 xwide, ytall
      real*8 drlat
c
c radius is the radius of the img mercator projection, that is, it is
c the number of pixels on the Equator per radian of longitude.
c twopi is 2 pi, and halfpi is pi/2; y is the argument to the 
c Gudermannian function, and pixelkm is the width of a pixel in km at
c the latitude drlat
c
      real*8 radius, twopi, halfpi, y, pixelkm
c
      halfpi = 2.0d+00 * datan(1.0d+00)
      twopi = 4.0d+00 * halfpi
      radius = nximg / twopi
      y = (nyimg/2 - jlatscl)/radius
      drlat = 2.0d+00 * datan(dexp(y)) - halfpi
      pixelkm = (360.0d+00 * 111.1949 * dcos(drlat))/nximg
c
      xwide = nximg * pixelkm
      ytall = nyfft * pixelkm
      return
      end
c
	subroutine getstuf2(filein,fileout,params,iwant,iouttype,outmul)
c
c Similar to getstuff, which is called by filter_grd, but set up for
c filter_img.  The prompts for filenames are a bit different, and there
c is a test to see if we are trying to write to an existing file, which
c is a no-no with this program.
c
	character*80 filein, fileout
	character*80 textin
	real params(13), outmul
	integer	iwant(6), iouttype
	logical exists
c	
	filein=textin('Enter input img v7.2 file name:', 
     +    '/geosat2/global_grav_2min/world_grav.img.7.2.high60')
	nc=index(filein,' ')
	filein(nc:nc) = '\0'
c
	fileout=textin('Enter output img filename:','filt.img.7.2')
	inquire(file=fileout, exist=exists)
	if (exists) then
	  write(*,888)
 888	  format('That name exists.  Remove it before continuing.')
        endif
	nc=index(fileout,' ')
	fileout(nc:nc) = '\0'
c
	outmul=realin('Enter multiplier to rescale data:', 1.0)
c
	write(*,311)
 311	format('-1 will do a lowpass cosine taper filter,')
	write(*,312)
 312	format('1 will do a highpass cosine taper filter,')
	write(*,313)
 313	format('0 will not do any cosine taper filter.')
	iwant(1)=intin('Enter -1,0,1 for the cosine filter',0)
	if (iwant(1) .ne. 0) then
	  write(*,314)
 314	  format('The cosine is ramped between h1 and h2,')
	  write(*,315)
 315	  format('which are harmonic degrees with h1 < h2')
	  params(1) = realin('Enter h1:', 50.0)
	  params(2) = realin('Enter h2:', 70.0)
	end if
c
	write(*,316)
 316	format('-1 will do a lowpass gaussian filter,')
	write(*,317)
 317	format('1 will do a highpass gaussian filter,')
	write(*,318)
 318	format('0 will not do any gaussian filter.')
	iwant(2)=intin('Enter -1,0,1 for the gaussian filter',0)
	if (iwant(2) .ne. 0) then
	  params(3) = realin('Enter sigma in km:', 30.0)
	end if
c
	iwant(3)=intin('Enter 1 for Wiener filter, else 0:',0)
	if (iwant(3) .ne. 0) then
	  params(4) = realin('Enter constant A in km**4:', 9500.0)
	end if
c
	write(*,319)
 319	format('-1 will do upward continuation,')
	write(*,320)
 320	format('1 will do downward continuation,')
	write(*,321)
 321	format('0 will not do any continuation.')
	iwant(4)=intin('Enter -1,0,1 for continuation',0)
	if (iwant(4) .ne. 0) then
	  params(5)=realin('Enter depth in km:', 0.0)
	end if
c
	iwant(5)=intin('Enter 1 to do flexural prediction:', 0)
	if (iwant(5) .ne. 0) then
	  params(6)=realin('Enter lambda in km:', 300.0)
	  params(7)=realin('Enter layer 2 thickness km:', 2.0)
	  params(8)=realin('Enter layer 3 thickness km:', 5.0)
	  params(9)=realin('Enter water density in SI:', 1030.0)
	  params(10)=realin('Enter layer 2 density, SI: ', 2650.0)
	  params(11)=realin('Enter layer 3 density, SI: ', 2950.0)
	  params(12)=realin('Enter mantle density, SI:', 3400.0)
	end if
c
	write(*,322)
 322	format('You may want to kill wavelengths longer than')
	write(*,323)
 323	format('the longest resolvable y wavelength.')
	iwant(6)=intin('Enter 1 to kill these, 0 not to:', 0)
	if (iwant(6) .ne. 0) then
	  params(13)=realin('Enter harmonic degree for kill:', 20.0)
	end if
c
	return
	end
C-*
C-*
      subroutine detrend(savex, nximg, nyfft, glmean)
      integer*4 nximg, nyfft
      real*4 savex(nximg,nyfft)
      real*8 glmean
C-*
      integer*4 i, j
C-*
      do 777 j = 1, nyfft
        do 776 i = 1, nximg
          savex(i,j) = savex(i,j) - glmean
 776    continue
 777  continue
      return
      end