SUBROUTINE MPC(UE,A,JDOF,MDOF,N,JTYPE,X,U,UINIT,MAXDOF,LMPC,
     1               KSTEP,KINC,TIME,NT,NF,TEMP,FIELD,LTRAN,TRAN)
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION UE(MDOF), A(MDOF,MDOF,N), JDOF(MDOF,N), X(6,N),
     1          U(MAXDOF,N), UINIT(MAXDOF,N), TIME(2), TEMP(NT,N),
     2          FIELD(NF,NT,N)
      PARAMETER( NTRIAL = 12, TOLU = 1.D-10, TOLF = TOLU )
C
      IF ( JTYPE .EQ. 1 ) THEN
C
C        INITIAL BEAM AXIS DIRECTORS ==> MUST BE RESET FOR DIFFERENT
C        INITIAL SECTION ORIENTATIONS.
C        FOR 4-2-1-2
C
         COSFI0 = 1.0
         SINFI0 = 0.0
         COSFIB = COS(U(6,3))
         SINFIB = SIN(U(6,3))
C
C        DEPENDENT SHELL DEGREES OF FREEDOM
C
         JDOF(1,1) = 1
         JDOF(2,1) = 5
         JDOF(3,1) = 6
         A(1,1,1)  = COSFI0*COSFIB - SINFI0*SINFIB
         A(2,2,1)  = A(1,1,1)
         A(3,3,1)  = 1.0
C
C        INDEPENDENT SHELL DEGREES OF FREEDOM
C
         JDOF(1,2) = 2
         JDOF(2,2) = 4
         A(1,1,2)  =  SINFI0*COSFIB + COSFI0*SINFIB
         A(2,2,2)  = -A(1,1,2)
C
C        INDEPENDENT BEAM DEGREES OF FREEDOM
C
         JDOF(1,3) = 1
         JDOF(2,3) = 2
         JDOF(3,3) = 6
         A(1,1,3)  = -A(1,1,1)
         A(1,2,3)  = -A(1,1,2)
         A(1,3,3)  = ( X(1,1) + U(1,1) - X(1,3) - U(1,3) )*A(2,2,2) + 
     1               ( X(2,2) + U(2,2) - X(2,3) - U(2,3) )*A(1,1,1)
         A(3,3,3)  = -1.0
C
C        RECOVERY OF DEPENDENT SHELL DEGREES OF FREEDOM
C
         UE(1) =  X(1,3) + U(1,3) - X(1,1) -
     1           (X(2,2) + U(2,2) - X(2,3) - U(2,3))*A(1,1,2)/A(1,1,1)
C
C        USE NEWTON LOOP TO SOLVE FOR ROTATION VARIABLES
C
         UE(2) = U(5,1)
         UE(3) = U(6,1)
         DO 100 K=1, NTRIAL
C
C          SUPPLY TRIGONOMETRIC FUNCTIONS
C
           PHI = SQRT( U(4,1)*U(4,1) + UE(2)*UE(2) + UE(3)*UE(3) )
           IF ( ABS(PHI) .GT. 0.01 ) THEN
              FPHI = ( PHI - SIN(PHI) ) / PHI**3
              GPHI = ( 1.0 - COS(PHI) ) / PHI**2
              HPHI = SIN(PHI) / PHI
           ELSE
              FPHI = 1.0/6.0 - (1.0/120.0)*PHI**2
              GPHI = 0.5 - (1.0/24.0)*PHI**2
              HPHI = 1.0 - (1.0/6.0)*PHI**2
           END IF
C
C          SUPPLY MATRIX COEFFICIENTS
C
           TMPO = U(4,1)*COSFI0 + UE(2)*SINFI0
           ABX  = COS(PHI)*COSFI0 - HPHI*UE(3)*SINFI0 +GPHI*U(4,1)*TMPO 
           ABY  = COS(PHI)*SINFI0 + HPHI*UE(3)*COSFI0 + GPHI*UE(2)*TMPO
           ABZ  = HPHI*(U(4,1)*SINFI0 - UE(2)*COSFI0) + GPHI*UE(3)*TMPO
           FY   = ABY - A(1,1,2)
           FZ   = ABZ
           DOT  = ABX*U(4,1) + ABY*UE(2) + ABZ*UE(3)
           TMPO = ABZ*U(4,1) - ABX*UE(3)
           BYY  = FPHI*UE(2)*TMPO + GPHI*(ABY*UE(2) - DOT)
           BYZ  = FPHI*UE(3)*TMPO - HPHI*ABX + GPHI*ABZ*UE(2)
           TMPO = ABX*UE(2) - ABY*U(4,1)
           BZY  = FPHI*UE(2)*TMPO + HPHI*ABX + GPHI*ABY*UE(3)
           BZZ  = FPHI*UE(3)*TMPO + GPHI*(ABZ*UE(3) - DOT)
C
C          SUMMED FUNCTION VALUES
C
           ERRF = ABS(FY) + ABS(FZ)
           IF ( ERRF .LE. TOLF ) RETURN
C
C          SOLVE LINEAR SYSTEM FOR CORRECTIONS (RETURNED IN FY AND FZ)
C
           DET  = BYY*BZZ - BYZ*BZY
           TMPO = ( BZZ*FY - BYZ*FZ ) / DET
           FZ   = ( BYY*FZ - BZY*FY ) / DET
           FY   = TMPO
C
C          UPDATE SOLUTION
C
           ERRU = ABS(FY) + ABS(FZ)
           UE(2) = UE(2) + FY
           UE(3) = UE(3) + FZ
           IF ( ERRU .LE. TOLU ) RETURN
 100     CONTINUE
C
      ELSE IF ( JTYPE .EQ. 2 ) THEN
C
C        COMPUTE NUMBER OF SHELL NODES AND NUMBER OF SHELLS
C
         NSHNOD = N - 1
         NSHELL = ( NSHNOD - 1 )
         UE(1)  = 0.
C
C        CONTRIBUTION FROM SHELL NODES
C
         DO 10 I=1, NSHNOD
            IF ( I .EQ. 1 .OR. I .EQ. NSHNOD ) THEN
               WEIGHT = 1.0
            ELSE
               WEIGHT = 2.0
            END IF
            JDOF(1,I) = 2
            A(1,1,I)  = WEIGHT
            IF ( I .NE. 1 ) UE(1) = UE(1) - WEIGHT*U(2,I)
 10      CONTINUE
C
C        CONTRIBUTION FROM BEAM NODE
C
         JDOF(1,N) = 2
         A(1,1,N)  = -2.0*NSHELL
         UE(1) = UE(1) - A(1,1,N)*U(2,N)
      END IF
C
      RETURN
      END