hm_read_fail_orthbiquad.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    hm_read_fail_orthbiquad    ../starter/source/materials/fail/orthbiquad/hm_read_fail_orthbiquad.F
!||--- called by ------------------------------------------------------
!||    hm_read_fail               ../starter/source/materials/fail/hm_read_fail.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                     ../starter/source/output/message/message.F
!||    hm_get_float_array_index   ../starter/source/devtools/hm_reader/hm_get_float_array_index.F
!||    hm_get_floatv              ../starter/source/devtools/hm_reader/hm_get_floatv.F
!||    hm_get_floatv_dim          ../starter/source/devtools/hm_reader/hm_get_floatv_dim.F
!||    hm_get_intv                ../starter/source/devtools/hm_reader/hm_get_intv.F
!||    hm_option_is_encrypted     ../starter/source/devtools/hm_reader/hm_option_is_encrypted.f
!||--- uses       -----------------------------------------------------
!||    hm_option_read_mod         ../starter/share/modules1/hm_option_read_mod.f
!||    message_mod                ../starter/share/message_module/message_mod.F
!||    submodel_mod               ../starter/share/modules1/submodel_mod.F
!||====================================================================
      SUBROUTINE hm_read_fail_orthbiquad(
     .           FAIL    ,MAT_ID   ,FAIL_ID  ,IRUPT    ,
     .           TITR    ,LSUBMODEL,UNITAB   )
C-----------------------------------------------
c    ROUTINE DESCRIPTION :
c    Orthotropic strain  failure model
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE fail_param_mod
      USE unitab_mod
      USE message_mod
      USE submodel_mod
      USE hm_option_read_mod
      USE names_and_titles_mod , ONLY : nchartitle
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C----------+---------+---+---+--------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER            ,INTENT(IN) :: FAIL_ID         ! failure model ID
      INTEGER            ,INTENT(IN) :: MAT_ID          ! material law ID
      INTEGER            ,INTENT(IN) :: IRUPT           ! failure model type number
      CHARACTER(LEN=NCHARTITLE) ,INTENT(IN) :: TITR     ! material model title
      TYPE(unit_type_)   ,INTENT(IN) :: UNITAB          ! table of input units
      TYPE(submodel_data),INTENT(IN) :: LSUBMODEL(*)    ! submodel table
      TYPE(fail_param_)  ,INTENT(INOUT) :: FAIL         ! failure model data structure
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: NANGLE,I,J,K,INFO,REG_FUNC,MFLAG,SFLAG,RATE_FUNC,NFUNC,NUPARAM,NUVAR
      INTEGER :: IPIV2(2),IPIV3(3)
      INTEGER ,PARAMETER :: NSIZE = 2
      INTEGER ,DIMENSION(NSIZE) :: IFUNC
      my_real :: pthk,ref_siz,ref_siz_unit,epsd0,cjc,rate_scale,ref_rate_unit,
     .    r1,r2,r4,r5,c5,c5_min,theta_myreal
      my_real, DIMENSION(:), ALLOCATABLE :: c1,c2,c3,c4,inst
      DOUBLE PRECISION A_1(2,2),B_1(2),A_2(3,3),B_2(3),
     .                 triax_1_lin,triax_2_lin,triax_3_lin,
     .                 triax_4_lin,triax_5_lin,triax_1_quad,
     .                 triax_2_quad,triax_3_quad,triax_4_quad,
     .                 triax_5_quad,cos2(10,10),xmin,ymin
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: THETA,THETA_RAD,Q_X11,Q_X12,Q_X13,
     .                                               Q_X21,Q_X22,Q_X23,Q_INST
      DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: X_1,X_2,AMAT,BVEC
      INTEGER, DIMENSION(:), ALLOCATABLE :: IPIV
      LOGICAL :: IS_AVAILABLE,IS_ENCRYPTED
      DATA triax_1_lin, triax_2_lin, triax_3_lin, triax_4_lin,
     .     triax_5_lin 
     .      / -0.33333333, 0.0, 0.33333333, 0.577350269, 0.66666667 /
      DATA triax_1_quad, triax_2_quad, triax_3_quad, 
     .     triax_4_quad, triax_5_quad 
     .      / 0.11111111, 0.0, 0.11111111, 0.33333333, 0.44444444 /
