      PROGRAM FIELDS

C     Compute pressure from Green's function.

C     The transform parameters are chosen to satisfy certain sampling
C     requirements.

C     It is assumed that the user chose the maximum possible 
C     spacing in computing g(k) and so deltak is not
C     allowed to be any larger.

C     The other part of the kernel is the exp(ikr) term.
C     this must be sampled on the k-axis finely enough to have
C     about 6 pts per wavelength.

C     DELTAR implies KMAX. KMAX is the upper limit of integration
C     and beyond the last point computed.  This assumes
C     the user had chosen KMAX as small as possible while still
C     covering the support of G(k).

C     Note 'WNO' is used for original data and 'K' for interpolated data

      INTEGER FLPFIL, PRTFIL
      PARAMETER ( FLPFIL = 5, PRTFIL = 6, 
     &   PI = 3.14159265, MAXNT = 500000, MAXNWN =  50000,
     &   MAXNSD = 501, MAXNRD = 101 )

      REAL WNO( MAXNWN ), R( MAXNT ),
     &   WNOIN( MAXNT ), RD( MAXNRD ), SD( MAXNSD ), ATTEN, FREQ
      COMPLEX G( MAXNWN ), GIN( MAXNT ), P( MAXNT ), TEMP( 2 * MAXNT )
      CHARACTER PLTITL*80, OPTION*80

C     *** Begin by reading in the data ***

      OPTION = 'RLPA'            ! default option
      READ( FLPFIL, * ) OPTION
      READ( FLPFIL, * ) RMINKM, RMAXKM, NR

      RMIN = 1000.0 * RMINKM
      RMAX = 1000.0 * RMAXKM

C     --- for a point source, take us a bit away from the source
      IF ( OPTION(1:1) .EQ. 'R' .AND. RMIN .LE. RMAX / 100000 ) THEN
         RMIN = RMAX / 100000
         WRITE( PRTFIL, * ) 'RMIN bumped to ', RMIN
      END IF

      DELTAR = ( RMAX - RMIN ) / ( NR - 1 )

C     *** Read the header records from GRNFIL ***

      CALL RDHEAD( PLTITL, SD, NSD, MAXNSD, RD, NRD, MAXNRD, 
     &   WNO, NWNO, MAXNWN, NKCORE, DELTAK, ATTEN, FREQ, NRECSR )

C     Set up for transform: need DELWNO, NT
C     DELTAK is what scooter used; DELWNO is for interpolation

C     --- If DELTAR is too big, take submultiple

      WNOMAX = WNO( 1 ) + 2.0 * PI / DELTAR

      IF ( WNOMAX .LT. WNO( NWNO ) ) THEN
         IRATIO = ( WNO( NWNO ) - WNO( 1 ) ) / ( WNOMAX - WNO( 1 ) ) + 1
         DELTAR = DELTAR / IRATIO
         NR     = IRATIO * ( NR - 1 ) + 1
         WNOMAX = WNO( 1 ) + 2.0 * PI / DELTAR
         WRITE( PRTFIL, * ) 'NR bumped to ', NR
      END IF

      NT = NR

C     Compute NT based on DELTAK (= 1/RMAX) sampling requirement

      NT2 = RMAX * ( WNOMAX - WNO( 1 ) )       ! DELTAK limit

      IF ( NT2 .GT. NT ) WRITE( PRTFIL, * ) 'NT bumped', NT, NT2
      NT = MAX( NT, NT2 )

C     Parameters when N must be a power of 2

      NT = 2 **( INT( LOG10( REAL( NT ) ) / 0.301 ) + 1 )
      DELWNO = 2.0 * PI / ( NT * DELTAR )

      WRITE( PRTFIL, * ) 'NT  used               = ', NT
      WRITE( PRTFIL, * ) 'NR  used               = ', NR
      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'DELTAK = ', DELTAK
      WRITE( PRTFIL, * ) 'DELWNO = ', DELWNO
      WRITE( PRTFIL, * ) 'DELTAR = ', DELTAR
      WRITE( PRTFIL, * )

      IF ( NT .GT. MAXNT ) THEN
         CALL ERROUT( PRTFIL, 'F', 'FIELDS',
     &             'Size of tranform exceeds array size' )
      ENDIF

C     Compute optimal stabilizing attenuation.
C     ATTEN is the value used by SCOOTER
C     ATTINT is the optimal value for the interpolated k-values

      ATTINT = ATTEN

      IF ( PLTITL(1:5) .EQ. 'SPARC' ) THEN
         WRITE( PRTFIL, * ) 'SPARC RUN: ATTINT SET TO ZERO'
         ATTINT = 0.0
         IF ( ATTEN .NE. 0.0 ) THEN
            CALL ERROUT( PRTFIL, 'F', 'FIELDS',
     &             'Stabilizing attenuation must vanish' )
         ENDIF
      ENDIF

