!ipk  last update sep 20 2013  add more output of progress and flushing of messages
      SUBROUTINE ADDTOMESH(IADDFIL,ISWT)

!     iswt = 0 ADD TO MESH
!     ISWT = 1 MERGE MESHES

      USE WINTERACTER
      USE BLK1MOD 
      USE BLK2MOD 

      INCLUDE 'D.INC'

!      INCLUDE 'BLK1.COM'
      INCLUDE 'BFILES.I90'

      IADD=IADDFIL+50
      CALL RDTOCLIP(IADD)

      IF(ISWT .EQ. 1) THEN
        CALL OUTLINES(1)  
        ISWT1=0
!        IF(NOUTLST(2) .EQ. 0) THEN
          ISWT2=1
!        ELSE
!          ISWT2=0
!        ENDIF
        CALL MERGEMESH1(ISWT1,ISWT2)
        write(90,*) 'finished mergemesh1'
        IF(ISWT2 .EQ. 0) CALL MERGEMESH
!         CALL MERGEMESH        
        write(90,*) 'finished mergemesh'
        flush(90)
      ENDIF

      CALL ADDMESH(0)
      write(90,*) 'finished addmesh'

      IF(ISWT .EQ. 1 ) THEN
        CALL WMessageBox(YesNo,QuestionIcon,CommonOK,'Do you wish to delete unused nodes?'//&
                           CHAR(13)//' ','Delete unused nodes?')
        !
!   If answer 'No', return
!
        IF (WInfoDialog(4).EQ.2) return
!                                                                       
!     Delete all unused nodes                                           
!                                                                       
        CALL DELETM(2) 
      ENDIF

      RETURN
      END


        SUBROUTINE RDTOCLIP(IUNIT)

        USE BLK1MOD 
!        INCLUDE 'BLK1.COM' 
        CHARACTER*80 ALINE 

        REWIND IUNIT
        READ(IUNIT) TITLE,NPSTO(1),NESTO(1) 
        WRITE(90,*) 'IN RDTOCLIP',IUNIT
		WRITE(90,*) TITLE,NPSTO(1),NESTO(1) 
        READ(IUNIT) ISLP,IPRT,IPNN,IPEN,IPO,IRO,IPP,IRFN              &
     &  ,IGEN,NXZL,NITST,ISCTXT,IFILL,IALTGM,NLAYD,xadded,yadded,ntempinc                      
        WRITE(90,*) ISLP,IPRT,IPNN,IPEN,IPO,IRO,IPP,IRFN              &
     &  ,IGEN,NXZL,NITST,ISCTXT,IFILL,IALTGM,NLAYD,xadded,yadded,ntempinc                      
        READ(IUNIT) HORIZ,VERT,XSALE,YSALE,XFACT,YFACT,AR,ANG 
        WRITE(90,*) HORIZ,VERT,XSALE,YSALE,XFACT,YFACT,AR,ANG 
        IF(IPP .GT. 0) READ(IIN) ALINE 

        READ(IUNIT) ((NOPSTO(J,K,1),K=1,8),IMATSTO(J,1),THTASTO(J,1),J=1,NESTO(1)) 

        READ(IUNIT)                                                   &
     &   (XUSRSTO(J,1),YUSRSTO(J,1),WDSTO(J,1),WIDTHSTO(J,1),SS1STO(J,1),SS2STO(J,1),WIDSSTO(J,1),         &
     &    WIDBSSTO(J,1),SSOSTO(J,1),BS1STO(J,1),J=1,NPSTO(1))                                       

        READ(IUNIT) NLSTSTO(1)
        IF(NLSTSTO(1) .GT. 0) THEN 
          READ(IUNIT) (LLISTSTO(J,1),J=1,NLSTSTO(1)),                            &
          ((ILISTSTO(J,I,1),I=1,LLISTSTO(J,1)),J=1,NLSTSTO(1))                        
        ENDIF 

        READ(IUNIT) NENTRYC,NLAYDC,NCLMSTO(1)
        IF(NENTRYC .GT. 0) THEN 
          READ(IUNIT) ((NEFC,J=1,3),I=1,NENTRYC) 
        ENDIF 
        IF(NLAYDC .GT. 0) THEN 
          READ(IUNIT) (LAYC,I=1,NPSTO(1)) 
        ENDIF 
        IF(NCLMSTO(1) .GT. 0) THEN 
          READ(IUNIT) ((ICCLNSTO(I,J,1),J=1,350),I=1,NCLMSTO(1)) 
        ENDIF 

        REWIND IUNIT
  	    RETURN
	    END

        SUBROUTINE ADDMESH(ISWT)

        USE BLK1MOD 
!        INCLUDE 'BLK1.COM' 

        
        ALLOCATABLE NODETRAN(:)
        DATA VDX9/-9.E9/

!       Loop through nodes assigning new number and adding to list

        IF(.NOT. ALLOCATED(NODETRAN)) ALLOCATE (NODETRAN(maxp))
        
        IF(ISWT .EQ. 0) THEN
          DO N=1,NPSTO(1)
            IF(XUSRSTO(N,1) .GT. VDX9) THEN
              CALL GETNOD(J)
              NODETRAN(N)=J
              XUSR(J)=XUSRSTO(N,1)
              YUSR(J)=YUSRSTO(N,1)
              WD(J)=WDSTO(N,1)
              WIDTH(J)=WIDTHSTO(N,1)
              SS1(J)=SS1STO(N,1)
              SS2(J)=SS2STO(N,1)
              WIDS(J)=WIDSSTO(N,1)
              WIDBS(J)=WIDBSSTO(N,1)
              SSO(J)=SSOSTO(N,1)
              BS1(J)=BS1STO(N,1)
              INSKP(J) = 0
              INEW(J) = 1
            ENDIF                                         
          ENDDO 
        ELSE
          DO N=1,NPSTO(1)
            NODETRAN(N)=N
          ENDDO
        ENDIF

!       Loop through elements assigning new number and adding to list

        DO N=1,NESTO(1)
          IF(IMATSTO(N,1) .GT. 0) THEN
            CALL GETELM(M)
            DO K=1,8
              IF(NOPSTO(N,K,1) .GT. 0) THEN
                J=NODETRAN(NOPSTO(N,K,1))
                NOP(M,K)=J
              ELSE
                NOP(M,K)=0
              ENDIF
            ENDDO
            IMAT(M)=IMATSTO(N,1)
            THTA(M)=THTASTO(N,1)
            IESKP(M)=0
            NCN = 2 
            IF (NOP(M,3) .NE. 0) NCN = 3 
            IF (NOP(M,4) .NE. 0) NCN = 4 
            IF (NOP(M,5) .NE. 0  .AND.  NOP(M,4) .NE. 0)  NCN = 5 
            IF (NOP(M,5) .NE. 0  .AND.  NOP(M,4) .EQ. 0)  NCN = 6 
            IF (NOP(M,6) .NE. 0) NCN = 6 
            IF (NOP(M,7) .NE. 0) NCN = 8 
            NCORN(M) = NCN 

          ENDIF
        ENDDO

        if(iswt .eq. 0) CALL RESCAL
  	    CALL HEDR

        RETURN
        END

        SUBROUTINE MERGEMESH1(ISWT1,ISWT2)

        USE BLK1MOD 
        USE BLK2MOD 
        USE WINTERACTER
        
!       INCLUDE 'BLK1.COM' 

        REAL*8 ELXMIN,ELXMAX,ELYMIN,ELYMAX,XLC,YLC,XXX,YYY
        LOGICAL LSTAT

        ALLOCATABLE ELXMIN(:),ELXMAX(:),ELYMIN(:),ELYMAX(:),KEY(:),NKEY(:)
        DIMENSION XOUT1(1000),YOUT1(1000)
        IF(.NOT. ALLOCATED(ELXMIN)) &
        ALLOCATE (ELXMIN(MAXE),ELXMAX(MAXE),ELYMIN(MAXE),ELYMAX(MAXE),KEY(MAXE),NKEY(MAXP))

        IF(ISWT2 .EQ. 0) GO TO 110
!       first eliminate any elements inside outline 
        CALL KCONST(0)
        NKEP=0
        DO K=1,10
          IF(NOUTLST(K) .LE. 0) THEN
            DO J=1,NPSTO(1)
              XXXX=XUSRSTO(J,1)
              YYYY=YUSRSTO(J,1)
              LSTAT=IGrInsidePolygon(XOUT(1,K),YOUT(1,K),-NOUTLST(K),XXXX,YYYY)
              IF(LSTAT) THEN
                NKEP(J)=1
              ENDIF
            ENDDO
          ENDIF
        ENDDO
        DO K=1,10
          IF(NOUTLST(K) .GT. 0) THEN
            DO J=1,NPSTO(1)
              IF(NKEP(J) .EQ. 1) CYCLE
              XXXX=XUSRSTO(J,1)
              YYYY=YUSRSTO(J,1)
!          WRITE(155,*) J,XXXX,YYYY
              LSTAT=IGrInsidePolygon(XOUT(1,K),YOUT(1,K),NOUTLST(K),XXXX,YYYY)
!          WRITE(155,*) J,LSTAT
              IF(LSTAT) THEN
                DO L=1,NDELM(J)
                  NCAN=NECON(J,L)
                  CALL DELEM(NCAN)
                ENDDO
              ENDIF
            ENDDO
          ENDIF
100       CONTINUE
        ENDDO
        IF(ISWT2 .EQ. 1) RETURN
!       First sort coordinates for min of element connection
      
!       List all limiting values
  110   CONTINUE
        DO N=1,NE
          IF(IMAT(N) .GT. 0) THEN
            ELXMIN(N)=XUSR(NOP(N,1))
            ELXMAX(N)=XUSR(NOP(N,1))
            ELYMIN(N)=YUSR(NOP(N,1))
            ELYMAX(N)=YUSR(NOP(N,1))
            DO M=2,8
              IF(NOP(N,M) .NE. 0) THEN
                ELXMIN(N)=MIN(ELXMIN(N),XUSR(NOP(N,M)))    
                ELXMAX(N)=MAX(ELXMAX(N),XUSR(NOP(N,M)))    
                ELYMIN(N)=MIN(ELYMIN(N),YUSR(NOP(N,M)))    
                ELYMAX(N)=MAX(ELYMAX(N),YUSR(NOP(N,M)))    
              ENDIF
            ENDDO
          ELSE
            ELXMIN(N)=VOID
            ELXMAX(N)=VOID
            ELYMIN(N)=VOID
            ELYMAX(N)=VOID
          ENDIF
        ENDDO

        CALL SORTDB(XUSRSTO,NKEY,NPSTO(1))

        CALL SORTDB(ELXMIN,KEY,NE)

!       Loop on elements to check for overlap

        
        DO KK=1,NESTO(1)
          IF (NOPSTO(KK,6,1) .EQ. 0) CYCLE
          IF(IMATSTO(KK,1) .GT. 0) THEN
            if(mod(kk,1000) .eq. 0) write(90,*) 'merged',kk
            flush(90)
            KL=1
  200       CONTINUE
            IF(ISWT1 .EQ. 0) THEN
              DO K=KL,8
                J=NOPSTO(KK,K,1) 
                IF(J .GT. 0) THEN
                  KLL=KL
                  XXX=XUSRSTO(J,1)
                  YYY=YUSRSTO(J,1)
                  GO TO 220
                ENDIF
              ENDDO
              KLL=8
              GO TO 400
  220         CONTINUE
            ELSE
              XXX=0.
              YYY=0.
              DO K=1,7,2
                JJ=NOPSTO(KK,K,1)
                IF(JJ .GT. 0) THEN
                  XXX=XXX+XUSRSTO(JJ,1)
                  YYY=YYY+YUSRSTO(JJ,1)
                ENDIF
              ENDDO
              IF(JJ .EQ. 0) THEN
                XXX=XXX/3.
                YYY=YYY/3.
              ELSE
                XXX=XXX/4.
                YYY=YYY/4.
              ENDIF 
            ENDIF          
!       Search on elements to find a startin point

            DO NN=1,NE

              N=KEY(NN)
              IF(IMAT(N) .GT. 0) THEN
!-                                                                      
!...... DETERMINE ELEMENT TYPE                                           
!-                                                                      
                NCN=8
                IT=1
                IF(NOP(N,7) .EQ. 0) THEN 
                  NCN=6
                  IT=2
                ENDIF
                IF(NOP(N,6) .EQ. 0) THEN
                  GOTO 350
                ENDIF 
!       Test for point inside an element        


!       Test for max and min within

                IF(XXX .GT. ELXMIN(N)) THEN
                  IF(XXX .GT. ELXMAX(N)) GO TO 350                      
                  IF(YYY .GT. ELYMIN(N)) THEN
                    IF(YYY .GT. ELYMAX(N)) GO TO 350

!       Now get local coordinate as final test

                    CALL GPTEV(N,XXX,YYY,XLC,YLC,IT,NCN)

                    IF(IT .EQ. 2) THEN
                      IF(XLC .LT. 0.  .OR.  YLC .LT. 0.  .OR.  XLC+YLC .GT. 1.) THEN
                        GO TO 350
                      ELSE
                        CALL DELEM(KK)
                        GO TO 400
                      ENDIF    
                    ELSE
                      IF(XLC .LT. -1.  .OR.  YLC .LT. -1.  .OR.  &
                         XLC .GT.  1.  .OR.  YLC .GT.  1.) THEN
                        GO TO 350
                      ELSE
                        CALL DELEM(KK)
                        GO TO 400
                      ENDIF                     
                    ENDIF

                  ENDIF
                ENDIF
              ENDIF 
  350         CONTINUE
            ENDDO
            KL=KLL+1
            IF(KL .LT. 8  .AND. ISWT1 .EQ. 0) GO TO 200
          ENDIF

!       Finished test 

  400     CONTINUE
        ENDDO
        RETURN
        END


      SUBROUTINE GPTEV(N,XSW,YSW,XG,YG,IT,NCN)
!-                                                                      
!......EVALUATE FUNCTION AT GRID POINTS                   
!-                                                                      
!-    N = ELEMENT NUMBER        
!_    XSW = X COORDINATE OF DESIRED POINT
!_    YSW = Y COORDINATE OF DESIRED POINT
!     XG  = X LOCAL COORDINATE
!     YG  = Y LOCAL COORDINATE
!     IT   = SWITCH FOR CHOICE BETWEEN LINEAR AND QUADRATIC WEIGHTING
!            = 1 FOR LINEAR
!            = 2 FOR QUADRATIC
!     FROM COMMON
!     NOP = LIST OF NODAL CONNECTIONS AROUND AN ELEMET
!     XUSR = REAL*8 ARRAY OF NODAL COORDINATES
!            

      USE BLK1MOD 
!      INCLUDE 'BLK1.COM'

      REAL*8 XN,DNX,DNY,XSW,YSW 
      DOUBLE PRECISION XG,YG,XK,YK,XP,YP
!-                                                                      
      DIMENSION X(9),Y(9),WGT(8)
!-                                                                      
      DATA TOL/0.01/ 
!-                                                                      

!-                                                                      
!......ESTABLISH LOCAL COORDINATES FOR EACH NODE POINT OF ELEMENT                 
!-                                                                      
      K1=NOP(N,1) 
      X(1)=0. 
      Y(1)=0. 
      DO 300 K=3,NCN,2
        K2=NOP(N,K) 
        X(K)=XUSR(K2)-XUSR(K1) 
        Y(K)=YUSR(K2)-YUSR(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.
      IF(IT .EQ. 2) THEN
        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))
      ELSE
        X(6)=(X(5)+X(7))/2.
        Y(6)=(Y(5)+Y(7))/2.
        X(8)=X(7)/2.
        Y(8)=Y(7)/2.

 	    xminl=min(x(1),x(3),x(5),x(7))
	    yminl=min(y(1),y(3),y(5),y(7))
	    xmaxl=max(x(1),x(3),x(5),x(7))
	    ymaxl=max(y(1),y(3),y(5),y(7))
      ENDIF


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

      XG=0. 
      YG=0. 
