      PROGRAM TORSIN9Q
C=======================================================================
C           2-DIM FEM PROGRAM FOR SOLVING TORSION EQUATION
C                  USING UPPER HALF BANDED MATRIX
C              EQUATION: DP/DX + DP/DY + ALPHA * P = 0.
C                ELEMENT : 4-NODED ISO-PARAMETERIC
C                ELEMENT : 9-NODED ISO-PARAMETERIC
C SEQUENCE: (-1,-1),(0,-1),(+1,-1),(+1,0),(+1,+1),(0,+1),
C           (-1,+1),(-1,0),(0,0)
C   ORIGINAL:1984 EIJI FUKUMORI BUFFALO NY & REVISED: 1994 ACHI
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( ND=9,MXE=20000,MXN=20000,MXB=3000,INTEPT=3,MXW=500 )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), IBND(MXB),
     * ITYPE(MXB), BVALUE(MXB),RK(MXN,MXW), RHS(MXN),SK(ND,ND), 
     * X(ND), Y(ND),BX(2,ND),SAI(INTEPT), W(INTEPT),SS(ND,ND),
     * BPP(2,ND,INTEPT,INTEPT), SF(ND,INTEPT,INTEPT),
     * DUMMY0(ND), DUMMY1(ND), XS(ND), YS(ND), BPPX(2,ND,ND)
      DIMENSION  SOURCE(MXN),TAUZX(MXN),TAUZY(MXN),NDREPEAT(MXN)
      CHARACTER INPFILE*12
      WRITE (*,*)' TORSIN9Q.FOR SOLVER'
C=======================================================================
      CALL GRULE ( INTEPT, SAI, W   )
      CALL DERIV ( ND, INTEPT, X, Y, SAI, BPP )
      CALL SHAPEF( ND, INTEPT, X, SAI, SF )
C=======================================================================
      INPFILE = 'TRS9DATA.QQQ'
      WRITE(*,*)' INPUT FILE NAME =', INPFILE
      CALL INPUT ( ND,MXE,MXN,MXB,NE,NNODE,NB,GM,THETA, NODEX,
     * XCOORD,YCOORD,IBND,ITYPE,BVALUE, INPFILE )
C=======================================================================
      CALL BANDWID ( MXE, ND, NE,NODEX, NBW )
      IF ( NBW .GT. MXW ) STOP 'Error ... NBW > MXW'
C=======================================================================
      CALL DAREA ( MXE,MXN,INTEPT,ND,BPP,W,NE,XCOORD,YCOORD,
     * NODEX, AREA )
C=======================================================================
      CALL GSM ( INTEPT,ND,NBW ,MXE,MXN,MXW,NE,NNODE,BPP,W,BX,SS,
     * NODEX,XCOORD,YCOORD,RK )
C=======================================================================
      CALL SRCTRM ( MXE,MXN,INTEPT,ND,NNODE,BPP,W,NE,GM,THETA,
     * SF,XCOORD,YCOORD,NODEX,SS, SOURCE )
C=======================================================================
      CALL FORM (MXN,MXB,MXW,NNODE,NB,NBW,RK,RHS,ITYPE,BVALUE,IBND,
     *                  SOURCE)
C=======================================================================
      CALL SYSTEM ( MXN, MXW, NNODE, NBW, RK, RHS )
C=======================================================================
      CALL MOMENT ( MXE,MXN,INTEPT,ND,BPP,W,NE,SF,
     * XCOORD,YCOORD,NODEX, RHS,TORQUE  )
C=======================================================================
      CALL STRESS ( ND,MXE,MXN,NE,NNODE,BPPX,NODEX,XCOORD,YCOORD,
     * XS, YS, DUMMY0, DUMMY1, RHS, NDREPEAT,TAUZX,TAUZY )
C=======================================================================
      CALL PRINT ( MXN, NNODE , RHS )
      CALL PRINTG ( ND,MXE,MXN, NNODE,NE,  XCOORD, YCOORD,
     *                      GM,THETA, NODEX,MXB, NB, IBND,ITYPE,BVALUE,
     *                      RHS,TORQUE,TAUZX,TAUZY, AREA )
C=======================================================================
      STOP
      END
