SUBROUTINE ABQMAIN
C====================================================================
C This program must be compiled and linked with the command:
C     abaqus make job=felbow
C Run the program using the command:
C     abaqus felbow
C==================================================================
C
C  PURPOSE: 
C 
C    This program reads results stored in an ABAQUS results (.fil) 
C    file and creates an input file for use in the Visualization 
C    module of ABAQUS/CAE (ABAQUS/Viewer) that permits
C    postprocessing and visualization of elbow element results.
C    
C
C  Input file names: `fname.fil', where fname is the root file name
C                                 of the input file.
C
C  Output file name:  output.dat
C
C==================================================================
C
C  VARIABLES USED BY THIS PROGRAM:
C
C   ndi    -- Number of direct components in stress/strain tensor.
C   nshr   -- Number of shear components in stress/strain tensor.
C   ndip1  -- ndi + 1
C   array  -- Real array containing values read from results file.
C               Equivalenced to jrray.
C   jrray  -- Integer array containing values read from results file.
C               Equivalenced to array.
C   fname  -- Root file name of input file (w/o .fil extension).
C   nru    -- Number of results files (.fil) to be read.
C   lrunit -- Array containing unit number and format of results files:
C               lrunit(1,*) --> Unit number of input file.
C               lrunit(2,*) --> Format of input file.
C   loutf  -- Format of output file:
C               0 --> Standard ASCII format.
C               1 --> ABAQUS results file ASCII format.
C               2 --> ABAQUS results file binary format.
C   junit  -- Unit number of file to be opened.
C   jrcd   -- Error check return code.
C               .EQ. 0 --> No errors.
C               .NE. 0 --> Errors detected.
C   key    -- Current record key identifier.
C   jelnum -- Current element number.
C   intpn  -- Integration point number.
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 
C  and sets the value of NPRECD to 1 or 2, depending on whether
C  the machine uses single or double precision. ABA_PARAM.INC is  
C  referenced from the \site (for NT) or /site (for Unix) 
C  ABAQUS release subdirectory when ABAQUS/Make is executed.
C
C==================================================================
C
      include 'aba_param.inc'
      parameter (melem=10000,mnode=10000,mkey=3000,mintpt=360)

      dimension array(513),jrray(nprecd,513),lrunit(2,1)
      equivalence (array(1),jrray(1,1))
C
C==================================================================
C
      character  fname*80
C
C SOME DIMENSION STATEMENTS MAY BE NEEDED IF YOU ARE DOING ADDITIONAL
C CALCULATIONS ON THE DATA
C
      dimension coords(3,mintpt)
      dimension iconn(3,melem),loc_el(melem)
      dimension naxipt(melem),iglb_nod(mnode),loc_nod(mnode)
      dimension islct(mkey),nintpt(melem),nodel(melem),iglb_el(melem)
      dimension xyz(3,mnode),xpos(mnode),tang1(3),tang2(3)
      dimension rot(3,3),xc(3),xpr(3,mintpt),scord(3,mintpt),xnorm(3)
      dimension isctchk(melem),icrdchk(melem),ivarchk(melem)

C The following line places array 'var' on the stack. This array is very
C large (~40MB in size). This will cause problems on Windows, where the size 
C of the stack is fixed and cannot be easily grown.
C
C      dimension var(513,melem)
C
C It is advisable to place large data structures like this into 
C dynamically-allocated memory, like so:

      ALLOCATABLE ::var(:,:)
      ALLOCATE(VAR(513,melem))

C
C==================================================================
C
C  GET THE NAME OF THE RESULTS FILE.
C
C==================================================================
C
      isect1 = 0
      intpt1 = 0
      icomp  = 1
C
      write(6,'(A$)') 'Enter the name of the input file (w/o .fil): '
      read(5,'(A)') fname
C
 100  write(6,*) 'Enter the postprocessing option:'
      write(6,*) ' 1 - variation along the riser '
      write(6,*) ' 2 - variation around the circumference'
      write(6,*) '     of the elbow'
      write(6,*) ' 3 - ovalization of elbow cross-section'
      read(5,*) ipost

      if (ipost .ne. 1 .and. ipost .ne. 2 .and.
     &    ipost .ne. 3) then
         write(6,*)'Invalid entry - try again'
         go to 100
      endif
C
      write(6,*)'Enter the first element number'
      read(5,*) jel_1
C
      if (ipost .ne. 1) then
        jel_2 = jel_1
      else
        write(6,*)'Enter the last element number'
        read(5,*) jel_2
      endif
      if (jel_2 .lt. jel_1) then
         write(6,*)'last elem number must be greater than first'
         return
      endif
      jdiff = jel_2 - jel_1 + 1
      if (jdiff .gt. melem) then
         write(6,*)'Too many elements - max is',melem
         return
      endif
C
      if (ipost .lt. 3) then
         write(6,*)'Enter the section point number'
         read(5,*) isect1
      endif
C
      if (ipost .eq. 1) then
        write(6,*)'Enter the integration point number'
        read(5,*) intpt1
      endif
C
      if (ipost .lt. 3) then

 998  write(6,999)
 999  format(' Enter the key for the variable you wish to process:')
C
      write(6,500)
 500  format(2x,'TEMP ......  2',/,
     &       2x,'COORD......  8',/,
     &       2x,'S.......... 11',/,
     &       2x,'SINV....... 12',/,
     &       2x,'SF......... 13',/,
     &       2x,'E.......... 21',/,
     &       2x,'PE......... 22',/,
     &       2x,'CE......... 23',/,
     &       2x,'IE......... 24',/,
     &       2x,'EE......... 25',/,
     &       2x,'THE........ 88',/,
     &       2x,'LE......... 89',/,
     &       2x,'NE......... 90',/,
     &       2x,'SP.........401',/,
     &       2x,'EP.........403',/,
     &       2x,'NEP........404',/,
     &       2x,'LEP........405',/,
     &       2x,'EEP........408',/,
     &       2x,'IEP........409',/,
     &       2x,'THEP.......410',/,
     &       2x,'PEP........411',/,
     &       2x,'CEP........412')
C
      read(5,*) key
C
      if (key .ne. 2 .and. key .ne. 8 .and.
     &    key .ne. 11 .and. key .ne. 12 .and.
     &    key .ne. 13 .and. key .ne. 21 .and. 
     &    key .ne. 22 .and. key .ne. 23 .and.
     &    key .ne. 24 .and. key .ne. 25 .and.
     &    key .ne. 88 .and. key .ne. 89 .and.
     &    key .ne. 90 .and. key .ne. 401 .and.
     &    key .ne. 403 .and. key .ne. 404 .and.
     &    key .ne. 405 .and. key .ne. 408 .and.
     &    key .ne. 409 .and. key .ne. 410 .and.
     &    key .ne. 411 .and. key .ne. 412) then
C
        write(6,*)'Invalid key - try again'
        go to 998
C
      endif
C
C     Integration point coordinates available only at middle surface
C     Overwrite the section point and set to the default
C
      if (key .eq. 8) isect1 = 0
C
      if (key .ge. 8) then
        write(6,*)'Enter the attribute number'
        read(5,*) icomp
      endif
C
      if (key .eq. 13) then
        write(6,*)'Section force data written only once per element'
        if (ipost .eq. 1) then
          write(6,*)'Will use default integration & section point'
          write(6,*)' '
          intpt1 = 1
          isect1 = 0
        elseif (ipost .gt. 1) then
          write(6,*)'Only use with option 1: Processing terminated'
          return
        endif
      endif
C          
      endif
C
      write(6,*)' Enter the step number'
      read(5,*) istep1
      write(6,*)' Enter the increment number'
      read(5,*) iinc1
C
C SELECT THE RECORDS TO BE PROCESSED
C
C IF THE ISLCT ARRAY IS SET TO 1 THE DATA WILL BE PROCESSED
C IF THE ISLCT ARRAY IS SET TO 0 THE DATA WILL NOT BE PROCESSED
C
      do ii = 1, mkey
        islct(ii) = 0
      end do
C
C ELEMENT DEFINITION
      islct(1900)=0
C NODE DEFINITION
      islct(1901)=0
C INCREMENT START
      islct(2000)=0
C ELEMENT HEADER
      islct(1)=0
C COORD
      islct(8)=0
C     
      if (ipost .lt. 3) islct(key) = 1
      keyprv=0
C
      nels = 0
      numnp = 0
C
      itimchk = 0
      ielchk  = 0
      do i = 1, melem
         loc_el(i) = 0
         isctchk(i) = 0
         icrdchk(i) = 0
         ivarchk(i) = 0
         do j = 1, 513
            var(j,i) = 0.0
         end do
      end do
C
C==================================================================
C
C  OPEN THE OUTPUT FILE.
C
C==================================================================
      open(unit=9,file='output.dat',status='unknown')
C
      nru = 1
      loutf = 0
      lrunit(1,1) = 8
      lrunit(2,1) = 2
C
      call initpf(fname,nru,lrunit,loutf)
C
      junit = 8
C
      call dbrnu(junit)
C
C REWIND FILE
C
      call dbfile(2,array,jrcd)
C
C
C==================================================================
C
C  Read records from the results (.fil) file and process the data.
C  Cover a maximum of 10 million records in the file. 
C
C==================================================================
C
      do 1000 k100 = 1, 100
      do 1000 k1 = 1, 99999
C
C READ RECORD FROM .FIL FILE AND PROCESS DATA.
C
      call dbfile(0, array, jrcd)

         if (jrcd .ne. 0 .and. keyprv .eq. 2001) then
           write(6,*)'End of file'
           close(junit)
           go to 1001
         else if (jrcd .ne. 0) then
           write(6,*)'Error reading input file'
           return
         endif
C
C LENGTH AND KEY RECORD NUMBER ARE ALWAYS NEEDED
C
         lenf= jrray(1,1)
         key = jrray(1,2)
C
C SPECIFIC RECORD NUMBERS FOLLOW...
C
C==================================================================
C
C ELEMENT DEFINITIONS
C
C==================================================================
        if ( key .eq. 1900 ) then
            jelnum=jrray(1,3)
            if (jelnum .gt. melem) then
               write(6,*)'Maximum global elem number too large'
               write(6,*)'Increase melem'
               return
            endif
            ctype=  array(4)
            if (jelnum .ge. jel_1 .and. 
     &          jelnum .le. jel_2) then
              ielchk = 1
              nels = nels + 1
              jel = nels
              iglb_el(jel) = jelnum
              loc_el(jelnum) = jel
              nodel(jel) = lenf - 4
              naxipt(jel) = 1
              if (ctype .ne. 'ELBOW31' .and.
     &            ctype .ne. 'ELBOW31B' .and.
     &            ctype .ne. 'ELBOW31C' .and.
     &            ctype .ne. 'ELBOW32') then
                 write(6,*)'Invalid element type encountered'
                 return
              endif
              if (ctype .eq. 'ELBOW32')naxipt(jel) = 2 
              do kk=5,lenf
                ii=kk-4
                iconn(ii,jel)=jrray(1,kk)
              end do
            endif
            keyprv = key
            goto 9099
C
C==================================================================
C
C NODE DEFINITIONS
C
C==================================================================
        else  if ( key .eq. 1901 ) then
            
            jnod = jrray(1,3)
            if (jnod .gt. mnode) then
               write(6,*)'Maximum global node number exceeded'
               write(6,*)'Increase mnode'
               return
            endif
            numnp = numnp + 1
            iglb_nod(numnp) = jnod
            loc_nod(jnod) = numnp
            xyz(1,numnp) =  array(4)
            xyz(2,numnp) =  array(5)
            xyz(3,numnp) =  array(6)
            keyprv = key
            goto 9099
C
C==================================================================
C
C CURRENT STEP AND INCREMENT NUMBER.
C
C==================================================================
        else if ( key .eq. 2000 ) then
            ttime=array(3)
            stime=array(4)
            sfreq=array(12)
            istep=jrray(1,8)
            iinc =jrray(1,9)
C
C           SET THE FLAG ITIME IF THIS IS THE REQUESTED STEP/INC
C
            if (istep .eq. istep1 .and.
     &         iinc  .eq. iinc1) then
              itime = 1
              itimchk = 1
            else
              itime = 0
            endif
c
            keyprv = key
            goto 9099
C
C==================================================================
C
C ELEMENT VARIABLES AT INTEGRATION POINTS WITHIN THE ELEMENT
C
C==================================================================
C
C ELEMENT DATA
C
         else if ( key .eq. 1 ) then
 
C
C    PROCESS THIS RECORD IF THIS IS THE CORRECT STEP/TIME
C
            if (itime .eq. 1) then
              jelnum = jrray(1,3)
              jel    = loc_el(jelnum)
              ielem  = 0
C
C    PROCESS THIS ELEMENT IF IT IS IN THE LIST OF DESIRED ELEMS
C
              if (jel .gt. 0) then
                ielem  =  1
                intpn  =  jrray(1,4)
                if (intpn .gt. mintpt) then
                   write(6,*)'Max number of int points exceeded'
                   write(6,*)'Increase mintpt'
                   return
                endif
                isect  =  jrray(1,5)
                ilocn  =  jrray(1,6)
                if (ilocn .ne. 0) then
                   write(6,*)'Element data must be at the int points'
                   return
                endif
                crbar  =  array(7)
                ndi    =  jrray(1,8)
                nshr   =  jrray(1,9)
                ndir   =  jrray(1,10)
                nsfc   =  jrray(1,11)
                if (isect .eq. isect1)then
                  isctflg = 1
                  isctchk(jel)  = 1
                else
                  isctflg = 0
                endif
                if ( (ipost .eq.1 .and. intpn .eq. intpt1) .or. 
     &                ipost .gt. 1) then
                  iptflg = 1
                else
                  iptflg = 0
                endif
              endif
            endif
            keyprv = key
            goto 9099
C
C==================================================================
C
C COORDINATES
C
C==================================================================
        else  if ( key .eq. 8 .and. ipost .gt. 1) then
C
C       PROCESS IF THIS IS THE CORRECT STEP/INC, ELEMENT, SECTION 
C       POINT, AND INTEGRATION POINT 
C
            icheck = itime*ielem*iptflg
            if (icheck .eq. 1 ) then
c            if (ipost .eq. 3) nintpt(jel) = intpn
               icrdchk(jel) = 1
               if (ipost .gt. 1) nintpt(jel) = intpn
               do kk=3,lenf
                  ii=kk-2
                  coords(ii,intpn)=array(kk)
               end do
C
            if (islct(key) .eq. 1 .and. ipost .lt. 3) then
              ivarchk(jel) = 1
              jvar = intpn
              nintpt(jel) = intpn
              do kk=3,lenf
                ii=kk-2
                var(ii,jvar)=array(kk)
              end do
              if (icomp .gt. ii) then
                 write(6,*)'Invalid component'
                 return
              endif
            endif
C
            endif
            keyprv = key
            goto 9099
C
C==================================================================
C
C VARIABLE RECORD
C
C==================================================================
C       PROCESS IF THIS IS THE DESIRED VARIABLE

        else if (islct(key) .eq. 1 .and. ipost .lt. 3) then
C
C       CONTINUE PROCESSING IF THIS IS THE CORRECT STEP/INC, ELEMENT, 
C       SECTION POINT, AND INTEGRATION POINT 
C
          icheck = itime*ielem*isctflg*iptflg
          if (icheck .eq. 1 ) then
            ivarchk(jel) = 1
            if (ipost .eq. 1) then
             jvar = jel
            else
             jvar = intpn
            endif
            nintpt(jel) = intpn
            do kk=3,lenf
              ii=kk-2
              var(ii,jvar)=array(kk)
            end do
            if (icomp .gt. ii) then
               write(6,*)'Invalid component'
               return
            endif
          endif
        keyprv = key
        goto 9099
C
C==================================================================
C
C    ANYTHING ELSE
C
C==================================================================
         else
c            write (6,9000) key
9000        format(I5,' *** NO CODING FOR THIS RECORD ***')
            keyprv = key
9099     end if
C
C==================================================================
C
C    END LOOPS
C
C==================================================================
 1000 continue
 1001 continue
C
C==================================================================
C     POSTPROCESS DATA HERE
C==================================================================
      do ii = 1, nels
         nintpt(ii) = nintpt(ii)/naxipt(ii)
      end do 
C
C==================================================================
C     VERIFY THAT APPROPRIATE DATA WERE WRITTEN TO RESULTS FILE
C==================================================================
c     iquit = 0
      if (ielchk .eq. 0) then
         write(6,*)'Desired element not found'
         return
      endif
C
      if (itimchk .eq. 0) then
         write(6,*)'Desired step/increment not found'
         return
      endif
C
      do iel = 1, nels
         if (isctchk(iel) .eq. 0 .and. ipost .lt. 3) then
            write(6,*)'Desired element and/or section point not found'
            return
         endif
         if (icrdchk(iel) .eq. 0 .and. ipost .gt. 1) then
            write(6,*)'Integration point coords not found'
            return
         endif
         if (ivarchk(iel) .eq. 0 .and. ipost .lt. 3) then
            write(6,*)'Desired element variable not found'
            return
         endif
      end do
C
C****************************************
C
C     PLOT VARIABLE ALONG ELBOW LENGTH
C
C****************************************
C
      if (ipost .eq. 1) then
C
         nn_old = iconn(nodel(1),1)
         do iel = 1, nels
            nn_new = iconn(1,iel)
            if (iel .gt. 1 .and. nn_new .ne. nn_old) then
               write(6,*)'Error in element connectivity'
               return
            endif
            nnum_b = iconn(1,iel)
            nnum_b = loc_nod(nnum_b)
            nnum_e = iconn(nodel(iel),iel)
            nnum_e = loc_nod(nnum_e)
            if (iel .eq. 1) xpos(nnum_b) = 0.0
            xlength = 0.d+0
            do jj = 1, 3
               delx = xyz(jj,nnum_e) - xyz(jj,nnum_b)
               xlength = xlength + delx*delx
            end do
            xlength = sqrt(xlength)
            xpos(nnum_e) = xpos(nnum_b) + xlength
            x_mid = 0.5*(xpos(nnum_e) + xpos(nnum_b))
C
C     EXTRACT VARIABLE INFORMATION AND WRITE TO FILE
C
            data = var(icomp,iel)
            write(9,2000)x_mid,data
 2000       format(5x,e17.9,5x,e17.9)
C
            nn_old=iconn(nodel(iel),iel)
C
         end do
         return
C
C
C****************************************
C
C     PLOT VARIABLE AROUND CIRCUMFERENCE OF ELBOW
C
C****************************************
C
       else if (ipost .eq. 2) then
C
         xpos(1) = 0.0
C
         locel = loc_el(jel_1)
         iaxial = 1
C
         if (naxipt(locel) .eq. 2) then
 3000       write(6,*)'Enter the axial station'
            read(5,*)iaxial
            if (iaxial .ne. 1 .and. iaxial .ne. 2) then
               write(6,*)'Try again - station must be 1 or 2'
               go to 3000
            endif
         endif
C
         do ipt = 1 , nintpt(locel)
            jpt = ipt + (iaxial - 1)*nintpt(locel)
            if (ipt .gt. 1) then
               xlength = 0.0
               do jj = 1, 3
                  delx = coords(jj,jpt) - coords(jj,jpt-1)
                  xlength = xlength + delx*delx
               end do
               xlength = sqrt(xlength)
               xpos(ipt) = xpos(ipt-1) + xlength
            endif
C
C     EXTRACT VARIABLE INFORMATION AND WRITE TO FILE
C
            data = var(icomp,jpt)
            write(9,2000)xpos(ipt),data
C
         end do

         ipt = 1
         jpt = ipt + (iaxial - 1)*nintpt(locel)
         data = var(icomp,jpt)
         xfinal = xpos(nintpt(locel)) + (xpos(2)- xpos(1)) 
         write(9,2000)xfinal,data

         return 
C
C****************************************
C
C     OVALIZATION PLOT
C
C****************************************
C
       else if (ipost .eq. 3) then
C
         locel = loc_el(jel_1)
         iaxial = 1
C
         if (naxipt(locel) .eq. 2) then
 3500       write(6,*)'Enter the axial station'
            read(5,*)iaxial
            if (iaxial .ne. 1 .and. iaxial .ne. 2) then
               write(6,*)'Try again - station must be 1 or 2'
               go to 3500
            endif
         endif
C
C     AVERAGE THE COORDS OF THE INTEGRATION POINTS TO GET THE CENTER
C            
         do i = 1, 3
            xc(i) = 0.0
         end do
C
         xnint = real(nintpt(locel))
C
         do ipt = 1 , nintpt(locel)
            jpt = ipt + (iaxial - 1)*nintpt(locel)
C
            do i = 1, 3
               xc(i) = xc(i) + coords(i,jpt)/xnint
            end do
C
         end do
C
C     SHIFT THE COORDINATES SO THAT CENTER OF SECTION IS AT ORIGIN
C
         do ipt = 1 , nintpt(locel)
            jpt = ipt + (iaxial - 1)*nintpt(locel)
            do i = 1, 3
               scord(i,jpt) = coords(i,jpt) - xc(i)
            end do
C
         end do
C
C     DETERMINE TANGENT VECTOR 1 (FROM CENTER OF CROSS-SECTION 
C                                 TO INT PT 1)
C
         jpt = 1 + (iaxial - 1)*nintpt(locel)
         xmag = 0.0
         do i = 1, 3
            tang1(i) = coords(i,jpt) - xc(i)
            xmag = xmag + tang1(i)*tang1(i)
         end do
         xmag = sqrt(xmag)
         do i = 1, 3
            tang1(i) =  tang1(i)/xmag
         end do
C
C     GUESS TANGENT VECTOR 2 (FROM CENTER OF CROSS-SECTION 
C                             TO INT PT 2)
C
         jpt = 2 + (iaxial - 1)*nintpt(locel)
         xmag = 0.0
         do i = 1, 3
            tang2(i) = coords(i,jpt) - xc(i)
            xmag = xmag + tang2(i)*tang2(i)
         end do
         xmag = sqrt(xmag)
         do i = 1, 3
            tang2(i) =  tang2(i)/xmag
         end do
C
C     DETERMINE THE NORMAL VECTOR, TANG1 X TANG2
C
         xnorm(1) =   tang1(2)*tang2(3) - tang1(3)*tang2(2)
         xnorm(2) = - tang1(1)*tang2(3) + tang1(3)*tang2(1)
         xnorm(3) =   tang1(1)*tang2(2) - tang1(2)*tang2(1)

         xmag = 0.0
         do i = 1, 3
            xmag = xmag + xnorm(i)*xnorm(i)
         end do
         xmag = sqrt(xmag)
         do i = 1, 3
            xnorm(i) = xnorm(i)/xmag
         end do
C
C     REDEFINE TANGENT VECTOR 2 
C
         tang2(1) = - tang1(2)*xnorm(3) + tang1(3)*xnorm(2)
         tang2(2) =   tang1(1)*xnorm(3) - tang1(3)*xnorm(1)
         tang2(3) = - tang1(1)*xnorm(2) + tang1(2)*xnorm(1)
C
C     DETERMINE THE ROTATION MATRIX
C
         do i = 1, 3
            rot(1,i) = tang1(i)
            rot(2,i) = tang2(i)
            rot(3,i) = xnorm(i)
         end do
C         
C     TRANSFORM COORDINATES INTO LOCAL SYSTEM
C
         do ipt = 1 , nintpt(locel)
            jpt = ipt + (iaxial - 1)*nintpt(locel)
            do i = 1, 3
               xpr(i,ipt) = 0.0
               do j = 1, 3
                  xpr(i,ipt) = xpr(i,ipt) + rot(i,j)*scord(j,jpt)
               end do
            end do
         end do
C
C     OUTPUT COORDS TO FILE
C
         do ipt = 1, nintpt(locel)
            write(9,4000)(xpr(i,ipt),i=1,2)
 4000       format(5x,e17.9,5x,e17.9,5x,e17.9)
         end do
C
         ipt = 1
         write(9,4000)(xpr(i,ipt),i=1,2)

       endif
C
      close (unit=9)
C
      IF( ALLOCATED(var) ) DEALLOCATE(var)

      return
      end