c.......................................................................
c
c     Program pfgrav3d conducts an interactive session with the user that
c  prepares input to and reports output from subroutine fgrav3d.  The
c  questions asked by pfgrav3d should be self-explanatory.  See the
c  comments in subroutine fgrav3d and Open-File Report 81-298 by R.Blakely
c  for more information.
c
c
c  This is version number 2.  The main differences between this version 
c  and the original version are as follows:
c
c  1.  The first question of version 2 asks the user, 'Do you want
c      verbose prompts?'  A yes answer causes conversation to be chatty;
c      a no answer causes terse conversation.
c
c  2.  Grids are checked for the presence of dvals and to be sure all
c      grids have the same dimensions.  Violations are punishable by
c      program termination.
c
c.......................................................................
c...
c
      parameter(isza=6400,iszs=80)
      common /input/ s(isza), ztop(isza)
      common /input1/ zbot(isza), nx, ny, naugx,
     &naugy, dx, dy, zh, nstop, err2, 
     &lulist, luterm
      common /output/ h(isza)
c      dimension xdummy(1), ydummy(1)
      character ztop_file*32, zbot_file*32, h_file*32, s_file*32
c     &print_file*32
      character id*56
      character answer*1,ansupr*1,ansbot*1,ansden*1
      parameter (dval = 1.e+37)
      logical*1 verbose
c......................................................................
c  --  Print a banner and determine whether verbose prompts are wanted
c......................................................................
      data nxdim / iszs /
      data nydim / iszs /
      index(i,j,ny) = ((i - 1) * ny) + j
      ngrid = 0
      write(unit=6, fmt=104)
  104 format(///,' Welcome to PFGRAV3D - Version 2...'
     &' Do you want verbose prompts ? (y or n)',$)
      read(unit=5, fmt=101) answer
  101 format(a1)
      verbose = .false.
      if (answer .eq. 'y') verbose = .true.
      if (verbose) write(unit=6, fmt=139) 
c................................................................
c  --  Open the formatted file to receive program output
c................................................................
  139 format(/,' PFGRAV3D is a conversational program for calculating'
     &' a rectangular',/,
     &' grid of gravity anomaly values from three rectangular'
     &' grids that',/,
     &' define the source:  the top surface, the bottom surface'
     &' and the',/,' density contrast.  Mathematical and'
     &' programming concepts',/,
     &' are described in USGS Open File 81-298.',//,
     &' PFGRAV3D expects all grids to be in standard format'
     &' with no flagged',/,
     &' values (dvals).  It also expects all grids to have the'
     &' same number',/,
     &' of rows and columns.  Grid distance units are assumed to be',/,
     &' kilometers and density units are assumed to be g/cc.')