C
C
      SUBROUTINE ISOPARA ( ND, E1 , E2 , F )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION F(ND)
      F(1) =  0.25*E1*(1.- E1)*E2*(1.- E2)
      F(2) = -0.50*(1.- E1*E1)*E2*(1.- E2)
      F(3) = -0.25*E1*(1.+ E1)*E2*(1.- E2)
      F(4) =  0.50*E1*(1.+ E1)*(1.- E2*E2)
      F(5) =  0.25*E1*(1.+ E1)*E2*(1.+ E2)
      F(6) =  0.50*(1.- E1*E1)*E2*(1.+ E2)
      F(7) = -0.25*E1*(1.- E1)*E2*(1.+ E2)
      F(8) = -0.50*E1*(1.- E1)*(1.- E2*E2)
      F(9) =       (1.- E1*E1)*(1.- E2*E2)     
      RETURN
      END
C
C
      SUBROUTINE GSM ( INTEPT,ND,NBW ,MXE,MXN,MXW,NE,NNODE,BPP,W,BX,SS,
     * NODEX,XCOORD,YCOORD,RK )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),RK(MXN,MXW),
     * SS(ND,ND), BX(2,ND),W(INTEPT), BPP(2,ND,INTEPT,INTEPT)
C====== RESET GLOBAL MATRIX RK ======
      DO I = 1, NNODE
      DO J = 1, NBW
      RK(I,J) = 0.
      END DO
      END DO
C===== INTEGRATION =======
      DO 500 IEL = 1 ,NE
      DO I = 1 , ND
      DO J = I , ND
      SS(I,J) = 0.
      END DO
      END DO
      DO 400 K = 1 , INTEPT
      DO 300 L = 1 , INTEPT
      WEIGHT = W(K) * W(L)
      YAC11 = 0.
      YAC12 = 0.
      YAC21 = 0.
      YAC22 = 0.
      DO I = 1 , ND
      YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODEX(IEL,I))
      YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODEX(IEL,I))
      END DO
      DETJ = YAC11 * YAC22 - YAC12 * YAC21
      Y11 =   YAC22 / DETJ
      Y12 = - YAC12 / DETJ
      Y21 = - YAC21 / DETJ
      Y22 =   YAC11 / DETJ
      DO J = 1 , ND
      BX(1,J) = Y11 * BPP(1,J,K,L) + Y12 * BPP(2,J,K,L)
      BX(2,J) = Y21 * BPP(1,J,K,L) + Y22 * BPP(2,J,K,L)
      END DO
      BETA = WEIGHT * DETJ
      DO I = 1 , ND
      DO J = I , ND
      SS(I,J) = SS(I,J) + ( BX(1,I)*BX(1,J)+BX(2,I)*BX(2,J) )*BETA
      END DO
      END DO
  300 CONTINUE
  400 CONTINUE
      DO I = 1 , ND
      DO J = I , ND
      SS(J,I) = SS(I,J)
      END DO
      END DO
C----ASSEMBLY
      DO K = 1 , ND
      I = NODEX(IEL,K)
      DO L = 1 , ND
      J = NODEX(IEL,L) - I + 1
      IF ( J .GT. 0 ) RK(I,J) = RK(I,J) - SS(K,L)
      END DO
      END DO
  500 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE INPUT ( ND,MXE,MXN,MXB,NE,NNODE,NB,GM,THETA,NODEX,
     * XCOORD,YCOORD,IBND,ITYPE,BVALUE, INPFILE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBND(MXB),
     * ITYPE(MXB), BVALUE(MXB)
      LOGICAL YES
      CHARACTER INPFILE*12
      IR = 1
      INQUIRE ( FILE=INPFILE, EXIST=YES )
      IF ( YES ) OPEN ( IR, FILE=INPFILE, STATUS='OLD' )
      IF ( .NOT. YES ) STOP' NO INPUT FILE'
      READ (IR,*) GM, THETA
      READ (IR,*) NE
      DO I = 1 , NE
      READ (IR,*) IEL,( NODEX(IEL,J),J=1,ND )
      END DO
      READ (IR,*) NNODE
      READ (IR,*) ( NODE,XCOORD(NODE),YCOORD(NODE),J=1,NNODE )
      READ (IR,*) NB
      READ (IR,*) (IBND(I), ITYPE(I), BVALUE(I),I=1,NB)
      CLOSE (IR)
      RETURN
      END
C
C
      SUBROUTINE PRINT ( MXN, NNODE, RHS )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN)
      CHARACTER OUTFILE*12, EXFILE*3
      LOGICAL YES
      OUTFILE ='SOLUTION.BIN'
      EXFILE ='NEW'
      IW = 1
      INQUIRE ( FILE=OUTFILE, EXIST=YES )
      IF ( YES ) EXFILE='OLD'
      OPEN ( IW, FILE=OUTFILE, STATUS=EXFILE,FORM='UNFORMATTED' )
      WRITE (IW) (RHS(I),I=1,NNODE)
      CLOSE (IW)
      RETURN
      END