C     --- Set up vector of range points
      DO J = 1, NR
         R( J ) = RMIN + ( J - 1 ) * DELTAR
      END DO

C     Construct shade file

      CALL SHADE( WNO, DELWNO, ATTEN, NWNO, G, NT, WNOIN, GIN,
     &      R, NR, RMIN, DELTAR, SD, NSD, RD, NRD, 
     &      PLTITL, FREQ, ATTINT, NRECSR, NKCORE,
     &      OPTION, P, TEMP )

      STOP
      END
C**********************************************************************C
      SUBROUTINE RDHEAD( PLTITL, SD, NSD, MAXNSD,
     &   RD, NRD, MAXNRD, 
     &   RK, NKPTS, MAXNWN, NKCORE, DELTAK, ATTEN, FREQ, NRECSR )

C     Routine to read header from disk file

      INTEGER GRNFIL, PRTFIL
      PARAMETER ( GRNFIL = 20, PRTFIL = 6 )

      REAL RK( * ), RD( * ), SD( * )
      CHARACTER PLTITL*80, PLTTYP*10

C     Double open to get record length

      OPEN( FILE = 'GRNFIL', UNIT = GRNFIL, STATUS = 'OLD',
     &      RECL = 20, ACCESS = 'DIRECT', FORM = 'UNFORMATTED' )
      READ( GRNFIL, REC = 1 ) LRECL
      CLOSE( GRNFIL )

      OPEN( FILE = 'GRNFIL', UNIT = GRNFIL, STATUS = 'OLD',
     &   RECL = 4 * LRECL, ACCESS = 'DIRECT', FORM = 'UNFORMATTED' )

C     Read data

      READ( GRNFIL, REC = 1 ) LRECL, PLTITL
      READ( GRNFIL, REC = 2 ) PLTTYP, XS, YS, THETA
      READ( GRNFIL, REC = 3 ) FREQ, NSD, NRD, NKPTS, NKCORE,
     &                        DELTAK, ATTEN, NRECSR

      IF ( NSD .GT. MAXNSD ) CALL ERROUT( PRTFIL, 'F', 'FIELDS:RDHEAD',
     &             'Too many source depths' )

      IF ( NRD .GT. MAXNRD ) CALL ERROUT( PRTFIL, 'F', 'FIELDS:RDHEAD',
     &             'Too many receiver depths' )
 
      READ( GRNFIL, REC = 4 ) ( SD( I ), I = 1, NSD )
      READ( GRNFIL, REC = 5 ) ( RD( J ), J = 1, NRD )

C     Read in the vector of k-space points

      IF ( NRECSR * NKCORE .GT. MAXNWN )
     &   CALL ERROUT( PRTFIL, 'F', 'FIELDS:RDHEAD',
     &             'Too many k-space points' )

      DO IREC = 1, NRECSR
         ISTART = ( IREC - 1 ) * NKCORE
         READ( GRNFIL, REC = 5 + IREC )
     &      ( RK( ISTART + I ), I = 1, NKCORE )
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SHADE( WNO, DELWNO, ATTEN, NWNO, G, NT, WNOIN, GIN,
     &   R, NR, RMIN, DELTAR, SD, NSD, RD, NRD, 
     &   PLTITL, FREQ, ATTINT, NRECSR, NKCORE,
     &   OPTION,  P, TEMP )

C     Performs the transforms to convert the Green's function file
C     to a shade file

C     Expects
C        WNO, NWNO: The wavenumber data
C        SD, RD: Source depth, receiver depth data
C     Returns
C        Nothing

      INTEGER PRTFIL, GRNFIL
      PARAMETER ( PRTFIL = 6, GRNFIL = 20, MAXNR = 2045 )

      LOGICAL INIFLG
      REAL R( * ), WNOIN( * ), SD( * ), RD( * ), WNO( * )
      COMPLEX G( * ), GIN( * ), P( * ), TEMP( * )
      CHARACTER*80 PLTITL, OPTION

C     Construct the shade file header

      INIFLG = .TRUE.
      IRECOU = 6

