C
C***********************************************************************
C*****************************************************************drlcgh
C
      Subroutine BuildWrho3(WrhoSprs,Wrho_t,IQrho,JQrho,w0,w1,w2,
     &                           nr,LRVS,x,y,nt,dt)
C
C  drl 18 Aug 99 -- from BuildWrho2
C        add computation of the diagonal matrix Wrho_t(i), where i is boundary node number.
C        this is needed to achieve temporal smoothing.
C        added three items to the argument list: 
C              Wrho_t -- computed herein.  One numver for each boundary node
C              w2       -- passed with value gotten from input stream
C              dt        -- the time step relevant to the inferred boundary conditions Rho 
C                              (not the time step used in the forward and reverse simulations).
C
C  drlcgh 9 July 99 -- from BuildWrho2 
C        use nt-1 instead of nt in denominator
C        related to changes in xWrho, xWrhox and xWrhoy
C
C  drl 6 July 99 -- from BuildWrho
C    purpose: 
c       a) to make the user interface to Wrho simpler, specifically the
C          term W1; and 
c       b) to account for variable mesh spacing in the diagonal term
C
c    a) The Laplacian (slope**2) part:  the average length of 
c         a boundary element is moved into the denominator of Wrho.  
C         This means that the new weight W1 for this matrix should be 
C         equal to the expected value of (boundary slope)**2
C   
C    b) The diagonal (size**2) part: this is fixed so that squared 
C         nodal values are weighted by their relative element sizes.  
C         The diagonal matrix is normalized by the average
C         element length.  So the weight for this matrix W1 retains its
C         original meaning, i.e. the expected size of (boundary
C         elevation)**2
C    In effect, Wrho now represents
C         { W0*Int[(z**2)ds] + W1*Int[(dz/ds)**2)ds] }/NT/Int[ds]
c    Integration of z**2 is performed by integral lumping.
C
C  drl 24 Nov 98 -- for casco BC weighting matrix W_rho
C  from BuildW in truxton
C  converted to sparse storage WrhoSprs (*) 24 nov 98
C  square, full matrix Wrho(nr, nr) commented out
C
C  Builds the Casco W_rho matrix WrhoSprs, IQrho, JQrho 
C    using standard sprspak storage. 
C  W0, W1 are the weights to be applied to the size and slope measures.  
C  W0 should be of order 1/(variance of rho)
C  W1 should be ???  
C  NR is the # of dirichlet BC's
C  LRVS(i) is the global node number of the ith Dirichlet node.
C  It is assumed that LRVS presents nodes in sequential, CCW order 
C    around a single simply-connected, closed contour. 
C  NT is the number of DoF in time belonging to Rho 
C    (# of 'j' steps for the BC's)
C  NR*NT is the total number of DoF in space-time for Rho
C  Wrho is divided by NR*NT to normalize w.r.t. the number of BC's
C
      INCLUDE 'CASCO.DIM'
c     INCLUDE 'Wrho.dim'
C       (need nrdim, nndim)
c     Parameter(NWdim=nrdim*3)
C
C Input:
C 
      Real X(nndim), Y(nndim)
      Integer LRVS(nndim)
      real   w0, w1, w2, dt
      integer nr, nt 

C Output:
C
C      Real Wrho(nrdim, nrdim)
      Real WrhoSprs(NWrhodim), Wrho_t(nrdim)
      Integer IQrho(nrdim), JQrho(NWrhodim)      
C
C
C Local Arrays:
C   INBE is the Bndry Element Incidence in global node numbers
C   InbeBN is the same but in boundary node numbers as ordered in LRVS      
C
      Integer INBE(NRdim,2), IBC(nndim)
      integer InbeBN(NRdim,2)
C
C local variables
      real factor
      integer i, i1,i2, i11, i22, i1test, i2test
      integer L
      integer nbe 
      real    ds, term, term2, term0
      real blength
      integer k11, k12, k21, k22
C 
C functions
      integer kay 
C
C Build the .bel file from the LRVS list 
C ASSUMING:
c nodes are ordered sequentially, counterclockwise.  Truxton subroutine
c BNodeOrder is intended to guarantee this; as of 11-24-98 that subr. 
c is not invoked; relying on default ordering from input.
c
      write(*,*)'Subr. BuildWrho: making boundary elements from the BC',
     &' node list LRVS.'
      write(*,*)'Sequential, CCW order is assumed.'
      write(*,*)'nr,nt ', nr,nt
      write(*,*)'w0,w1,w2,dt', w0,w1,w2,dt
      nbe=nr
      do L=1, nr-1
         inbe(L,1)=lrvs(L)
         inbe(L,2)=lrvs(L+1)
C
         inbeBN(L,1)=L
         inbeBN(L,2)=L+1
      enddo
C
      inbe(nr,1)=lrvs(nr)
      inbe(nr,2)=lrvs(1)
C
      inbeBN(nr,1)=nr
      inbeBN(nr,2)=1
C
C Compute IBC(I), i.e. the inverse of LRVS
C IBC(i) is the boundary node number of global node number i. 
C IBC(i) = 0 if global node i is an interior node
C IBC(i) is between 1 and NR if global node i is on the boundary
C Note 11-24-98: New construction inbeBN above (added to get to sparse
C    storage mode) should make ibc unneeded.  
C    I am leaving it in with a consistency check below,
C    in the assembly loop.  If there is a need to clean house, 
C    ibc could be tossed.
C
      do i=1,nndim
         IBC(i)=0
      enddo
      do i=1,nr
         IBC(LRVS(i))=i
      enddo
C
C Check for flaws in IBC, INBE
C
c      write(2,*)'Subr. BuildW: Opening BuildW.out for .bel diagnosis'
c      open(unit=19,file='BuildW.out',status='unknown')
c      write(19,*) ' INBE:'
c      Do L=1,NBE
c      I1=INBE(L,1)
c      I2=INBE(L,2)
c      write(19,*) I1, I2
c      enddo
c      write(19,*) ' IBC:'
c      do I=1,nn
c      if(IBC(i).ne.0) then
c      write(19,*) IBC(i)
c      endif
c      enddo
c      write(19,*) ' LRVS:'
c      do I=1,nr
c      write(19,*) LRVS(i)
c      enddo   
C      
C Initialize Wrho, WrhoSprs, Wrho_t
C
C      do i=1, nr
C      do j=1, nr
C      Wrho(i,j)=0.0      
C      enddo
C      enddo
      do i=1,NWrhodim
         WrhoSprs(i)=0.
      enddo
      do i=1, nr
      Wrho_t(i)=0.
      enddo
C
C Build the sparse pointers for WrhoSprs
C
      call sprsbldw1(inbebn, nrdim, 2, iqrho, jqrho, nr, nr, 3, 1)
C
C Make the Laplacian portion, weight w1
C   and the diagonal portion, weight w0
C
      blength=0.0
      do L=1, nbe
C   Boundary Node Numbers (location in BC list LRVS)
        I1=IBC( INBE(L,1) )
        I2=IBC( INBE(L,2) )
C (temporary test)      
        i1test=inbebn(l,1)
        i2test=inbebn(l,2)
        if (i1test.ne.i1)write(*,*) 'problem in bel area'
        if (i2test.ne.i2)write(*,*) 'problem in bel area'
C (end test)      
C   Global Node Numbers:      
        I11= INBE(L,1) 
        I22= INBE(L,2) 
C
        if(i1.le.0.or.i2.le.0) then
        write(*,*) 'error: element #', L
        stop
      endif
      if(i1.gt.nr.or.i2.gt.nr) then  
      write(*,*) 'error: element #', L
      stop
      endif
C
      ds= sqrt( (x(I22)-x(I11))**2 + (y(I22)-y(I11))**2 )
      term = w1/ds
      term0= w0*ds/2.
      term2= w2*ds/2.
      blength=blength+ds
c      Wrho(i1, i1) = Wrho(i1, i1) + term
c      Wrho(i1, i2) = Wrho(i1, i2) - term
c      Wrho(i2, i1) = Wrho(i2, i1) - term
c      Wrho(i2, i2) = Wrho(i2, i2) + term
      k11=kay(i1,i1,iqrho, jqrho)
      k12=kay(i1,i2,iqrho, jqrho)
      k21=kay(i2,i1,iqrho, jqrho)
      k22=kay(i2,i2,iqrho, jqrho)
      WrhoSprs(k11)= WrhoSprs(k11) + term+term0
      WrhoSprs(k12)= WrhoSprs(k12) - term
      WrhoSprs(k21)= WrhoSprs(k21) - term
      WrhoSprs(k22)= WrhoSprs(k22) + term+term0
c
      Wrho_t(i1)=Wrho_t(i1)+term2
      Wrho_t(i2)=Wrho_t(i2)+term2
C
      enddo
C
C Add the diagonal term w0
C 
c      do i=1,nr
C      Wrho(i,i) = Wrho(i,i) + w0
c      k=kay(i,i,iqrho, jqrho)
c      WrhoSprs(k)= WrhoSprs(k) + w0
c      enddo
C
C Finally, Normalize Wrho by Blength*(NT-1)
C 
      factor=1.0/(1.0*(nt-1)*blength)
      do i=1,nr*3
      WrhoSprs(i)= WrhoSprs(i)*factor
      enddo
C
C and Normalize Wrho_t by Blength*(NT-1)*dt*dt
C
      factor=1.0/( 1.0*(nt-1)*blength*dt*dt )
      do i=1, nr
      Wrho_t(i)=Wrho_t(i)*factor
      enddo
C
C Done
C
      return
      end             
      
C
C***********************************************************************
C*****************************************************************drlcgh
C