C
C
      SUBROUTINE PRINTG ( ND,MXE,MXN, NNODE,NE,  XCOORD, YCOORD,
     *                      GM,THETA, NODEX,MXB, NB, IBND,ITYPE,BVALUE,
     *                      RHS,TORQUE,TAUZX,TAUZY, AREA )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN), XCOORD(MXN), YCOORD(MXN),NODEX(MXE,ND),
     *  IBND(MXB),ITYPE(MXB), BVALUE(MXB), TAUZX(MXN),TAUZY(MXN)
      CHARACTER OUTFILE*12, EXFILE*3
      LOGICAL YES
      OUTFILE ='SOLUTION.FEM'
      EXFILE ='NEW'
      IW = 1
      INQUIRE ( FILE=OUTFILE, EXIST=YES )
      IF ( YES ) EXFILE='OLD'
      OPEN ( IW, FILE=OUTFILE, STATUS=EXFILE,FORM='FORMATTED' )
      WRITE (IW,*) 'GM,THETA=', GM,THETA
      WRITE (IW,*)
      WRITE (IW,*) 'ELEMENTS  ', NE
      DO I = 1, NE
      WRITE (IW,*) I, (NODEX(I,J),J=1,ND)
      END DO
      WRITE (IW,*)
      WRITE (IW,*) 'NODES  ',NNODE
      DO I = 1, NNODE
      WRITE (IW,*) I,  XCOORD(I), YCOORD(I)
      END DO
      WRITE (IW,*)
      WRITE (IW,*) 'RESULTS:  AREA-OF-DOMAIN =', AREA
      WRITE (IW,*)
      WRITE (IW,*) 'BOUNDARYS  ',NB
      DO I = 1, NB
      WRITE (IW,*) IBND(I),ITYPE(I), BVALUE(I)
      END DO
      WRITE (IW,*)
      WRITE (IW,*) 'RESULTS: STRESS FUNCTION'
      DO I = 1, NNODE
      WRITE (IW,*) I,  RHS(I)
      END DO
      WRITE (IW,*)
      WRITE (IW,*) 'RESULTS:  TORQUE =', TORQUE
      WRITE (IW,*)
      WRITE (IW,*) 'RESULTS: NODE, TAUZX,  TAUZY'
      DO I = 1, NNODE
      WRITE (IW,*) I,  TAUZX(I),  TAUZY(I)
      END DO
      CLOSE (IW)
      RETURN
      END
C
C
      SUBROUTINE BANDWID ( MXE, ND, NE, NODEX , NBW  )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND)
      NBW = 0
      DO 30 I = 1 , NE
      DO 20 J = 1 , ND-1
      DO 10 K = J+1 , ND
      NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
   10 CONTINUE
   20 CONTINUE
   30 CONTINUE
      NBW = NBW + 1
      write(6,*) ' bandwidth =', nbw
      RETURN
      END
C
C
      SUBROUTINE SYSTEM ( MXN, MXW, NUMNP, MBAND, A, B )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION A(MXN,MXW) , B(MXN)
C---------- ELIMINATION ------------------
      DO N = 1 , NUMNP 
      DO L = 2 , MBAND 
      C = A(N,L) / A(N,1) 
      I = N + L - 1 
      IF ( I .LE. NUMNP ) THEN
      J = 0 
      DO K = L , MBAND 
      J = J + 1 
      A(I,J) = A(I,J) - C * A(N,K)
      END DO
      A(N,L) = C
      B(I) = B(I) - A(N,L) * B(N)
      ENDIF
      END DO
      B(N) = B(N) / A(N,1)
      END DO
C---------- BACKSUBSTITUTION -------------
      N = NUMNP
      DO WHILE ( N .GT. 0 )
      DO K = 2 , MBAND
      L = N + K - 1 
      IF ( L .LE. NUMNP ) B(N) = B(N)-A(N,K)*B(L) 
      END DO
      N = N - 1 
      END DO
      RETURN
      END 
C
C
      SUBROUTINE FORM (MXN,MXB,MXW,NNODE,NB,NBW,A,RHS,IBC,BV,IBND,
     * SOURCE)
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION RHS(MXN),A(MXN,MXW),IBC(MXB),BV(MXB),IBND(MXB),
     * SOURCE(MXN)
