!IPK  LAST UPDATE FEB 11 2002 ADD LOCK AS VARIABLE
!ipk  last update Feb 10 1997                                           
      SUBROUTINE GRIDSB(ISWTIN)
!                                                                       
!  Routines to control interpolation of nodal elevations                
!                                                                       
      USE WINTERACTER
      USE BLKMAP
      USE BLK1MOD

      include 'd.inc'

!      INCLUDE 'BLK1.COM' 

      INCLUDE 'TXFRM.COM'

!iPK APR94
      COMMON /RECOD/ IRECD,TSPC


!IPK MAY02       COMMON  /TXFRM/  XS, YS, TXSCAL 
!                                                                       
!IPKJAN94      INTEGER*2 LISTM                                          
!            INTEGER LISTM 
!ipk feb94 add ARF then remove may97                                    

!      DIMENSION LISTM(1000),listt(1600,4),nlf(4),icomp(4),xnear(4) 
!ipk feb03      common /mapc/imap(maxpl),NCRS(MAXPL) 
!ipk sep97 add NCRS above                                               
!                                                                       
!ipknov93      CHARACTER*1 ANS,ANSW(10)                                 
      CHARACTER*1 ANS,ANSW(10),IFLAG 
      CHARACTER*63  STRELS 
      DATA STRELS/' You have tried set to set elevation with no mapfile"'/ 
!                                                                       
      DATA ANSW/'m','a','f','s','k','u','t','w','h','q'/ 
!JUN08	  DATA ISWTAGN/0/
!ipk feb94 add DATA and FUNCTION below                                  
!      DATA ARF/-180.,-90.,0.,90.,180./                                 
!      ANGN(K,L)=                                                       
!     +  ATAN2((CMAP(K,2)-CORD(L,2)),(CMAP(K,1)-CORD(L,1)))*57.296      
!                                                                       
!  Draw box around selections                                           
!                                                                       
!IPK SEP97                                                              
100   CONTINUE 
      IDONET=0
      NHTP = 9 
      NMESS = 0 
      NBRR = 0 
      IERREL=0
      IF(ISWTIN .EQ. -1) GO TO 190
      CALL HEDR 
!                                                                       
!  Get answer                                                           
!                                                                       
  110 call xyloc(XPT,YPT,ANS,IBOX) 
      IF(IRMAIN .EQ. 1) RETURN 
      IF(ANS .EQ. 'c') THEN
        if(ibox .eq. 0) go to 110 
        ANS=ANSW(IBOX) 
      ENDIF 
      IF(ANS .EQ. 'm') THEN 
!                                                                       
!  This option allows changes to bottom elevations                      
!                                                                       
        CALL ADDPTH 
        IF(IRMAIN .EQ. 1) RETURN 
        GO TO 100 
                                                                        
      ELSEIF (ANS .EQ. 'a')  THEN 
!                                                                       
!  All nodes                                                            
!                                                                       
        ISWT = -1 
        DO N=1,NP
          IF(INEW(N) .EQ. 1) WD(N)=-9999.
        ENDDO
      ELSEIF(ANS .EQ. 'f') THEN 
!                                                                       
!  Fill nodes                                                           
!                                                                       
        ISWT = 0 
      ELSEIF(ANS .EQ. 's') THEN 
!                                                                       
!  Single node at a time                                                
!                                                                       
        ISWT = 1 

!ipk feb02 add lock/unlock and remove cdata
                                                                        
      ELSEIF(ANS .EQ. 'k') THEN 
!                                                                       
!   lock node
!                                                                       
!  Get M from mouse                       
!                                                                       
  115     CONTINUE
          NHTP=0 
          NMESS=21 
          NBRR=3 
          CALL HEDR 
          IBOX=1 
          CALL PROX(CORD(1,1),CORD(1,2),NP,XX,YY,M,IFLAG,INSKP,IBOX) 
          IF(IRMAIN .EQ. 1) RETURN 
          if(iflag .eq. 'q') go to 100
          lock(m)=1
          go to 115
      ELSEIF(ANS .EQ. 'u') THEN 