!-                                                                      
!......ITERATE TO FIND LOCAL COORDINATE                                 
!-                                                                      
      DO  ITER=1,10 
        DXKDX=0. 
        DXKDY=0. 
        DYKDX=0. 
        DYKDY=0. 
        XK=-XP 
        YK=-YP 
        DO  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) 
        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 
      END DO 
!-                                                                      
!......NOW GET WEIGHTING FUNCTIONS FOR QUAD FUNCTION
!-                                                                      
  420 CONTINUE 


      RETURN
      END

      SUBROUTINE DELEM(J)
!                                                                       
      USE BLK1MOD 
!      INCLUDE 'BLK1.COM' 
!                                                                       
!-                                                                      
!......DELETE ELEMENT 
!                                                                       
!  Search for elements that attach to node J and remove them            
!                      

      IMATSTO(J,1)=0 
      DO  KK=1,8 
        NOPSTO(J,KK,1)=0 
      ENDDO
!                                                                       

      RETURN 
      END                                           
                   

      SUBROUTINE MERGEMESH

      USE BLK1MOD 
      LOGICAL LSTAT
!      INCLUDE 'BLK1.COM' 

!     Loop on element to be added

      DO N=1,NESTO(1)
        
        IF(IMATSTO(N,1) .NE. 0) THEN
            if(mod(n,1000) .eq. 0) write(90,*) 'adding',n,nesto(1)
            flush(90)
          IF(IMATSTO(N,1) .GT. 900  .AND.  IMATSTO(N,1) .LT. 904) THEN
            X1=XUSRSTO(NOPSTO(N,1,1),1)
            Y1=YUSRSTO(NOPSTO(N,1,1),1)
            CALL CHECKIN(X1,Y1,LSTAT)
            IF(ISTATUS .EQ. 5) THEN
              CALL DELEM(N)
              GO TO 400
            ENDIF
            GO TO 400
          ENDIF    
          
