      PROGRAM FIELDMO

C     Generate file of replica vectors

      INTEGER PRTFIL, FLPFIL
      PARAMETER ( FLPFIL = 5, PRTFIL = 6, 
     &   MXM = 1500, MXNW = 4094, MXNWZ = 210000, 
     &   MXNSD = 51, MXNRD = 201, MXNOC = 9 )

      INTEGER CTR( MXNOC )
      REAL      SD( MXNSD ), RD( MXNRD ),
     &          W( MXNW ), ALPHA( MXNOC )
      COMPLEX   CK( MXNOC, MXM ), CKTEMP( MXM), PHITEMP( MXM, MXNRD ),
     &          CKA( MXM ), CKB( MXM ),
     &          PHIS( MXNOC, MXM, MXNSD ),
     &          PHIR( MXNOC, MXM, MXNRD ),
     &          PHISA( MXM, MXNRD ), PHISB( MXM, MXNRD ),
     &          PHIRA( MXM, MXNRD ), PHIRB( MXM, MXNRD ),
     &          P( MXNWZ )
      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
C      CALL STR$UPCASE( OPT, OPT )
      COMP = OPT( 3:3 )

      READ( FLPFIL, * ) MLIMIT
      MLIMIT = MIN( MXM, MLIMIT )

      READ( FLPFIL, * ) NOC

C     *** Read ocean weights ***

      NW = 0
      W( 3 ) = -999.9
      READ( FLPFIL, * ) NW, ( W( IA ), IA = 1, NW )

      CALL SUBTAB( W, NW )
      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Number of ocean weights = ', NW
      IF ( NW .GE. 1 )
     &   WRITE( PRTFIL, * ) ( W( IA ), IA = 1, MIN( NW, 51 ) )

      WRITE( PRTFIL, * )

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

      ZMIN = -1.0E10
      ZMAX = +1.0E10

      CALL SDRD( FLPFIL, PRTFIL, ZMIN, ZMAX,
     &   MXNSD, NSD, SD,
     &   MXNRD, NRD, RD )

      WRITE( PRTFIL, * )

      IF ( NW * NRD .GT. MXNWZ )
     &   CALL ERROUT( PRTFIL, 'F', 'FIELD',
     &                'Number of oceans x Depths too large' )

      READ( FLPFIL, * ) R
      R = 1000.0 * R

C     *** Read in modes ***

      M = MLIMIT

      DO 2000 IOC = 1, NOC

         CALL GETMOD( IOC, ' ', MXM, SD, NSD, 'N',
     &      CKTEMP, PHITEMP, MNEW, FREQ, TITLE )

C        --- copy into the big matrix
         DO 1000 MODE = 1, MNEW
            DO 500 ISD = 1, NSD
               PHIS( IOC, MODE, ISD ) = PHITEMP( MODE, ISD )
  500       CONTINUE
 1000    CONTINUE

         CALL GETMOD( IOC, ' ', MXM, RD, NRD, COMP,
     &      CKTEMP, PHITEMP, MNEW, FREQ, TITLE )

C        --- copy into the big matrix
         DO 1800 MODE = 1, MNEW
            DO 1500 IRD = 1, NRD
               PHIR( IOC, MODE, IRD ) = PHITEMP( MODE, IRD )
 1500       CONTINUE
            CK( IOC, MODE ) = CKTEMP( MODE )
 1800    CONTINUE

         M = MIN( M, MNEW )
 2000 CONTINUE   ! Next ocean

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

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

      DO 8000 J = 1, NSD

C        --- Initialize counter
         DO 3700 IOC = 1, NOC
            CTR( IOC ) = 0
 3700    CONTINUE

         DO 7000 I = 1, NW ** ( NOC - 2 )

C           --- increment the counter and assign the weight
            IF ( NOC .GE. 3 ) THEN
            DO 6500 IOC = 3, NOC
               IF ( CTR( IOC ) .GT. NW ) THEN
                  CTR( IOC ) = 1
               ELSE
                  CTR( IOC ) = CTR( IOC ) + 1
               ENDIF
               ALPHA( IOC ) = W( CTR( IOC ) )
 6500       CONTINUE
            ENDIF

C           ****** Interpolate the modes ******

            DO 6000 MODE = 1, M

C              --- initialize
               CKA( MODE ) = CK( 1, MODE )
               CKB( MODE ) = CK( 2, MODE )

               DO 2600 ISD = 1, NSD
                  PHISA( MODE, ISD ) = PHIS( 1, MODE, ISD )
                  PHISB( MODE, ISD ) = PHIS( 2, MODE, ISD )
 2600          CONTINUE

               DO 2700 IRD = 1, NRD
                  PHIRA( MODE, IRD ) = PHIR( 1, MODE, IRD )
                  PHIRB( MODE, IRD ) = PHIR( 2, MODE, IRD )

 2700          CONTINUE

C              --- add in weighted data from other oceans

               IF ( NOC .GE. 3 ) THEN
               DO 4000 IOC = 3, NOC
C                 --- following can be done more efficiently but ...
                  CKA( MODE ) = CKA( MODE )
     &                  + ALPHA( IOC ) * CK( IOC, MODE )
                  CKB( MODE ) = CKB( MODE )
     &                  + ALPHA( IOC ) * CK( IOC, MODE )

                  DO 3600 ISD = 1, NSD
                     PHISA( MODE, ISD ) = PHISA( MODE, ISD )
     &                  + ALPHA( IOC ) * PHIS( IOC, MODE, ISD )
                     PHISB( MODE, ISD ) = PHISB( MODE, ISD )
     &                  + ALPHA( IOC ) * PHIS( IOC, MODE, ISD )
 3600             CONTINUE

                  DO 3800 IRD = 1, NRD
                     PHIRA( MODE, IRD ) = PHIRA( MODE, IRD )
     &                  + ALPHA( IOC ) * PHIR( IOC, MODE, IRD )
                     PHIRB( MODE, IRD ) = PHIRB( MODE, IRD )
     &                  + ALPHA( IOC ) * PHIR( IOC, MODE, IRD )
 3800             CONTINUE

 4000          CONTINUE  ! Next ocean
               ENDIF

 6000       CONTINUE  ! Next mode

            CALL EVALMO( PHISA( 1, J ), PHISB( 1, J ),
     &         CKA, CKB, PHIRA, PHIRB, M,
     &         RD, NRD, W, NW,
     &         R, OPT, P )

            CALL WRTFLD( P, NRD, NW, IREC )
 7000    CONTINUE   ! Next weight

 8000 CONTINUE  ! Next source depth

      STOP
      END