C
      DATA  cos2/
     1 1.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,   
     2 0.   ,1.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     3 -1.  ,0.   ,2.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     4 0.   ,-3.  ,0.   ,4.   ,0.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     5 1.   ,0.   ,-8.  ,0.   ,8.   ,0.   ,0.   ,0.   ,0.   ,0.   ,     
     6 0.   ,5.   ,0.   ,-20. ,0.   ,16.  ,0.   ,0.   ,0.   ,0.   ,     
     7 -1.  ,0.   ,18.  ,0.   ,-48. ,0.   ,32.  ,0.   ,0.   ,0.   ,     
     8 0.   ,-7.  ,0.   ,56.  ,0.   ,-112.,0.   ,64.  ,0.   ,0.   ,     
     9 1.   ,0.   ,-32. ,0.   ,160. ,0.   ,-256.,0.   ,128. ,0.   ,     
     a 0.   ,9.   ,0.   ,-120.,0.   ,432. ,0.   ,-576 ,0.   ,256. /
C=======================================================================
      is_encrypted = .false.
      is_available = .false.
C--------------------------------------------------
C     (IS OPTION CRYPTED)
C--------------------------------------------------
 
      CALL hm_option_is_encrypted(is_encrypted)
 
C======================================================================================
C     EXTRACT DATA
!---------------
! -> Card1
!---------------
! Percentage of integration failure, flags and size regularization
      CALL hm_get_floatv ('Pthk'           ,pthk   ,is_available,lsubmodel,unitab)
      CALL hm_get_intv   ('MAT_MFLAG'      ,mflag  ,is_available,lsubmodel)
      CALL hm_get_intv   ('MAT_SFLAG'      ,sflag  ,is_available,lsubmodel)
      CALL hm_get_intv   ('MAT_refanglemax',nangle ,is_available,lsubmodel)
      ! Error message
      IF (nangle > 10) THEN 
        CALL ancmsg(msgid=2015,msgtype=msgerror,
     .              anmode=aninfo_blind_1,i1=mat_id,c1=titr)      
      ENDIF 
      CALL hm_get_intv   ('fct_IDel' ,reg_func    ,is_available,lsubmodel)
      CALL hm_get_floatv ('EI_ref'   ,ref_siz     ,is_available,lsubmodel,unitab)
! Default values      
      IF (pthk == zero) pthk = one - em06
      pthk = min(pthk, one)
      pthk = max(pthk,-one)
      IF (sflag == 0)    sflag = 2
! Units
      IF ((ref_siz == zero).AND.(reg_func > 0)) THEN
        CALL hm_get_floatv_dim('EI_ref' ,ref_siz_unit,is_available, lsubmodel, unitab)
        ref_siz = one*ref_siz_unit
      ENDIF
!---------------
! -> Card2
!---------------
! Size regularization and percentage of integration failure
      CALL hm_get_floatv ('MAT_C5'    ,c5          ,is_available,lsubmodel,unitab)
      CALL hm_get_floatv ('MAT_EPSD0' ,epsd0       ,is_available,lsubmodel,unitab)
      CALL hm_get_floatv ('MAT_CJC'   ,cjc         ,is_available,lsubmodel,unitab) 
      CALL hm_get_intv   ('fct_IDrate',rate_func   ,is_available,lsubmodel)
      CALL hm_get_floatv ('RATE_scale',rate_scale  ,is_available,lsubmodel,unitab)