!     loop on sides

          DO M=1,7,2
            N1=NOPSTO(N,M,1)
            IF(M .EQ. 3  .AND.  NOPSTO(N,5,1) .EQ. 0) GO TO 400
            IF(N1 .GT. 0) THEN
              IF((M .EQ. 5  .AND.  NOPSTO(N,7,1) .EQ. 0)  .OR. (M .EQ. 7)) THEN
                N2=NOPSTO(N,1,1)
              ELSE  
                N2=NOPSTO(N,M+2,1)
              ENDIF
              IF(NKEP(N1) .EQ. 1  .AND.  NKEP(N2) .EQ. 1) GO TO 380

!     Now loop trough existing elements

              DO I=1,NE
                IF(IMAT(I) .NE. 0) THEN
                  DO J=1,7,2
                    M1=NOP(I,J)
                    IF(J .EQ. 3  .AND.  NOP(I,5) .EQ. 0) GO TO 360
                    IF(M1 .GT. 0) THEN
                      IF((J .EQ. 5  .AND.  NOP(I,7) .EQ. 0)  .OR. (J .EQ. 7)) THEN
                        M2=NOP(I,1)
                      ELSE  
                        M2=NOP(I,J+2)
                      ENDIF
                      if(m2 .eq. 0) cycle
                      X1=XUSRSTO(N1,1)
                      X2=XUSRSTO(N2,1)
                      Y1=YUSRSTO(N1,1)
                      Y2=YUSRSTO(N2,1)
                      X3=XUSR(M1)
                      X4=XUSR(M2)
                      Y3=YUSR(M1)
                      Y4=YUSR(M2)
                      CALL IGrIntersectLine(X1,Y1,X2,Y2,X3,Y3,X4,Y4,XINTER,YINTER,ISTATUS)
                      IF(ISTATUS .EQ. 5) THEN
                        CALL DELEM(N)
                        GO TO 400
                      ENDIF
                    ENDIF
                  ENDDO
                ENDIF
