C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  Routines and Functions for W_rho multiplications
C  subroutine xWRHO(b, x, Wrho, Wrho_t, IQrho, JQrho, NR, NT)
c  real function xWRHOy(x, Wrho, Wrho_t, y, IQrho, JQrho, NR, NT) 
c  real function xWRHOx(x, Wrho, Wrho_t, IQrho, JQrho, NR, NT) 
c  real function d1WDELTAd2(d1 , Wdelta, d2, Nobs) 
c  real function dWDELTAd(d , Wdelta, Nobs) 
c
c  drl 18 aug 99
c      adding temporal smoothing (weight w2) to xWRHO, xWRHOy, and xWRHOx.
C      this requires a two-matrix representation of the covariance matrix: 
C           Wrho     -- the spatial smoothing as before, 
C                                      incorportating (w0,w1) (unchanged)
C           Wrho_t -- the temporal smoothing matrix.  This is a 
C                     node-indexed diagonal matrix which multiplies 
C                     the second-temporal-difference matrix.  
C                     It incorporates w2.
C
C  drl and cgh 9 July 1999
C      changes to xWRHO, xWRHOy, xWRHOx to implement trapezoidal
C      rule at the end points of the time integral. 
C
C drl and cgh 2 Dec '98 
C
C The xWRHO... functions and subroutine are not true library 
C functions.  They contain the statement 'include 'CASCO.DIM''
C Which is used to allocate space for scratch arrays and 
C dimension the 2d arrays. 
C  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      subroutine xWRHO(b, x, Wrho, Wrho_t, IQrho, JQrho, NR, NT) 
C
C  drl 18 aug 99: added Wrho_t, the time-smoothing feature (embeds w2)
C
C  drl and cgh 9 July 1999
C changes to implement trapezoidal rule at the end points of the time integral. 
C
C 2 Dec 98, drl
C For Casco computation of gradient of omega
C
C Computes b=x*Wrho
C
C b(i,j) and x(i,j) are dimensioned the same: 
C   i is spatial index; i=1,nr
C   j is time index; j=1,nt
C Wrho:
C   square (nr, nr), symmetric positive definite weight matrix for rho
C   independent of time
C   stored in standard Sprspak form, with pointers IQrho, JQrho
C 
C For any fixed time j:
C   b(i,j)=x(i,j)*Wrho*scale  
C Do this for all values of j=1,nt
C
      include 'CASCO.DIM'
C Arrays passed:
C
      Real Wrho(NWrhoDIM), x(nrdim,ntbcdim), b(nrdim, ntbcdim)
      Real Wrho_t(nrdim)
      Integer IQrho(nrdim), JQrho(nrdim)
      Integer nr, nt
C
C Local arrays:
C
      real xj(nrdim), bj(nrdim)
      integer i, j 
      real    scale 
C
      do j=1,nt
	scale = 1.0 
	if((j.eq.1).or.(j.eq.nt)) scale = 0.5 
      	do i=1,nr
      	   xj(i)=x(i,j)
      	enddo
      	call SPRSpreMLT(Wrho,IQrho,JQrho,xj,bj,NR)
      	do i=1,nr 
      	   b(i,j) = scale*bj(i)
      	enddo
      enddo
C
C  18 aug 99, drl: the temporal smoothing part.
C
      do i=1,nr
         j=1
         b(i,j)=b(i,j)+Wrho_t(i)*( -x(i,j+1)+x(i,j)         )
         do j=2,nt-1
      b(i,j)=b(i,j)+Wrho_t(i)*( -x(i,j+1)+2*x(i,j)-x(i,j-1) )
         enddo
         j=nt
         b(i,j)=b(i,j)+Wrho_t(i)*(          x(i,j)-x(i,j-1) )
      enddo
C
      Return
      End
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      real function xWRHOy(x, Wrho, Wrho_t, y, IQrho, JQrho, NR, NT) 
C
C  drl 18 aug 99: added Wrho_t, the time-smoothing feature (embeds w2)
C
C  drl and cgh 9 July 1999
C changes to implement trapezoidal rule at the end points of the time integral. 
C
C 2 Dec 98, drl
C For Casco computation of quadratic forms using Wrho
C Rho contributions to Theta use this
C
C Computes z=x*Wrho*y
C
C x(i,j) and y(i,j) are dimensioned the same: 
C   i is spatial index; i=1,nr
C   j is time index; j=1,nt
C Wrho:
C   square (nr, nr), symmetric positive definite weight matrix for rho
C   independent of time
C   stored in standard Sprspak form, with pointers IQrho, JQrho
C 
C For any fixed time j:
C   zj=x(i,j)*Wrho*y(i,j) 
C Sum the zj over all values of j=1,nt
C
      include 'CASCO.DIM'
C Arrays passed:
C
      Real Wrho(NWrhoDIM), x(nrdim,ntbcdim), y(nrdim, ntbcdim)
      Real Wrho_t(nrdim)
      Integer IQrho(nrdim), JQrho(nrdim)
      Integer nr, nt