! Units and default values
      IF ((rate_scale == zero).AND.(rate_func > 0)) THEN 
        CALL hm_get_floatv_dim('RATE_scale' ,ref_rate_unit  ,is_available, lsubmodel, unitab)
        rate_scale = ref_rate_unit*one
      ENDIF
      IF (rate_func > 0) THEN 
        cjc   = zero
        epsd0 = zero
      ELSE
        rate_scale = zero
      ENDIF
      IF (cjc == zero .OR. epsd0 == zero) THEN 
        cjc   = zero
        epsd0 = zero
      ENDIF
!---------------
! -> Card3
!---------------
! Biquad parameter by angle
      IF (.NOT.ALLOCATED(c1))   ALLOCATE(c1(nangle))
      IF (.NOT.ALLOCATED(c2))   ALLOCATE(c2(nangle))
      IF (.NOT.ALLOCATED(c3))   ALLOCATE(c3(nangle)) 
      IF (.NOT.ALLOCATED(c4))   ALLOCATE(c4(nangle))
      IF (sflag == 3) THEN 
        IF (.NOT.ALLOCATED(inst)) ALLOCATE(inst(nangle))
        inst = zero
      ENDIF
      c5_min = infinity
      ! Material selector
      IF (mflag == 0) THEN 
        ! Loop over angles (must be equally distributed between 0 and pi/2)
        DO j = 1, nangle
          CALL hm_get_float_array_index('MAT_C1',c1(j),j,is_available,lsubmodel,unitab)
          CALL hm_get_float_array_index('MAT_C2',c2(j),j,is_available,lsubmodel,unitab)
          CALL hm_get_float_array_index('MAT_C3',c3(j),j,is_available,lsubmodel,unitab)
          CALL hm_get_float_array_index('MAT_C4',c4(j),j,is_available,lsubmodel,unitab)
          ! Default values
          IF (c3(j) == zero) c3(j) = 0.6d0
          IF (c1(j) == zero .AND. c2(j) == zero .AND. c4(j) == zero .AND. c5 == zero) THEN
            c1(j)  = 3.5d0*c3(j)
            c2(j)  = 1.6d0*c3(j)
            c4(j)  = 0.6d0*c3(j)
            c5_min = min(c5_min,1.5d0*c3(j))
          ENDIF
          ! If necking instability is activated
          IF (sflag == 3) THEN
            CALL hm_get_float_array_index('MAT_INST',inst(j),j,is_available,lsubmodel,unitab)
            IF (inst(j) <= zero)   THEN
              CALL ancmsg(msgid=2016,msgtype=msgwarning,
     .                    anmode=aninfo_blind_1,i1=mat_id,c1=titr)    
              sflag = 2
            ELSEIF (inst(j) >= c4(j)) THEN
              CALL ancmsg(msgid=2017,msgtype=msgwarning,
     .                    anmode=aninfo_blind_1,i1=mat_id,c1=titr)
              sflag = 2
            ENDIF           
          ENDIF
        ENDDO
      ELSE
        ! RX parameters
        IF (mflag == 1) THEN ! Mild Seel (c3 = 0.6)
          r1 = 3.5d0 
          r2 = 1.6d0 
          r4 = 0.6d0 
          r5 = 1.5d0
        ELSEIF (mflag == 2) THEN ! DP600 (c3 = 0.5)
          r1 = 4.3d0 
          r2 = 1.4d0
          r4 = 0.6d0
          r5 = 1.6d0
        ELSEIF (mflag == 3) THEN ! Boron (c3 = 0.12)
          r1 = 5.2d0 
          r2 = 3.1d0 
          r4 = 0.8d0 
          r5 = 3.5d0 
        ELSEIF (mflag == 4) THEN ! Aluminium AA5182 (c3 = 0.3)
          r1 = 5.0d0
          r2 = 1.0d0 
          r4 = 0.4d0
          r5 = 0.8d0
        ELSEIF (mflag == 5) THEN ! Aluminium AA6082-T6 (c3 = 0.17)
          r1 = 7.8d0
          r2 = 3.5d0   
          r4 = 0.6d0
          r5 = 2.8d0
        ELSEIF (mflag == 6) THEN ! Plastic light_eBody PA6GF30 (c3 = 0.1)
          r1 = 3.6d0
          r2 = 0.6d0
          r4 = 0.5d0
          r5 = 0.6d0
        ELSEIF (mflag == 7) THEN ! Plastic light_eBody PP T40 ( c3=0.11 )
          r1 = 10.0d0
          r2 = 2.7d0
          r4 = 0.6d0
          r5 = 0.7d0
        ELSEIF (mflag == 99) THEN ! user scalling factors
          CALL hm_get_floatv ('MAT_R1' ,r1    ,is_available,lsubmodel,unitab)
          CALL hm_get_floatv ('MAT_R2' ,r2    ,is_available,lsubmodel,unitab)
          CALL hm_get_floatv ('MAT_R4' ,r4    ,is_available,lsubmodel,unitab)
          CALL hm_get_floatv ('MAT_R5' ,r5    ,is_available,lsubmodel,unitab)
        ELSE ! ELSE --> Mild Seel
          r1 = 3.5d0
          r2 = 1.6d0
          r4 = 0.6d0
          r5 = 1.5d0
        ENDIF   
        ! Loop over angles (must be equally distributed between 0 and pi/2)
        DO j = 1, nangle
          CALL hm_get_float_array_index('MAT_C3'  ,c3(j)  ,j,is_available,lsubmodel,unitab)
          ! Default values 
          IF (c3(j) == zero) THEN 
            IF (mflag == 1) THEN 
              c3(j) = 0.6d0
            ELSEIF (mflag == 2) THEN 
              c3(j) = 0.5d0
            ELSEIF (mflag == 3) THEN
              c3(j) = 0.12d0
            ELSEIF (mflag == 4) THEN
              c3(j) = 0.3d0
            ELSEIF (mflag == 5) THEN
              c3(j) = 0.17d0
            ELSEIF (mflag == 6) THEN
              c3(j) = 0.1d0
            ELSEIF (mflag == 7) THEN
              c3(j) = 0.11d0
            ENDIF
          ENDIF
          ! Computation of C1,C2,C4,C5
          c1(j)  = r1*c3(j)
          c2(j)  = r2*c3(j)
          c4(j)  = r4*c3(j)
          c5_min = min(c5_min,r5*c3(j))
          ! If necking instability is activated
          IF (sflag == 3) THEN
            CALL hm_get_float_array_index('MAT_INST',inst(j),j,is_available,lsubmodel,unitab)
            IF (inst(j) <= zero)   THEN
              CALL ancmsg(msgid=2016,msgtype=msgwarning,
     .                    anmode=aninfo_blind_1,i1=mat_id,c1=titr)    
              sflag = 2
            ELSEIF (inst(j) >= c4(j)) THEN
              CALL ancmsg(msgid=2017,msgtype=msgwarning,
     .                    anmode=aninfo_blind_1,i1=mat_id,c1=titr)
              sflag = 2
            ENDIF           
          ENDIF
        ENDDO
      ENDIF
      ! Default value
      IF (c5 == zero) c5 = c5_min