C     If necessary, reduce range sampling to fit in the record

      NRSMPL = NR

      IF ( NR .GT. MAXNR ) THEN

         WRITE( PRTFIL, * )
     &      'Range axis subsampled to fit max record size'

         IRSKIP = NR / ( MAXNR + 1 ) + 1
         NRSMPL = NR / IRSKIP

         DO IR = 1, NRSMPL
            R( IR ) = R( ( IR - 1 ) * IRSKIP + 1 )
         END DO
         WRITE( PRTFIL, * ) 'DELTAR = ', IRSKIP * DELTAR
         WRITE( PRTFIL, * )
      ENDIF
         
      CALL HEADER( PLTITL, SD, NSD, RD, NRD, 
     &   R, NRSMPL, FREQ, ATTEN, '     ', 0.0, 0.0, 0.0 )

C     *** Loop over all source/rcvr combos ***

      DO ISD = 1, NSD

         DO IRD  = 1, NRD

            ISR    = ( ISD - 1 ) * NRD + IRD   ! Index of source/receiver
            IRECIN = 5 + ISR * NRECSR

            DO IREC = 1, NRECSR
               ISTART = ( IREC - 1 ) * NKCORE
               READ( GRNFIL, REC = IRECIN + IREC )
     &                       ( G( ISTART + I ), I = 1, NKCORE )
            END DO

            WRITE( PRTFIL, * ) 'SD, RD = ', SD( ISD ), RD( IRD )

C           *** Interpolate ***

            CALL INTERP( WNO, DELWNO, ATTEN, ATTINT, G, NWNO,
     &                   WNOIN, GIN, NT, OPTION )

C           *** Evaluate either Fourier or Hankel transform ***

            IF ( OPTION(1:1) .EQ. 'X' ) THEN
               CALL FTS( NT, WNOIN( 1 ), DELWNO, ATTINT, 
     &            RMIN, NR, DELTAR, GIN, TEMP, P, INIFLG )
            ELSE
               CALL HTS( NT, WNOIN( 1 ), DELWNO, ATTINT, 
     &            RMIN, NR, DELTAR, GIN, TEMP, P, INIFLG, OPTION )
            ENDIF

            INIFLG = .FALSE.

C           If necessary, reduce range sampling to fit in the record

            IF ( NR .GT. MAXNR ) THEN
               DO IR = 1, NRSMPL
                  P( IR ) = P( IRSKIP * IR + 1 )
               END DO
            ENDIF

C           *** Write out the field ***

            CALL WRTFLD( P, 1, NRSMPL, IRECOU )

         END DO   ! next receiver depth
      END DO   ! next source depth

      RETURN
      END
C**********************************************************************C
      SUBROUTINE INTERP( WNO, DELWNO, ATTEN, ATTINT, G, NWNO,
     &   WNOIN, GIN, NT, OPTION )

C     Produces an evenly sampled kernel from input G

      INTEGER PRTFIL
      COMPLEX CI
      PARAMETER ( PRTFIL = 6, CI = ( 0.0, 1.0 ) )

      REAL WNO( * ), WNOIN( * )
      COMPLEX G( * ), GIN( * ), WNOT2, X( 9 ), F( 9 ), PADE
      CHARACTER ERRMSG*80, OPTION*80

C     ------ Initialize interpolation data
      ISIZE  = 3
      ICTR   = 2
      IRIGHT = 3
      NEW    = ISIZE

      DO I = 1, ISIZE
         X( I ) = ( WNO( I ) + CI * ATTEN )**2
         F( I ) = G( I )
      END DO

C     *** Main loop: take steps of delwno along the k-axis ***

      DO I = 1, NT
         WNOIN( I ) = WNO( 1 ) + ( I - 1 ) * DELWNO

C        --- Desirable/possible to advance interpolation window?
 5000    IF ( WNOIN( I ) .GT. WNO( ICTR ) ) THEN
            IF ( IRIGHT .LT. NWNO ) THEN
               ICTR   = ICTR   + 1
               IRIGHT = IRIGHT + 1
               NEW    = MOD( NEW, ISIZE ) + 1

               X( NEW ) = ( WNO( IRIGHT ) + CI * ATTEN )**2
               F( NEW ) = G( IRIGHT )
C              ------ Check again
               GOTO 5000
            ENDIF
         ENDIF

C        --- Interpolate or zero fill
         IF ( WNOIN( I ) .LE. WNO( NWNO ) ) THEN
            WNOT2 = ( WNOIN( I ) + CI * ATTINT )**2
            IF ( OPTION(4:4) .EQ. 'O' ) THEN
               GIN( I ) = PC(   WNOT2, X, F, ISIZE, ERRMSG )   ! polynomial
            ELSE
               GIN( I ) = PADE( WNOT2, X, F, ISIZE, ERRMSG )   ! Pade
            ENDIF

            IF ( ERRMSG .NE. ' ' ) WRITE( PRTFIL, * ) ERRMSG
         ELSE
            GIN( I ) = 0.0
         ENDIF

      END DO

      RETURN
      END
