      SUBROUTINE ABQMAIN
C
C     PROGRAM  FMAXIMP
C====================================================================
C
C   INPUT FILE --  `FNAME'.fil, max_round_input.dat 
C
C   OUTPUT FILE -- `IMPFNAME`
C
C====================================================================
C
C  The use of ABA_PARAM.INC eliminates the need to have different 
C  versions of the code for single and double precision.  
C  ABA_PARAM.INC defines an appropriate IMPLICIT REAL statement and 
C  sets the value of NPRECD to 1 or 2, depending on whether the 
C  machine uses single or double precision.
C
C====================================================================
C  ARRAY   = Real array containing values read from results file
C             (.fil). Equivalenced to JRRAY.
C  JRRAY   = Integer array containing values read from results file
C             (.fil). Equivalenced to ARRAY.
C  LRUNIT  = Array containing unit number and format of results file
C              LRUNIT(1,*) --> Unit number of input file.
C              LRUNIT(2,*) --> Format of input file.
C  DISP    = Contains the eigenvector for a particular eigenmode
C  COORD   = Original coordinate data
C  INODE   = Original node label
C  IDOF    = DOF for the element
C  JEIGNO  = Array of mode shapes used for calculating the perturbed
C             mesh
C  FNAME   = Name of the results file
C  NODEMAX = Number of nodes in the model
C  IELMAX  = Number of elements in the model
C====================================================================
C
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION  ARRAY(513), JRRAY(NPRECD,513), LRUNIT(2,1)
      EQUIVALENCE (ARRAY(1), JRRAY(1,1))
C===================================================================
C  ITOTAL must be greater than or equal to the number of nodes in the
C  model; MEIGEN must be greater than or equal to the number of eigenmodes.
C===================================================================
      PARAMETER (ITOTAL = 8000,MEIGEN=20)
C
C====================================================================
C
      DIMENSION DISP(3,ITOTAL,MEIGEN), COORD(3,ITOTAL),FCTR(MEIGEN)
      DIMENSION ALPHA(MEIGEN), NPAIR(MEIGEN)
      DIMENSION OLD_COORD(3,ITOTAL)
      DIMENSION INODE(ITOTAL), IDOF(30), JEIGNO(MEIGEN),EIGVAL(MEIGEN)
      CHARACTER FNAME*80,OUTFILE*80,INPFILE*80
      PARAMETER (INPFILE = 'max_round_input.dat')
C
C====================================================================
C  Define flags and counters.
C  
C====================================================================
      ICYCLE = 0
      I1901  = 0
      I101   = 0
      I      = 1
      K      = 1
C
C====================================================================
C  Define file access variables.
C
C====================================================================
      NRU = 1
      LRUNIT(1,NRU) = 8
      LRUNIT(2,NRU) = 2
      LOUTF = 0
C
C====================================================================
C  Open input and output files.
C
C====================================================================
      OPEN(UNIT=16,FILE=INPFILE,STATUS='OLD')
C
      READ(16,'(A)', IOSTAT = J ) OUTFILE
      IF (J .NE. 0)  GOTO 950
C
      OPEN(UNIT=15,FILE=OUTFILE,STATUS='UNKNOWN',IOSTAT = J)
      IF (J .NE. 0) THEN
        WRITE(*,900) OUTFILE
        GOTO 950
      ENDIF
C
C====================================================================
C  Get the name of the results (.fil) file.
C
C====================================================================
      WRITE(*,2000)
      WRITE(6,*) ' Enter the name of the results file (w/o .fil):'
      READ(16,'(A)', IOSTAT = J ) FNAME
      IF (J .NE. 0) GOTO 950
      READ(16,*) N_CIRCUM
      READ(16,*) N_LENGTH
C
C====================================================================
C  Access ABAQUS libraries to set up input file.
C
C====================================================================
      CALL  INITPF (FNAME, NRU, LRUNIT, LOUTF)
C                                                             
      JUNIT = LRUNIT(1,NRU)
C
      CALL  DBRNU (JUNIT)
C
C====================================================================
C  Read a record from the input file.
C  On the first pass through the file obtain the number of nodes for
C  a diagnostic check.
C====================================================================
      CALL DBFILE (0, ARRAY, JRCD)
      DO WHILE (JRCD .EQ. 0)     
        IF (JRRAY(1,2) .EQ. 1980) IEIGNO = JRRAY(1,3)
        IF (JRRAY(1,2) .EQ. 1921 ) THEN
          NODEMAX = JRRAY(1,8)
          IELMAX  = JRRAY(1,7)
          ICYCLE = ICYCLE +1
        ENDIF
        CALL DBFILE (0, ARRAY, JRCD)
      ENDDO
C
      CALL DBFILE (2, ARRAY, JRCD)
C===================================================================
C  User is given a choice of eigenmodes for calculating the perturbed
C  mesh.
C
C===================================================================
      READ(16,*)ACT_IMP
C
      WRITE(*,2010) NODEMAX, IELMAX
C
 5    READ(16,*,ERR = 950) JEIGNO(I)
      IF (JEIGNO(I) .EQ. 0) GOTO 10
      I=I+1
      GOTO 5
C
   10 CONTINUE
      NMODES = I-1
      IF (NMODES.GT.MEIGEN)THEN
         WRITE (6,*) 'INCREASE MEIGEN TO' ,NMODES
         STOP
      END IF
C
      CALL DBFILE (0, ARRAY, JRCD)
C
      DO WHILE (JRCD .EQ. 0)
C
C====================================================================
C  If this is the first pass through the file and the current record
C  is the nodal coordinate record, then read the original nodal 
C  coordinates and the node numbers. Make sure that the third 
C  coordinate exists before saving it.
C
C====================================================================
      IF (JRRAY(1,2) .EQ. 1901 .AND. ICYCLE .LE. 1) THEN
         I1901 = I1901 + 1
         INODE(I1901)  = JRRAY(1,3)
         COORD(1,I1901) = ARRAY(4)
         COORD(2,I1901) = ARRAY(5)
         COORD(3,I1901) = 0.0D0         
         IF (JRRAY(1,1) .GE. 6) COORD(3,I1901) = ARRAY(6)
C
         DO II = 1, 3
           OLD_COORD(II,I1901) = COORD(II,I1901)
         END DO         
C
C====================================================================
C  If this is the first pass through the file and the current record
C  is the active degree of freedom record, save the active d.o.f.  
C  If the d.o.f. is active in the model, IDOF(XX) equals the 
C  position of d.o.f. XX in the output arrays. If the d.o.f. is not 
C  active, IDOF(XX) is zero for d.o.f. XX (i.e., for planar models 
C  IDOF(1) = 1, IDOF(2) = 2, IDOF(3) = 0, IDOF(4) = 0, IDOF(5) = 0,
C  IDOF(6) = 3, etc.). ITRANS equals the number of translational 
C  d.o.f.'s in the model.
C
C====================================================================
      ELSE IF (JRRAY(1,2) .EQ. 1902) THEN
         DO 15 IXX = 1, JRRAY(1,1)-2
            IDOF(IXX) = JRRAY(1,IXX+2)
   15    CONTINUE
         ITRANS = 3
         IF (IDOF(3) .EQ. 0) ITRANS = 2
C
C====================================================================
C  If the current record is the modal record, save the current 
C  eigenvalue number and eigenvalue.
C
C====================================================================
C
      ELSE IF (JRRAY(1,2) .EQ. 1980) THEN
         IEIGNO = JRRAY(1,3)
         DO J = 1, NMODES
           IF (JEIGNO(J) .EQ. IEIGNO) THEN
              K = J
              EIGVAL(K) = ARRAY(4)
           END IF
         ENDDO          
C
C====================================================================
C  If the current record is the displacement record and the current 
C  eigenvalue was requested, read the displacement data.  The data 
C  will be in the coordinate system specified in the 
C  `*NODE FILE,GLOBAL=' option.  If nodal transformations were 
C  performed and GLOBAL=NO was used, the displacements will be in 
C  the local system.  If nodal transformations were used and 
C  GLOBAL=YES, the results will be in the global system.  In all 
C  other cases the results will be in the global system.  Also, 
C  make sure that degrees of freedom are active in the model before
C  saving them in the appropriate array location.
C  
C====================================================================
C
      ELSE IF (JRRAY(1,2) .EQ. 101 .AND. IEIGNO .EQ. JEIGNO(K)) THEN 
         I101 = I101 + 1
         DISP(1,I101,K) = ARRAY(4)
         DISP(2,I101,K) = ARRAY(5)
         IF (IDOF(3) .NE. 0) DISP(3,I101,K) = ARRAY(IDOF(3)+3)
         IF (INODE(I101) .EQ. 0) INODE(I101) = JRRAY(1,3)
           IF( I101 .EQ. NODEMAX ) THEN
             WRITE(6,16) JEIGNO(K)
  16         FORMAT(/,2X,'Nodal coordinate data being computed for',
     1             ' eigenvalue . . .',I5)