C======================================================================================
C    COMPUTING FITTING PARAMETERS
C======================================================================================
      IF (.NOT.ALLOCATED(x_1))  ALLOCATE(x_1(nangle,3))
      IF (.NOT.ALLOCATED(x_2))  ALLOCATE(x_2(nangle,3))
      !================================================
      ! Loop over the angle
      !================================================
      DO j = 1,nangle
C
        ! Coefficient for the first parabole
        ! ---------------------------------------
        a_1(1,1) = triax_1_lin
        a_1(1,2) = triax_1_quad
        a_1(2,1) = triax_3_lin
        a_1(2,2) = triax_3_quad
        b_1(1)   = c1(j) - c2(j)
        b_1(2)   = c3(j) - c2(j)
C
        ! Fitting the first quadratic function
#ifndef WITHOUT_LINALG
        CALL dgesv(2, 1, a_1, 2, ipiv2, b_1, 2, info)  
#else
        WRITE(6,*) "Error: Blas/Lapack required"
#endif
        x_1(j,1)   = c2(j)
        x_1(j,2:3) = b_1(1:2)
C      
        ! Coefficient for the second parabole
        !----------------------------------------
        a_2(1,1) = one
        a_2(1,2) = triax_3_lin
        a_2(1,3) = triax_3_quad
        a_2(2,1) = one
        a_2(2,2) = triax_4_lin
        a_2(2,3) = triax_4_quad
        a_2(3,1) = one
        a_2(3,2) = triax_5_lin
        a_2(3,3) = triax_5_quad
        b_2(1)   = c3(j)
        b_2(2)   = c4(j)
        b_2(3)   = c5
