      SUBROUTINE TRIANINT(NMAP,M,ISWT,ITIME)

      USE BLKMAP
      USE BLK1MOD
      SAVE
!      INCLUDE 'BLK1.COM' 

	  DIMENSION WGT(8)
      REAL*8 XMINL,YMINL,XMAXL,YMAXL
!	  data itime/0/

!     LOOK FOR MATCHING POINTS
      ITIMESKP=0
      DO K=1,MAXPTS
        DISQ=(XUSR(M)-XMAP(K))**2+(YUSR(M)-YMAP(K))**2
        IF(DISQ .LT. 1.) THEN
          WD(M)=VAL(K)
          FPN = WD(M)*10. 
          X = CORD(M,1) 
          Y = CORD(M,2) - .11 
          IF(X .GT. 0.  .AND.  X .LT. HSIZE  .AND.                         &
     &     Y .GT. 0.  .AND.  Y .LT. 7.5) THEN                           
            CALL RRED 
            CALL NUMBR(X,Y,0.1,FPN,0.0,-1) 
          endif
          ITIMESKP=1
		  GO TO 300
        ENDIF
      ENDDO

!     Search for element that has circumcircle around the node

      IF(ISWT .NE. 0) THEN
	    IF(ITIME .EQ. 0) NSTART=1
	  ELSE
	    NSTART=1
	  ENDIF
      DO N=NSTART,NELTS
	    IF(NOPEL(N,1) .EQ. 0) GO TO 200
	    if(RADS(N) .eq. 0.) then
          CALL CCENTRE(XMAP(NOPEL(N,1)),XMAP(NOPEL(N,2)),XMAP(NOPEL(N,3))      &
          &,YMAP(NOPEL(N,1)),YMAP(NOPEL(N,2)),YMAP(NOPEL(N,3))                  &
          &,XCEN(N),YCEN(N),RADS(N)) 
        endif		  

        IF(RADS(N)+XCEN(N) .GE. XUSR(M)) THEN
		  NSTART=N
		  GO TO 210
		ENDIF
  200   CONTINUE
      ENDDO
  210 CONTINUE
  220 continue        
!      WRITE(155,*) M,NSTART
	  DO N=NSTART,NELTS
	    IF(NOPEL(N,1) .EQ. 0) GO TO 250
	    if(RADS(N) .eq. 0.) then
          CALL CCENTRE(XMAP(NOPEL(N,1)),XMAP(NOPEL(N,2)),XMAP(NOPEL(N,3))      &
          &,YMAP(NOPEL(N,1)),YMAP(NOPEL(N,2)),YMAP(NOPEL(N,3))                  &
          &,XCEN(N),YCEN(N),RADS(N)) 
        endif		  
        xminl=min(XMAP(NOPEL(N,1)),XMAP(NOPEL(N,2)),XMAP(NOPEL(N,3)))
        xmaxl=max(XMAP(NOPEL(N,1)),XMAP(NOPEL(N,2)),XMAP(NOPEL(N,3)))
        yminl=min(YMAP(NOPEL(N,1)),YMAP(NOPEL(N,2)),YMAP(NOPEL(N,3)))
        ymaxl=max(YMAP(NOPEL(N,1)),YMAP(NOPEL(N,2)),YMAP(NOPEL(N,3)))
!        IF(M .EQ. 6316) THEN
!          WRITE(156,'(2I6,6F15.2)') M,N,XUSR(M),XMINL,XMAXL,YUSR(M),YMINL,YMAXL
!        ENDIF
        if(xusr(m) .lt. xminl-0.01  .or.  xusr(m)  .gt. xmaxl+0.01) then
          go to 250
        elseif(yusr(m) .lt. yminl-0.01  .or.  yusr(m)  .gt. ymaxl+0.01) then
		  go to 250
	    endif
!        IF(M .EQ. 6316) WRITE(156,*) 'PASSED X AND Y TEST',N

        DISQ=(XUSR(M)-XCEN(N))**2+(YUSR(M)-YCEN(N))**2

		IF(DISQ .LE. RADS(N)**2*1.0001) THEN

!            IF(M .EQ. 6316) write(156,*) m,n,disq,rads(n)**2,xusr(m),xcen(n)

!     We have a candidate

          CALL GETWT(N,XUSR(M),YUSR(M),WGT,1)
          DO K=1,3
		    IF(WGT(K) .LT. -1E-4  .OR. WGT(K) .GT. 1.0001) THEN
			  WRITE(142,*) 'REJECT',m,n,disq,rads(n)**2,wgt(1),wgt(2),wgt(3)
              if(nstart .gt. 1) then
                nstart=1
                go to 220
              endif
			  GO TO 250
			ENDIF
		  ENDDO
          WD(M)=WGT(1)*VAL(NOPEL(N,1))+WGT(2)*VAL(NOPEL(N,2))+WGT(3)*VAL(NOPEL(N,3))
          FPN = WD(M)*10. 
          X = CORD(M,1) 
          Y = CORD(M,2) - .11 
          IF(X .GT. 0.  .AND.  X .LT. HSIZE  .AND.                         &
     &     Y .GT. 0.  .AND.  Y .LT. 7.5) THEN                           
            CALL RRED 
            CALL NUMBR(X,Y,0.1,FPN,0.0,-1) 
		  endif
		  GO TO 300
		ENDIF
  250   CONTINUE
	  ENDDO
  300 CONTINUE
	  IF(ITIMESKP .EQ. 0) ITIME=1

      RETURN
	  END


      SUBROUTINE GETWT(N,XSW,YSW,WGT,ISWT)