C
C====================================================================
C     FACTOR should be entered as a perturbation factor in terms of a
C     percentage value multiplied by a geometric parameter 
C     (e.g., shell thickness).
C====================================================================
C
             WRITE(6,2020)
             READ(16,*,IOSTAT = J) FCTR(K)
             FACTOR = FCTR(K)
             IF (J .NE. 0) GOTO 950
             ICYCLE = ICYCLE + 1
             I101 = 0
C
C====================================================================
         ENDIF
      ENDIF
C
      CALL DBFILE(0, ARRAY, JRCD)
      ENDDO
C
C====================================================================
C  Next line is added for diagnostics.
C
C====================================================================
      IF (NODEMAX .NE. I1901) THEN
       WRITE(*,2025) NODEMAX,I1901
       GOTO 950
      ENDIF
C
      IF (ICYCLE .LE. 1) THEN
        WRITE(*,30)
 30     FORMAT(2X,'. . . NO EIGENVECTORS WERE FOUND . . .',/
     1      2X,'The input file for the buckling analysis must contain:',
     2    /,2X,'--> *NODE FILE <--',
     3    /,2X,'--> U          <--')
        GOTO 950
      ENDIF   
     
C====================================================================
C  Determine if there are repeated eigenvalues.
C
C====================================================================

      TOL = 1D-6
      NP = 0
      DO I = 1, NMODES
         NPAIR(I) = 0
      ENDDO