C
        ! Fitting the second quadratic function
#ifndef WITHOUT_LINALG
        CALL dgesv(3, 1, a_2, 3, ipiv3, b_2, 3, info)
#endif
        x_2(j,1:3) = b_2(1:3)
C
      ENDDO
C
C======================================================================================
C    COMPUTING COSINE INTERPOLATION
C======================================================================================
c      
      IF (.NOT.ALLOCATED(theta))     ALLOCATE(theta(nangle))
      IF (.NOT.ALLOCATED(theta_rad)) ALLOCATE(theta_rad(nangle))
c
      ! Computation of angles and check curves
      DO j = 1, nangle
        IF (nangle > 1) THEN 
          theta(j)     = (j-1)*(90.0d0/(nangle-1))
          theta_rad(j) = theta(j)*(pi/180.0d0)
        ELSE
          theta(j)     = zero
          theta_rad(j) = zero
        ENDIF
c
        ! Check if minimum of first parabola is negative 
        xmin = -x_1(j,2)/(two*x_1(j,3))
        ymin = x_1(j,3)*(xmin**2) + x_1(j,2)*xmin + x_1(j,1)
        IF (ymin < zero) THEN 
          theta_myreal = theta(j)
          CALL ancmsg(msgid=3002,
     .                msgtype=msgwarning,
     .                anmode=aninfo_blind_1,
     .                i1=mat_id,
     .                c1=titr,
     .                i2=j,
     .                r1=theta_myreal)          
        ENDIF
c
        ! Check if minimum of second parabola is negative 
        IF (sflag == 1) THEN 
          xmin = -x_2(j,2)/(two*x_2(j,3))
          ymin = x_2(j,3)*(xmin**2) + x_2(j,2)*xmin + x_2(j,1)
          IF (ymin < zero) THEN 
            theta_myreal = theta(j)
            CALL ancmsg(msgid=3003,
     .                  msgtype=msgwarning,
     .                  anmode=aninfo_blind_1,
     .                  i1=mat_id,
     .                  c1=titr,
     .                  i2=j,
     .                  r1=theta_myreal)              
          ENDIF
        ENDIF
c
      ENDDO
c
      ! Allocation of the A-MATRIX and the Pivot vector
      ALLOCATE (amat(nangle,nangle),ipiv(nangle))
c
      ! Filling the A-MATRIX
      DO j = 1,nangle
        DO i = 1,nangle
          amat(j,i) = zero
          DO k = 1,i
            amat(j,i) = amat(j,i) + cos2(k,i)*(cos(two*theta_rad(j)))**(k-1)
          ENDDO
        ENDDO
      ENDDO
c
      ! Allocation of factors
      ALLOCATE(q_x11(nangle),q_x12(nangle),q_x13(nangle),
     .         q_x21(nangle),q_x22(nangle),q_x23(nangle))