C  BOUNDARY CONDITIONS ARE COUPLED HERE.
C  ITYPE(I) = 1 ---> DIRICHLET, ITYPE(I) = 2 ---> NEUMANN
C----- REFORMATION OF VECTOR RHS DUE TO FIRST KIND B.C. 
C------ CLEAR RIGHT HAND SIDE.
      DO I = 1 , NNODE
      RHS (I) = -SOURCE(I)
      END DO
      DO 50 K = 1 , NB
      IF ( IBC(K) .EQ. 2 ) GO TO 50 
      I = IBND(K) 
      DO 20 J = 2 , NBW 
      I = I - 1 
      IF ( I.LE. 0 ) GO TO 30 
      RHS(I) = RHS(I) - BV(K) * A(I,J)
   20 CONTINUE
   30 I = IBND(K) 
      DO 40 J = 2 , NBW 
      I = I + 1 
      IF ( I .GT. NNODE ) GO TO 50
      RHS(I) = RHS(I) - BV(K) * A(IBND(K),J)
   40 CONTINUE
   50 CONTINUE
C-----REFORMATION OF MATRIX A.
      DO 70 K = 1 , NB
      I = IBND (K)
      IF ( IBC(K) .EQ. 2 ) GO TO 65 
      RHS(I) = BV(K) 
      A(I,1) = 1. 
      DO 60 J = 2 , NBW 
      L = I - J + 1 
      A(I,J) = 0. 
      IF ( L  .LE. 0 ) GO TO 60 
      A( L ,J) = 0. 
   60 CONTINUE
      GO TO 70
   65 RHS (I) = RHS(I) - BV(K)
   70 CONTINUE
      RETURN
      END 
C
C
      SUBROUTINE GRULE ( N , SAI , W )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(N) , W(N)
      IF ( N .LT. 2 ) STOP'N<2'
      IF ( N .GT. 6 ) STOP'N>6'
      GO TO ( 99, 20, 30, 40, 50, 60 ) , N
   99 STOP
   20 SAI(1) = 0.5773502692
      W(1) = 1.
      GO TO 88
   30 SAI(1) = 0.7745966692
      SAI(2) = 0.
      W(1) = 5./ 9.
      W(2) = 8./ 9.
      GO TO 88
   40 SAI(1) = 0.3399810436
      SAI(2) = 0.8611363116
        W(1) = 0.6521451549
        W(2) = 0.3478548451
      GO TO 88
   50 SAI(1) = 0.9061798459
      SAI(2) = 0.5384693101
      SAI(3) = 0.
        W(1) = 0.2369268850
        W(2) = 0.4786286705
        W(3) = 5.12 / 9.
      GO TO 88
   60 SAI(1) = 0.2386191861
      SAI(2) = 0.6612093865
      SAI(3) = 0.9324695142
        W(1) = 0.4679139346
        W(2) = 0.3607615730
        W(3) = 0.1713244924
   88 NN = N / 2
      DO 11 I = 1 , NN
      J = N - I + 1
      SAI(J) = - SAI(I)
      W(J) = W(I)
   11 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE DERIV ( ND, INTEPT, F0, F1, SAI, BPP )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION SAI(INTEPT),BPP(2,ND,INTEPT,INTEPT),F0(ND),F1(ND)
      DO 40 K = 1 , INTEPT
      E1 = SAI (K)
      DO 30 L = 1 , INTEPT
      E2 = SAI (L)
COMPUTATION OF BP(J) = DN(J) / DETA1
      CALL ISOPARA ( ND, E1+0.5 , E2 , F1 )
      CALL ISOPARA ( ND, E1-0.5 , E2 , F0 )
      DO I = 1 , ND
      BPP(1,I,K,L) = F1(I) - F0(I)
      END DO
COMPUTATION OF BP(J)= DN(J)/DETA2
      CALL ISOPARA ( ND, E1 , E2+0.5 , F1 )
      CALL ISOPARA ( ND, E1 , E2-0.5 , F0 )
      DO I = 1 , ND
      BPP(2,I,K,L) = F1(I) - F0(I)
      END DO
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE SHAPEF ( ND, INTEPT, F, SAI, SF )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION F(ND) , SAI(INTEPT) , SF(ND,INTEPT,INTEPT)
      DO 40 K = 1 , INTEPT
      E1 = SAI (K)
      DO 30 L = 1 , INTEPT
      E2 = SAI( L )
      CALL ISOPARA ( ND, E1 , E2 , F )
      DO 20 I = 1 , ND
      SF(I,K,L) = F(I)
   20 CONTINUE
   30 CONTINUE
   40 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE SRCTRM ( MXE,MXN,INTEPT,ND,NNODE,BPP,W,NE,GM,THETA,
     * SF,XCOORD,YCOORD,NODEX,SS, SOURCE )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
     * SOURCE(MXN),
     * SS(ND,ND),BPP(2,ND,INTEPT,INTEPT),W(INTEPT),SF(ND,INTEPT,INTEPT)