360             CONTINUE                
              ENDDO
            ENDIF
380         CONTINUE              
          ENDDO
        ENDIF
  400   CONTINUE
      ENDDO

      RETURN
      END

      SUBROUTINE CHECKIN(X1,Y1,LSTAT)
      USE BLK1MOD
      LOGICAL LSTAT
      DIMENSION XP(4),YP(4)
!     Now loop trough existing elements

      DO I=1,NE
        IF(IMAT(I) .NE. 0) THEN
          JJ=0
          DO J=1,7,2
            INODE=NOP(I,J)
            IF(INODE .GT. 0) THEN
              JJ=JJ+1
              XP(JJ)=XUSR(INODE)
              YP(JJ)=YUSR(INODE)
            ENDIF
          ENDDO
          LSTAT=IGrInsidePolygon(XP,YP,JJ,X1,Y1)
          IF(LSTAT) RETURN
        ENDIF
      ENDDO
      RETURN
      END
      SUBROUTINE KCONST(isw1) 
!                                                                       
!     ESTABLISH ELEMENT CONNECTED TO ELEMENT TABLE                      
!                                                                       
      USE BLK1MOD
      USE BLK2MOD 
!      INCLUDE 'BLK1.COM' 
!      INCLUDE 'BLK2.COM' 
!                                                                       
!     INITIALIZE                                                        
!                   
      NCM=11
      DO 200 J=1,NCM 
        DO 200 N=1,NPSTO(1) 
  200   NECON(N,J)=0 
      DO 230 N=1,NPSTO(1) 
  230 NDELM(N)=0 
!                                                                       
!     FORM TABLE OF ELEMENTS CONNECTED TO EACH NODE                     
!                                                                       
      DO 300 M=1,NESTO(1)
        IF(IMATSTO(M,1) .EQ. 0) GO TO 300 
		if(isw1 .eq. 1) then
		  if(imat(m) .eq. 999) go to 300
		endif
        DO 280 K=1,8 
          N=NOPSTO(M,K,1) 
          IF (N .GT. 0) THEN 
            NDELM(N)=NDELM(N)+1 
            J=NDELM(N) 
            NECON(N,J)=M 
!ipkoct93          ELSE                                                 
!ipkoct93            GO TO 300                                          
          ENDIF 
  280   CONTINUE 
  300 END DO 
      RETURN
      END
    