c
      ! Initialization of tables
      q_x11(1:nangle) = zero
      q_x12(1:nangle) = zero
      q_x13(1:nangle) = zero
      q_x21(1:nangle) = zero
      q_x22(1:nangle) = zero
      q_x23(1:nangle) = zero
      IF (sflag == 3) THEN 
        ALLOCATE(q_inst(nangle))
        q_inst(1:nangle) = zero
      ENDIF
c
      ! Filling the B vector with experimental points
      IF (sflag == 3) THEN 
        ALLOCATE (bvec(nangle,7))
      ELSE
        ALLOCATE (bvec(nangle,6))
      ENDIF
      bvec(1:nangle,1) = x_1(1:nangle,1)
      bvec(1:nangle,2) = x_1(1:nangle,2)
      bvec(1:nangle,3) = x_1(1:nangle,3)
      bvec(1:nangle,4) = x_2(1:nangle,1)
      bvec(1:nangle,5) = x_2(1:nangle,2)
      bvec(1:nangle,6) = x_2(1:nangle,3)
      IF (sflag == 3) bvec(1:nangle,7) = inst(1:nangle)
c        
      ! Initialization of the Pivot vector
      ipiv(1:nangle) = 0
c
      ! Solving the A-MATRIX * x = B vector system
#ifndef WITHOUT_LINALG
      IF (sflag == 3) THEN 
        CALL dgesv(nangle, 7, amat, nangle, ipiv, bvec, nangle, info)
      ELSE
        CALL dgesv(nangle, 6, amat, nangle, ipiv, bvec, nangle, info)
      ENDIF
#else
      WRITE(6,*) "Error: Blas/Lapack required"
#endif
c        
      ! Recovering the solution
      q_x11(1:nangle)  = bvec(1:nangle,1)
      q_x12(1:nangle)  = bvec(1:nangle,2)
      q_x13(1:nangle)  = bvec(1:nangle,3)
      q_x21(1:nangle)  = bvec(1:nangle,4)
      q_x22(1:nangle)  = bvec(1:nangle,5)
      q_x23(1:nangle)  = bvec(1:nangle,6)
      IF (sflag == 3) q_inst(1:nangle) = bvec(1:nangle,7)
c
c----------------------------------
      ! -> Number of parameters
      nuparam = 7
      IF (sflag == 3) THEN 
        nuparam = nuparam + 7*nangle
      ELSE
        nuparam = nuparam + 6*nangle
      ENDIF
      ! -> Number of functions
      nfunc   = 0
      IF (rate_func /= 0) THEN 
        nfunc = nfunc + 1
        ifunc(nfunc) = rate_func
      ENDIF
      IF (reg_func /= 0) THEN 
        nfunc = nfunc + 1 
        ifunc(nfunc) = reg_func
      ENDIF
      ! -> Number of user variables
      nuvar = 3   
C======================================================================================
c     Filling buffer tables
C======================================================================================
      fail%KEYWORD = 'ORTH-BIQUAD' 
      fail%IRUPT   = irupt 
      fail%FAIL_ID = fail_id 
      fail%NUPARAM = nuparam
      fail%NIPARAM = 0
      fail%NUVAR   = nuvar
      fail%NFUNC   = nfunc
      fail%NTABLE  = 0
      fail%NMOD    = 0
      fail%PTHK    = pthk
c            
      ALLOCATE (fail%UPARAM(fail%NUPARAM))
      ALLOCATE (fail%IPARAM(fail%NIPARAM))
      ALLOCATE (fail%IFUNC (fail%NFUNC))
      ALLOCATE (fail%TABLE (fail%NTABLE))
c
      fail%IFUNC(1:nfunc) = ifunc(1:nfunc)