C
      DO I=1, NMODES-1
        IF (NPAIR(I) .EQ. 0) THEN
         EIG1 = EIGVAL(I)
         DO J = I+1, NMODES
            EIG2 = EIGVAL(J)
            DIFF = ABS ((EIG1 - EIG2)/EIG1)
            IF (DIFF.GT.TOL) GO TO 100
            NP = NP+1
            NPAIR(I) =  NP
            NPAIR(J) = -NP
            GO TO 101           
 100     END DO
        ENDIF
 101  END DO

C====================================================================
C  Determine the scale factors for the eigenvectors.
C
C  Each of the eigenvectors corresponding to the repeated 
C  eigenvalues will be scaled in such a way that the maximum radial 
C  imperfection is introduced along the X-axis when both the eigenvectors 
C  are used to introduce the imperfection.
C
C  The maximum radial imperfection for the unique eigenvectors and the 
C  linear combination of the repeated eigenvectors is scaled to one.
C====================================================================
C      
      DO I = 1, NMODES
C
            IF (NPAIR(I).GE.0) THEN
C
               RAD_MAX = 0.D0
               ALPHA(I)= 1.D0
               TAN_MAX = 0.D0
C
C===================================================================
C  Find the node with the maximum radial displacement.
C===================================================================
C
               DO K = 1, NODEMAX