!                                                                       
!   unlock node
!                                                                       
!  Get M from mouse                       
!                                                                       
  120     CONTINUE
          NHTP=0 
          NMESS=21 
          NBRR=3 
          CALL HEDR 
          IBOX=1 
          CALL PROX(CORD(1,1),CORD(1,2),NP,XX,YY,M,IFLAG,INSKP,IBOX) 
          IF(IRMAIN .EQ. 1) RETURN 
          if(iflag .eq. 'q') go to 100
          lock(m)=0
          go to 120
      ELSEIF(ANS .EQ. 't') THEN 
!                                                                       
!  Create data for layers                                               
!                                                                       
        CALL ADDLAY 
        IF(IRMAIN .EQ. 1) RETURN 
        GO TO 100 
      ELSEIF(ANS .EQ. 'w') THEN 
!                                                                       
!  This option allows changes to nodal widths                           
!                                                                       
        CALL ADDWID 
        IF(IRMAIN .EQ. 1) RETURN 
        GO TO 100 
!                                                                       
!  Call to help screen                                                  
!                                                                       
      ELSEIF(ANS .EQ. 'h') THEN 
        CALL HELPS(4) 
        IF(IRMAIN .EQ. 1) RETURN 
        GO TO 100 
!                                                                       
      ELSEIF(ANS .EQ. 'q') THEN 
!                                                                       
!  Writeout and return                                                  
!                                                                       
        CALL WRTOUT(0) 
        RETURN 
      ENDIF 
      
190   CONTINUE
      
      IF(IMP .EQ. 0) THEN
        CALL SYMBL(0.,7.25,0.20,STRELS,0.,63) 
        go to 100
      endif
!                                                                       
!  Establish size for range                                             
!               
      call setrng(xnears,nmap)

      ITIME=0
      ICOUNTF=0
      MM=0 
  200 MM=MM+1 
!      write(90,*) 'gridsb-111',mm,np,iswt,inew(mm)
      IF(MM .LE. NP) THEN 
!                                                                       
!  Decode which alternative we are processing
! ipk feb 03 determine interpolation method                           
!
        IF(MM .EQ. 1  .AND.  ISWTAGN .EQ. 0) THEN

          IF(IRECD .EQ. 2) THEN
	        iswtintp=0
	        iswtagn=0
            go to 210		    
          ENDIF
          IF(IGUNIT .EQ. 203) THEN
            ISWTINTP=1
	        iswtagn=1
            GO TO 210
          ENDIF
          CALL WMessageBox(YesNo,QuestionIcon,CommonOK,'Do you wish to interpolate '//&
          CHAR(13)//'from the triangulated map file?'//&
           'or from a grid file',&
                 'Select Interpolation method?')
!   If answer 'Yes'  set interpolate switch to 1
!
          IF (WInfoDialog(4) .EQ. 2) then
	        iswtintp=0
	      ELSE
	        iswtintp=1
          ENDIF    

          CALL WMessageBox(YesNo,QuestionIcon,CommonOK,'Ask this question again?'//&
          CHAR(13)//'  ' ,&
                 'Ask again?')
!   If answer 'Yes'  set again switch to 0
!
          IF (WInfoDialog(4) .EQ. 2) then
	        iswtagn=1
          ELSE
	        iswtagn=0
          ENDIF    
	    ENDIF

  210   CONTINUE

        IF(iswtintp .eq. 1  .AND. IGUNIT .NE. 203) then
	      if(iswt .ne. 1) then

!     Sort the x-coordinates

            call sortdb(xusr,ncrs,np)
		  else
		    ncrs(mm)=mm
		  endif
	      m=ncrs(mm)
        else
	      m=mm
	    endif 
! IPK OCT 2 1991                                                        
        IF(ISWT .EQ. 1) THEN
!  Single node at a time         ISWT = 1 
!                                                                       
!  Get M from mouse and set MM to NP                                    
!                                                                       
          NHTP=0 
          NMESS=21 