c
      fail%UPARAM(1) = pthk
      fail%UPARAM(2) = sflag
      fail%UPARAM(3) = ref_siz
      fail%UPARAM(4) = epsd0
      fail%UPARAM(5) = cjc
      fail%UPARAM(6) = rate_scale
      fail%UPARAM(7) = nangle
      IF (sflag == 3) THEN 
        DO j = 1,nangle
          fail%UPARAM(8  + 7*(j-1)) = q_x11(j)
          fail%UPARAM(9  + 7*(j-1)) = q_x12(j)
          fail%UPARAM(10 + 7*(j-1)) = q_x13(j)
          fail%UPARAM(11 + 7*(j-1)) = q_x21(j)
          fail%UPARAM(12 + 7*(j-1)) = q_x22(j)
          fail%UPARAM(13 + 7*(j-1)) = q_x23(j)
          fail%UPARAM(14 + 7*(j-1)) = q_inst(j)
        ENDDO    
      ELSE
        DO j = 1,nangle
          fail%UPARAM(8  + 6*(j-1)) = q_x11(j)
          fail%UPARAM(9  + 6*(j-1)) = q_x12(j)
          fail%UPARAM(10 + 6*(j-1)) = q_x13(j)
          fail%UPARAM(11 + 6*(j-1)) = q_x21(j)
          fail%UPARAM(12 + 6*(j-1)) = q_x22(j)
          fail%UPARAM(13 + 6*(j-1)) = q_x23(j)
        ENDDO          
      ENDIF
c--------------------------
c     Printout data
c-------------------------- 
      IF (is_encrypted) THEN
        WRITE(iout,'(5X,A,//)')'CONFIDENTIAL DATA'
      ELSE
        WRITE(iout,1000)
        IF (mflag /= 0) WRITE(iout, 1100) mflag
        DO j = 1,nangle
          WRITE(iout,1200) theta(j),c1(j),c2(j),c3(j),c4(j),c5,
     &    x_1(j,3),x_1(j,2),x_1(j,1),x_2(j,3),x_2(j,2),x_2(j,1)           
          IF (sflag == 3) WRITE(iout, 1900) inst(j)
        ENDDO
        WRITE(iout,1300) pthk,sflag
        IF (reg_func > 0) WRITE(iout,1400) reg_func,ref_siz
        IF (epsd0 > zero) THEN 
          WRITE(iout,1500) epsd0,cjc
        ELSEIF (rate_func > 0) THEN 
          WRITE(iout,1600) rate_func,rate_scale
        ENDIF
        WRITE(iout,2000)
      ENDIF
c--------------------------
c     Deallocation
c--------------------------  
      IF (ALLOCATED(c1))        DEALLOCATE(c1)
      IF (ALLOCATED(c2))        DEALLOCATE(c2)
      IF (ALLOCATED(c3))        DEALLOCATE(c3)
      IF (ALLOCATED(c4))        DEALLOCATE(c4)
      IF (ALLOCATED(inst))      DEALLOCATE(inst)
      IF (ALLOCATED(x_1))       DEALLOCATE(x_1)
      IF (ALLOCATED(x_2))       DEALLOCATE(x_2)
      IF (ALLOCATED(theta))     DEALLOCATE(theta)
      IF (ALLOCATED(theta_rad)) DEALLOCATE(theta_rad)
      IF (ALLOCATED(q_x11))     DEALLOCATE(q_x11)
      IF (ALLOCATED(q_x12))     DEALLOCATE(q_x12)
      IF (ALLOCATED(q_x13))     DEALLOCATE(q_x13)
      IF (ALLOCATED(q_x21))     DEALLOCATE(q_x21)
      IF (ALLOCATED(q_x22))     DEALLOCATE(q_x22)
      IF (ALLOCATED(q_x23))     DEALLOCATE(q_x23)
      IF (ALLOCATED(q_inst))    DEALLOCATE(q_inst)
      IF (ALLOCATED(amat))      DEALLOCATE(amat)
      IF (ALLOCATED(ipiv))      DEALLOCATE(ipiv)