C===================================================================
C  Determine the local radial and tangential directions of node K.
C===================================================================
C
                  RADIUS = SQRT((OLD_COORD(1,K)**2)+(OLD_COORD(2,K)**2))
C
                  UNIT_R1 = OLD_COORD(1,K)/RADIUS
                  UNIT_R2 = OLD_COORD(2,K)/RADIUS
C
                  UNIT_T1 = -UNIT_R2
                  UNIT_T2 = UNIT_R1
C
C==================================================================
C  Project the displacement vector onto the local directions to 
C  obtain the radial and tangential displacements.
C==================================================================
C
                  RAD_DISP=DISP(1,K,I)*(UNIT_R1)+
     *                     DISP(2,K,I)*(UNIT_R2)
C
                  TAN_DISP=DISP(1,K,I)*(UNIT_T1)+
     *                     DISP(2,K,I)*(UNIT_T2)
C       
                  IF (ABS(RAD_DISP).GT.ABS(RAD_MAX))THEN
                     RAD_MAX=RAD_DISP
                     TAN_MAX=TAN_DISP
                     N_RMAX= K
                  END IF
C
               END DO
C
C==================================================================
C  Determine the scale factors.
C
C  For unique eigenvalues, scale so that max UR=1.0.
C==================================================================
C
               IF (NPAIR(I).EQ.0) THEN
                  IF (ABS(RAD_MAX).GT.0) THEN
                     ALPHA(I)=ALPHA(I)/RAD_MAX
                  END IF
C==================================================================
C  For repeated eigenvalues...
C==================================================================
               ELSE
C
C==================================================================
C  Calculate the corresponding node on the x-axis (NodeX) and save U.
C==================================================================
C
                  MAXROW = (N_RMAX/N_CIRCUM)
C
                  IF (MOD(N_RMAX,N_CIRCUM) .EQ. 0) THEN
                     NODEX = (N_CIRCUM*(MAXROW-1))+1
                  ELSE
                     NODEX = N_CIRCUM*MAXROW+1
                  ENDIF
C
                  X1 = DISP(1,NODEX,I)
                  Y1 = DISP(2,NODEX,I)
C
C==================================================================
C  Save U for this node in the corresponding eigenvector 
C  for the repeated eigenvalue.
C==================================================================
C
                  DO J=I+1, NMODES
                     IF (NPAIR(J) .EQ. -1*NPAIR(I))THEN 
                        X2 = DISP(1, NODEX, J)
                        Y2 = DISP(2, NODEX, J)
                        ALPHA(J)=1.D0
                        GO TO 150
                     ENDIF
                  END DO
 150              CONTINUE

C==================================================================
C  Check the determinant (>0) and calculate the scale factors.
C==================================================================

                  DETERM = (X1*Y2-Y1*X2)
                  IF (ABS(DETERM).GT.0.D0) THEN
                     ALPHA(I) = (Y2*RAD_MAX-X2*TAN_MAX)/DETERM
                     ALPHA(J) = (X1*TAN_MAX-Y1*RAD_MAX)/DETERM
                  ELSE
                     WRITE(6,*) 'DETERMINANT IS ZERO; WILL NOT SCALE'
                  END IF
                  IF (ABS(RAD_MAX).GT.0.D0) THEN
                     ALPHA(I)=ALPHA(I)/RAD_MAX
                     ALPHA(J)=ALPHA(J)/RAD_MAX
                  END IF
               END IF 
            ENDIF
      END DO

C====================================================================
C  Compute the perturbed mesh.  This section assumes that nodal 
C  displacements are in the GLOBAL coordinate system. If they are
C  not, the correct transformations should be applied prior to 
C  perturbing the mesh.  The user should supply this coding. Also, 
C  only the translational degrees of freedom should be used for the 
C  perturbation (X = Xo + u).
C====================================================================
C
      DO I = 1, NMODES
         FACTOR = FCTR(I)*ALPHA(I)
         CALL NODEGEN(COORD,DISP(1,1,I),I1901,FACTOR)
      END DO 