C
      TRSCONST = 2.*GM*THETA
      DO I = 1, NNODE
      SOURCE(I) = 0.
      END DO
      DO 500 IEL = 1 ,NE
      DO I = 1 , ND
      DO J = 1 , ND
      SS(I,J) = 0.
      END DO
      END DO
      DO 400 K = 1 , INTEPT
      DO 300 L = 1 , INTEPT
      WEIGHT = W(K) * W(L) * TRSCONST
      YAC11 = 0.
      YAC12 = 0.
      YAC21 = 0.
      YAC22 = 0.
      DO I = 1 , ND
      YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODEX(IEL,I))
      YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODEX(IEL,I))
      END DO
      DETJ = YAC11 * YAC22 - YAC12 * YAC21
      BETA = WEIGHT * DETJ
      DO I = 1 , ND
      DO J = 1 , ND
      SS(I,J) = SS(I,J) + SF(I,K,L) * SF(J,K,L) * BETA
      END DO
      END DO
  300 CONTINUE
  400 CONTINUE
C
      DO K = 1 , ND
      SUM = 0
      DO L = 1 , ND
      SUM = SUM + SS(K,L)
      END DO
      SOURCE(NODEX(IEL,K)) = SOURCE(NODEX(IEL,K)) + SUM
      END DO
  500 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE MOMENT ( MXE,MXN,INTEPT,ND,BPP,W,NE,SF,
     * XCOORD,YCOORD,NODEX, RHS,TORQUE  )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),RHS(MXN),
     * SS(ND,ND),BPP(2,ND,INTEPT,INTEPT),W(INTEPT),SF(ND,INTEPT,INTEPT)
      TORQUE = 0.
      DO 500 IEL = 1 ,NE
      DO 400 K = 1 , INTEPT
      DO 300 L = 1 , INTEPT
      WEIGHT = W(K) * W(L)
      YAC11 = 0.
      YAC12 = 0.
      YAC21 = 0.
      YAC22 = 0.
      DO I = 1 , ND
      YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODEX(IEL,I))
      YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODEX(IEL,I))
      END DO
      DETJ = YAC11 * YAC22 - YAC12 * YAC21
      BETA = WEIGHT * DETJ
      DO I = 1 , ND
      TORQUE  = TORQUE  + SF(I,K,L) * BETA*RHS(NODEX(IEL,I))
      END DO
  300 CONTINUE
  400 CONTINUE
  500 CONTINUE
      TORQUE  = 2.*TORQUE
      RETURN
      END
C
C
      SUBROUTINE ELNDPT ( ND, XS, YS )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION XS(ND), YS(ND)
      IF ( ND .LT. 4 ) STOP'ND<4'
      IF ( ND .GT. 9 ) STOP'ND>9'
      IF ( ND .EQ. 4) THEN
      XS(1) = -1.
      XS(2) = 1.
      XS(3) = 1.
      XS(4) = -1.
      YS(1) = -1.
      YS(2) = -1.
      YS(3) = 1.
      YS(4) = 1.
      END IF
      IF ( ND .EQ. 9) THEN
      XS(1) = -1.
      XS(2) = 0.
      XS(3) = 1.
      XS(4) = 1.
      XS(5) = 1.
      XS(6) = 0.
      XS(7) = -1.
      XS(8) = -1.
      XS(9) = 0.
      YS(1) = -1.
      YS(2) = -1.
      YS(3) = -1.
      YS(4) = 0.
      YS(5) = 1.
      YS(6) = 1.
      YS(7) = 1.
      YS(8) = 0.
      YS(9) = 0.
      END IF
      RETURN
      END
C
C
      SUBROUTINE STRESS ( ND ,MXE,MXN,NE,NNODE, BPPX,NODEX,
     * XCOORD,YCOORD,
     * XS, YS, DUMMY0, DUMMY1, RHS, NDREPEAT,TAUZX,TAUZY )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),RHS(MXN),
     *  BPPX(2,ND,ND),XS(ND),YS(ND),DUMMY0(ND),DUMMY1(ND),
     *  TAUZX(MXN),TAUZY(MXN), NDREPEAT(MXN)
