THIS IS A CONSTITUTIVE MODEL FOR PREDICTION OF THE NON LINEAR 
MATERIAL BEHAVIOR OF METAL ANISOTROPIC MATERIALS(PLASTICITY)
USING FREDRIK-ARMESTRONG KINEMATIC CRITERION/ VOYIADJIS KINEMATIC 
HARDENING CRITERION
THE FOLLOWING INDIVIDUALS HELPED IN DEVELOPING THIS PROGRAM:
G.Z VOYIADJIS -P.I.KATAN-I.N BASUROCHOWDHURY
MODIFIED BY RASHID K.ABU AL-RUB 2001

USING RADIAL RETURN ALGORITHM



C
      PROGRAM MAIN
C
      IMPLICIT NONE
      INTEGER MAX_MAT_TYPE,INCREM,NIT,NDIVER,I_OUT,I_IN,MATNUM
       INTEGER STRS_STRN_REL,ICOUNT,IOCNT,PLANE_STRAIN
      INTEGER EVAL_STIFF_OR_EVAL_STRESS,EVAL_STIFF,EVAL_STRESS
      INTEGER ITERATIONS,K1,K2,K3,K4,DIVER_STOP
       INTEGER ISTART,IFINAL,RESTART
      INTEGER IYIEL,IEND,I,K,ITEST,AXISYMMETRIC,PLANE_STRESS
      INTEGER MAT_ELAS,MAT_PLAS,MAT_ELAS_DAM,MAT_PLAS_DAM,J,LAST
      INTEGER MATTYPE,INCREMENTS,OUTPUT_INTR,ELEM_TYPE,P2X
      LOGICAL INITIAL_CORRECTION,IYIELD
      REAL*8 STRESS_IN(3,3),STRAIN,STRESS_INC(6),EDOTEL,STRESS_ITR(6)
      REAL*8 STRESS_VEC(6),STRN,STRS,SDOTV,STRN1(6),STRS1(6),DE,STRESS
      REAL*8 NUX,NUY,NUZ,DLINC,P2Z,ONE,DEPINV(6,6)
      REAL*8 POISS,SYIELD,YOUNG,AD,STRPLA,STRELA
      REAL*8 DEP,DEPM,EX,EY,EZ,P1X,P1Y,P1Z,P2Y
      REAL*8 P6X,P6Y,P6Z,P7X,P7Y,P7Z
       REAL*8 CONV_FAC,ENRG,ENRG1
      REAL*8 CENTER,BETA_CONST,C_CONST,GAMA_CONST,Q_ISOTROPIC
      REAL*8 R_ISOTROPIC,ISOTROPIC_CONST
      CHARACTER*12 INP_FILE,OUT_FILE
      COMMON/INPUT8/INCREMENTS,ITERATIONS
      PARAMETER(MAX_MAT_TYPE=10)
      COMMON/INPUT5/NUX(MAX_MAT_TYPE),NUY(MAX_MAT_TYPE),
          NUZ(MAX_MAT_TYPE),EX(MAX_MAT_TYPE),
          EY(MAX_MAT_TYPE),EZ(MAX_MAT_TYPE),
          P1X(MAX_MAT_TYPE),P1Y(MAX_MAT_TYPE),
          P1Z(MAX_MAT_TYPE),P2X(MAX_MAT_TYPE),
          P2Y(MAX_MAT_TYPE),P2Z(MAX_MAT_TYPE)
      COMMON/PLASTICITY/P6X(MAX_MAT_TYPE),P6Y(MAX_MAT_TYPE),
          P6Z(MAX_MAT_TYPE),P7X(MAX_MAT_TYPE)
          P7Y(MAX_MAT_TYPE),P7Z(MAX_MAT_TYPE)
      COMMON/XXX16/SYIELD
      COMMON/INPUTF/MATYPE(MAX_MAT_TYPE),
      COMMON/INPUTB/CONV_FAC,ENRG1,NDIVER,DIVER,STOP
      COMMON/CONTR1/INCREM,NIT
      COMMON/ELSTR1/STRN(6)
      COMMON/ELSTR2/STRS(6)
       COMMON/ADMAT1/AD(3,3,3,3)
       COMMON/MATER1/DEP(6,6)
      COMMON/IN_IO/I_OUT,I_IN
       COMMON/STRAIN_INCR/DE(6),EDOTEL(3,3),SDOTV(6)
       COMMON/OUT1/STRESS(6),STRAIN(6),STRELA(6),STRPLA(6)
       COMMON/OUT2/CENTER(6),INITIAL_CORRECTION_IYIELD

       COMMON/MAT_CONT/BETA_CONST,C_CONST,GAMA_CONST,Q_ISOTROPIC,
           R_ISOTROPIC,ISOTROPIC_CONST
       PARAMETER (EVAL_STIFF=0,EVAL_STRESS=1)
       PARAMETER (PLANE_STRESS=1,PLAIN_STRAIN=2,AXISYMMETRIC=3)
       PARAMETER
       (MAT_ELAS=1,MAT_PLAS=2,ELAS_DAM=3,MAT_PLAS_DAM=4)
       DATA ONE/1.0D0/
       C
       C====OPEN INPUT AND OUTPUT FILES
       C
       WRITE(*,10)
       10 FORMAT(2X,'PLEASE ENTER THE  INPUT FILE NAME (12-CHARACTER MAX):'
         ,/)
       READ(*,'(12A)')INP_FILE
       WRITE (*,20)
       20 FORMAT(2/,2X,'PLEASE WRITE THE OUTPUTFILE NAME (12-CHARACTER MAX):'
         ,/)
       I_IN=11
       I_OUT=13
       OPEN(UNIT=I_IN,FILE=INP_FILE,STATUS='UNKNOWN')
       OPEN(UNIT=I_OUT,FILE=OUT_FILE,STATUS='UNKNOWN')
       OPEN(UNIT=2,FILE='TEST.OUT',STATUS='UNKNOWN')
       OPEN(UNIT=7,FILE='RESTF.DAT',STATUS='UNKNOWN')
       C
       C====READ THE INPUT
       C
       DO I=1,3
         READ(I_IN,*)(STRESS_IN(I,J),J=1,3)
       END DO
       READ(I_IN,*)MATNUM,MATYPE(MATNUM)
       READ(I_IN,*)YOUNG,POISS
       NUX(MATINUM)=POISS
       EX(MATINUM)=YOUNG
       READ(I_IN,*)INCREMENTS,ITERATIONS
       READ(I_IN,*)CONV_FAC,DIVER_STOP
       READ(I_IN,*)OUTPUT_INTR
       READ(I_IN,*)STRS_STRN_REL,ELEM_TYPE
       READ(I_IN,*)RESTART
       C
       C====PRINTING THE INPUT DATA
       C
       WRITE(I_OUT,30)
       30 FORMAT(1X,'THE STRESS TENSOR:',/)
       DO I=1,3
         WRITE(I_OUT,35)(STRESS_IN(I,J),J=1,3)
         END DO
      35 FORMAT (3(2X,E12.5))
        WRITE(I_OUT,40)YOUNG,POISS
      40 FORMAT(/,2X,'E=',E12.5,5X,'v=',F5.3)
      WRITE(*,*)'INCREMENTS',INCREMENTS,'ITERATIONS',ITERATIONS
      WRITE(*,*)CONV_FAC,DIVER_STOP
      WRITE(*,*)OUTPUT_INTR,STRS,STRN_REL,ELEM_TYPE
      WRITE(*,*)'RESTART='RESTART
      C
      C INCREMENTS=NUMBER OF STRESS INCREMENTS
      C ITERATIONS=NUMBER OF ITERATIONS
      C OUTPUT_INTR=NUMBER OF INCREMENT AT WHICH OUTPUT IS REQUIRED
      C RESTART=START THE RUN FROM THE LAST CONVERGED INCREMENT
      C
      IF(ELEM_TYPE.GT.300)THEN
        IEEND=6
        ELSE
          IEEND=4
          ENDIF
      C
      C INITIALIZATION
      C
        DO K=1, IEND
          STRESS(K)=0.0D0
          STRAIN(K)=0.0D0
          STRELA(K)=0.0D0
          STRPLA(K)=0.0D0
          CENTER(K)=0.0D0
          END DO
         Q_ISOTROPIC=0.0D0
         R_ISOTROPIC=0.0D0
         INITIAL_CORRECTION=FALSE
         IYIELD=FALSE
      C
      C
      C
      C
      C
      IF(RESTART.EQ.1) THEN
        DO K=1,IEND
          READ(7,*)STRESS(K),STRAIN(K),STRELA(K),STRPLA(K),CENTER(K)
          END DO
          READ(7,*)ISTART

          READ(7,*)Q_ISOTROPIC,R_ISOTROPIC,INITIAL_CORRECTION,IYIELD
          REWIND 7
          CALL VECTOR(ELEM_TYPE,STRESS_IN,STRESS_VEC,ONE)
          DLINC=DFLOAT(INCREMENTS)
          DO K=1,IEND
            STRESS_INCR(K)=(STRESS_VEC(K)-STRESS(K))/DLINC
            END DO
            IFINAL=ISTART+INCREMENTS
            ISTART=ISTART+1
            ELSE
              C
              C ICOUNT=ITERATION COUNT FOR THE RUN
      C IOCNT =INCREMENT COUNT FROM THE START OR SINCE THE LAST
      C    OUTPUT. WHEN'IOCNT' IS EQUAL TO 'OUTPUT_INTR' A
      COMPLETE
      C   OUTPUT WILL BE GENERATED.
      C   INCR=INCREMENT NUMBER
      C  NIT=ITERATION NUMBER
      C  ITERATIONS=MAXIMUM NUMBER OF ITERATIONS ALLOWED
      C
        CALL VECTOR(ELEM_TYPE,STRESS_IN,STRESS_VEC,ONE)
        DLINC=DFLOAT(INCREMENTS)
        DO K=1,IEND
          STRESS_INCR(K)=STRESS_VEC(K)/DLINC
          END DO
          ISTART=1
          IFINAL=INCREMENTS
          ENDIF
          ICOUNT=0
          IOCNT=0
          IF(RESTART.EQ.1)THEN
            DO I=1,IEND
               STRS(I)=STRESS(I)
               STRS1(I)=STRESS(I)
               STRN(I)=STRAIN(I)
               STRN1(I)=STRAIN(I)
               END DO
               ELSE
                 DO I=1, IEND
                   STRS(I)=0.0D0
                   STRS1(I)=0.0D0
                   STRN1(I)=0.0D0
                   END DO
                   ENDIF
       C
       C    START OF
       C    INCREMENT LOOP
       C
         DO INCREMENT=ISTART,IFINAL
           IOCNT=IOCNT+1
           write(2,*)'INCREMENT=',increm
       C
       C STRESS_ITR=TOTAL APPLIED STRESS AT THE END OF THE INCREMENT
       C
        DO I=1,IEND
          STRESS_ITR(I)=STRESS_INCR(I)+STRS(I)
          END DO
          write(2,*)'TOTAL STRESS=',stress_itr(1)
       C
       C   START OF
       C  ITERATION LOOP
       C
          DO NIT=1,ITERATIONS
       C
       C  CALCULATION OF THE STRESS INCREMENT
       C
           DO I=1,IEND
             SDOTV(I)=STRESS_ITR(I)-STRS(I)
             ENDDO
              write(2,*)'STRESS INCREMENT=',sdotv(1)
       C
       C     CALCULATION OF THE STRESS INCREMENT
       C
          CALL MATMOD(ELEM_TYPE,MATNUM,STRS_STRN_REL,
          I_OUT,EVAL_STIFF)
          CALL DINV(DEP,IEND,DEPINV)
          DO K1=1,IEND
            DE(K1)=0.0
            DO K2=1,IEND
              DE(K1)=DE(K1)+DEPINV(K1,K2)*SDOTV(K2)
              ENDDO
              DO I=1, IEND
                STRS1(I)=STRS(1)
                STRN1(I)=STRN(I)
                ENDDO      
       C
       C    UPDATING THE STRESS INCREMENT
       C
             CALL MATMOD(ELEM_TYPE,MATNUM,STRS_STRN_REL,
             I_OUT,EVAL_STRESS)
       C
       C   ITERATION CONVERGENCE
       C
          CALL CHECK(STRS1,STRN1,IEND,ITEST,I_OUT)
          IF (ITEST.EQ.1) THEN
            WRITE(*,*)
            GOTO 600
            ELSE IF (ITEST.EQ.1)THEN
              WRITE(*,*)
              GO TO 590
              END IF
              END DO
              WRITE(*,*)
       C
       C    END OF
       C      ITERATION LOOP
       C
           IF(ITERATION.EQ.1)GO TO 600
             WRITE(I_OUT,1003)INCREM,INCREM-1
             PRINT*,'MAXIMUM NUMBER OF ITERATIONS EXCEEDED.'//
             'PROGRAM TERMINATED'
          590 IF INCREM.LE.1) GO TO 800
             WRITE(*,*)'WRITING OUTPUT FOR LOAD INCREMENT#'
             WRITE(*,*)INCREM
             WRITE(I_OUT_1004)INCREM
             CALL OUTPUT(I_OUT,ELEM_TYPE,MATNUM,STRS_STRN_REL)
             GO TO 800
           600 ICOUNT=ICOUNT+NIT
             IF(OUTPUT_INTR.GT.0)THEN
               IF(MOD(IOCNT,OUTPUT_INTR).EQ.0)THEN
                 WRITE(*,*)'WRITING OUTPUT FOR LOAD INCREMENT#'
                 WRITE(*,*)INCREM
                 CALL OUTPUT(I_OUT,ELEM_TYPE,MATNUM,STRS_STRN_REL)
                 END IF
                 END IF
          C
          C SAVING THE RESTART NECESSARY RESULTS
          C
            DO K=1,IEND
              WRITE(7,*)STRESS(K),STRELA(K),STRPLA(K),CENTER(K)
              END DO
              WRITE(7,*)INCREM
              WRITE(7,*)INCREM

              WRITE(7,*)S_ISOTROPIC,R_ISOTROPIC,INITIAL_CORRECTION,IYIELD
              REWIND 7
          C
          END DO
          800 WRITE(I_OUT,1002)ICOUNT
          1002 FORMAT(//1X,'>>>>>>>TOTAL NUMBER OF ITERATIONS FOR THIS 
          RUN IS'
             .,'=',15)
          1003 FORMAT(1/X,'>>>>>>>PROGRAM TERMINATED DUE TO EXCEEDING
          THE '/
          .9X,'ALLOWABLE NUMBER OF ITERATIONS AT LOAD INCREMENT',I4//
          .1X,'>>>>>>>OUTPUTS ARE FOR THE LAST CONVERGED INCREMENT',I4)
          1004 FORMAT(///1X,'>>>>>>>OUTPUTS AT INCREMENT',I4)
          C
          END