C
C====================================================================
C    Compute the max out-of-round imperfection and close the file.
C====================================================================
C   
      DMAX_IMP=0.D0
C
      DO K = 1, NODEMAX
C
         R_OLD = SQRT(OLD_COORD(1,K)**2 + OLD_COORD(2,K)**2)
         R_NEW = SQRT(    COORD(1,K)**2 +     COORD(2,K)**2)
         TEMP=ABS(R_NEW - R_OLD)
         IF (TEMP.GT.DMAX_IMP) THEN
           DMAX_IMP = TEMP
           NN=INODE(K)
         END IF
      ENDDO
C
C===================================================================
C     Final result.
C===================================================================
C
      WRITE(*,*) 'MAX IMPERFECTION = ',DMAX_IMP,' AT NODE ',NN
      DO K = 1, NMODES
        FACTOR = FCTR(K)*ALPHA(K)
        WRITE(15,120)JEIGNO(K), FACTOR*(ACT_IMP/DMAX_IMP)
      END DO
      
 120  FORMAT(I6,',',E16.9)
C
      CLOSE (15)
      CLOSE (16)
      STOP ' . . . PROGRAM FINISHED SUCCESSFULLY . . . '
C
C====================================================================
C
 900  FORMAT(//,
     1        /,2X,'TROUBLE OPENING FILE',1X,A)
 950  WRITE(*,1000)
 1000 FORMAT(//,
     1        /,2X,' . . . TROUBLE READING DATA . . .   ',
     2        /,2X,' . . .   PROGRAM STOPPED    . . .   ',/)
 2000 FORMAT(//,
     1 /,'   +----------------------------------------------+',
     2 /,'   |                                              |',
     3 /,'   |    P R O G R A M ---  F M A X I M P          |',
     4 /,'   |                                              |',
     5 /,'   +----------------------------------------------+',//)
 2010 FORMAT(//,
     1        /,2X,'Number of nodes in model       . . . . . . . ',I5,
     2        /,2X,'Number of elements in model    . . . . . . . ',I5)
 2015 FORMAT(//,
     1        /,2X,'Number of mode shapes available . . . . . . .',I5,
     2       //,2X,'Enter the mode shape(s) to be used in calculating',
     3        /,2x,'the perturbed mesh (zero when finished):')
 2020 FORMAT(//,
     1       /,2X,'Enter the imperfection factor to be introduced ',
     2       /,2X,'into the geometry for this eigenmode:          ')
 2025 FORMAT(//,
     1        /,2X,'. . . TROUBLE READING COORDINATE DATA . . .  ',
     2        /,2X,'Number of coordinates in model . . . . . .',I5,
     3        /,2X,'Number of coordinates read . . . . . . . .',I5)                
      STOP
      END
C====================================================================
C
      SUBROUTINE NODEGEN(COORD,DISP,I1901,FACTOR)
C
C====================================================================
C PURPOSE:  Defines new coordinate data based upon a fraction of the 
C           eigenvector obtained in a buckling analysis.
C               
C INPUT:
C
C    COORD   = Current coordinate data
C    DISP    = Displacement data (eigenvector)
C    I1901   = Total number of nodes
C    FACTOR  = Imperfection factor (e.g., percentage of shell 
C              thickness)
C====================================================================
C
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION COORD(3,*),DISP(3,*)
C
      DO I = 1, I1901
        COORD(1,I) = COORD(1,I) + FACTOR*DISP(1,I)
        COORD(2,I) = COORD(2,I) + FACTOR*DISP(2,I)
        COORD(3,I) = COORD(3,I) + FACTOR*DISP(3,I)
      ENDDO
      RETURN
      END