C====================================
      CALL ELNDPT  ( ND, XS, YS )
      CALL DERIVX ( ND, DUMMY0, DUMMY1, XS, YS, BPPX )
C===== INTEGRATION =======
      DO I = 1, NNODE
      NDREPEAT(I) = 0
      TAUZX(I) = 0.
      TAUZY(I) = 0.
      END DO
      DO 500 IEL = 1 ,NE
      DO 400 K = 1 , ND
      YAC11 = 0.
      YAC12 = 0.
      YAC21 = 0.
      YAC22 = 0.
      DO I = 1 , ND
      YAC11 = YAC11 + BPPX(1,I,K) * XCOORD(NODEX(IEL,I))
      YAC12 = YAC12 + BPPX(1,I,K) * YCOORD(NODEX(IEL,I))
      YAC21 = YAC21 + BPPX(2,I,K) * XCOORD(NODEX(IEL,I))
      YAC22 = YAC22 + BPPX(2,I,K) * YCOORD(NODEX(IEL,I))
      END DO
      DETJ = YAC11 * YAC22 - YAC12 * YAC21
      Y11 =   YAC22 / DETJ
      Y12 = - YAC12 / DETJ
      Y21 = - YAC21 / DETJ
      Y22 =   YAC11 / DETJ
      DTDX = 0.
      DTDY = 0.
      DO J = 1 , ND
      DNDX = Y11*BPPX(1,J,K)+Y12*BPPX(2,J,K)
      DNDY = Y21*BPPX(1,J,K)+Y22*BPPX(2,J,K)
      DTDX = DTDX + DNDX*RHS(NODEX(IEL,J))
      DTDY = DTDY + DNDY*RHS(NODEX(IEL,J))
      END DO
      NODE = NODEX(IEL,K)
      TAUZX(NODE) = TAUZX(NODE) + DTDY
      TAUZY(NODE) = TAUZY(NODE) - DTDX
      NDREPEAT(NODE) = NDREPEAT(NODE) + 1
  400 CONTINUE
  500 CONTINUE
      DO I = 1 , NNODE
      TAUZX(I) = TAUZX(I)/NDREPEAT(I)
      TAUZY(I) = TAUZY(I)/NDREPEAT(I)
      END DO
      RETURN
      END
C
C
      SUBROUTINE DERIVX ( ND, DUMMY0, DUMMY1, XS, YS, BPPX )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION XS(ND),YS(ND),BPPX(2,ND,ND),DUMMY0(ND),DUMMY1(ND)
      DO 40 K = 1 , ND
      E1 = XS (K)
      E2 = YS (K)
COMPUTATION OF BP(J) = DN(J) / DETA1
      CALL ISOPARA ( ND, E1+0.5 , E2 , DUMMY1 )
      CALL ISOPARA ( ND, E1-0.5 , E2 , DUMMY0 )
      DO I = 1 , ND
      BPPX(1,I,K) = DUMMY1(I) - DUMMY0(I)
      END DO
COMPUTATION OF BP(J)= DN(J) / DETA2
      CALL ISOPARA ( ND, E1 , E2+0.5 , DUMMY1 )
      CALL ISOPARA ( ND, E1 , E2-0.5 , DUMMY0 )
      DO I = 1 , ND
      BPPX(2,I,K) = DUMMY1(I) - DUMMY0(I)
      END DO
   40 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE DAREA ( MXE,MXN,INTEPT,ND,BPP,W,NE,XCOORD,YCOORD,
     * NODEX, AREA )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
     * BPP(2,ND,INTEPT,INTEPT),W(INTEPT)
      AREA = 0.
      DO IEL = 1 ,NE
      DO K = 1 , INTEPT
      DO L = 1 , INTEPT
      WEIGHT = W(K) * W(L)
      YAC11 = 0.
      YAC12 = 0.
      YAC21 = 0.
      YAC22 = 0.
      DO I = 1 , ND
      YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODEX(IEL,I))
      YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODEX(IEL,I))
      YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODEX(IEL,I))
      END DO
      DETJ = YAC11 * YAC22 - YAC12 * YAC21
      AREA = AREA + WEIGHT * DETJ
      END DO
      END DO
      END DO
      RETURN
      END
