c
c       Modelling of multiples frequency slice by frequency slice  
c      
c       inputs: ifrmax = number of frequency samples. 
c               [If n is the number of samples in time used in fft,
c                frmax is the maximal frequency, and dt is the sample interval 
c                in time of the data, then ifrmax = 1+int(framx*n*dt)]
c               nrecv = number of receivers per shot 
c                    (we assume in this subroutine that nrecv is equal to 
c                     the number of shots)
c               cpres = array containing pressure data. These data are 
c                       assumed to be in split spread. cpres is physically
c                       stored in nfreqxmxm array.
c               cvz   = array containing the vertical component of the 
c                       particle velocity. These data are assumed to be 
c                       in split spread. For the cases in which the 
c                       vertical component of the particle velocity is not 
c                       available used the formula in Box 10-5 to compute 
c                       this component from pressure data.
c                       cvz is physically stored in nfreqxmxm array.
c
c       returns: cm1: second term of the Kirchhoff series
c                cm2: third term of the Kirchhoff series
c
        subroutine predmult(ifrmax,nrecv,nfreq,m,cvz,cpres,cm1,cm2)
c
        complex cvz(nfreq,m,m), cpres(nfreq,m,m)
        complex cm1(nfreq,m,m), cm2(nfreq,m,m)
c
        do iom=1,ifrmax
           print *,' frequency= ',iom
c 
c       *********************************************************
c       * 1. compute of the second term of the Kirchhoff series *
c       *********************************************************
c
           do ixs=1,nrecv
              do ixr=1,nrecv
                 cm1(iom,ixr,ixs) = cmplx(0.,0.)
                 do ix=1,nrecv
                    cm1(iom,ixr,ixs) = cm1(iom,ixr,ixs) +
     +                      (cvz(iom,ix,ixs)*cpres(iom,ixr,ix))
                 enddo
              enddo
           enddo
c
c       *********************************************************
c       * 2. compute of the second term of the Kirchhoff series *
c       *********************************************************
c
           do ixs=1,nrecv
              do ixr=1,nrecv
                 cm2(iom,ixr,ixs) = cmplx(0.,0.)
                 do ix=1,nrecv
                    cm2(iom,ixr,ixs) = cm2(iom,ixr,ixs) +
     +                      (cvz(iom,ix,ixs)*cm1(iom,ixr,ix))
                 enddo
              enddo
           enddo
c
        enddo
c
        return 
        end