C
C Local arrays:
C
      Real xj(NRdim), bj(NRdim)
      real z
      integer i , j
      real    scale 
C
C summing over all time j:
C   z = sum_j {x(i,j)*Wrho *y(i,j)}
C
      z=0.0
      do j=1,NT
	 scale = 1.0 
	 if((j.eq.1).or.(j.eq.nt)) scale = 0.5 
         do i=1,nr
         xj(i)=x(i,j)
         enddo
         call SPRSpreMLT(Wrho,IQrho,JQrho,xj,bj,NR)
         do i=1,nr 
         z = z + scale*bj(i)*y(i,j)
         enddo
      enddo

      xWRHOy = z
C
C  18 aug 99, drl: the temporal smoothing part.
C
      z=0.0
      do i=1,nr
         j=1
         z=z+Wrho_t(i)*y(i,j)*( -x(i,j+1)+x(i,j)            )
         do j=2,nt-1
            z=z+Wrho_t(i)*y(i,j)*( -x(i,j+1)+2*x(i,j)-x(i,j-1) )
         enddo
         j=nt
         z=z+Wrho_t(i)*y(i,j)*(                x(i,j)-x(i,j-1) )
      enddo
C
      xWRHOy=xWRHOy+z
C
      Return
      End
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      real function xWRHOx(x, Wrho, Wrho_t, IQrho, JQrho, NR, NT) 
C
C  drl 18 aug 99: added Wrho_t, the time-smoothing feature (embeds w2)
C
C  drl and cgh 9 July 1999
C changes to implement trapezoidal rule at the end points of the time integral. 
C
C 2 Dec 98, drl
C For Casco computation of quadratic forms using Wrho
C Rho contributions to Omega and to Theta use this
C
C Computes z=x*Wrho*x
C
C x(i,j): 
C   i is spatial index; i=1,nr
C   j is time index; j=1,nt
C Wrho:
C   square (nr, nr), symmetric positive definite weight matrix for rho
C   independent of time
C   stored in standard Sprspak form, with pointers IQrho, JQrho
C 
C For any fixed time j:
C   zj=x(i,j)*Wrho*x(i,j) 
C Sum the zj over all values of j=1,nt
C
      include 'CASCO.DIM'
C Arrays passed:
C
      Real Wrho(NWrhoDIM), x(nrdim,ntbcdim)
      Real Wrho_t(nrdim)
      Integer IQrho(nrdim), JQrho(nrdim)
      Integer nr, nt
C
C Local arrays:
C
      Real xj(NRdim), bj(NRdim)
      real z
      integer i , j
      real    scale 
C
C summing over all time j:
C   z = sum_j {x(i,j)*Wrho*x(i,j)}
C
      z=0.0
      do j=1,NT
	 scale = 1.0 
	 if((j.eq.1).or.(j.eq.nt)) scale = 0.5 
         do i=1,nr
         xj(i)=x(i,j)
         enddo
         call SPRSpreMLT(Wrho,IQrho,JQrho,xj,bj,NR)
         do i=1,nr 
         z = z + scale*bj(i)*x(i,j)
         enddo
      enddo
      xWRHOx = z
C
C  18 aug 99, drl: the temporal smoothing part.
C
      z=0.0
      do i=1,nr
         j=1
         z=z+Wrho_t(i)*x(i,j)*( -x(i,j+1)+x(i,j)               )
         do j=2,nt-1
            z=z+Wrho_t(i)*x(i,j)*( -x(i,j+1)+2*x(i,j)-x(i,j-1) )
         enddo
         j=nt
         z=z+Wrho_t(i)*x(i,j)*(                x(i,j)-x(i,j-1) )
      enddo
C
      xWRHOx=xWRHOx+z
C

      Return
      End
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  Routines for W_Delta multiplications
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      real function d1WDELTAd2(d1 , Wdelta, d2, Nobs) 
C
C 2 Dec 98, drl
C For Casco computation of quadratic forms using Wdelta
C Computes z=d1*Wdelta*d2 
C
C Arrays passed
      Real d1(*), d2(*), Wdelta(*)
      Integer Nobs  
C
C local
      real z 
      integer i
C  
      z=0.0
      do i=1,nobs
         z=z+d1(i)*d2(i)*wdelta(i)
      enddo
      d1WDELTAd2 = z
C
      Return
      End
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      real function dWDELTAd(d , Wdelta, Nobs) 
C
C 2 Dec 98, drl
C For Casco computation of quadratic forms using Wdelta
C Computes z=d*Wdelta*d
C
C Arrays passed
        Real d(nobs),Wdelta(nobs)
        Integer Nobs  

C local 
        real z
	integer i 
C  
        z=0.0
        do i=1,nobs
           z=z+d(i)*d(i)*wdelta(i)
        enddo
        dWDELTAd = z
C
      Return
      End
     