c-----------------------------------------------------
 1000 FORMAT(
     & 5x,'  ------------------------------------------   ',/
     & 5x,'    FAILURE CRITERION : ORTHOTROPIC BIQUAD     ',/,
     & 5x,'  ------------------------------------------   ',/)
 1100 FORMAT(
     & 5x,'MATERIAL PARAMETER SELECTOR  . . . . . . . . . . .=',i10)
 1200 FORMAT(     
     & 5x,'|| FAILURE STRAINS FOR ANGLE',f5.1,' DEG',/,
     & 5x,'  -------------------------------------------------',/,
     & 5x,'  SIMPLE COMPRESSION C1  . . . . . . . . . . . . .=',1pg20.13/
     & 5x,'  SHEAR C2 . . . . . . . . . . . . . . . . . . . .=',1pg20.13/
     & 5x,'  SIMPLE TENSION C3  . . . . . . . . . . . . . . .=',1pg20.13/
     & 5x,'  PLANE STRAIN C4  . . . . . . . . . . . . . . . .=',1pg20.13/
     & 5x,'  BIAXIAL TENSION C5 . . . . . . . . . . . . . . .=',1pg20.13/
     & 5x,'                                                   ',/
     & 5x,'  LOW STRESS TRIAXIALITY PARABOLA PARAMETER A. . .=',1pg20.13/ 
     & 5x,'  LOW STRESS TRIAXIALITY PARABOLA PARAMETER B. . .=',1pg20.13/ 
     & 5x,'  LOW STRESS TRIAXIALITY PARABOLA PARAMETER C. . .=',1pg20.13/ 
     & 5x,'                                                   ',/
     & 5x,'  HIGH STRESS TRIAXIALITY PARABOLA PARAMETER D . .=',1pg20.13/ 
     & 5x,'  HIGH STRESS TRIAXIALITY PARABOLA PARAMETER E . .=',1pg20.13/ 
     & 5x,'  HIGH STRESS TRIAXIALITY PARABOLA PARAMETER F . .=',1pg20.13/)
 1300 FORMAT(
     & 5x,'ELEMENT DELETION :',/,
     & 5x,'SHELL ELEMENT DELETION PARAMETER PTHICKFAIL. . . .=',1pg20.13,/,
     & 5x,'  > 0.0 : FRACTION OF FAILED THICKNESS             ',/,
     & 5x,'  < 0.0 : FRACTION OF FAILED INTG. POINTS OR LAYERS',/,
     & 5x,'S-FLAG . . . . . . . . . . . . . . . . . . . . . .=',i10,/)
 1400 FORMAT(
     & 5x,'ELEMENT LENGTH REGULARIZATION USED:',/,
     & 5x,'REGULARIZATION FUNCTION ID . . . . . . . . . . . .=',i10,/,
     & 5x,'REFERENCE ELEMENT LENGTH . . . . . . . . . . . . .=',1pg20.13,/)
 1500 FORMAT(
     & 5x,'JOHNSON-COOK STRAIN-RATE DEPENDENCY:',/,
     & 5x,'REFERENCE STRAIN-RATE  . . . . . . . . . . . . . .=',1pg20.13,/,
     & 5x,'JOHNSON-COOK PARAMETER . . . . . . . . . . . . . .=',1pg20.13,/)
 1600 FORMAT(
     & 5x,'TABULATED STRAIN-RATE DEPENDENCY:',/,
     & 5x,'STRAIN-RATE DEPENDENCY FUNCTION ID . . . . . . . .=',i10,/,
     & 5x,'STRAIN-RATE SCALE FACTOR . . . . . . . . . . . . .=',1pg20.13,/)
 1900 FORMAT(
     & 5x,'  INSTABILITY STRAIN . . . . . . . . . . . . . . .=',1pg20.13,//)
 2000 FORMAT(
     & 5x,' --------------------------------------------------',//)
c----------- 
      END