!ipk jun08          NBRR=0 
          NBRR=1 
          CALL HEDR 
          IBOX=1 
          CALL PROX(CORD(1,1),CORD(1,2),NP,XX,YY,INODE,IFLAG,INSKP,IBOX) 
          IF(IRMAIN .EQ. 1) RETURN 
          if(iflag .eq. 'q') go to 100
          M=INODE 
          MM=NP 
	    endif
        IF(INEW(M) .EQ. 0) GO TO 200 
! IPK END OCT 2 1991                                                    
                                     
        IF(ISWT .EQ. -1) THEN 
!  All nodes           ISWT = -1 
!ipk feb02
          if(lock(m) .eq. 1) go to 200

        ELSEIF(ISWT .EQ. 0) THEN 
!  Fill nodes          ISWT = 0 
!ipk feb02
          IF(WD(M) .GT. -9000.  .or.  lock(m) .eq. 1) go to 200

        ENDIF
!        write(90,*) 'gridsb-138', m,mm,iswt,wd(m),xnears

        IF(ISWTINTP .EQ. 0) THEN
          if(lock(m) .eq. 0) CALL SETELV(XNEARS,NMAP,M,ISWT)
		ELSE
          if(nelts .eq. 0  .and.  igunit .ne. 203) then
            CALL WMessageBox(YesNo,QuestionIcon,CommonOK,'No triangulated exists'//&
            CHAR(13)//'Do you wish to triangulate now?' ,&
                 'NO TRIANGULATION AVAILABLE?')
!   If answer 'Yes'  set triangulate now
!
            IF (WInfoDialog(4) .EQ. 2) then
	          return
	        ELSE
	          call triang
              IDONET=1
            ENDIF    

          endif
          if(lock(m) .eq. 0) then
            if(igunit .ne. 203) then
              CALL TRIANINT(NMAP,M,ISWT,ITIME)
            else
              call GETGRDELEV(M,IERREL)
            endif  
          endif
        ENDIF