!-                                                                      
!......SUBROUTINE TO EVALUATE FUNCTION AT GRID POINTS                   
!-                                                                      
!-    N = ELEMENT NUMBER        
!_    XSW = X COORDINATE OF DESIRED POINT
!_    YSW = Y COORDINATE OF DESIRED POINT
!     WGT(8) = ARRAY OF WEIGHTING FUNCTIONS
!     ISWT   = SWITCH FOR CHOICE BETWEEN LINEAR AND QUADRATIC WEIGHTING
!            = 1 FOR LINEAR
!            = 2 FOR QUADRATIC
!     FROM COMMON
!     NOP = LIST OF NODAL CONNECTIONS AROUND AN ELEMET
!     CORD = REAL*8 ARRAY OF NODAL COORDINATES
!                                                                       
      USE BLKMAP
      USE BLK1MOD
      REAL*8 XN,DNX,DNY,XSW,YSW 
      DOUBLE PRECISION XG,YG,XK,YK,XP,YP
!      INCLUDE 'BLK1.COM'
!-                                                                      
      DIMENSION X(9),Y(9),WGT(8)
!-                                                                      
      DATA TOL/0.01/ 
!-                                                                      

!-                                                                      
!......DETERMINE ELEMENT TYPE                                           
!-                                                                      
!IPKOCT93 ADD                                                           
      if(n .eq. 1910) then
         aaa=0
      endif
      NCN=6 
      IT=2 
!-                                                                      
!......ESTABLISH LOCAL COORDINATES FOR EACH NODE POINT OF ELEMENT                 
!-                                                                      
      K1=NOPEL(N,1) 
      X(1)=0. 
      Y(1)=0. 
      DO 300 K=3,NCN,2
        K2=NOPEL(N,K/2+1) 
        X(K)=XMAP(K2)-XMAP(K1) 
        Y(K)=YMAP(K2)-YMAP(K1) 
  300 END DO 
      X(2)=X(3)/2.
      Y(2)=Y(3)/2.
      X(4)=(X(3)+X(5))/2.
      Y(4)=(Y(3)+Y(5))/2.
      X(6)=X(5)/2.
      Y(6)=Y(5)/2.
	  xminl=min(x(1),x(3),x(5))
	  yminl=min(y(1),y(3),y(5))
	  xmaxl=max(x(1),x(3),x(5))
	  ymaxl=max(y(1),y(3),y(5))


!-                                                                      
!......ESTABLISH LOCAL COORDINATES OF DESIRED POINT                     
!-                                                                      
      XP=XSW-XMAP(K1) 
      YP=YSW-YMAP(K1) 

      if(xp .lt. xminl  .or.  xp  .gt. xmaxl) then
	    wgt(1)=2.0
        return
      elseif(yp .lt. yminl  .or.  yp  .gt. ymaxl) then
        wgt(1)=2.0
		return
	  endif
      XG=0. 
      YG=0. 
!-                                                                      
!......ITERATE TO FIND LOCAL COORDINATE                                 
!-                                                                      
      DO 400 ITER=1,10 
        DXKDX=0. 
        DXKDY=0. 
        DYKDX=0. 
        DYKDY=0. 
        XK=-XP 
        YK=-YP 
        DO 350 K=2,NCN 
          XK=XK+XN(IT,K,XG,YG)*X(K) 
          YK=YK+XN(IT,K,XG,YG)*Y(K) 
          DXKDX=DXKDX+DNX(IT,K,XG,YG)*X(K) 
          DYKDX=DYKDX+DNX(IT,K,XG,YG)*Y(K) 
          DXKDY=DXKDY+DNY(IT,K,XG,YG)*X(K) 
          DYKDY=DYKDY+DNY(IT,K,XG,YG)*Y(K) 
  350   END DO 
        DET=DXKDX*DYKDY-DXKDY*DYKDX 
        DX=(-DYKDY*XK+DXKDY*YK)/DET 
        DY=( DYKDX*XK-DXKDX*YK)/DET 
        XG=XG+DX 
        YG=YG+DY 
        IF(ABS(DX).LT.TOL .AND. ABS(DY).LT.TOL) GO TO 420 
  400 END DO 
!-                                                                      
!......NOW GET WEIGHTING FUNCTIONS FOR QUAD FUNCTION
!-                                                                      
  420 CONTINUE 
      DO  K=1,NCN 
        WGT(K)=XN(IT,K,XG,YG)
      END DO 

      IF(ISWT .EQ. 1) THEN
!-
!-     REDUCE TO LINEAR FUNCTION BY ADDING TERMS
!-
        DO K=2,NCN,2
          WGT(K-1)=WGT(K-1)+WGT(K)/2.
          IF(K .LT. NCN) THEN
            WGT(K+1)=WGT(K+1)+WGT(K)/2.
          ELSE
            WGT(1)=WGT(1)+WGT(K)/2.
          ENDIF
        ENDDO
!-    
!-     THEN COMPACT ARRAY
!-
        DO K=1,NCN/2
          WGT(K)=WGT(2*K-1)
        ENDDO

	  ENDIF

      RETURN 
      END                                           