c      write(unit=6, fmt=125)
c  125 format(/,' Name of print-output file = ?' $)
c      read(unit=5, fmt=901) print_file
c      open(unit=20, file=print_file, form='formatted', status='unknown')
c      write(unit=20, fmt=105)
c......................................................................
c  --  Determine if user wants a flat, horizontal upper surface.  If so,
c      enter the constant level
c......................................................................
c  105 format(///,' PROGRAM PFGRAV3D - VERSION 2....',//)
      if (verbose) write(unit=6, fmt=140) 
      if (.not. verbose) write(unit=6, fmt=240) 
  140 format(/,' You will have a choice for both the upper and'
     &' lower surface of the',/,
     &' body:  It can be flat and horizontal [in which'
     &' case you will be',/,
     &' asked for that single level] or it can be bumpy'
     &' and/or tilted [in',/,
     &' which case you will be asked for the name of a'
     &' file containing the',/,' gridded surface].',//,
     &' Is the upper surface of the body flat and horizontal ? ',$)
  240 format(/,' Is upper surface flat ? '$)
      read(unit=5, fmt=101) ansupr
c......................................................................
c  --  Open and read the standard file containing the gridded upper surf
cace
c......................................................................
      if (ansupr .eq. 'y') goto 47
      write(unit=6, fmt=103) 
  103 format(/,' Standard file containing gridded upper surface = ? '$)
      read(unit=5, fmt=901) ztop_file
  901 format(a32)
c      write(unit=20, fmt=128) ztop_file
  128 format(/,' FILE OF UPPER SURFACE = ',a32)
      open(unit=11, file=ztop_file, form='unformatted', status='old') 
      call read_write(0, 11, id, ny, nx, y0, dy, x0, dx, ztop)
      write(unit=6, fmt=102) id, nx, x0, dx, ny, y0, dy
c      write(unit=20, fmt=102) id, nx, x0, dx, ny, y0, dy
      ngrid = ngrid + 1
      nxtest = nx
      nytest = ny
c......................................................................
c  -- Correct upper surface units
c......................................................................
      ntotal = nx * ny
      if (verbose) write(unit=6, fmt=122) 
      if (.not. verbose) write(unit=6, fmt=222) 
  122 format(/,' What single number should these values be multiplied'
     &' by so that',/,
     &' they will increase vertically down and have the same'
     &' units as dx',/,
     &' and dy [For example, if the topographic values are'
     &' in meters and',/,
     &' increase upwards and if dx and dy are in'
     &' kilometers, then',/,' -.001 is the appropriate response] = ?'
     & ,$)
  222 format(/,' Conversion factor = ? ',$)
      read(unit=5, fmt=*) topscale
      tmax = -.1e20
      tmin = .1e20
      do 18 i = 1, ntotal
      if (ztop(i) .eq. dval) call gerror(1, *10)
      ztop(i) = ztop(i) * topscale
      tmax = amax1(tmax,ztop(i))
   18 tmin = amin1(tmin,ztop(i))
      write(unit=6, fmt=117) tmax, tmin
c......................................................................
c  --  Printer contour the upper surface
c......................................................................
  117 format(/,' The maximum and minimum of the upper surface are '
     &,2g10.3/)
c      write(unit=6, fmt=116)
c  116 format(/,' Printer contour the upper surface ? '$)
c      read(unit=5, fmt=101) answer
c      if (answer .eq. 'n') goto 13
c      if ((nx .lt. 41) .and. (ny .lt. 41)) goto 2
c      write(unit=6, fmt=133)
c  133 format(/,' Grid dimensions are too great for contouring')
c      goto 13
c    2 continue
c      write(unit=6, fmt=120)
c  120 format(/,' Contour interval = ? '$)
c      read(unit=5, fmt=*) pc
c      call contp(pc, ztop, nx, ny, dx, dy, 0., 0., 20, 0., 0, 0, xdummy
c     &, ydummy)
   13 continue
c......................................................................
c  --  User wants a flat upper surface, so read that level
c......................................................................
      goto 49
   47 continue
      ntotal = nxdim * nydim
      if (verbose) write(unit=6, fmt=141) 
      if (.not. verbose) write(unit=6, fmt=241) 
  141 format(/,' The elevation of the flat, horizontal upper'
     &' surface (same units as',/,
     &' dx and dy; remember, z is positive down) = ? '$)
  241 format(/,' Elevation = ? '$)
      read(unit=5, fmt=*) zconst
      do 48 i = 1, ntotal
   48 ztop(i) = zconst
c      write(unit=20, fmt=142) zconst
  142 format(/,' TOP OF SOURCE IS FLAT AND HORIZONTAL AT = ',g10.3/)
c......................................................................
c  --  See if user wants a flat, horizontal bottom; if so, read that lev
cel
c......................................................................
   49 continue
      if (verbose) write(unit=6, fmt=106) 
      if (.not. verbose) write(unit=6, fmt=206) 
  106 format(/,' Is the bottom surface of the body flat'
     &' and horizontal ? '$)
  206 format(/,' Is bottom surface flat ? '$)
      read(unit=5, fmt=101) ansbot
c......................................................................
c  --  Open and read the standard file containing the gridded bottom sur
cface
c......................................................................
      if (ansbot .eq. 'y') goto 3
      write(unit=6, fmt=107) 
  107 format(/,' Standard file containing gridded bottom'
     &' surface = ? '$)
      read(unit=5, fmt=901) zbot_file
c      write(unit=20, fmt=129) zbot_file
  129 format(/,' FILE OF BOTTOM SURFACE = ',a32)
      open(unit=12, file=zbot_file, form='unformatted', status='old') 
      call read_write(0, 12, id, ny, nx, y0, dy, x0, dx, zbot)
      write(unit=6, fmt=102) id, nx, x0, dx, ny, y0, dy
c...............................................................
c  --  Test this grid for grid errors
c...............................................................
c      write(unit=20, fmt=102) id, nx, x0, dx, ny, y0, dy
      ngrid = ngrid + 1
      if (ngrid .gt. 1) goto 1
      nxtest = nx
      nytest = ny
    1 if ((nx .ne. nxtest) .or. (ny .ne. nytest)) call gerror(2, *10)
c......................................................................
c  --  Correct bottom surface units
c......................................................................
      ntotal = nx * ny
      if (verbose) write(unit=6, fmt=122) 
      if (.not. verbose) write(unit=6, fmt=222) 
      read(unit=5, fmt=*) topscale
      do 20 i = 1, ntotal
      if (zbot(i) .eq. dval) call gerror(1, *10)
   20 zbot(i) = zbot(i) * topscale
      goto 5
c......................................................................
c  --  User wants a flat bottom so create zbot
c......................................................................
    3 continue
      if (verbose) write(unit=6, fmt=108) 
      if (.not. verbose) write(unit=6, fmt=241) 
  108 format(/,' The elevation of the flat, horizontal bottom'
     &' surface (same units',/,
     &' as dx and dy; remember, z is positive down) = ? '$)
      read(unit=5, fmt=*) zconst
      do 8 i = 1, ntotal
      zbot(i) = zconst
    8 continue
c      write(unit=20, fmt=109) zconst
  109 format(/,' BOTTOM OF SOURCE HAS A FLAT HORIZONTAL BOTTOM = '
     &,g10.3/)
c......................................................................
c  --  See if user wants constant density; if not, open and read
c      the standard file containing the magnetization
c......................................................................
    5 continue
      if (verbose) write(unit=6, fmt=134) 
      if (.not. verbose) write(unit=6, fmt=234) 
  134 format(/,' Do you want the density contrast'
     &' to be uniform',/,
     &' everywhere between the top and bottom surfaces ? '$)
  234 format(/,' Use a constant density ? '$)
      read(unit=5, fmt=101) ansden
      if (ansden .eq. 'n') goto 37
      if(ansupr.eq.'y'.and.ansbot.eq.'y'.and.ansden.eq.'y')then
         write(6,121)
  121    format(' ** ERROR...At least one standard grid is required')
         stop
         end if
      write(unit=6, fmt=135) 
  135 format(/,' What is the value of the density = ? '$)
      read(unit=5, fmt=*) c_mag
c      write(unit=20, fmt=136) c_mag
  136 format(/,' BODY HAS UNIFORM DENSITY = ',g10.3/)
      do 38 i = 1, ntotal
      s(i) = c_mag
   38 continue
      goto 11
   37 continue
      write(unit=6, fmt=100) 
  100 format(/,' Standard file containing gridded density = ? '$)
      read(unit=5, fmt=901) s_file
c      write(unit=20, fmt=127) s_file
  127 format(/,' FILE OF DENSITY VALUES = ',a32)
      open(unit=10, file=s_file, form='unformatted', status='old') 
      call read_write(0, 10, id, ny, nx, y0, dy, x0, dx, s)
      write(unit=6, fmt=102) id, nx, x0, dx, ny, y0, dy
c      write(unit=20, fmt=102) id, nx, x0, dx, ny, y0, dy
c............................................................
c  --  Test this grid for grid errors
c............................................................
  102 format(/,15h HEADER:  id = ,a56,/,10x,5hnx = ,i5,5x,5hx0 = 
     &,g10.3,5x,5hdx = ,g10.3,/,10x,5hny = ,i5,5x,5hy0 = ,g10.3,5x,
     &5hdy = ,g10.3//)
      ngrid = ngrid + 1
      if (ngrid .gt. 1) goto 4
      nxtest = nx
      nytest = ny
    4 if ((nx .ne. nxtest) .or. (ny .ne. nytest)) call gerror(2, *10)
      do 14 i = 1, ntotal
c......................................................................
c  --  Printer contour the density
c......................................................................
      if (s(i) .eq. dval) call gerror(1, *10)
   14 continue
c      write(unit=6, fmt=114)
c  114 format(/,' Printer contour the density ? '$)
c      read(unit=5, fmt=101) answer
c      if (answer .eq. 'n') goto 11
c      if ((nx .lt. 41) .and. (ny .lt. 41)) goto 6
c      write(unit=6, fmt=133)
c      goto 11
c    6 continue
c      smax = -.1e20
c      smin = .1e20
c      do 12 i = 1, ntotal
c      smax = amax1(smax,s(i))
c      smin = amin1(smin,s(i))
c   12 continue
c      write(unit=6, fmt=115) smax, smin
c  115 format(/,' The maximum and minimum density contrast are '
c     &,2g10.3,//,'  Contour interval = ? '$)
c      read(unit=5, fmt=*) pc
c      call contp(pc, s, nx, ny, dx, dy, 0., 0., 20, 0., 0, 0, xdummy,
c     &ydummy)
c......................................................................
c  --  Open the standard file to receive the gridded anomaly
c......................................................................
   11 continue
      write(unit=6, fmt=111) 
  111 format(/,' Standard file to receive gridded anomaly'
     &' output = ? '$)
      read(unit=5, fmt=901) h_file
c      write(unit=20, fmt=130) h_file
  130 format(/,' FILE OF ANOMALY = ',a32)
c......................................................................
c  --  Request a title for the output
c......................................................................
      open(unit=13, file=h_file, form='unformatted', status='unknown') 
      if (verbose) write(unit=6, fmt=131) 
      if (.not. verbose) write(unit=6, fmt=231) 
  131 format(/,' Type an identifying label for this job, 56'
     &' characters or less.',/,
     &' This becomes part of the header record of the'
     &' output standard grid')
  231 format(/,' ID = ? '$)
      read(unit=5, fmt=138) id
c......................................................................
c  --  Request the various input parameters
c......................................................................
  138 format(a56)
      if (verbose) call fftdims(ny, nx)
   24 continue
      if (verbose) write(unit=6, fmt=110) 
      if (.not. verbose) write(unit=6, fmt=210) 
  110 format(/,' Enter the following five parameters on one line.'
     &' Separate each',/,' number by a comma or a blank...',//,
c     &49h sincl = Inclination of magnetization in degrees,
c     &15h positive below,/,23h         the horizontal,/,
c     &49h sdecl = Declination of magnetization in degrees,
c     &17h positive east of,/,22h         of true north,/,
c     &51h fincl = Inclination of regional field [in degrees],/,
c     &51h fdecl = Declination of regional field [in degrees],/,
c     &56h azim  = Azimuth of positive x axis in degrees, positive
c     &8h east of,/,19h         true north,/,
     &' zh    = Survey elevation in same units as dx and dy,'
     &' z positive',/,'         down',/,
     &' nstop = Limit on number of summations, typically 20',/,
     &' err2  = Convergence criterion for summations, typically 0.05'
     &,/,' naugx = Number of rows in expanded grids [Must be'
     &' greater than or',/,
     &'         equal to nx; use above list of fft dimensions'
     &' as a guide]',/,
     &' naugy = Number of columns in expanded grids [Must be'
     &' greater than',/,'         or equal to ny]',/)
  210 format(/,' Enter:  zh,nstop,'
     &'err2,naugx,naugy',/)
c......................................................................
c  --  check for various errors in input parameters
c......................................................................
      read(unit=5, fmt=*) zh, nstop, 
     &err2, naugx, naugy
      if ((((naugx .le. nxdim) .and. (naugy .le. nydim)) .and. (naugx
     & .ge. nx)) .and. (naugy .ge. ny)) goto 17
      write(unit=6, fmt=132) 
  132 format(/,44h ERROR...naugx or naugy either exceed array ,/,
     &53hdimensions or they are less than nx or ny...try again/)
      goto 24
   17 continue
c      write(unit=20, fmt=123) zh,
c     &nstop, err2, naugx, naugy
  123 format(//,29h THE INPUT PARAMETERS WERE...,//,
c     &50h  sincl = Inclination of magnetization          = ,g10.3,/,
c     &50h  sdecl = Declination of magnetization          = ,g10.3,/,
c     &50h  fincl = Inclination of regional field         = ,g10.3,/,
c     &50h  fdecl = Declination of regional field         = ,g10.3,/,
c     &50h  azim  = Azimuth of positive x axis            = ,g10.3,/,
     &50h  zh    = Elevation of survey                   = ,g10.3,/,
     &50h  nstop = Limit on number of summations         = ,i5,/,
     &50h  err2  = Convergence criterion for summations  = ,g10.3,/,
     &50h  naugx = Rows added to arrays                  = ,i5,/,
     &50h  naugy = Columns added to arrays               = ,i5//)
      luterm = 6
      lulist = 20
      naugx = naugx - nx
c......................................................................
c  --  Call the subroutine fgrav3d
c......................................................................
      naugy = naugy - ny
c......................................................................
c  --  Print the gridded anomaly
c......................................................................
      call fgrav3d
c      write(unit=6, fmt=126) nx * ((ny / 10) + 1)
c  126 format(/,46h Include anomaly values in print output (list ,/,
c     &8hwill be ,i7,15h lines long) ? $)
c      read(unit=5, fmt=101) answer
c      if (answer .eq. 'n') goto 19
c      write(unit=20, fmt=112)
c  112 format(//,2x,11hX DECREASES,2x,1h|,10x,19hCALCULATED ANOMALY
c     &12h Y INCREASES,14h TO THE RIGHT),/,6x,4hDOWN,5x,1h|,/,1x,14(1h-),
c     &1h|,110(1h-),/,15x,1h|)
c      i = nx
c      do 9 l = 1, nx
c      i1 = index(i,1,ny)
c      i2 = i1 + 9
c      x = x0 + ((i - 1) * dx)
c      write(unit=20, fmt=113) x, (h(m),m = i1, i2)
c  113 format(5x,g10.3,1h|,10g11.4)
c      if (ny .gt. 10) goto 9
c      i3 = i2 + 1
c      i4 = index(i,ny,ny)
c      write(unit=20, fmt=124) (h(m),m = i3, i4)
c  124 format(15x,1h|,10g11.4)
c    9 i = i - 1
c......................................................................
c  --  Send the anomaly output to the standard file
c......................................................................
   19 continue
c......................................................................
c  --  Calculate and report certain statistics about the anomaly
c......................................................................
      call read_write(1, 13, id, ny, nx, y0, dy, x0, dx, h)
      sumh = 0.
      hmax = -.1e20
      hmin = .1e20
      do 16 i = 1, ntotal
      sumh = sumh + h(i)
      hmax = amax1(hmax,h(i))
      hmin = amin1(hmin,h(i))
   16 continue
      sumh = sumh / ntotal
      write(unit=6, fmt=137) hmax, hmin, sumh
c......................................................................
c  --  Contour the anomaly if desired
c......................................................................
  137 format(/,48h THE MAX, MIN, AND AVERAGE VALUES OF THE ANOMALY
     &5h ARE ,g11.4,/,53x,g11.4,/,53x,g11.4/)
c      write(unit=6, fmt=118)
c  118 format(/,31h Printer contour the anomaly ? $)
c      read(unit=5, fmt=101) answer
c      if (answer .eq. 'n') goto 15
c      if ((nx .lt. 41) .and. (ny .lt. 41)) goto 7
c      write(unit=6, fmt=133)
c      goto 15
c    7 continue
c      write(unit=6, fmt=119)
c  119 format(/,22h Contour interval = ? $)
c      read(unit=5, fmt=*) pc
c      call contp(pc, h, nx, ny, dx, dy, 0., 0., 20, 0., 0, 0, xdummy,
c     &ydummy)
c......................................................................
c  --  Termination of program
c......................................................................
   15 continue
   10 continue
      close(unit=10) 
      close(unit=11) 
      close(unit=12) 
      close(unit=13) 
c      close(unit=20)
      close(unit=6) 
      close(unit=5) 
      stop 
      end
      subroutine gerror(iflag, *)
      goto (10, 20), iflag
   10 write(unit=6, fmt=100) 
  100 format(/,44h ERROR...a dval was encountered in the grid.)
      return 1
   20 write(unit=6, fmt=101) 
  101 format(/,46h ERROR...Grid dimensions do not match previous
     &9h grid(s).)
      return 1
      end
c     Prints out numbers less than 1024 with no prime factor greater tha
cn 5.    
c     The FFT algorithm works fastest for these numbers.                
c        
c     Subroutine prints out good numbers between                        
c        
c       min of (ncol,nrow) and max of (ncol+100,nrow+100).              
c        
c                                                                       
c        
      subroutine fftdims(ncol, nrow)
c     first 2 cols from table III of Singleton's paper on fft....       
c        
c       (Singleton,R.C., 1969, An algotithm for computing the mixed     
c        
c          radix fast Fourier transform:  IEEE Trans., Audio and Electro
c-       
c          acoustics: v. AU-17, p. 93-103.                              
c        
c                                                                       
c        
      dimension nlist(104), ngood(100)
c                                                                       
c        
      data nlist / 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 
     &27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96
     &, 100, 108, 120, 125, 128, 135, 144, 150, 160, 162, 180, 192, 200
     &, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375
     &, 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625
     &, 640, 648, 675, 720, 729, 750, 768, 800, 810, 864, 900, 960, 972
     &, 1000, 1024, 1080, 1125, 1152, 1200, 1215, 1250, 1280, 1296, 1350
     &, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875 /
      min = min0(ncol,nrow)
      if (min .eq. 1) min = max0(ncol,nrow)
      max = max0(ncol + 100,nrow + 100)
      n = 0
      do 10 i = 1, 104
      if (nlist(i) .lt. min) goto 10
      if (nlist(i) .gt. max) goto 20
      n = n + 1
      ngood(n) = nlist(i)
c                                                                       
c        
   10 continue
   20 write(unit=6, fmt=902) 
  902 format(38h Suggested values for naugx and naugy:)
      write(unit=6, fmt=101) (ngood(i),i = 1, n)
c                                                                       
c        
  101 format(6h      ,14i5)
      return 
      end
      subroutine read_write(mode, lun, id, ncol, nrow, col1, delcol, 
     &row1, delrow, array)
      dimension array(1), work(1000)
      character id*56
c
c.......................................................................
c
c     This subroutine reads a standard file into a one-dimensional array
c.
c  Inputs are mode (0 for read, 1 for write) and lun (logical unit numbe
cr).
c  Outputs are id,ncol,nrow,col1,delcol,row1,delrow, and array.
c
c.......................................................................
c
      character pgm*8
      if (mode .eq. 1) goto 5
      read(unit=lun) id, pgm, ncol, nrow, nz, col1, delcol, row1, delrow
      goto 6
    5 nz = 1
      pgm = 'PFMAG3D '
      write(unit=lun) id, pgm, ncol, nrow, nz, col1, delcol, row1, 
     &delrow
c
    6 continue
      do 1 j = 1, nrow
      l = 0
      i1 = ((j - 1) * ncol) + 1
      i2 = (i1 + ncol) - 1
      if (mode .eq. 1) goto 3
      call row_read(lun, work, ncol)
      do 2 i = i1, i2
      l = l + 1
    2 array(i) = work(l)
      goto 1
    3 continue
      do 4 i = i1, i2
      l = l + 1
    4 work(l) = array(i)
      call row_write(lun, work, ncol)
    1 continue
      return 
      end
      subroutine row_read(lun, work, ncol)
      dimension work(ncol)
      read(unit=lun) dummy, work
      return 
      end
      subroutine row_write(lun, work, ncol)
      dimension work(ncol)
      dummy = 999.
      write(unit=lun) dummy, work
      return 
      end
c      donald plouff printer contour routine
c          pc,dc---the increment in the contour values (read in)
c        a(y,x)---the input array containing the data to be contoured.
c index for equivalent 1d array is n*(x-1)+y
c the one-dimensinal a-array is filled (without voids) in rows, starting
c   from minimum x (xmn) toward maximum x. within each row the progres-
c  sion is from minimum y (ymn) to maximum y.
c  m---the  final  value of x-index in a-array, starting at 1
c  n---the  final  value of y-index in a-array, starting at 1
c x increases to right and y increases downward
c dx,dy are respective grid intervals in distance units
c iout is unit number of printer
c scl is scale in distance units per inch. default fills page.
c lb is number of (xb,yb) locations to be superimposed on plot. if this
c   facility is not needed, then lb=0 (also xmn,ymn not needed).
c   note that this subroutine changes the values of xb and yb.
c nx is number of intervals in x-direction (right) less than 41.
c ny is number of intervals in y-direction (down)  less than 41.
c      subroutine contp(pc, a, m, n, dx, dy, xmn, ymn, iout, scl, laxcol
c     &, lb, xb, yb)
c      integer pnt(131), sym(35), c, iy(41), jx(41)
cc     dimension a(41,41),xb(1),yb(1)
c      dimension a(1), xb(1), yb(1)
cc needs carriage with 132 columns available or change naxcol
c      data sym / '1', '2', '3', '4', '5', '6', '7', '8', '9', 'a', 'b',
c     &'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o',
c     &'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z' /
c      data ihor / '-' /
c      data ivert / '|' /
c      data jplus / '+' /
c      data iblank / ' ' /
c      data jex / '*' /
c      naxcol = laxcol - 1
c      if (naxcol .eq. (-1)) naxcol = 131
c      if (m .lt. 2) goto 150
c      if (m .gt. 41) goto 150
c      if (n .lt. 2) goto 150
c      if (n .gt. 41) goto 150
c      dc = pc
c      nb = lb
c      nx = m - 1
c      ny = n - 1
c      nt = n * m
c      xmx = xmn + (dx * nx)
c      ymx = ymn + (dy * ny)
c      if (scl .le. 0.0) goto 24
c      maxcol = (((10.0 * dx) * nx) / scl) + 1.5
c      if (maxcol .gt. naxcol) goto 24
c      goto 38
c   24 maxcol = naxcol
c      write(unit=iout, fmt=203)
c  203 format(5x,39hprinter scale is expanded to fill page.)
c   38 write(unit=iout, fmt=603) m, n, nt, xmn, xmx, ymn, ymx, dc
c  603 format(5x,22hexpect printer plot of,i3,3h by,i3,1h(,i3,8h) array.,
c     &3h x=,f8.1,3h to,f8.1,3h y=,f8.1,3h to,f8.1,13h contoured at,f7.1,
c     &6h units)
c      fmin = 1.0e20
c      fmax = -1.0e20
c      do 4 i = 1, nt
c      aij = a(i)
c      if (aij .gt. fmax) fmax = aij
c      if (aij .lt. fmin) fmin = aij
c    4 continue
c   29 ic = 1.0 + (fmin / dc)
c      if (ic .lt. 0) ic = ic - 1
c      jc = fmax / dc
c      if (jc .lt. 0) jc = jc - 1
c      ncont = (jc - ic) + 1
c      if (ncont .lt. 36) goto 35
c      write(unit=iout, fmt=600) dc
c  600 format(5x,29hrequested contour interval of,f9.4,11h is doubled,
c     &44h because more than 35 contours are produced.)
c      dc = 2.0 * dc
c      goto 29
c   35 cd = 1.0 / dc
cc     cf=dc*jc
c      ci = dc * ic
c      width = maxcol - 1
cc number of columns in one grid interval
c      ux = nx / width
c      fjinc = 1.0 / ux
c      fjincp = 0.5 + fjinc
cc maintains true proportions with x expanded to maxcol-1 columns
c      colefm = 1.0 - fjinc
c      uy = (dx * ux) / (dy * 0.6)
c      if (uy .lt. 1.0) goto 11
c      uy = 1.0 / uy
c      write(unit=iout, fmt=601) uy
c  601 format(5x,30hy-distances are exaggerated by,f9.4,
c     &15h for contouring)
cc number of columns per unit distance
c      uy = 1.0
c   11 dlx = fjinc / dx
c      fiinc = 1.0 / uy
c      dly = fiinc / dy
c      row = 1.0
c      do 5 k = 1, ny
c      row = row + fiinc
c    5 iy(k) = row + 0.5
c      lastrw = iy(ny)
c      col = 1.0
c      do 6 k = 1, m
c      jx(k) = col + 0.5
c    6 col = col + fjinc
c      if (nb .eq. 0) goto 39
c      l = 0
c      nf = 1
c    7 do 8 k = nf, nb
c      kk = k
c      dum = xb(k)
c      if (dum .lt. xmn) goto 9
c      if (dum .gt. xmx) goto 9
c      kx = 1.5001 + (dlx * (dum - xmn))
c      dum = yb(k)
c      if (dum .lt. ymn) goto 9
c      if (dum .gt. ymx) goto 9
c      ky = 1.5001 + (dly * (dum - ymn))
cc charges xb/yb. provide another array, if values needed.
c      l = l + 1
c      xb(l) = kx
c    8 yb(l) = ky
c      goto 12
c    9 nf = kk + 1
c      write(unit=iout, fmt=602) xb(kk), yb(kk)
c  602 format(5x,13hbody point x=,f10.3,3h y=,f10.3,11h is outside,
c     &14h contour plot.)
c      if (kk .lt. nb) goto 7
c   12 nb = l
c   39 cm = ci - dc
cc changes values inputted as 'a' (restored later).
c      do 26 i = 1, nt
c   26 a(i) = cd * (a(i) - cm)
cc loop of rows from top to bottom
c      jr = 1
c      do 27 jrow = 1, lastrw
cc major grid index (i) increases once each fiinc rows
c      row = 1.0 + (uy * (jrow - 1))
c      i = row
c      ful = row - i
c      ind = 0
c      kt = i - n
c      do 25 l = 1, maxcol
c   25 pnt(l) = iblank
cc loop of grid intervals from left to right
c      coleft = colefm
c      do 28 j = 1, nx
c      coleft = coleft + fjinc
c      if (jrow .eq. lastrw) goto 36
c      kt = kt + n
c      kb = kt + n
c      dum = a(kt)
c      zl = dum + (ful * (a(kt + 1) - dum))
c      dum = a(kb)
c      zr = dum + (ful * (a(kb + 1) - dum))
c      goto 37
c   36 ind = ind + n
c      zl = a(ind)
c      kb = ind + n
c      zr = a(kb)
c   37 den = zr - zl
c      if (abs(den) .gt. 0.000001) goto 30
c      izc = zl
c      fl = izc
c      if (fl .ne. zl) goto 28
c      left = coleft + 0.5
cc loop of contours in this grid interval
c      jrt = coleft + fjincp
c      do 31 l = left, jrt
c   31 pnt(l) = sym(izc)
c      goto 28
c   30 if (zl .gt. zr) goto 32
c      fmin = zl
c      fmax = zr
c      goto 33
c   32 fmin = zr
c      fmax = zl
c   33 lmin = fmin + 0.9999
c      lmax = fmax
c      if (lmax .eq. 0) goto 28
c      if (lmax .lt. lmin) goto 28
c      ratio = fjinc / den
c      do 34 izc = lmin, lmax
c      c = coleft + (ratio * (izc - zl))
c   34 pnt(c) = sym(izc)
c   28 continue
c      if (nb .eq. 0) goto 15
c      row = jrow
c      nf = 1
c   13 do 14 k = nf, nb
c      kk = k
c      if (yb(k) .eq. row) goto 16
c   14 continue
c      goto 15
c   16 nb = nb - 1
c      kx = xb(kk) + 0.5
c      pnt(kx) = jex
c      if (kk .gt. nb) goto 15
c      nf = kk
c      do 17 k = nf, nb
c      kp = k + 1
c      xb(k) = xb(kp)
c   17 yb(k) = yb(kp)
c      goto 13
c   15 itest = 0
c      if (jrow .eq. lastrw) goto 18
c      if (jrow .eq. 1) goto 18
c      if (jrow .eq. iy(jr)) goto 19
c      if (pnt(maxcol) .eq. iblank) pnt(maxcol) = ivert
c      if (pnt(1) .eq. iblank) pnt(1) = ivert
c   20 write(unit=iout, fmt=501) (pnt(k),k = 1, maxcol)
c  501 format(1x,131a1)
c      goto 27
c   19 itest = 1
c      jr = 1 + jr
c   18 do 21 k = 1, m
c      kx = jx(k)
c      if (pnt(kx) .eq. iblank) pnt(kx) = jplus
c   21 continue
c      if (itest .eq. 1) goto 20
c      do 23 k = 1, maxcol
c      if (pnt(k) .eq. iblank) pnt(k) = ihor
c   23 continue
c      if (jrow .eq. lastrw) goto 9000
c      write(unit=iout, fmt=502) (pnt(k),k = 1, maxcol)
c  502 format(1h1,131a1)
c   27 continue
c 9000 write(unit=iout, fmt=501) (pnt(k),k = 1, maxcol)
c      write(unit=iout, fmt=505)
c  505 format(//,7h legend,/,7h symbol,8x,5hvalue,12x,11hx increases
c     &34h to right and y increases downward)
c      zc = ci - dc
c      do 120 i = 1, ncont
c      zc = zc + dc
c  120 write(unit=iout, fmt=503) sym(i), zc
c  503 format(1h ,3x,a1,6x,f14.4)
c      write(unit=iout, fmt=504)
c  504 format(4x,1h+,6x,10hgrid point,14x,1h*,6x,13hlocation mark)
c      do 10 i = 1, nt
c   10 a(i) = cm + (dc * a(i))
c      return
c  150 write(unit=iout, fmt=900) m, n
c  900 format(1h0,7hgrid of,i4,3h by,i4,19h (x,y) is too large9h (41) or
c     &,30htoo small for printer contour.)
c      return
c      end
c
c.......................................................................
c..
c
      subroutine fgrav3d()
c
c.........................................................................
c
c     Subroutine fgrav3d computes a rectangular grid of anomaly values, h(x,y),
c  from three rectangular grids which define the distribution of mass:
c  the top surface, ztop(x,y); the bottom surface, zbot(x,y); and
c  the density, s(x,y).  The corners and sample intervals of each
c  of these four grids are assumed to be the same.
c
c  Algorithm
c
c    The method is described by R. L. Parker (Geophysical Journal, v. 31,
c  p. 447, 1973).  It solves the summation
c
c
c                     m
c         h = IF[B * SUM( A * F[s * (zbot**n - ztop**n)])]
c                    n=0
c
c
c  where F and IF denote two-dimensional Fourier transforms and inverse
c  transforms, respectively.  The summation is continued until sufficient
c  convergence occurs.  The rate of convergence depends on the elevation of
c  the origin with respect to the magnetic body and the level of optimum
c  convergence changes at different parts of the algorithm.  Thus upward
c  and downward continuation occurs many times in the course of a run.
c
c  Input
c
c     Although the arrays h, ztop, zbot, and s are singly dimensioned, they
c  nevertheless represent two-dimensional arrays.  These increase first along
c  the y axis and then along the x axis.  For example,
c
c
c      x axis
c
c        ^
c        |
c        |
c      s((nx-1)*ny+1)..........................s(nx*ny)
c        .                                       .
c        .                                       .
c        .                                       .
c        .                                       .
c      s(ny+1)       s(ny+2)       s(ny+3) ... s(ny+nx)
c      s(1)          s(2)          s(3)    ... s(ny)     --->  y axis
c
c
c
c  All of the input arrays and parameters are supplied to fgrav3d through
c  labeled common.  An itemized list follows:
c
c     s      - The array of density values
c     ztop   - The array representing the top surface
c     zbot   - The array representing the bottom surface
c     nx     - The number of array elements in the x direction
c     ny     - The number of array elements in the y direction
c     naug   - The number of rows and columns to add to the input arrays
c     dx     - the sample interval of the arrays in the x direction
c     dy     - the sample interval of the arrays in the y direction
c     zh     - the elevation of the survey
c     nstop  - the limit of summations, if convergence does not occur
c     err2   - the definition of convergence (see note below)
c     lulist - the logical unit number to receive printed output
c
c
c  Units
c
c     The assumed units of distance for zbot, ztop, dx, dy, and zh are
c  kilometers.  The unit of density for s is g/cc.
c
c  Output
c
c     The only output is the array h, the gravity anomaly in units of mgals.
c
c
c  Notes
c
c     1 - Convergence is assumed to occur when the absolute value of the
c         nth term is err2 times the absolute value of the sum of all previous 
c         terms.  A typical value for err2 is .001
c     2 - Subroutine fgrav3d calls a fast-Fourier transform subroutine called
c         fft written by R.C. Singleton (IEEE Trans. Audio Electroac., v.
c         AU-17, p. 93, 1969).  Subroutine fft is faster if the size of the
c         array being transformed can be factored into many small prime
c         numbers, and in fact certain dimensions are not permitted by fft
c         (see the documentation in fft for details).  Because the arrays
c         being sent to fgrav3d are not necessarily going to have optimal
c         dimensions, a facility is available to augment rows and columns.
c         The augmentation also helps reduce edge effects produced by the
c         Fourier transform.  These additional rows and columns are
c         discarded by fgrav3d before sending results back to the calling
c         program.
c
c     3 - For further information, see the paper, 'A program for rapidly
c         computine the magnetic anomaly over digital topography', USGS
c         open file report, 1981, by R. J. Blakely.
c
c.........................................................................
c
      parameter(isza=6400,iszs=80)
      common /input/ s(isza), ztop(isza)
      common /input1/ zbot(isza), nx, ny,
     &naugx, naugy, dx, dy, zh, nstop, err2, lulist, luterm
      common /output/ h(isza)
      common /inter1/ work1(iszs, iszs), work2(iszs, iszs)
      common /inter3/ fcmplx(isza)
      common /inter2/ sum(isza), k(isza), kx(iszs), ky(iszs), fact(31)
      common /inter4/ hcmplx(isza), expbot(isza), exptop(isza)
      complex sum, hcmplx, fcmplx
      real kx, ky, k
c      character mode*3
c.......................................................................
c..
c  --  Augment the arrays to reduce edge effects and speed up fft
c.......................................................................
c..
      data twopi / 6.2831853 /
      data epslon / .1e-10 /
      data g / 6.670e-8 /
      data gal_to_mgal / 1.e3 /
      data akm_to_cm / 1.e5 /
      index(i,j,ny) = ((i - 1) * ny) + j
      if ((naugx .eq. 0) .and. (naugy .eq. 0)) goto 17
      call augment(s, ny, nx, naugy, naugx, 1)
      call augment(ztop, ny, nx, naugy, naugx, 1)
      call augment(zbot, ny, nx, naugy, naugx, 1)
      nx = nx + naugx
      ny = ny + naugy
c......................................................................
c  --  Initialize a few parameters
c......................................................................
   17 continue
      do 22 i = 1, nstop
      fact(i) = 1. / fac(i,ier)
   22 continue
      dkx = twopi / (nx * dx)
      dky = twopi / (ny * dy)
      ntotal = nx * ny
      nnx = (nx / 2) + 1
c......................................................................
c  --  Initialize the wavenumber coordinates
c......................................................................
      nny = (ny / 2) + 1
      do 3 i = 1, nx
      if (i - nnx) 4, 4, 5
    4 kx(i) = (i - 1) * dkx
      goto 3
    5 kx(i) = ((i - nx) - 1) * dkx
    3 continue
      do 6 j = 1, ny
      if (j - nny) 7, 7, 8
    7 ky(j) = (j - 1) * dky
      goto 12
    8 ky(j) = ((j - ny) - 1) * dky
   12 continue
      sky = ky(j) ** 2
      do 6 i = 1, nx
      ij = index(i,j,ny)
      k(ij) = sqrt((kx(i) ** 2) + sky)
c......................................................................
c  --  Find maxima and minima of ztop and zbot arrays; fixup bad values
c      of ztop and zbot
c......................................................................
    6 continue
      ztpmax = -.1e20
      ztpmin = .1e20
      zbtmax = -.1e20
      zbtmin = .1e20
      nerr1 = 0
      nerr2 = 0
      do 1 i = 1, ntotal
      if (ztop(i) .gt. (zh + epslon)) goto 10
      nerr2 = nerr2 + 1
      ztop(i) = zh + epslon
   10 continue
      if (ztop(i) .le. zbot(i)) goto 2
      zbot(i) = ztop(i)
      nerr1 = nerr1 + 1
    2 continue
      ztpmax = amax1(ztpmax,ztop(i))
      ztpmin = amin1(ztpmin,ztop(i))
      zbtmax = amax1(zbtmax,zbot(i))
      zbtmin = amin1(zbtmin,zbot(i))
    1 continue
      if (nerr1 .gt. 0) write(unit=luterm, fmt=100) nerr1
c      if (nerr1 .gt. 0) write(unit=lulist, fmt=100) nerr1
  100 format(/,49h **WARNING  - THE BODY HAS NEGATIVE THICKNESS AT 
     &,i6,/,52h              POINTS.  BODY WILL HAVE ZERO THICKNESS,/,
     &33h              AT THESE LOCATIONS./)
      if (nerr2 .gt. 0) write(unit=luterm, fmt=106) nerr2
c      if (nerr2 .gt. 0) write(unit=lulist, fmt=106) nerr2
c......................................................................
c  --  Find median values for ztop and zbot; determine if either are fla
ct
c......................................................................
  106 format(/,51h **WARNING  - BODY RISES ABOVE SURVEY ELEVATION AT 
     &,i6,/,45h              POINTS.  BODY WILL BE TRUNCATED,/,
     &26h              ACCORDINGLY./)
      iflttp = 0
      ifltbt = 0
      if (ztpmax .eq. ztpmin) iflttp = 1
      if (zbtmax .eq. zbtmin) ifltbt = 1
      ztpmed = ztpmin + ((ztpmax - ztpmin) / 2.)
      zbtmed = zbtmin + ((zbtmax - zbtmin) / 2.)
      write(unit=luterm, fmt=108) ztpmin, ztpmax, ztpmed, zbtmin, zbtmax
     &, zbtmed
c      write(unit=lulist, fmt=108) ztpmin, ztpmax, ztpmed, zbtmin, zbtmax
c     &, zbtmed
c......................................................................
c  --  Find the upward continuation arrays, adjust surfaces to median le
cvels,
c      and make the magnetization complex
c......................................................................
  108 format(/,48h HIGH, LOW, AND MEDIAN POINTS OF TOPOGRAPHY ARE 
     &,//,3g11.4,///,47h THE HIGH, LOW, AND MEDIAN POINTS OF THE BOTTOM,
     &13h SURFACE ARE ,//,3g11.4/)
      iflag1 = 0
      iflag2 = 0
      do 11 i = 1, ntotal
      sum(i) = 0.
      arg1 = - (k(i) * (ztpmed - zh))
      arg2 = - (k(i) * (zbtmed - zh))
      if (arg1 + 88.028) 30, 31, 31
   30 iflag1 = iflag1 + 1
      exptop(i) = 0.
      goto 32
   31 exptop(i) = exp(arg1)
   32 continue
      if (arg2 + 88.028) 33, 34, 34
   33 iflag2 = iflag2 + 1
      expbot(i) = 0.
      goto 35
   34 expbot(i) = exp(arg2)
   35 continue
      ztop(i) = ztop(i) - ztpmed
      zbot(i) = zbot(i) - zbtmed
      fcmplx(i) = s(i)
c......................................................................
c  --  Report problems with upward continuation
c......................................................................
   11 continue
      if (iflag1 .gt. 0) write(unit=luterm, fmt=111) iflag1
c      if (iflag1 .gt. 0) write(unit=lulist, fmt=111) iflag1
  111 format(/,48h **WARNING - UPWARD CONTINUATION EXPONENTIAL FOR,/,
     &38h             UPPER SURFACE UNDERFLOWS ,i6,6h TIMES,/,
     &41h             AND THESE WILL BE MADE ZERO./)
      if (iflag2 .gt. 0) write(unit=luterm, fmt=112) iflag2
c      if (iflag2 .gt. 0) write(unit=lulist, fmt=112) iflag2
c......................................................................
c  --  Do the n=0 term because its simple
c......................................................................
  112 format(/,48h **WARNING - UPWARD CONTINUATION EXPONENTIAL FOR,/,
     &38h             LOWER SURFACE UNDERFLOWS ,i6,6h TIMES,/,
     &41h             AND THESE WILL BE MADE ZERO./)
      call fftsub(work1, work2, nx, ny, -1)
      do 29 i = 2, ntotal
      hcmplx(i) = - ((fcmplx(i) * (expbot(i) - exptop(i))) / k(i))
   29 continue
c......................................................................
c  --  Get ready for summation
c......................................................................
      hcmplx(1) = - (fcmplx(1) * (ztpmed - zbtmed))
      write(unit=luterm, fmt=110) err2
c      write(unit=lulist, fmt=110) err2
  110 format(/,25h PERFORMANCE OF SUMMATION,/,17h TRYING TO REACH ,f5.3,
     &10h OF STOTAL/)
      write(unit=luterm, fmt=101)
c      write(unit=lulist, fmt=101)
c......................................................................
c  --  Begin summations
c......................................................................
  101 format(/,4x,1hN,4x,6hSTOTAL,5x,5hSLAST,4x,6hSRATIO/)
      n = 0
c......................................................................
c  --  Combine upper surface with density and Fourier transform; if uppe
cr
c      is flat, skip this step
c......................................................................
   14 n = n + 1
      if (iflttp .eq. 1) goto 13
      do 21 i = 1, ntotal
      fcmplx(i) = s(i) * (ztop(i) ** n)
   21 continue
c......................................................................
c  --  Multiply by upward continuation array
c......................................................................
      call fftsub(work1, work2, nx, ny, -1)
      do 15 i = 1, ntotal
      sum(i) = - (exptop(i) * fcmplx(i))
   15 continue
c......................................................................
c  --  Combine lower surface with density and Fourier transform; if lowe
cr
c      surface is flat, skip this step
c......................................................................
   13 continue
      if (ifltbt .eq. 1) goto 16
      do 24 i = 1, ntotal
      fcmplx(i) = s(i) * (zbot(i) ** n)
   24 continue
c......................................................................
c  --  Multiply by upward continuation array and add to previous result 
cfor
c      upper surface
c......................................................................
      call fftsub(work1, work2, nx, ny, -1)
      do 23 i = 1, ntotal
      sum(i) = (expbot(i) * fcmplx(i)) + sum(i)
   23 continue
c......................................................................
c  --  Multiply this term by -k**n/n! and find its contribution; find th
ce
c      contribution of the sum of all previous terms; add the new term
c......................................................................
   16 continue
      slast = 0.
      stotal = 0.
      do 26 i = 1, ntotal
      if ((n .ne. 1) .or. (i .ne. 1)) goto 20
      arg = 1.
      goto 26
   20 continue
      arg = ((- k(i)) ** (n - 1)) * fact(n)
      sum(i) = arg * sum(i)
      slast = slast + cabs(sum(i))
      stotal = stotal + cabs(hcmplx(i))
      hcmplx(i) = hcmplx(i) + sum(i)
   26 continue
      slast = slast / ntotal
      stotal = stotal / ntotal
      stest = err2 * stotal
c......................................................................
c  --  check for convergence of sum
c......................................................................
      sratio = slast / (stotal + .1e-20)
      write(unit=luterm, fmt=104) n, stotal, slast, sratio
c      write(unit=lulist, fmt=104) n, stotal, slast, sratio
  104 format(i5,3g10.3)
c......................................................................
c  --  Summation has converged; multiply by G and inverse Fourier transf
corm
c......................................................................
      if ((n .lt. nstop) .and. (slast .ge. stest)) goto 14
      factor = (twopi * g) * gal_to_mgal * akm_to_cm
      do 9 j = 1, ny
      do 9 i = 1, nx
      ij = index(i,j,ny)
      fcmplx(ij) = hcmplx(ij) * factor
    9 continue
c......................................................................
c  --  Put real part of complex anomaly back into real array
c......................................................................
      call fftsub(work1, work2, nx, ny, 1)
      do 19 i = 1, ntotal
      h(i) = fcmplx(i) / ntotal
c......................................................................
c  --  Unaugment the arrays
c......................................................................
   19 continue
      if ((naugx .eq. 0) .and. (naugy .eq. 0)) goto 18
      nx = nx - naugx
      ny = ny - naugy
      call augment(h, ny, nx, naugy, naugx, 0)
   18 continue
      return 
      end
c
      function fac(n, ier)
      ier = 0
      if (n .gt. 1) goto 10
      if (n .ge. 0) goto 20
      ier = 1
   20 fac = 1.
      return 
   10 fac = n
      fac2 = fac
   30 fac2 = fac2 - 1.
      if (fac2 .eq. 1.) return 
      fac = fac * fac2
      goto 30
      end
c
      subroutine augment(f, ncol, nrow, naugcol, naugrow, iadd)
      parameter(isza=6400,iszs=80)
      dimension f(isza)
c
c.......................................................................
c...
c
c     Subroutine augment adds or subtracts naugcol columns and naugrow
c  rows to or from array f.  If iadd=1, the new array f has dimensions
c  (ncol+naugcol) X (nrow+naugrow) and the new
c  values added to f are calculated to be straight-line slopes so as to
c  make f periodic.  For example, if f(2,1)=4, f(2,nrow)=0, and naugrow=
c3, then
c  f(2,nrow+1)=3, f(2,nrow+2)=2, and f(2,nrow+3)=1.
c     if iadd=0, subroutine augment takes an array f of dimensions
c  (ncol+naugcol) X (nrow+naugrow) and strips off the last naugcol colum
cns
c  and naugrow rows to leave f with dimensions ncol X nrow.
c
c.......................................................................
c...
c
      common /inter5/ temp(isza)
      index(i,j,n) = ((i - 1) * n) + j
      nrp1 = nrow + 1
      ncp1 = ncol + 1
      nrpna = nrow + naugrow
      ncpna = ncol + naugcol
      if (iadd .eq. 0) goto 9
      do 1 i = 1, nrow
      ij1 = index(i,1,ncol)
      ij2 = index(i,ncol,ncol)
      discon = (f(ij2) - f(ij1)) / (naugcol + 1)
      do 2 j = 1, ncol
      ij3 = index(i,j,ncol)
      ij4 = index(i,j,ncpna)
    2 temp(ij4) = f(ij3)
      do 3 j = ncp1, ncpna
      ij4 = index(i,j,ncpna)
    3 temp(ij4) = f(ij2) - (discon * (j - ncol))
    1 continue
      do 4 j = 1, ncpna
      ij1 = index(1,j,ncpna)
      ij2 = index(nrow,j,ncpna)
      discon = (temp(ij2) - temp(ij1)) / (naugrow + 1)
      do 5 i = nrp1, nrpna
      ij = index(i,j,ncpna)
    5 temp(ij) = temp(ij2) - (discon * (i - nrow))
    4 continue
      do 6 i = 1, nrpna
      do 6 j = 1, ncpna
      ij = index(i,j,ncpna)
    6 f(ij) = temp(ij)
c
      return 
    9 continue
      do 7 i = 1, nrow
      do 7 j = 1, ncol
      ij1 = index(i,j,ncol)
      ij2 = index(i,j,ncpna)
    7 temp(ij1) = f(ij2)
      do 8 i = 1, nrow
      do 8 j = 1, ncol
      ij = index(i,j,ncol)
    8 f(ij) = temp(ij)
      return 
      end
      subroutine fftsub(work1, work2, nx, ny, isign)
      parameter(isza=6400,iszs=80)
      dimension work1(ny, nx), work2(ny, nx)
      common /inter3/ fcmplx(isza)
c
c.......................................................................
c
c  --     This routine splits a complex array fcmplx into its real and
c      imaginary parts, calls Singleton's subroutine fft, and puts the
c      results back into fcmplx
c
c.......................................................................
c.
c
      complex fcmplx, cmplx
      index(i,j,n) = ((i - 1) * n) + j
      ntotal = nx * ny
      do 1 i = 1, nx
      do 1 j = 1, ny
      ij = index(i,j,ny)
      work1(j,i) = real(fcmplx(ij))
      work2(j,i) = aimag(fcmplx(ij))
    1 continue
      call fft(work1, work2, ntotal, ny, ny, isign)
      call fft(work1, work2, ntotal, nx, ntotal, isign)
      do 2 i = 1, nx
      do 2 j = 1, ny
      ij = index(i,j,ny)
      fcmplx(ij) = cmplx(work1(j,i),work2(j,i))
    2 continue
      return 
      end
c ********************* start of fft *************************
c                                  up to 2/4/73
c
c  ** for real data, use together with realtr, a separate routine
c      -see realtr for details
c
c  ** does not run under watfiv
c
c  ** as an example of how to call this program,  a test driver
c          is available in a separate library
c
c  multivariate complex fourier transform, computed in place
c    using mixed-radix fast fourier transform algorithm.
c  by r. c. singleton, stanford research institute, sept. 1968
c  arrays a and b originally hold the real and imaginary
c    components of the data, and return the real and
c    imaginary components of the resulting fourier coefficients.
c  multivariate data is indexed according to the fortran
c    array element successor function, without limit
c    on the number of implied multiple subscripts.
c    the subroutine is called once for each variate.
c    the calls for a multivariate transform may be in any order.
c  ntot is the total number of complex data values.
c  n is the dimension of the current variable.
c  nspan/n is the spacing of consecutive data values
c    while indexing the current variable.
c  the sign of isn determines the sign of the complex
c    exponential, and the magnitude of isn is normally one.
c  a tri-variate transform with a(n1,n2,n3), b(n1,n2,n3)
c    is computed by
c      call fft(a,b,n1*n2*n3,n1,n1,1)
c      call fft(a,b,n1*n2*n3,n2,n1*n2,1)
c      call fft(a,b,n1*n2*n3,n3,n1*n2*n3,1)
c  for a single-variate transform,
c    ntot = n = nspan = (number of complex data values), e.g.
c      call fft(a,b,n,n,n,1)
c  the data can alternatively be stored in a single complex array c
c   in standard fortran fashion, i.e., alternating real and
c   imaginary parts.  then with most fortran compilers, the
c   complex array c can be equivalenced to a real array a, the
c   magnitude of isn changed to two to give the correct indexing
c   increment, and a and a(2) used to pass the initial addresses
c   for the sequences of real and imaginary values, e.g.,
c      complex c(ntot)
c      real a(2*ntot)
c      equivalence (c(1),a(1))
c      call fft(a,a(2),ntot,n,nspan,2)
c  arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
c    are used for temporary storage.  if the available storage
c    is insufficient, the program is terminated by a stop.
c    maxf must be .ge. the maximum prime factor of n.
c    maxp must be .gt. the number of prime factors of n.
c    in addition, if the square-free portion k of n has two or
c    more prime factors, then maxp must be .ge. k-1.
      subroutine fft(a, b, ntot, n, nspan, isn)
c  array storage in nfac for a maximum of 15 prime factors of n.
c  if n has more than one square-free factor, the product of the
c    square-free factors must be .le. 210
      dimension a(1), b(1)
c  array storage for maximum prime factor of 91
      dimension nfac(11), np(209)
      dimension at(91), ck(91), bt(91), sk(91)
c  the following two constants should agree with the array dimensions.
      equivalence (ii, i)
      maxf = 91
      maxp = 209
      if (n .lt. 2) return 
      inc = isn
      c72 = 0.30901699437494742d0
      s72 = 0.95105651629515357d0
      s120 = 0.86602540378443865d0
      rad = 6.2831853071796d0
      if (isn .ge. 0) goto 10
      s72 = - s72
      s120 = - s120
      rad = - rad
      inc = - inc
   10 nt = inc * ntot
      ks = inc * nspan
      kspan = ks
      nn = nt - inc
      jc = ks / n
      radf = (rad * float(jc)) * 0.5
      i = 0
c  determine the factors of n
      jf = 0
      m = 0
      k = n
      goto 20
   15 m = m + 1
      nfac(m) = 4
      k = k / 16
   20 if ((k - ((k / 16) * 16)) .eq. 0) goto 15
      j = 3
      jj = 9
      goto 30
   25 m = m + 1
      nfac(m) = j
      k = k / jj
   30 if (mod(k,jj) .eq. 0) goto 25
      j = j + 2
      jj = j ** 2
      if (jj .le. k) goto 30
      if (k .gt. 4) goto 40
      kt = m
      nfac(m + 1) = k
      if (k .ne. 1) m = m + 1
      goto 80
   40 if ((k - ((k / 4) * 4)) .ne. 0) goto 50
      m = m + 1
      nfac(m) = 2
      k = k / 4
   50 kt = m
      j = 2
   60 if (mod(k,j) .ne. 0) goto 70
      m = m + 1
      nfac(m) = j
      k = k / j
   70 j = (((j + 1) / 2) * 2) + 1
      if (j .le. k) goto 60
   80 if (kt .eq. 0) goto 100
      j = kt
   90 m = m + 1
      nfac(m) = nfac(j)
      j = j - 1
c  compute fourier transform
      if (j .ne. 0) goto 90
  100 sd = radf / float(kspan)
      cd = 2.0 * (sin(sd) ** 2)
      sd = sin(sd + sd)
      kk = 1
      i = i + 1
c  transform for factor of 2 (including rotation factor)
      if (nfac(i) .ne. 2) goto 400
      kspan = kspan / 2
      k1 = kspan + 2
  210 k2 = kk + kspan
      ak = a(k2)
      bk = b(k2)
      a(k2) = a(kk) - ak
      b(k2) = b(kk) - bk
      a(kk) = a(kk) + ak
      b(kk) = b(kk) + bk
      kk = k2 + kspan
      if (kk .le. nn) goto 210
      kk = kk - nn
      if (kk .le. jc) goto 210
      if (kk .gt. kspan) goto 800
  220 c1 = 1.0 - cd
      s1 = sd
  230 k2 = kk + kspan
      ak = a(kk) - a(k2)
      bk = b(kk) - b(k2)
      a(kk) = a(kk) + a(k2)
      b(kk) = b(kk) + b(k2)
      a(k2) = (c1 * ak) - (s1 * bk)
      b(k2) = (s1 * ak) + (c1 * bk)
      kk = k2 + kspan
      if (kk .lt. nt) goto 230
      k2 = kk - nt
      c1 = - c1
      kk = k1 - k2
      if (kk .gt. k2) goto 230
      ak = c1 - ((cd * c1) + (sd * s1))
      s1 = ((sd * c1) - (cd * s1)) + s1
      c1 = 2.0 - ((ak ** 2) + (s1 ** 2))
      s1 = c1 * s1
      c1 = c1 * ak
      kk = kk + jc
      if (kk .lt. k2) goto 230
      k1 = (k1 + inc) + inc
      kk = ((k1 - kspan) / 2) + jc
      if (kk .le. (jc + jc)) goto 220
c  transform for factor of 3 (optional code)
      goto 100
  320 k1 = kk + kspan
      k2 = k1 + kspan
      ak = a(kk)
      bk = b(kk)
      aj = a(k1) + a(k2)
      bj = b(k1) + b(k2)
      a(kk) = ak + aj
      b(kk) = bk + bj
      ak = (- (0.5 * aj)) + ak
      bk = (- (0.5 * bj)) + bk
      aj = (a(k1) - a(k2)) * s120
      bj = (b(k1) - b(k2)) * s120
      a(k1) = ak - bj
      b(k1) = bk + aj
      a(k2) = ak + bj
      b(k2) = bk - aj
      kk = k2 + kspan
      if (kk .lt. nn) goto 320
      kk = kk - nn
      if (kk .le. kspan) goto 320
c  transform for factor of 4
      goto 700
  400 if (nfac(i) .ne. 4) goto 600
      kspnn = kspan
      kspan = kspan / 4
  410 c1 = 1.0
      s1 = 0
  420 k1 = kk + kspan
      k2 = k1 + kspan
      k3 = k2 + kspan
      akp = a(kk) + a(k2)
      akm = a(kk) - a(k2)
      ajp = a(k1) + a(k3)
      ajm = a(k1) - a(k3)
      a(kk) = akp + ajp
      ajp = akp - ajp
      bkp = b(kk) + b(k2)
      bkm = b(kk) - b(k2)
      bjp = b(k1) + b(k3)
      bjm = b(k1) - b(k3)
      b(kk) = bkp + bjp
      bjp = bkp - bjp
      if (isn .lt. 0) goto 450
      akp = akm - bjm
      akm = akm + bjm
      bkp = bkm + ajm
      bkm = bkm - ajm
      if (s1 .eq. 0) goto 460
  430 a(k1) = (akp * c1) - (bkp * s1)
      b(k1) = (akp * s1) + (bkp * c1)
      a(k2) = (ajp * c2) - (bjp * s2)
      b(k2) = (ajp * s2) + (bjp * c2)
      a(k3) = (akm * c3) - (bkm * s3)
      b(k3) = (akm * s3) + (bkm * c3)
      kk = k3 + kspan
      if (kk .le. nt) goto 420
  440 c2 = c1 - ((cd * c1) + (sd * s1))
      s1 = ((sd * c1) - (cd * s1)) + s1
      c1 = 2.0 - ((c2 ** 2) + (s1 ** 2))
      s1 = c1 * s1
      c1 = c1 * c2
      c2 = (c1 ** 2) - (s1 ** 2)
      s2 = (2.0 * c1) * s1
      c3 = (c2 * c1) - (s2 * s1)
      s3 = (c2 * s1) + (s2 * c1)
      kk = (kk - nt) + jc
      if (kk .le. kspan) goto 420
      kk = (kk - kspan) + inc
      if (kk .le. jc) goto 410
      if (kspan .eq. jc) goto 800
      goto 100
  450 akp = akm + bjm
      akm = akm - bjm
      bkp = bkm - ajm
      bkm = bkm + ajm
      if (s1 .ne. 0) goto 430
  460 a(k1) = akp
      b(k1) = bkp
      a(k2) = ajp
      b(k2) = bjp
      a(k3) = akm
      b(k3) = bkm
      kk = k3 + kspan
      if (kk .le. nt) goto 420
c  transform for factor of 5 (optional code)
      goto 440
  510 c2 = (c72 ** 2) - (s72 ** 2)
      s2 = (2.0 * c72) * s72
  520 k1 = kk + kspan
      k2 = k1 + kspan
      k3 = k2 + kspan
      k4 = k3 + kspan
      akp = a(k1) + a(k4)
      akm = a(k1) - a(k4)
      bkp = b(k1) + b(k4)
      bkm = b(k1) - b(k4)
      ajp = a(k2) + a(k3)
      ajm = a(k2) - a(k3)
      bjp = b(k2) + b(k3)
      bjm = b(k2) - b(k3)
      aa = a(kk)
      bb = b(kk)
      a(kk) = (aa + akp) + ajp
      b(kk) = (bb + bkp) + bjp
      ak = ((akp * c72) + (ajp * c2)) + aa
      bk = ((bkp * c72) + (bjp * c2)) + bb
      aj = (akm * s72) + (ajm * s2)
      bj = (bkm * s72) + (bjm * s2)
      a(k1) = ak - bj
      a(k4) = ak + bj
      b(k1) = bk + aj
      b(k4) = bk - aj
      ak = ((akp * c2) + (ajp * c72)) + aa
      bk = ((bkp * c2) + (bjp * c72)) + bb
      aj = (akm * s2) - (ajm * s72)
      bj = (bkm * s2) - (bjm * s72)
      a(k2) = ak - bj
      a(k3) = ak + bj
      b(k2) = bk + aj
      b(k3) = bk - aj
      kk = k4 + kspan
      if (kk .lt. nn) goto 520
      kk = kk - nn
      if (kk .le. kspan) goto 520
c  transform for odd factors
      goto 700
  600 k = nfac(i)
      kspnn = kspan
      kspan = kspan / k
      if (k .eq. 3) goto 320
      if (k .eq. 5) goto 510
      if (k .eq. jf) goto 640
      jf = k
      s1 = rad / float(k)
      c1 = cos(s1)
      s1 = sin(s1)
      if (jf .gt. maxf) goto 998
      ck(jf) = 1.0
      sk(jf) = 0.0
      j = 1
  630 ck(j) = (ck(k) * c1) + (sk(k) * s1)
      sk(j) = (ck(k) * s1) - (sk(k) * c1)
      k = k - 1
      ck(k) = ck(j)
      sk(k) = - sk(j)
      j = j + 1
      if (j .lt. k) goto 630
  640 k1 = kk
      k2 = kk + kspnn
      aa = a(kk)
      bb = b(kk)
      ak = aa
      bk = bb
      j = 1
      k1 = k1 + kspan
  650 k2 = k2 - kspan
      j = j + 1
      at(j) = a(k1) + a(k2)
      ak = at(j) + ak
      bt(j) = b(k1) + b(k2)
      bk = bt(j) + bk
      j = j + 1
      at(j) = a(k1) - a(k2)
      bt(j) = b(k1) - b(k2)
      k1 = k1 + kspan
      if (k1 .lt. k2) goto 650
      a(kk) = ak
      b(kk) = bk
      k1 = kk
      k2 = kk + kspnn
      j = 1
  660 k1 = k1 + kspan
      k2 = k2 - kspan
      jj = j
      ak = aa
      bk = bb
      aj = 0.0
      bj = 0.0
      k = 1
  670 k = k + 1
      ak = (at(k) * ck(jj)) + ak
      bk = (bt(k) * ck(jj)) + bk
      k = k + 1
      aj = (at(k) * sk(jj)) + aj
      bj = (bt(k) * sk(jj)) + bj
      jj = jj + j
      if (jj .gt. jf) jj = jj - jf
      if (k .lt. jf) goto 670
      k = jf - j
      a(k1) = ak - bj
      b(k1) = bk + aj
      a(k2) = ak + bj
      b(k2) = bk - aj
      j = j + 1
      if (j .lt. k) goto 660
      kk = kk + kspnn
      if (kk .le. nn) goto 640
      kk = kk - nn
c  multiply by rotation factor (except for factors of 2 and 4)
      if (kk .le. kspan) goto 640
  700 if (i .eq. m) goto 800
      kk = jc + 1
  710 c2 = 1.0 - cd
      s1 = sd
  720 c1 = c2
      s2 = s1
      kk = kk + kspan
  730 ak = a(kk)
      a(kk) = (c2 * ak) - (s2 * b(kk))
      b(kk) = (s2 * ak) + (c2 * b(kk))
      kk = kk + kspnn
      if (kk .le. nt) goto 730
      ak = s1 * s2
      s2 = (s1 * c2) + (c1 * s2)
      c2 = (c1 * c2) - ak
      kk = (kk - nt) + kspan
      if (kk .le. kspnn) goto 730
      c2 = c1 - ((cd * c1) + (sd * s1))
      s1 = s1 + ((sd * c1) - (cd * s1))
      c1 = 2.0 - ((c2 ** 2) + (s1 ** 2))
      s1 = c1 * s1
      c2 = c1 * c2
      kk = (kk - kspnn) + jc
      if (kk .le. kspan) goto 720
      kk = ((kk - kspan) + jc) + inc
      if (kk .le. (jc + jc)) goto 710
c  permute the results to normal order---done in two stages
c  permutation for square factors of n
      goto 100
  800 np(1) = ks
      if (kt .eq. 0) goto 890
      k = (kt + kt) + 1
      if (m .lt. k) k = k - 1
      j = 1
      np(k + 1) = jc
  810 np(j + 1) = np(j) / nfac(j)
      np(k) = np(k + 1) * nfac(j)
      j = j + 1
      k = k - 1
      if (j .lt. k) goto 810
      k3 = np(k + 1)
      kspan = np(2)
      kk = jc + 1
      k2 = kspan + 1
      j = 1
c  permutation for single-variate transform (optional code)
      if (n .ne. ntot) goto 850
  820 ak = a(kk)
      a(kk) = a(k2)
      a(k2) = ak
      bk = b(kk)
      b(kk) = b(k2)
      b(k2) = bk
      kk = kk + inc
      k2 = kspan + k2
      if (k2 .lt. ks) goto 820
  830 k2 = k2 - np(j)
      j = j + 1
      k2 = np(j + 1) + k2
      if (k2 .gt. np(j)) goto 830
      j = 1
  840 if (kk .lt. k2) goto 820
      kk = kk + inc
      k2 = kspan + k2
      if (k2 .lt. ks) goto 840
      if (kk .lt. ks) goto 830
      jc = k3
c  permutation for multivariate transform
      goto 890
  850 k = kk + jc
  860 ak = a(kk)
      a(kk) = a(k2)
      a(k2) = ak
      bk = b(kk)
      b(kk) = b(k2)
      b(k2) = bk
      kk = kk + inc
      k2 = k2 + inc
      if (kk .lt. k) goto 860
      kk = (kk + ks) - jc
      k2 = (k2 + ks) - jc
      if (kk .lt. nt) goto 850
      k2 = (k2 - nt) + kspan
      kk = (kk - nt) + jc
      if (k2 .lt. ks) goto 850
  870 k2 = k2 - np(j)
      j = j + 1
      k2 = np(j + 1) + k2
      if (k2 .gt. np(j)) goto 870
      j = 1
  880 if (kk .lt. k2) goto 850
      kk = kk + jc
      k2 = kspan + k2
      if (k2 .lt. ks) goto 880
      if (kk .lt. ks) goto 870
      jc = k3
  890 if (((2 * kt) + 1) .ge. m) return 
c  permutation for square-free factors of n
      kspnn = np(kt + 1)
      j = m - kt
      nfac(j + 1) = 1
  900 nfac(j) = nfac(j) * nfac(j + 1)
      j = j - 1
      if (j .ne. kt) goto 900
      kt = kt + 1
      nn = nfac(kt) - 1
      if (nn .gt. maxp) goto 998
      jj = 0
      j = 0
      goto 906
  902 jj = jj - k2
      k2 = kk
      k = k + 1
      kk = nfac(k)
  904 jj = kk + jj
      if (jj .ge. k2) goto 902
      np(j) = jj
  906 k2 = nfac(kt)
      k = kt + 1
      kk = nfac(k)
      j = j + 1
c  determine the permutation cycles of length greater than 1
      if (j .le. nn) goto 904
      j = 0
      goto 914
  910 k = kk
      kk = np(k)
      np(k) = - kk
      if (kk .ne. j) goto 910
      k3 = kk
  914 j = j + 1
      kk = np(j)
      if (kk .lt. 0) goto 914
      if (kk .ne. j) goto 910
      np(j) = - j
      if (j .ne. nn) goto 914
c  reorder a and b, following the permutation cycles
      maxf = inc * maxf
      goto 950
  924 j = j - 1
      if (np(j) .lt. 0) goto 924
      jj = jc
  926 kspan = jj
      if (jj .gt. maxf) kspan = maxf
      jj = jj - kspan
      k = np(j)
      kk = ((jc * k) + ii) + jj
      k1 = kk + kspan
      k2 = 0
  928 k2 = k2 + 1
      at(k2) = a(k1)
      bt(k2) = b(k1)
      k1 = k1 - inc
      if (k1 .ne. kk) goto 928
  932 k1 = kk + kspan
      k2 = k1 - (jc * (k + np(k)))
      k = - np(k)
  936 a(k1) = a(k2)
      b(k1) = b(k2)
      k1 = k1 - inc
      k2 = k2 - inc
      if (k1 .ne. kk) goto 936
      kk = k2
      if (k .ne. j) goto 932
      k1 = kk + kspan
      k2 = 0
  940 k2 = k2 + 1
      a(k1) = at(k2)
      b(k1) = bt(k2)
      k1 = k1 - inc
      if (k1 .ne. kk) goto 940
      if (jj .ne. 0) goto 926
      if (j .ne. 1) goto 924
  950 j = k3 + 1
      nt = nt - kspnn
      ii = (nt - inc) + 1
      if (nt .ge. 0) goto 924
c  error finish, insufficient array storage
      return 
  998 isn = 0
      write(unit=*, fmt=999) 
      stop 
  999 format(44h0ARRAY BOUNDS EXCEEDED WITHIN SUBROUTINE FFT)
      end