!        write(90,*) 'gridsb-141', m,iswt,wd(m)
        if(wd(m) .lt. -9997.) THEN
          icountf=icountf+1
          WD(M)=-9998.
        ENDIF
        GO TO 200 
      ENDIF 
      IF(IDONET .EQ. 1) THEN
        CALL RDMAP(2,99,0,0) ! XXXXX 
        CLOSE(99)
      ENDIF
      
      CALL WMessageBox(YesNo,QuestionIcon,CommonOK,'Do you wish to plot contours?'//&
        CHAR(13)//' ','PLOT CONTOURS?')
!
!   If answer 'No', return
!
      IF (WInfoDialog(4).EQ.2) THEN
        go to 220     
      ENDIF
215   menus=13
      call conout(menus)
      MENUS=12
      CALL CONOUT(MENUS)

!ipkjan94      IF(ISWT .EQ. -1) GO TO 210                               
220 if(icountf .gt. 0) then
        
        CALL FMESS(ICOUNTF,ISWTT)
!
!   If answer 'Yes', use search for adjacent nodes
!
        IF(ISWTT .EQ. 1) then
          call fillin(icountf)
          CALL WMessageBox(YesNo,QuestionIcon,CommonOK,'Do you wish to plot contours again?'//&
          CHAR(13)//' ','PLOT CONTOURS?')
!
!   If answer 'No', return
!
          IF (WInfoDialog(4).EQ.2) THEN
            IF(ICOUNTF .GT. 0) GO TO 220
          ELSE
            GO TO 215
          ENDIF
        END IF
      endif
      IF(ISWTIN .EQ. -1) RETURN 
      IF(ISWT .EQ. 1) THEN
!ipk jun08        CALL PROX(CORD(1,1),CORD(1,2),NP,XX,YY,INODE,IFLAG,INSKP,IBOX) 
        IF(IRMAIN .EQ. 1) RETURN 
        if(iflag .eq. 'q') go to 100
        M=INODE 
        MM=NP-1
		GO TO 200
	  ELSEIF(ISWT .EQ. -1) THEN
	    GO TO 100
	  ENDIF
!ipk jun08      IF(ABS(ISWT) .EQ. 1) GO TO 100 
      RETURN 
      END                                           
!ipk jul98 revise call                                                  
!IPK      SUBROUTINE GRIDIN(I,SOLN,LISTM,NLG)                           
      SUBROUTINE GRIDIN(XZ,YZ,SOLN,LISTM,NLG) 
!                                                                       
!  Routine to interpolate values from map to node points                
!                                                                       
!       I is the location in the CORD array to be interpolated          
!       SOLN is the interpolated value developed                        
!       NLG is the number of entries in the map array                   
                                                                        
      USE BLKMAP
      USE BLK1MOD
!      INCLUDE 'BLK1.COM' 

!IPK JAN94      INTEGER*2 LISTM,NLIST,NKEY                              
!IPK DEC09      INTEGER LISTM,NLIST,NKEY 
      INTEGER LISTM,NLIST 
      DIMENSION LISTM(*) 
      DIMENSION NLIST(1000),ADIST(1000),WT(1000) 
!                                                                       
!  Function statements                                                  
!                                                                       
      DIST(K,X,Y)=SQRT((CMAP(K,1)-X)**2+(CMAP(K,2)-Y)**2) 
!IPK FEB97      ANG(K,X,Y)=ATAN2((CMAP(K,2)-Y),(CMAP(K,1)-X))*57.296    
      ANG(K,X,Y)=ATAN2((CMAP(K,1)-X),(CMAP(K,2)-Y))*57.296 
!                                                                       
!  Initialize                                                           
!                                                                       
      TOL=120.
!      DO KK=1,NLG
!        WRITE(90,*) 'LISTM',KK,LISTM(KK) 
!      ENDDO
!                                                                       
!  Form list of distances from I to data locations                      
!                                                                       
      NTMP=0 
      IPAS=1 
!IPK JUL98      X=CORD(I,1)                                             
!PK  JUL98      Y=CORD(I,2)                                             
      X=XZ 
      Y=YZ 
  140 CONTINUE 
      DO 150 KK=1,NLG 
        K=LISTM(KK) 
        IF(K .EQ. NTMP) THEN 
          ADIST(KK)=-VOID 
        ELSE 
          ADIST(KK)=DIST(K,X,Y) 
        ENDIF 
  150 END DO 
!                                                                       
!  Sort order for nearest points                                        
!                                                                       
      CALL SORT(ADIST,NKEY,NLG) 
       do nnnn=1,nlg                                                    
         nn=nkey(nnnn)                                                  
         nzz=listm(nn)
         WRITE(90,*) NZZ,X,Y,cmap(nzz,1),cmap(nzz,2)
         ATMP=ANG(NZZ,X,Y)                                                 
         write(90,*) nnnn,adist(nnnn),val(nzz),ATMP  
       enddo                                                            
!      read(*,*) aaa                                                    
      KK=1 
!                                                                       
!  Search through sorted list                                           
!                                                                       
      INIT=1 
  180 DO 260 K=1,NLG 
        NN=NKEY(K) 
        IF(NN .EQ. 0) GO TO 260 
        N=LISTM(NN) 
!                                                                       
!  Initialize                                                           
!                                                                       
        IF(N .EQ. NTMP) GO TO 260 
        IF(INIT .EQ. 1) THEN 
          NLIST(1)=N 
          YY=(CMAP(N,2)-Y) 
          XX=(CMAP(N,1)-X) 
          IF(YY .EQ. 0.  .AND.  XX .EQ. 0.) THEN 
            RANGEF=0. 
            RANGEB=360. 
          ELSE 
            RANGEF=ANG(N,X,Y) 
            RANGEB=ANG(N,X,Y)+360. 
          ENDIF 
          INIT=2 
          GO TO 260 
        ENDIF 
!                                                                       
!  Skip out if already processed                                        
!                                                                       
        YY=(CMAP(N,2)-Y) 
        XX=(CMAP(N,1)-X) 
        IF(YY .EQ. 0.  .AND.  XX .EQ. 0.) THEN 
          ANGLE=0. 
        ELSE 
          ANGLE=ANG(N,X,Y) 
        ENDIF 
  200   CONTINUE 
!        write(90,*) 'angls',n,angle,rangef,rangeb,val(n)               
!                                                                       
!  Test if angle greater than high value                                
!                                                                       
        IF(ANGLE .GT. RANGEF) THEN 
!                                                                       
!  Test if the new point is within the allowable range                  
!                                                                       
          IF(ANGLE .LT. RANGEF+TOL) THEN 
!                                                                       
!  Process this point                                                   
!                                                                       
            RANGEF=ANGLE 
            KK=KK+1 
            NLIST(KK)=N 
            NKEY(K)=0 
!                                                                       
!  Test if we now have enough points to exit                            
!                                                                       
            IF(RANGEF+TOL .GT. RANGEB) THEN 
              GO TO 300 
            ELSE 
              GO TO 180 
            ENDIF 
!                                                                       
!  Test if angle lies inside the already spanned area                   
!  If so it cannot be used                                              
!                                                                       
          ELSEIF(ANGLE .GT. RANGEB) THEN 
            NKEY(K)=0 
!                                                                       
!  Test if it close enough to the low value                             
!                                                                       
          ELSEIF(ANGLE .GT. RANGEB-TOL) THEN 
!                                                                       
!  Process this point                                                   
!                                                                       
            RANGEB=ANGLE 
            KK=KK+1 
            NLIST(KK)=N 
            NKEY(K)=0 
!                                                                       
!  Test if we have enough points to exit                                
!                                                                       
            IF(RANGEF+TOL .GT. RANGEB) THEN 
              GO TO 300 
            ELSE 
              GO TO 180 
            ENDIF 
!                                                                       
!  Not a usable point at this time,  move on to next point              
!                                                                       
          ELSE 
            GO TO 260 
!                                                                       
!  Increase angle by 360 and test again                                 
!                                                                       
          ENDIF 
        ELSE 
          ANGLE=ANGLE+360. 
          GO TO 200 
        ENDIF 
  260 END DO 
!                                                                       
!  We have finished loop without completing polygon                     
!                                                                       
      GO TO 500 
!                                                                       
!  Process least squares fit on this list                               
!                                                                       
  300 CONTINUE 
!        WRITE(90,*) 'least squares list',KK,x,y                        
!        WRITE(90,*) (NLIST(N),N=1,KK)                                  
!        do n=1,kk                                                      
!        write(90,*) nlist(n),cmap(nlist(n),1),cmap(nlist(n),2)         
!     +            ,val(nlist(n)),dist(nlist(n),x,y)                    
!        enddo                                                          
!ipk feb97 changes to refine processing                                 
!                                                                       
!    Check if points are close together relative to the centre point    
!                                                                       
!      write(90,*) kk,x,y,nlg
      do n=1,kk 
        l=nlist(n) 
        dc=dist(l,x,y) 
        xx=cmap(l,1) 
        yy=cmap(l,2) 
        if(n .lt. kk) then 
          do m=n+1,kk 
            ll=nlist(m) 
            dr=dist(ll,xx,yy) 
            if(dr .lt. 0.1*dc) then 
              if(kk .gt. 3) then 
                ds=dist(ll,x,y) 
                if(ds .lt. dc) then 
                  ndrp=n 
                else 
                  ndrp=m 
                endif 
!                                                                       
!     drop this point                                                   
!                                                                       
                do mm=ndrp,kk-1 
                  nlist(mm)=nlist(mm+1) 
                enddo 
                kk=kk-1 
                go to 300 
              else 
                go to 310 
              endif 
            endif 
          enddo 
        endif 
      enddo 
  310 continue 
!ipk feb97 end changes for processing                                   
!        WRITE(90,*) '310',kk                                           
!        WRITE(90,*) (NLIST(N),N=1,KK)                                  
!ipk feb97 chnage to add weighting                                      
       do n=1,kk 
!ipk jul98         if(dist(nlist(n),CORD(I,1),CORD(I,2)) .gt. 0.) then  
!ipk jul98           wt(n)=1./dist(nlist(n),CORD(I,1),CORD(I,2))        
         if(dist(nlist(n),XZ,YZ) .gt. 0.) then 
           wt(n)=1./dist(nlist(n),XZ,YZ) 
         else 
           soln=val(nlist(n)) 
           return 
         endif 
       enddo 
!IPK JUL98      CALL ALSQ(KK,NLIST,I,SOLN,WT)                           
      CALL ALSQ(KK,NLIST,XZ,YZ,SOLN,WT) 
!ipk feb97 end changes                                                  
!                                                                       
!  final value is SOLN                                                  
!                                                                       
      RETURN 
  500 TOL=TOL+25. 
      IF(TOL .GT. 180.) GO TO 550 
      IF(RANGEF+TOL .GT. RANGEB) THEN 
        GO TO 300 
      ENDIF 
      GO TO 180 
  550 CONTINUE 
!c      write(90,*) ' in trouble split',rangef,rangeb                   
      SPLIT=(RANGEF+RANGEB)/2.-180. 
      AMIN=180. 
      DO 600 N=1,KK 
        IF(NLIST(N) .EQ. NTMP) GO TO 600 
        YY=(CMAP(NLIST(N),2)-Y) 
        XX=(CMAP(NLIST(N),1)-X) 
        IF(YY .EQ. 0.  .AND.  XX .EQ. 0.) THEN 
          ANGL=0. 
        ELSE 
          ANGL=ANG(NLIST(N),X,Y) 
        ENDIF 
        IF(IPAS .EQ. 2) GO TO 600 
!                                                                       
!  Find line closest to split                                           
!                                                                       
        IF(ABS(SPLIT-ANGL) .LT. AMIN) THEN 
          AMIN=ABS(SPLIT-ANGL) 
!          write(90,*) 'ntmp reset',ntmp,nlist(n),amin,split            
          NTMP=NLIST(N) 
        ENDIF 
          ANGLP=ANGL-360. 
        IF(ABS(SPLIT-ANGLP) .LT. AMIN) THEN 
          AMIN=ABS(SPLIT-ANGLP) 
         NTMP=NLIST(N) 
        ENDIF 
!  580   WRITE(90,*) NLIST(N),ANGL                                      
  600 END DO 
      IF(IPAS .EQ. 1) THEN 
        IPAS=2 
        X=CMAP(NTMP,1) 
        Y=CMAP(NTMP,2) 
        TOL=120. 
        write(90,*) 'INTERP FOR ',xz,yz,' MOVED TO',x,y,ntmp 
        GO TO  140 
      ENDIF 
      WRITE(90,*) 'ERROR NO POLYGON RANGEF,RANGEB',RANGEF,RANGEB,SPLIT 
      WRITE(90,*) 'OPPOSITE NODE AND ANGULAR DIFF',NTMP,AMIN 
      SOLN=-9998. 
      RETURN 
      END                                           
!                                                                       
!      FUNCTION ANG(K,X,Y)                                              
!                                                                       
!      INCLUDE 'BLK1.COM'                                               
!                                                                       
!      YY=(CMAP(K,2)-Y)                                                 
!      XX=(CMAP(K,1)-X)                                                 
!      IF(YY .EQ. 0.  .AND.  XX .EQ. 0.) THEN                           
!        ANG=0.                                                         
!      ELSE                                                             
!        ANG=ATAN2(YY,XX)*57.296                                        
!      ENDIF                                                            
!      RETURN                                                           
!      END                                                              
                                                                        
      SUBROUTINE SORT(A,NKEY,N) 
!*********************************** .....SORT.....                     
!-                                                                      
!......SORT IS A SIMPLE SHELL SORT ROUTINE                              
!-                                                                      
!     SHELL SORT                                                        
      SAVE 
!                                                                       
!IPK JAN94      INTEGER*2 NKEY                                          
      DIMENSION A(*),NKEY(1) 
      IF(N.LT.2) RETURN 
      DO 90 J=1,N 
        NKEY(J)=J 
   90 END DO 
      ID = N 
  100 ID = ID / 2 
  110 IB = 1 
  120 GO TO 200 
  130 IB = IB + 1 
      IF( IB .LE. ID ) GO TO 200 
      IF( ID .GT. 1 ) GO TO 100 
      RETURN 
  200 I = IB 
  210 K = I + ID 
  220 IF( A(NKEY(I)) .LE. A(NKEY(K)) ) GO TO 250 
      NKT = NKEY(K) 
      NKEY(K) = NKEY(I) 
      J = I 
  230 K = J - ID 
      IF( K .LT. 1 ) GO TO 240 
      IF( A(NKT) .GT. A(NKEY(K)) ) GO TO 240 
      NKEY(J) = NKEY(K) 
      J = K 
      GO TO 230 
  240 NKEY(J) = NKT 
  250 I = I + ID 
      IF( I + ID .LE. N ) GO TO 210 
      GO TO 130 
      END                                           
!ipk feb97 add weighting                                                
!iok jul98      SUBROUTINE ALSQ(NPTS,NLIST,I,SOLN,WT)                   
      SUBROUTINE ALSQ(NPTS,NLIST,xx,yy,SOLN,WT) 
!                                                                       
!  Least squares routine                                                
!                                                                       
!      INCLUDE 'PARAM.COM'                                              
      USE BLKMAP
      USE BLK1MOD
!      INCLUDE 'BLK1.COM' 
!IPK JAN94      INTEGER*2 NLIST                                         
      REAL*8 A,R,B,S,X,Y,ATR,ATR2,BTR,C,T,X3,X2,X1
      DIMENSION A(3,3),R(3),B(2,2),S(2),wt(*) 
      DIMENSION NLIST(*) 
!                                                                       
!  Initialize matrices                                                  
!                                                                       
!ipk jul98      X=CORD(I,1)                                             
!ipk jul98      Y=CORD(I,2)                                             
       X=XX 
       Y=YY 
!      write(*,*) (nnn,cmap(nnn,1),cmap(nnn,2),nnn=1,16)                
!      write(*,*) (nlist(n),n=1,npts)                                   
      DO 160 K=1,3 
        R(K)=0. 
        DO 150 J=1,3 
          A(J,K)=0. 
  150   CONTINUE 
  160 END DO 
!                                                                       
!  Form A and R matrices                                                
!                                                                       
      DO 200 N=1,NPTS 
        KK=NLIST(N) 
!        write(*,*) cmap(kk,1),cmap(kk,2),val(kk)                       
!ipk feb97 add weighting                                                
        A(1,1)=A(1,1)+1.0*wt(n) 
        A(1,2)=A(1,2)+CMAP(KK,1)*wt(n) 
        A(1,3)=A(1,3)+CMAP(KK,2)*wt(n) 
        A(2,2)=A(2,2)+CMAP(KK,1)**2*wt(n) 
        A(2,3)=A(2,3)+CMAP(KK,1)*CMAP(KK,2)*wt(n) 
        A(3,3)=A(3,3)+CMAP(KK,2)**2*wt(n) 
        R(1)=R(1)+VAL(KK)*wt(n) 
        R(2)=R(2)+CMAP(KK,1)*VAL(KK)*wt(n) 
        R(3)=R(3)+CMAP(KK,2)*VAL(KK)*wt(n) 
!ipk feb97 end addition of weighting                                    
  200 END DO 
!      read(*,*) al                                                     
!                                                                       
!  Solve equations                                                      
!                                                                       
      ATR=A(1,2)/A(1,1) 
      ATR2=A(1,3)/A(1,1) 
      B(1,1)=A(2,2)-ATR*A(1,2) 
      B(1,2)=A(2,3)-ATR*A(1,3) 
      S(1)=R(2)-ATR*R(1) 
      B(2,2)=A(3,3)-ATR2*A(1,3) 
      S(2)=R(3)-ATR2*R(1) 
      BTR=B(1,2)/B(1,1) 
      C=B(2,2)-BTR*B(1,2) 
      T=S(2)-BTR*S(1) 
      X3=T/C 
      X2=S(1)/B(1,1)-BTR*X3 
      X1=R(1)/A(1,1)-ATR*X2-ATR2*X3 
!                                                                       
!  Substitute to get interpolated value                                 
!                                                                       
      SOLN=X1+X2*X+X3*Y 
      RETURN 
      END                                           
!                                                                       
!ipksep97 new routine for soring map lines                              
!                                                                       
      SUBROUTINE SORTMAP(A,NKEY,N,IMAP) 
!*********************************** .....SORT.....                     
!-                                                                      
!......SORT IS A SIMPLE SHELL SORT ROUTINE                              
!-                                                                      
!     SHELL SORT                                                        
      SAVE 
!                                                                       
!IPK JAN94      INTEGER*2 NKEY                                          
      DIMENSION A(*),NKEY(1),IMAP(*) 
      DATA VOID/1.E35/ 
      IF(N.LT.2) RETURN 
      DO 90 J=1,N 
        NKEY(J)=J 
        IF(IMAP(J) .LT. 0) A(J)=VOID 
   90 END DO 
      ID = N 
  100 ID = ID / 2 
  110 IB = 1 
  120 GO TO 200 
  130 IB = IB + 1 
      IF( IB .LE. ID ) GO TO 200 
      IF( ID .GT. 1 ) GO TO 100 
      RETURN 
  200 I = IB 
  210 K = I + ID 
  220 IF( A(NKEY(I)) .LE. A(NKEY(K)) ) GO TO 250 
      NKT = NKEY(K) 
      NKEY(K) = NKEY(I) 
      J = I 
  230 K = J - ID 
      IF( K .LT. 1 ) GO TO 240 
      IF( A(NKT) .GT. A(NKEY(K)) ) GO TO 240 
      NKEY(J) = NKEY(K) 
      J = K 
      GO TO 230 
  240 NKEY(J) = NKT 
  250 I = I + ID 
      IF( I + ID .LE. N ) GO TO 210 
      GO TO 130 
      END                                           
      subroutine fillin(icountf)
      USE BLKMAP
      USE BLK1MOD
      USE BLK2MOD
      DIST(N,M)=(cord(n,1)-cord(m,1))**2+(cord(n,2)-cord(m,2))**2      
      CALL KCON(0)
      MCOUNT=0
      MCOUNTF=0
      DO N=1,NP
        IF(WD(N) .LT. -9997. .and. WD(N) .GT. -9998.5) THEN
          MCOUNT=MCOUNT+1
          DISTCUR=1.E20
          NADJCT=0
          DO K=1,NDELM(N)
            J=NECON(N,K)
            DO I=1,NCORN(J)
              NC=NOP(J,I)
              IF(NC .NE. 0  .AND.  NC .NE. N) THEN
                IF(WD(NC) .GT. -9997.) THEN
                  distance=dist(n,nc)
                  if(distance .lt. distcur) then
                    distcur=distance
                    nadjct=nc
                  endif
                ENDIF
              ENDIF
            ENDDO
          ENDDO
          if(nadjct .gt. 0) then
            wd(n)=wd(nadjct)
          else
            mcounfT=mcountf+1
          ENDIF
          if(mcount .eq. icountf) THEN
            ICOUNTF=MCOUNTF
            return
          ENDIF
        endif
      enddo        
      ICOUNTF=MCOUNTF
      return
      end
      
      SUBROUTINE FMESS(N1,N2)
      use winteracter

      implicit none

      include 'D.inc'
      INCLUDE 'BFILES.I90'

!
! Declare window-type and message variables
!
      TYPE(WIN_STYLE)    :: WINDOW

      TYPE(WIN_MESSAGE) :: MESSAGE

      integer  :: N1,N2,IERR
!      real     :: 
      character*3 :: sub

      call wdialogload(IDD_FBED)
      ierr=infoerror(1)

      CALL WDialogPutInteger(idf_integer1,n1)


      CALL WDialogSelect(IDD_FBED)
      ierr=infoerror(1)

      CALL WDialogShow(-1,-1,0,Modal)
      ierr=infoerror(1)
      
      DO
        IF (WInfoDialog(ExitButton) .EQ. IDOK) THEN
          N2=1
          RETURN
        ELSEIF(WInfoDialog(ExitButton) .EQ. IDCANCEL) THEN
          N2=0
          RETURN
        ENDIF
      ENDDO
      RETURN
      END
    