      PROGRAM FIELD

C     Generate file of replica vectors

      INTEGER PRTFIL, FLPFIL
      PARAMETER ( FLPFIL = 5, PRTFIL = 6,
     &   MAXM = 1500, MAXNR = 9000, MAXNRZ = 210000, 
     &   MAXNSD = 51, MAXNRD = 201, MAXNPROF = 502 )

      REAL      SD( MAXNSD ), RD( MAXNRD ), RR( MAXNRD ),
     &          R( MAXNR ), RPROF( MAXNPROF )
      COMPLEX   CK( MAXM ), PHIS( MAXM, MAXNSD ),
     &          PHIR( MAXM, MAXNRD ), P( MAXNRZ )
      CHARACTER OPT*50, SHTITL*80, TITLE*80, COMP*1

C     Read source, rcvr depths and rcvr displacements

      SHTITL(1:1) = '$'

      READ( FLPFIL, * ) SHTITL
      READ( FLPFIL, * ) OPT
      READ( FLPFIL, * ) MLIMIT

      COMP   = OPT( 3:3 )
      MLIMIT = MIN( MAXM, MLIMIT )

C     *** Read profile ranges ***

      RPROF( 3 ) = -999.9
      READ( FLPFIL, * ) NPROF, ( RPROF( I ), I = 1, NPROF )

C     ( an extra profile is used at `infinite' range )
      IF ( NPROF .GT. MAXNPROF - 1 ) THEN
         WRITE( PRTFIL, * ) 'NPROF, MAXNPROF = ', NPROF, MAXNPROF
         CALL ERROUT( PRTFIL, 'F', 'FIELD', 'Too many profiles' )
      ENDIF

C     EVALAD needs a profile at zero range

      IF ( RPROF( 1 ) .NE. 0.0 )
     &   CALL ERROUT( PRTFIL, 'F', 'FIELD',
     &   'The first profile must be at 0 km' )

      CALL SUBTAB( RPROF, NPROF )

C     *** Read receiver ranges ***

      CALL RANGES( FLPFIL, PRTFIL, MAXNR, NR, R )

C     *** Read source/receiver depths ***

      ZMIN = -1.0E10
      ZMAX = +1.0E10

      CALL SDRD( FLPFIL, PRTFIL, ZMIN, ZMAX,
     &   MAXNSD, NSD, SD,
     &   MAXNRD, NRD, RD )

C     *** Read receiver ranges ***

      NRR     = 0
      RR( 2 ) = 0.0
      RR( 3 ) = -999.9
      READ( FLPFIL, * ) NRR, ( RR( IR ), IR = 1, NRR )

      IF ( NRR .NE. NRD ) THEN
         WRITE( PRTFIL, * ) 'NRR, NRD = ', NRR, NRD
         CALL ERROUT( PRTFIL, 'W', 'FIELD', 'NRR being set to NRD' )
         NRR = NRD
      ENDIF

      CALL SUBTAB( RR, NRD )
      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Number of receivers = ', NRD
      IF ( NRD .GE. 1 )
     &   WRITE( PRTFIL, * ) ( RR( IR ), IR = 1, MIN( NRD, 51 ) )

      WRITE( PRTFIL, * )

      IF ( NR * NRD .GT. MAXNRZ )
     &   CALL ERROUT( PRTFIL, 'F', 'FIELD',
     &                'Number of ranges x Depths too large' )

C     *** Read in modes ***

      IPROF = 1
      CALL GETMOD( IPROF, ' ', MAXM, SD, NSD, 'N',
     &   CK, PHIS, MSRC, FREQ, TITLE )
      CALL GETMOD( IPROF, ' ', MAXM, RD, NRD, COMP,
     &   CK, PHIR, MSRC, FREQ, TITLE )

C     ------ Generate header
      IF ( SHTITL(1:1) .EQ. '$' ) SHTITL = TITLE
      CALL HEADER( SHTITL, SD, NSD, RD, NRD,
     &   R, NR, FREQ, ATTEN, '     ', 0.0, 0.0, 0.0 )
      IREC = 6

C     *** MAIN LOOP: For each source evaluate and write the field ***

      DO J = 1, NSD
         M = MIN( MLIMIT, MSRC )   ! Set number of propagating modes

         IF    ( NPROF .EQ. 1      ) THEN   ! Range-independent case
               CALL EVAL( PHIS( 1, J ), PHIR, NRD, 
     &            R, NR, RR, CK, M, OPT, P )
         ELSE
            IF ( OPT(2:2) .EQ. 'C' ) THEN   ! Range-dependent case
C              --- Coupled mode theory
               CALL EVALCM( RPROF, NPROF, PHIS( 1, J ), PHIR, RD, NRD, 
     &            R, NR, CK, M, OPT, P )
            ELSE
C              --- Adiabatic mode theory
               CALL EVALAD( RPROF, NPROF, PHIS( 1, J ), PHIR, RD, NRD, 
     &            R, NR, M, OPT, P )
            ENDIF
         ENDIF

         CALL WRTFLD( P, NRD, NR, IREC )

      END DO   ! next source depth

      STOP
      END
