text:
      SUBROUTINE ESCAPE
*
*
*       Escaper removal.
*       ----------------
*
      INCLUDE 'common6.h'
      COMMON/BINARY/  CM(4,MMAX),XREL(3,MMAX),VREL(3,MMAX),
     &                HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
     &                NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAG(MMAX)
      CHARACTER*11  WHICH1
*
*
*       Adopt twice the tidal radius as escape condition.
      RESC2 = 4.0*RTIDE**2
*       For tidal cutoff check only energy
      IF(KZ(23).GE.3)THEN
          RESC2 = 0.D0
          ETID = ZMASS/RTIDE
          ZMOLD = ZMASS
      END IF
*
      RTIDE2 = RTIDE**2
      NCORR = 0
      NCRIT1 = 0
      NCRIT2 = 0
      I3HI = 0
      DO 1 K = 1,3
          CMR(K) = 0.0D0
          CMRDOT(K) = 0.0D0
    1 CONTINUE
*
*       Set the distance (squared) with respect to the density centre.
      I = IFIRST
    5 RI2 = (X(1,I) - RDENS(1))**2 + (X(2,I) - RDENS(2))**2 +
     &                               (X(3,I) - RDENS(3))**2
      IF (RI2.LT.RTIDE2) NCRIT1 = NCRIT1 + 1
*
      IF(KZ(23).LE.2) THEN
*       See whether escape is indicated (retain ghost particles).
      IF (RI2.GT.RESC2.AND.RI2.LT.1.0D+10) THEN
          RJMIN2 = 1.0D+10
*       Find distance to the nearest neighbour and calculate potential.
          POTI = 0.0D0
          DO 8 J = IFIRST,NTOT
              IF (J.EQ.I) GO TO 8
              RIJ2 = (X(1,I) - X(1,J))**2 + (X(2,I) - X(2,J))**2 +
     &                                      (X(3,I) - X(3,J))**2
              POTI = POTI + BODY(J)/SQRT(RIJ2)
              IF (RIJ2.LT.RJMIN2) THEN
                  RJMIN2 = RIJ2
                  JMIN = J
              END IF
    8     CONTINUE
          VI2 = XDOT(1,I)**2 + XDOT(2,I)**2 + XDOT(3,I)**2
*       Check escape criterion for tidal case or isolated system.
          EI = 0.5*VI2 - POTI
          IF (KZ(14).GT.0.OR.EI.GT.0.0) GO TO 30
          IF (NAME(I).LT.0) GO TO 30
      END IF
*
      ELSE
*       Check tidal cutoff criterion - use last computed phi-value
*       Find distance to the nearest neighbour
          POTI = 0.0D0
          NNB = LIST(1,I)
          RJMIN2 = 1.0D+10
          DO 9 L = 1,NNB
              J = LIST(L+1,I)
              RIJ2 = (X(1,I) - X(1,J))**2 + (X(2,I) - X(2,J))**2 +
     &                                      (X(3,I) - X(3,J))**2
              IF (RIJ2.LT.RJMIN2) THEN
                  RJMIN2 = RIJ2
                  JMIN = J
              END IF
    9     CONTINUE
          VI2 = XDOT(1,I)**2 + XDOT(2,I)**2 + XDOT(3,I)**2
          POTI = -PHIDBL(I)
*       Check escape criterion for tidal case or isolated system.
          EI = 0.5*VI2 - POTI + ETID
*
          IF (EI.GT.0.0.AND.RI2.GT.RTIDE2) THEN
          NCRIT2 = NCRIT2 + 1
          GO TO 30
          END IF
      END IF
*
   10 I = I + 1
   12 IF (N.EQ.2) GO TO 13
*
      IF (I.LE.NTOT) GO TO 5
      IF (NCORR.EQ.0) GO TO 25
*
*       Form centre of mass terms.
   13 DO 15 I = 1,N
          DO 14 K = 1,3
              CMR(K) = CMR(K) + BODY(I)*X(K,I)/ZMASS
              CMRDOT(K) = CMRDOT(K) + BODY(I)*XDOT(K,I)/ZMASS
   14     CONTINUE
   15 CONTINUE
      CMR(4) = SQRT(CMR(1)**2 + CMR(2)**2 + CMR(3)**2)
*
      JLAST = MIN(2*NCORR,NMAX)
      if(rank.eq.0)
     &WRITE (6,18)  N, NSESC, NBESC, ZMASS, BE(3), CMR(4), RESC, STEPI,
     &              RSI, ZMASS/FLOAT(N), NCRIT1, (JLIST(J),J=1,JLAST)
   18 FORMAT (/,' ESCAPE    N =',2I6,I4,F8.4,F12.6,2F7.2,F7.3,F7.2,F9.5,
     &                                        I6,2X,6I6,/,5(10X,20I6,/))
*
      IF (KZ(23).EQ.4)THEN
          JLAST = MIN(2*NCORR,NMAX)
          if(rank.eq.0) WRITE (6,20)  (ILIST(J),J=1,JLAST)
   20     FORMAT (/,' ESCAPE ANGLES ',11(2I4,2X),9(/,15X,11(2I4,2X)))
      END IF
*
*       Check updating of global index for chain c.m.
      IF (NCH.GT.0) THEN
          CALL CHFIND
      END IF
*
*       Set phase indicator < 0 for new NLIST in routine INTGRT (Hermite).
      IPHASE = -1
*
   25 RETURN
*
   30 A2 = (X(1,JMIN) - RDENS(1))**2 + (X(2,JMIN) - RDENS(2))**2 +
     &                                 (X(3,JMIN) - RDENS(3))**2
*       See whether nearest body satisfies escape condition or RIJ > 10*RMIN.
      IF (A2.GT.RESC2.AND.KZ(14).GT.0) GO TO 40
      IF (RJMIN2.GT.100.0*RMIN2) GO TO 40
*
      IF(KZ(23).GE.3)THEN
          VI2 = XDOT(1,JMIN)**2 + XDOT(2,JMIN)**2 + XDOT(3,JMIN)**2
          EI = 0.5*VI2 + PHIDBL(JMIN) + ETID
          IF(EI.GT.0.D0) GO TO 40
      END IF
*
      A3 = XDOT(1,JMIN)**2 + XDOT(2,JMIN)**2 + XDOT(3,JMIN)**2
      A4 = (XDOT(1,I) - XDOT(1,JMIN))**2 +
     &     (XDOT(2,I) - XDOT(2,JMIN))**2 +
     &     (XDOT(3,I) - XDOT(3,JMIN))**2
      A5 = (BODY(I) + BODY(JMIN))/SQRT(RJMIN2)
*       Check velocity of binary component in case of bound pair.
      A6 = 2.0*A5/A4
      IF (A6.GT.1.0.AND.A3.GT.2.0*ZMASS/SQRT(A2)) GO TO 40
*       Retain escaper if dynamical effect on neighbour is significant.
      IF (A6.GT.0.01) GO TO 10
*
*       Form optional output diagnostics.
   40 X1 = X(1,I) - RDENS(1)
      Y1 = X(2,I) - RDENS(2)
      Z1 = X(3,I) - RDENS(3)
      A3 = ABS(Y1/X1)
      A4 = ABS(Z1)/SQRT(X1**2 + Y1**2)
      ILIST(2*NCORR+1) = DATAN(A3)*180.0/3.14159
      ILIST(2*NCORR+2) = DATAN(A4)*180.0/3.14159
*       Escape angles with respect to the X-axis and XY-plane.
      RESC = SQRT(RI2)
      STEPI = STEP(I)
      RSI = RS(I)
*
*       Accumulate escaper names and save current name in case I > N.
      NCORR = NCORR + 1
      JLIST(NCORR) = NAME(I)
      NAMEI = NAME(I)
      KSTARI = KSTAR(I)
      IF (NAME(I).LT.0) KSTARI = KSTARI + 30
      IF (BODY(I).LE.0.0D0) GO TO 50
*
*       Obtain binding energy of body #I and update optional tidal radius.
      ZK = 0.5D0*BODY(I)*VI2
      IF (KZ(14).GT.0) THEN
          CALL XTRNLV(I,I)
          ZK = ZK + HT
          RTIDE = (ZMASS/TIDAL(1))**0.3333
      END IF
      EI = ZK - BODY(I)*POTI
*
*       Update tidal radius in case of tidal cutoff (R.Sp.)
      IF(KZ(23).GE.3)THEN
          RTOLD = RTIDE
          RTIDE = RTIDE*(ZMASS/ZMOLD)**0.3333
      END IF
*
*       Correct total energy.
      BE(3) = BE(3) - EI
*
*       Update total mass and save energy & number of single escapers.
      IF (I.LE.N) THEN
          ZMASS = ZMASS - BODY(I)
          E(4) = E(4) + EI
          NPOP(4) = NPOP(4) + 1
          NSESC = NSESC + 1
      END IF
*
*       Include optional escape output on unit 11.
      IF (KZ(23).EQ.2.OR.KZ(23).EQ.4) THEN
          TESC = TSCALE*TTOT
          EESC = EI/BODY(I)
          VB2 = 2.0*ZKIN/ZMASS
*       Distinguish between tidal field and isolated system (ECRIT at RTIDE).
          IF (KZ(14).GT.0) THEN
              ECRIT = -1.5*(TIDAL(1)*ZMASS**2)**0.3333
              EESC = 2.0*(EESC - ECRIT)/VB2
              EESC = MIN(EESC,999.0D0)
              BESC = ZMBAR*BODY(I)
          ELSE
              EESC = 2.0*EESC/VB2
              BESC = BODY(I)/BODYM
          END IF
          VKM = SQRT(VI2)*VSTAR
      if(rank.eq.0)then
          WRITE (11,45)  TESC, BESC, EESC, VKM, KSTARI
   45     FORMAT (F8.1,3F6.1,I4)
      end if
      END IF
*
*       Reduce particle number & total membership and check NNBMAX.
   50 N = N - 1
      NTOT = NTOT - 1
      NNBMAX = MIN(N/2,NNBMAX)
      ZNBMAX = 0.9*FLOAT(NNBMAX)
      ZNBMIN = MAX(0.01*FLOAT(NNBMAX),1.0)
*       Set indicator for removal of c.m. or KS components.
      KCASE = 1
*
*       Update COMMON arrays to remove the escaper and correct F & FDOT.
      CALL REMOVE(I,1)
*
*       Delete escaper from neighbour lists and reduce higher locations.
   60 DO 150 J = 1,NTOT
          NNB = LIST(1,J)
          IF (NNB.EQ.0) GO TO 150
          L = 2
   70     IF (LIST(L,J).NE.I) GO TO 130
*
*       Move up the remaining list members and reduce membership by one.
          DO 80 K = L,NNB
              LIST(K,J) = LIST(K+1,J)
   80     CONTINUE
          LIST(1,J) = LIST(1,J) - 1
*       Reduce the steps to minimize error effect (do not allow DT < 0).
*         STEP(J) = MAX(0.5D0*STEP(J),TIME - T0(J))
*         STEPR(J) = MAX(0.5D0*STEPR(J),TIME - T0R(J))
          IF (LIST(1,J).GT.0) GO TO 130
*
*       Add a distant body as neighbour if list only contains escaper.
          K = IFIRST - 1
  100     K = K + 1
          RJK2 = (X(1,J) - X(1,K))**2 + (X(2,J) - X(2,K))**2 +
     &                                  (X(3,J) - X(3,K))**2
          IF (RJK2.LT.RESC2.AND.K.LT.N) GO TO 100
          LIST(1,J) = 1
          LIST(2,J) = K
          GO TO 150
*
*       Reduce higher particle locations by one.
  130     IF (LIST(L,J).GT.I) LIST(L,J) = LIST(L,J) - 1
          L = L + 1
          IF (L.LE.LIST(1,J) + 1) GO TO 70
  150 CONTINUE
*
*       Update list of old KS components (remove #I and rename > I).
      NNB = LISTR(1)
      DO 170 L = 2,NNB+1
          IF (LISTR(L).EQ.I) THEN
*       Remove both components of pair and reduce membership by two.
              J = 2*KVEC(L-1)
              DO 165 K = J,NNB-1
                  LISTR(K) = LISTR(K+2)
  165         CONTINUE
              LISTR(1) = LISTR(1) - 2
          END IF
  170 CONTINUE
*
*       Reduce higher particle locations by one (separate loop for pairs).
      DO 180 L = 2,NNB+1
          IF (LISTR(L).GT.I) LISTR(L) = LISTR(L) - 1
  180 CONTINUE
*
*       Update list of high velocity particles (remove #I and rename > I).
      NNB = LISTV(1)
      DO 190 L = 2,NNB+1
          IF (LISTV(L).EQ.I) THEN
*       Remove escaper and reduce the membership.
              DO 185 K = L,NNB
                  LISTV(K) = LISTV(K+1)
  185         CONTINUE
              LISTV(1) = LISTV(1) - 1
          END IF
*       Reduce higher particle locations by one (or three for c.m.).
          IF (LISTV(L).GT.I) THEN
              LISTV(L) = LISTV(L) - 1
              IF (I.GT.N + 1) LISTV(L) = LISTV(L) - 2
          END IF
  190 CONTINUE
*
*       Check special case of second KS component removal.
      IF (KCASE.GT.1) GO TO 205
*
*       See whether the escaper is a single particle or c.m.
      IF (I.LE.N + 1) GO TO 12
*
*       Prepare removal of regularized pair.
      IPAIR = I - N - 1
*
*       Skip correction if ghost is also merged binary (NAMEI = 0 below).
      IF (NAMEI.NE.0.AND.BODY(2*IPAIR-1).GT.0.0D0) THEN
          ZMB = BODY(2*IPAIR-1) + BODY(2*IPAIR)
          ZM2 = BODY(2*IPAIR)
          EB = H(IPAIR)*BODY(2*IPAIR-1)*BODY(2*IPAIR)/ZMB
          SEMI = -0.5*ZMB/H(IPAIR)
          ECC2 = (1.0 - R(IPAIR)/SEMI)**2 + TDOT2(IPAIR)**2/(SEMI*ZMB)
          ECC = SQRT(ECC2)
          PMIN  = SEMI*(1.0 - ECC)
          RATIO = MAX(RADIUS(2*IPAIR-1),RADIUS(2*IPAIR))/PMIN
          NAME1 = NAME(2*IPAIR-1)
          KSTAR1 = KSTAR(2*IPAIR-1)
          PB = DAYS*SEMI*SQRT(SEMI/ZMB)
      ELSE
*       Obtain two-body elements of ghost binary and update energies.
          EB = 0.0D0
          RATIO = 0.0
          BODYCM = CM(3,IM) + CM(4,IM)
          SEMI = -0.5*BODYCM/H(JPAIR)
          ZMU = CM(3,IM)*CM(4,IM)/BODYCM
          PB = DAYS*SEMI*SQRT(SEMI/BODYCM)
          ECC2 = (1.0-R(JPAIR)/SEMI)**2 + TDOT2(JPAIR)**2/(BODYCM*SEMI)
          ECC = SQRT(ECC2)
          EB1 = ZMU*H(JPAIR)
          BE(3) = BE(3) - EB1
          EMERGE = EMERGE - EB1
          DEB = DEB + EB1
      END IF
*
*       Update total energy (ECOLL with EB < -1 & #27 > 0 affects BINOUT).
      BE(3) = BE(3) - EB
*
*       Check optional diagnostics for hierarchical systems.
      IF (NAMEI.LE.0.AND.(KZ(11).EQ.1.OR.KZ(11).EQ.3)) THEN
          IPHASE = 7
          CALL HIARCH(IPAIR)
      END IF
*
*       Specify binary type (0: standard; -1: merged binary).
      M = 0
      IF (NAMEI.LE.0) M = -1
*
*       Distinguish between actual and ghost binary (mergers come later).
      IF (BODY(2*IPAIR-1).GT.0.0D0) THEN
          NAME2 = NAME(2*IPAIR)
*
*       Include rare case of higher-order system (4, 5 or 6 members).
          IF (NAMEI.LT.-2*NZERO) THEN
              JM = 1
              DO 195 K = 1,NMERGE
                  IF (NAMEM(K).EQ.NAMEI) JM = K
  195         CONTINUE
              I3HJ = N
              DO 196 J = IFIRST,NTOT
                  IF (NAME(J).EQ.NAMEG(JM)) I3HJ = J
  196         CONTINUE
*       Identify quartet, quintuplet or sextuplet.
              IF (NAME2.LE.NZERO.AND.NAME(I3HJ).LE.NZERO) THEN
                  WHICH1 = '  QUARTET  '
*       Include both types of quintuplet: [[B,S],B] and [[B,B],S].
              ELSE IF (NAME2.LE.NZERO.OR.NAME(I3HJ).LE.NZERO) THEN
                  WHICH1 = ' QUINTUPLET'
              ELSE
                  WHICH1 = ' SEXTUPLET '
              END IF
              ZM1 = CM(1,JM)*ZMBAR
              ZM2 = CM(2,JM)*ZMBAR
              EB3 = CM(1,JM)*CM(2,JM)*HM(JM)/(CM(1,JM) + CM(2,JM))
              SEMI3 = -0.5*(CM(1,JM) + CM(2,JM))/HM(JM)
              PB3 = DAYS*SEMI3*SQRT(SEMI3/(CM(1,JM) + CM(2,JM)))
              if(rank.eq.0)
     &        WRITE (6,199)  WHICH1, NAME(2*IPAIR-1), NAME(2*IPAIR),
     &                       NAME(I3HJ), ZM1, ZM2, EB3, SEMI3, PB3
  199         FORMAT (/,A11,' ESCAPE    NM =',3I6,'  M =',2F5.1,
     &                    '  EB3 =',F8.4,'  A3 =',1P,E8.1,'  P3 =',E8.1)
          END IF
*
          if(rank.eq.0)then
          VI = SQRT(0.5*VI2*ZMASS/ZKIN)
          WRITE (6,200)  IPAIR, NAME(2*IPAIR-1), NAME2,
     &                   KSTAR(2*IPAIR-1), KSTAR(2*IPAIR), KSTARI,
     &                   LIST(2,2*IPAIR), BODY(2*IPAIR-1)*ZMBAR,
     &                   BODY(2*IPAIR)*ZMBAR, EB, RATIO, VI, ECC, EI, PB
  200     FORMAT (/,' BINARY ESCAPE    KS =',I5,'  NM =',2I6,
     &                '  K* =',4I3,'  M =',2F5.1,'  EB =',F8.4,
     &                '  R*/PM =',F5.2,'  V/<V> =',F5.2,'  E =',F5.2,
     &                '  EI =',F8.5,'  P =',1P,E8.1)
          end if
      ELSE
          if(rank.eq.0)
     &    WRITE (6,202)  IPAIR, NAME(2*IPAIR-1), NAME(2*IPAIR),
     &                   LIST(2,2*IPAIR), KSTAR(JM), CM(3,IM)*ZMBAR,
     &                   CM(4,IM)*ZMBAR, EB1, ECC, SEMI, PB
  202     FORMAT (/,' QUAD ESCAPE    KS =',I5,'  NM =',2I6,
     &                '  K* =',2I3,'  M =',2F5.1,'  EB =',F8.4,
     &                '  E =',F7.3,'  A =',1P,E8.1,'  P =',E8.1)
          EB = EB + EB1
          EI = 0.0
      END IF
*
*       Accumulate escaping binary energies and increase the counter.
      IF (LIST(2,2*IPAIR).EQ.-1) THEN
          E(5) = E(5) + EB
          E(6) = E(6) + EI
          NPOP(5) = NPOP(5) + 1
      ELSE
          E(7) = E(7) + EB
          E(8) = E(8) + EI
          NPOP(6) = NPOP(6) + 1
      END IF
*
      IF (M.EQ.0) THEN
          NBESC = NBESC + 1
      ELSE
          NMESC = NMESC + 1
      END IF
*
*       Reduce particle number, pair index & single particle index.
      N = N - 1
      NPAIRS = NPAIRS - 1
      IFIRST = 2*NPAIRS + 1
*
*       Move up all tables of regularized pairs below IPAIR.
      IF (IPAIR.LE.NPAIRS) THEN
          CALL REMOVE(IPAIR,2)
      END IF
*
*       Increase index for removing KS components.
  205 KCASE = KCASE + 1
      IF (KCASE.LE.3) THEN
*       Remove COMMON arrays of the second component before the first.
          I = 2*IPAIR + 2 - KCASE
          NTOT = NTOT - 1
*       Reduce NTOT by 3 and N by 2 when KS pair escapes.
          CALL REMOVE(I,3)
          GO TO 60
      END IF
*
*       Check selection of possible ghost in higher-order system.
      IF (I3HI.GT.0) THEN
          I = 0
          DO 208 J = IFIRST,NTOT
              IF (NAME(J).EQ.NAMEG(JM)) THEN
                  I = J
              END IF
  208     CONTINUE
          IF (I.GT.0) THEN
              I3HI = 0
              GO TO 50
          END IF
      END IF
*
*       Include the case of escaping merger.
      IF (NAMEI.GE.0) GO TO 250
*
*       Locate current position in the merger table (standard case).
      IM = 0
      DO 210 K = 1,NMERGE
          IF (NAMEM(K).EQ.NAMEI) IM = K
  210 CONTINUE
*       Skip on failed detection just in case.
      IF (IM.EQ.0) GO TO 250
*
*       Include case of higher-order system (outer single or binary).
      DEB = 0.0
      IF (NAMEI.LT.0) THEN
*       Determine the ghost index.
          I3HI = 0
          DO 215 J = IFIRST,NTOT
              IF (NAME(J).EQ.NAMEG(IM)) I3HI = J
  215     CONTINUE
*       Specify nominal escape distance for ghost removal.
          X(1,I3HI) = 1.0D+04
*       Define possible KS pair index for quadruple system correction.
          JPAIR = I3HI - N
      END IF
*
*       Consider current ghost unless deeper hierarchy is present.
      JM = IM
      IF (I3HI.GT.0) THEN
*       Search for c.m. name one level below current (NAMEI < 0).
          IF (NAMEI.LT.-2*NZERO) THEN
              DO 220 K = 1,NMERGE
                  IF (NAMEM(K).EQ.NAMEI + 2*NZERO) JM = K
  220         CONTINUE
*       Use previous merger index to look for binary ghost at earlier level.
              I30 = I3HI
              DO 225 J = IFIRST,NTOT
                  IF (NAME(J).EQ.NAMEG(JM)) I3HI = J
  225         CONTINUE
              IF (I3HI.GT.N) THEN
                  JPAIR = I3HI - N
              ELSE
*       Adopt original solution on failure to identify binary.
                  I3HI = I30
              END IF
*       Set nominal escape distance of any new ghost (I3HI <= N possible).
              X(1,I3HI) = 1.0D+04
          END IF
      END IF
*
*       Copy merger energy to respective energy bins (primordial or new).
      ZMU = CM(1,JM)*CM(2,JM)/(CM(1,JM) + CM(2,JM))
      IF (MIN(NAME1,NAMEG(JM)).LE.2*NBIN0) THEN
          E(5) = E(5) + ZMU*HM(JM) + DEB
      ELSE
          E(7) = E(7) + ZMU*HM(JM) + DEB
      END IF
*
*       Reduce membership if IM is last (otherwise done in RESET).
      IF (IM.EQ.NMERGE) THEN
          NMERGE = NMERGE - 1
      END IF
*
*       Identify merged ghost particle (single body or binary c.m.).
      JCOMP = -1
      DO 230 J = IFIRST,NTOT
          IF (BODY(J).EQ.0.0D0.AND.NAME(J).EQ.NAMEG(IM)) JCOMP = J
  230 CONTINUE
*       Include search over lower level on failed identification.
      IF (JCOMP.EQ.-1.AND.I3HI.GT.0) THEN
          DO 232 J = IFIRST,NTOT
              IF (BODY(J).EQ.0.0D0.AND.NAME(J).EQ.NAMEG(JM)) JCOMP = J
  232     CONTINUE
      END IF
*
*       Skip if correct ghost not identified (note I3HI # JCOMP if JM # IM).
      IF (JCOMP.GT.0) THEN
          I = I3HI
*       Form two-body elements and period of inner binary.
          NAMEI = 0
          NPOP(7) = NPOP(7) + 1
          BODYCM = CM(1,JM) + CM(2,JM)
          SEMI0 = -0.5*BODYCM/HM(JM)
          ZMU = CM(1,JM)*CM(2,JM)/BODYCM
          EB = ZMU*HM(JM)
          PB = DAYS*SEMI0*SQRT(SEMI0/BODYCM)
          BE(3) = BE(3) - EB
          EMERGE = EMERGE - EB
*       Include extra correction for mergers between binary pairs.
          IF (JM.LT.IM) THEN
              ZMU2 = CM(1,IM)*CM(2,IM)/(CM(1,IM) + CM(2,IM))
              EB2 = ZMU2*HM(IM)
              BE(3) = BE(3) - EB2
              EMERGE = EMERGE - EB2
*       Copy smaller index for QUAD binary output.
              IM = JM
          END IF
          RJ = 0.0
          TD2 = 0.0
          DO 235 K = 1,4
              RJ = RJ + UM(K,JM)**2
              TD2 = TD2 + 2.0*UM(K,JM)*UMDOT(K,JM)
  235     CONTINUE
          ECC2 = (1.0 - RJ/SEMI0)**2 + TD2**2/(BODYCM*SEMI0)
          ECC0 = SQRT(ECC2)
          Q = ZM2/BODYCM
          XFAC = (1.0 + Q)*(1.0 + ECC)/SQRT(1.0 - ECC)
          PCRIT = 2.8*XFAC**0.4*SEMI0
*
          if(rank.eq.0)
     &    WRITE (6,240)  NAME1, NAMEG(JM), KSTAR1, KSTAR(JCOMP),
     &                   KSTARM(JM), CM(1,JM)*ZMBAR, CM(2,JM)*ZMBAR,
     &                   EB, ECC0, PMIN/PCRIT, SEMI0, PB
  240     FORMAT (/,' HIARCH ESCAPE    NM =',2I6,'  K* =',3I3,
     &              '  M =',2F5.1,'  EB =',F8.4,'  E =',F7.3,
     &              '  PM/PC =',1P,E8.1,'  A =',E8.1,'  P =',E8.1)
*
*       Remove the ghost particle (NAME = 0 & EB = 0 for second binary).
          GO TO 50
      END IF
*
  250 I = IPAIR + N
      GO TO 12
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!>    @brief domain of validity for module const
!>    @param[in] ismpl local sample index
!>    @details
!>    Checking whether pres/temp/(conc) are in domain of props
!>    validity. Version for property module const. \n
!>    \n
!>    For concentration, an error is thrown and the execution is
!>    stopped if the concentration is outside the physical values.
      subroutine check_domain(ismpl)
        use arrays, only: pres, temp, conc
        use mod_genrl, only: i0, j0, k0
        use mod_genrlc, only: def_props
        use mod_conc, only: ntrac
        use mod_linfos, only: linfos

        implicit none

        ! Sample index
        integer :: ismpl

        ! Iteration counters
        integer :: i, j, k, l

        ! counters for the values outside domain of validity
        ! pres
        integer :: icountp
        ! temp
        integer :: icountt
        ! conc
        integer :: icountc

        ! min/max boundaries of the domain of validity
        ! pres
        double precision, parameter ::  pmin = 0.01d6
        double precision, parameter ::  pmax = 110.0d6
        ! temp
        double precision, parameter ::  tmin = 0.0d0
        double precision, parameter ::  tmax = 350.0d0
        ! conc
        double precision, parameter ::  cmin = 0.0d0
        double precision, parameter ::  cmax = 1.0d5
        ! numerical boundary
        double precision, parameter ::  csmin = 1.0d-30

        ! records the overall min/max of values if they are outside
        ! domain of validity
        double precision :: dpmax, dtmax, dcmax, dhmax
        double precision :: dpmin, dtmin, dcmin, dhmin

        intrinsic trim

        ! Set counters to zero
        icountp = 0
        icountt = 0
        icountc = 0

        ! Set overall min/max to boundaries of the domain of validity
        dpmax = pmax
        dpmin = pmin
        dtmax = tmax
        dtmin = tmin
        dcmax = cmax
        dcmin = cmin

        ! Check pres
        do k = 1, k0
          do j = 1, j0
            do i = 1, i0
              if (pres(i,j,k,ismpl)<pmin) then
                ! Set min counter
                icountp = icountp + 1
                ! Set new overall minimum
                dpmin = min(dpmin,pres(i,j,k,ismpl))
                ! Change pres value to minimum of the domain of validity
                pres(i,j,k,ismpl) = pmin
              end if
              if (pres(i,j,k,ismpl)>pmax) then
                ! Set max counter
                icountp = icountp + 1
                ! Set new overall maximum
                dpmax = max(dpmax,pres(i,j,k,ismpl))
                ! Change pres value to maximum of the domain of validity
                pres(i,j,k,ismpl) = pmax
              end if
            end do
          end do
        end do

        ! Check temp
        do k = 1, k0
          do j = 1, j0
            do i = 1, i0
              if (temp(i,j,k,ismpl)<tmin) then
                icountt = icountt + 1
                dtmin = min(dtmin,temp(i,j,k,ismpl))
                temp(i,j,k,ismpl) = tmin
              end if
              if (temp(i,j,k,ismpl)>tmax) then
                icountt = icountt + 1
                dtmax = max(dtmax,temp(i,j,k,ismpl))
                temp(i,j,k,ismpl) = tmax
              end if
            end do
          end do
        end do

        ! Check conc
        do k = 1, k0
          do j = 1, j0
            do i = 1, i0
              do l = 1, ntrac
                if (conc(i,j,k,l,ismpl).gt.cmax) then
                  icountc = icountc +1
                  dcmax = max(dcmax, conc(i,j,k,l,ismpl))
                  conc(i,j,k,l,ismpl) = cmax
                end if
                if (conc(i,j,k,l,ismpl)<cmin .and. &
                  conc(i,j,k,l,ismpl)<-csmin) then
                  icountc = icountc + 1
                  dcmin = min(dcmin,conc(i,j,k,l,ismpl))
                  conc(i,j,k,l,ismpl) = cmin
                end if
                if (conc(i,j,k,l,ismpl)<csmin) then
                  ! very small conc values set to zero to avoid
                  ! numerically instabilities
                  conc(i,j,k,l,ismpl) = cmin
                end if
              end do
            end do
          end do
        end do

!       disable the warning output for linfos(3)==-1
        if (linfos(3)>=0) then
          if (icountp/=0) write(*,'(3A,1I8,1A,1e16.7,1A,1e16.7,1A)') &
            'warning: pres not in domain of validity of module <', &
            trim(def_props), '> at ', icountp, ' points (min', dpmin, &
            ', max', dpmax, ')!'
          if (icountt/=0) write(*,'(3A,1I8,1A,1e16.7,1A,1e16.7,1A)') &
            'warning: temp not in domain of validity of module <', &
            trim(def_props), '> at ', icountt, ' points (min', dtmin, &
            ', max', dtmax, ')!'
          if (icountc/=0) write(*,'(3A,1I8,1A,1e16.7,1A,1e16.7,1A)') &
            'warning: conc not in domain of validity of module <', &
            trim(def_props), '> at ', icountc, ' points (min', dcmin, &
            ', max', dcmax, ')!'

          ! error outputs for hard physical concentration boundaries
          if (dcmax > cmax) then
            write(unit = *, fmt = *) "[E1] Error in check_domain.f90:", &
                "  maximum concentration dcmax= ", dcmax, &
                " larger than allowed maximum value cmax=", cmax
            stop
          end if
          if (dcmin > cmin) then
            write(unit = *, fmt = *) "[E2] Error in check_domain.f90:", &
                "  minimum concentration dcmin= ", dcmin, &
                " smaller than allowed minimum value cmin=", cmin
            stop
          end if
        end if

        return

      end subroutine check_domain


! { dg-do compile }
!
! Check whether MODULE PROCEDUREs are properly treated
! They need to be contained in a procedure, i.e. an
! interface in another procedure is invalid; they may, however,
! come from a use-associated procedure.
! (The PROCEDURE statement allows also for non-module procedures
!  if there is an explicit interface.)
!
! PR fortran/33228
!
module inclmod
  implicit none
  interface
    subroutine wrong1(a)
      integer :: a
    end subroutine wrong1
  end interface
  interface gen_incl
    module procedure ok1
  end interface gen_incl
  external wrong2
  external wrong3
  real wrong3
contains
  subroutine ok1(f)
    character :: f
  end subroutine ok1
end module inclmod

module a
  use inclmod
  implicit none
  interface gen
    subroutine ok1_a(a,b)
      integer :: a,b
    end subroutine ok1_a
    module procedure ok1, ok2_a
  end interface gen
contains
  subroutine ok2_a(a,b,c)
     integer :: a,b,c
  end subroutine ok2_a
end module a

module b
  use inclmod
  interface gen_wrong_0
    module procedure gen_incl  ! { dg-error "Cannot change attributes" }
  end interface gen_wrong_0
end module b

module c
  use inclmod
  interface gen_wrong_1
    module procedure wrong1  ! { dg-error "is not a module procedure" }
  end interface gen_wrong_1
end module c

module d
  use inclmod
  interface gen_wrong_2
    module procedure wrong2  ! { dg-error "Cannot change attributes" }
  end interface gen_wrong_2
end module d

module e
  use inclmod
  interface gen_wrong_3
    module procedure wrong3  ! { dg-error "Cannot change attributes" }
  end interface gen_wrong_3
end module e

module f
  implicit none
  interface
    subroutine wrong_a(a)
      integer :: a
    end subroutine wrong_a
  end interface
  interface gen_wrong_4
    module procedure wrong_a  ! { dg-error "is not a module procedure" }
  end interface gen_wrong_4
end module f

module g
  implicit none
  external wrong_b
  interface gen_wrong_5
    module procedure wrong_b  ! { dg-error "has no explicit interface" }
  end interface gen_wrong_5
end module g

module h
  implicit none
  external wrong_c
  real wrong_c
  interface gen_wrong_6
    module procedure wrong_c  ! { dg-error "has no explicit interface" }
  end interface gen_wrong_6
end module h

end


text:
!> \file exponential_cov.f90
!! Set of tools to build exponential covariance matrices.
!! By now, this toolbox support only the uniformly spaced grid: fixed discretization step.
!!  \author Innocent Souopgui
!!  \date 12/2013
!!  \version 1.0
!<

!> Procedures to build exponential covariance matrices
!!
!<
module exponential_cov
  use general_constant
  use general_tools
  use debug_tools

implicit none

  private
  public exp_cov_param, se_cov_matrix, load_ecp

  !> Derived type for 1D exponential covariance matrix parameters
  type exp_cov_param
    integer:: n !< Number of grid point
    integer:: k !< Number of significant eigenvalues/singular values to be kept
    real(dp):: dx !< discretization step
    real(dp):: s !< smoothing parameter of the covariance function (standard deviation at a single point)
    real(dp):: l !< decorrelation length
  end type exp_cov_param

  !> Computes a square exponential covariance matrix
  interface se_cov_matrix
    module procedure se_cov_matrix_scalar
    module procedure se_cov_matrix_type
  end interface se_cov_matrix

contains

  !> \brief Loads the parameters of an exponential covariance matrix from namelist file
  !! \param[in] ecp data structure containing the exponential covariance matrix parameters
  !! \param[in] fName name of the namelist file
  !! \param[in] dName name associated with physical space dimension where the covariance matrix will be applied.
  !! \param[in] n size of the covariance matrix (number of columns/rows)
  !! \param[in] delta discretization step
  !!
  !! \details The namelist is supposed to holds
  !! - the dimension name to make the difference with other covariance parameters in the namelist file
  !! - the number of significant singular values
  !! - the smoothing parameter
  !! - and the decorrelation length
  !! By default, dName picks values in ["X", "Y", "Z"], but can be anything meaningful to the program writer
  !<
  subroutine load_ecp(ecp, fName, dName, n, delta)
    INTEGER, PARAMETER:: fid = 68
    type(exp_cov_param), intent(out) :: ecp
    CHARACTER(len=*), INTENT(IN) :: fName, dName
    INTEGER, INTENT(IN) :: n
    REAL(cp), INTENT(IN) :: delta
    !local var
    CHARACTER(len=ip_snl) :: dimName, ala_dimName
    INTEGER :: nsv
    REAL(cp):: smoothing, decorrelation
    namelist/NAM_EXPONENTIAL_COV_PARAM/&
      dimName,&
      nsv,&
      smoothing,&
      decorrelation

    open(fid,file=fName,form='FORMATTED',status='OLD')
    dimName = "{#}@"
    ala_dimName = trim(dName)
    call uppercase(ala_dimName)
    call debug(ala_dimName, "In load_ecp, looking for dimName = ", tag=dALLWAYS)
    do while(dimName.ne.ala_dimName)
      read(fid, NAM_EXPONENTIAL_COV_PARAM)!reading the block
      call uppercase(dimName)
      call debug(dimName, "  found dimName = ", tag=dALLWAYS)
    end do
    close(fid)

    ecp = exp_cov_param(n, nsv, delta, smoothing, decorrelation)
  end subroutine load_ecp

  !> @brief Computes the square exponential covariance matrix for a 1D uniform grid, using exp_cov_param derived type
  !! @param [in, out] C square matrix of covariance
  !! @param [in] xcp covariance parameters in the x dimension
  !!
  !! @see se_cov_matrix_scalar
  !>
  subroutine se_cov_matrix_type(C, xcp)
    real(dp), dimension(:,:), intent(in out) :: C
    type(exp_cov_param), intent(in) :: xcp

    call se_cov_matrix_scalar(C, xcp%dx, xcp%s, xcp%l)
  end subroutine se_cov_matrix_type

  !> @brief Computes the square exponential covariance matrix for a 1D uniform grid, using scalar parameters
  !! @param [in, out] C square matrix of covariance
  !! @param [in] dx discretization step
  !! @param [in] s smoothing parameter of the covariance function (standard deviation)
  !! @param [in] lx correlation length
  !!
  !! The covariance betwen the grid points i and j is given by
  !! \f$ C(i,j) = s^2\exp(\frac{-(i\Delta_x-j\Delta_x)^2}{2l_x^2}) \f$
  !>
  subroutine se_cov_matrix_scalar(C, dx, s, lx)
    real(dp), dimension(:,:), intent(in out) :: C
    real(dp), intent(in) :: dx, s, lx
    !local variables
    real(dp) :: x, s2, lx2
    integer :: ii, jj

    s2 = s*s
    lx2 = lx*lx

    !compute the upper triangular part
    do jj=2, size(C,2)
      do ii=1, jj-1
        x = (ii-jj)*dx
        C(ii, jj) = s2*exp(-x*x/(2*lx2))
      end do
    end do

    !The diagonal element are all identical
    do jj=1,size(C,2)
      C(jj, jj) = s2
    end do

    !computes the lower triangular part by symmetry
    do jj = 1, size(C,2)
      do ii = jj+1, size(C,1)
        C(ii, jj) = C(jj, ii)
      end do
    end do

    !I strongly believe that the next (commented line) was not necessary
    !C = dx*C
  end subroutine se_cov_matrix_scalar

  !> @brief Covariance function using exponential of the square
  !! @param [in] x1, x2 points between which the covariance is evaluated
  !! @param [in] s smoothing parameter of the covariance function
  !! @param [in] lx decorrelation length
  !!
  !<
  function se_cov_function(x1, x2, s, lx)result(f)
    real(dp), intent(in) :: x1, x2, s, lx
    !local variables
    real(dp) :: f, x

    x = x1-x2
    f = s*s*exp(-x*x/(2*lx*lx))
  end function se_cov_function

  !> @brief Covariance function using exponential of absolute value
  !! @param [in] x1, x2 points between which the covariance is evaluated
  !! @param [in] s smoothing parameter of the covariance function; this is the varaince
  !! @param [in] lx decorrelation length
  !!
  !<
  function ae_cov_function(x1, x2, s, lx)result(f)
    real(dp), intent(in) :: x1, x2, s, lx
    real(dp) :: f

    f = s*s*exp(-abs(x1-x2)/lx);
  end function ae_cov_function
end module exponential_cov

max_stars_repo_name: souopgui/da-toolbox

max_issues_repo_name: souopgui/da-toolbox

max_forks_repo_name: souopgui/da-toolbox


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

! Integer function to open a series of files whose names are stored
! in a file caled FSOPEN.FIL. The current program version can handle
! only three tokens:
! 1. The unit number,
! 2. The file name,
! 3. The format of the file (optional).
!
      INTEGER FUNCTION FSOPEN ()
      CHARACTER*80  LINE
      CHARACTER*40  NAME, OPTIONS
      INTEGER FSUNIT, STOI, GETS, POSIT
      LOGICAL LEXIST
!
! Check if the file with the filedef information exists.
!
      FSOPEN = 0
      INQUIRE (FILE = 'FSOPEN.FIL', EXIST = LEXIST)
      IF ( .NOT. LEXIST) THEN
         FSOPEN = 1
         GO TO 100
      END IF
!
! Read a record form the filedefinition file. Tokenize it.
!
      OPEN (99, FILE = 'FSOPEN.FIL')
5     CONTINUE
      POSIT = 1
      READ (99, 10, END=100) LINE
10    FORMAT (A80)
!
! Get unit number.
!
      IF (STOI (LINE,POSIT,80,FSUNIT) .NE. 0) THEN
         FSOPEN = 2
         GO TO 100
      END IF
!
! Get file name.
!
      IF (GETS (LINE,POSIT,80,40,NAME,LENNAM) .NE. 0) THEN
         FSOPEN = 3
         GO TO 100
      END IF
!
! Get format (optional)
! Open the file.
!
      IF (GETS (LINE,POSIT,80,40,OPTIONS,LENOP) .EQ. 0) THEN
         OPEN (UNIT=FSUNIT, FILE= NAME, FORM=OPTIONS)
      ELSE
         OPEN (UNIT=FSUNIT, FILE= NAME)
      END IF
      GO TO 5
100   CONTINUE
      CLOSE (99, STATUS = 'DELETE')
      RETURN
      END

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      SUBROUTINE PLA32
C
C     THIS ROUTINE READS THE ESTNLS DATA BLOCK CREATED IN SUBROUTINE
C     PLA31, AND CALLS THE PROPER ELEMENT ROUTINE TO COMPUTE ELEMENT
C     STRESSES.
C     ELEMENT STRESS INFORMATION IS STORED BY THE ELEMENT ROUTINE IN
C     /STROUT/.  THE ELEMENT ROUTINE ALSO UPDATES THE EST ENTRY WHICH
C     HAS BEEN COMMUNICATED TO IT VIA /PLA32E/.  NOTE THAT THIS UPDATED
C     EST ENTRY DOES NOT CONTAIN DISPLACEMENT VECTOR INFORMATION.
C
      INTEGER         BUFSZ,BUFR1,BUFR2,BUFR3,CSTM,DIT,ESTNLS,CASECC,
     1                ONLES,ESTNL1,EOR,CLSRW,FILE,IZ(1),IESTBK(100),
     9                ESTWDS(40),ELTYPE,OUTRW,PLACNT,PLSETN,PLANOS(2),
     2                OSTRT(7),ESTT(7),SETNO
      DIMENSION       NAME(2),NSTWDS(40),NWDSP2(40),P(4),IP(4),DUM2(2),
     1                TUBSAV(20),ICHAR(9),ITITLE(3),IY(30)
      COMMON /BLANK / PLACNT,PLSETN
      COMMON /SYSTEM/ BUFSZ
      COMMON /CONDAS/ PI,TWOPI,RADEG,DEGRA,S4PISQ
      COMMON /ZZZZZZ/ Z(1)
      COMMON /PLA32C/ GAMMA,GAMMAS,IPASS
      COMMON /PLA32E/ ESTBK(100)
C
C     SCRATCH BLOCK USED BY ELEMENT ROUTINES (325 SINGLE PRECISION
C     CELLS)  AND OUTPUT BLOCK FOR ELEMENT STRESSES
C
      COMMON /PLA32S/ XXXXXX(325)
      COMMON /SOUT  / YYYYYY(30)
      EQUIVALENCE     (Z(1),IZ(1)),(ESTBK(1),IESTBK(1)),(P(1),IP(1)),
     1                (YYYYYY(1),IY(1))
      DATA    NAME  / 4HPLA3, 4H2          /
      DATA    ITITLE/ 4HLOAD, 4H FAC,4HTOR /
      DATA    CSTM  , MPT,DIT,ESTNLS,CASECC/ 101,102,103,301,106 /
      DATA    ONLES , ESTNL1  / 201,202    /
      DATA    INRW  , OUTRW,EOR,NEOR,CLSRW / 0,1,1,0,1 /
      DATA    PLANOS/ 1103,11 /
C
C    1        ROD       BEAM      TUBE      SHEAR     TWIST
C    2        TRIA1     TRBSC     TRPLT     TRMEM     CONROD
C    3        ELAS1     ELAS2     ELAS3     ELAS4     QDPLT
C    4        QDMEM     TRIA2     QUAD2     QUAD1     DAMP1
C    5        DAMP2     DAMP3     DAMP4     VISC      MASS1
C    6        MASS2     MASS3     MASS4     CONM1     CONM2
C    7        PLOTEL    REACT     QUAD3     BAR       CONE
C    8        TRIARG    TRAPRG    TORDRG    CORE      CAP
C
      DATA    ESTWDS/
     1           21,        0,       20,        0,        0,
     2           38,        0,        0,       27,       21,
     3            0,        0,        0,        0,        0,
     4           32,       32,       37,       43,        0,
     5            0,        0,        0,        0,        0,
     6            0,        0,        0,        0,        0,
     7            0,        0,        0,       50,        0,
     8            0,        0,        0,        0,        0 /
      DATA    NSTWDS/
     1            5,        0,        5,        0,        0,
     2           17,        0,        0,        8,        5,
     3            0,        0,        0,        0,        0,
     4            8,       17,       17,       17,        0,
     5            0,        0,        0,        0,        0,
     6            0,        0,        0,        0,        0,
     7            0,        0,        0,       16,        0,
     8            0,        0,        0,        0,        0 /
      DATA    NWDSP2/
     1           33,        0,       32,        0,        0,
     2           56,        0,        0,       36,       33,
     3            0,        0,        0,        0,        0,
     4           44,       50,       61,       67,        0,
     5            0,        0,        0,        0,        0,
     6            0,        0,        0,        0,        0,
     7            0,        0,        0,       62,        0,
     8            0,        0,        0,        0,        0 /
C
C     DEFINE POSITION IN CASECC RECORD OF DESTINATION (PRINTER, PUNCH,
C     ETC.) OF ELEMENT STRESSES.
C
      DATA          IDEST /24/
C
C
C     DETERMINE SIZE OF CORE, DEFINE BUFFERS AND INITIALIZE CORE POINTER
C     AND COUNTERS
C
      IZMAX = KORSZ(Z)
      BUFR1 = IZMAX - BUFSZ
      BUFR2 = BUFR1 - BUFSZ
      BUFR3 = BUFR2 - BUFSZ
      LEFT  = BUFR3 - 1
      IPASS = PLACNT- 1
      ICSTM = 0
      NCSTM = 0
      DO 5 I = 1,7
      OSTRT(I) = 0
    5 ESTT(I)  = 0
C
C     ATTEMPT TO READ CSTM INTO CORE
C
      FILE = CSTM
      CALL OPEN (*20,CSTM,Z(BUFR1),INRW)
      CALL FWDREC (*9020,CSTM)
      CALL READ (*9020,*10,CSTM,Z(ICSTM+1),LEFT,EOR,NCSTM)
      CALL MESAGE (-8,0,NAME)
   10 LEFT = LEFT - NCSTM
      CALL CLOSE (CSTM,CLSRW)
      CALL PRETRS (Z(ICSTM+1),NCSTM)
   20 IMAT = NCSTM
C
C     COMPUTE GAMMA AND GAMMAS FROM THE PROPER PLFACT CARD
C
      FILE = MPT
      CALL PRELOC (*9010,Z(BUFR1-3),MPT)
      CALL LOCATE (*9040,Z(BUFR1-3),PLANOS,IFLAG)
   30 CALL READ (*9020,*9030,MPT,SETNO,1,NEOR,IFLAG)
      IF (SETNO .EQ. PLSETN) GO TO 40
   35 CALL READ (*9020,*9030,MPT,NN,1,NEOR,IFLAG)
      IF (NN .EQ. -1) GO TO 30
      GO TO 35
   40 IF (PLACNT .LE. 4) GO TO 45
      CALL READ (*9020,*9030,MPT,0,-(PLACNT-4),NEOR,IFLAG)
   45 NWDSRD = 4
      IF (PLACNT .LT. 4) NWDSRD = PLACNT
      CALL READ (*9020,*9030,MPT,P,NWDSRD,NEOR,IFLAG)
      IF (IP(NWDSRD) .NE. -1) GO TO 48
      IF (PLACNT-3) 42,43,44
  42  GAMMAS = 1.0
      GO TO 46
  43  GAMMAS = (P(2) - P(1))/P(1)
      GO TO 46
  44  GAMMAS = (P(3) - P(2))/(P(2) - P(1))
  46  GAMMA = 1.0
      GO TO 65
  48  A = P(2) - P(1)
      IF (PLACNT-3) 50,55,60
   50 GAMMAS = 0.0
      GAMMA  = A/P(1)
      GO TO 65
  55  GAMMAS = A/P(1)
      GAMMA = (P(3) - P(2))/A
      GO TO 65
   60 WORD = P(3) - P(2)
      GAMMAS = WORD/A
      GAMMA  = (P(4) - P(3))/WORD
   65 CALL CLOSE (MPT,CLSRW)
C
C     READ MPT AND DIT FILES.  NOTE MINUS SIGN ON DIT TO TRIGGER PLA
C     FLAG.
C
      CALL PREMAT (IZ(IMAT+1),Z(IMAT+1),Z(BUFR1-3),LEFT,MUSED,MPT,-DIT)
      LEFT = LEFT  - MUSED
      ICC  = NCSTM + MUSED
C
C     READ CASECC INTO OPEN CORE
C
      FILE = CASECC
      CALL OPEN (*9010,CASECC,Z(BUFR1),INRW)
      CALL FWDREC (*9020,CASECC)
      CALL READ (*9020,*68,CASECC,Z(ICC+1),LEFT,EOR,NCC)
      CALL MESAGE (-8,0,NAME)
   68 LEFT = LEFT - NCC
      CALL CLOSE (CASECC,CLSRW)
C
C OPEN INPUT FILE
C
      FILE = ESTNLS
      CALL OPEN (*9010,ESTNLS,Z(BUFR1),INRW)
      CALL FWDREC (*9020,ESTNLS)
C
C     OPEN THE ELEMENT STRESS FILE FOR OUTPUT AND BUILD HEADER WHICH IS
C     NON-CHANGING.
C
      FILE = ONLES
      CALL OPEN (*9010,ONLES,Z(BUFR2),OUTRW)
      CALL FNAME (ONLES,DUM2)
      CALL WRITE (ONLES,DUM2,2,EOR)
C
C     THE FOLLOWING INDICES HAVE TO CHANGE  WHEN THERE ARE CHANGES IN
C     THE FORMAT OF THE CASECC DATA BLOCK
C
      IONLES = ICC + NCC
      IZ(IONLES+1) = IZ(ICC+18) + 100
      IZ(IONLES+2) = 5
      IZ(IONLES+4) = IZ(ICC+1)
      IZ(IONLES+5) = IZ(ICC+4)
      IZ(IONLES+6) = 0
      IZ(IONLES+7) = 0
      IZ(IONLES+8) = 0
      IZ(IONLES+9) = 0
      ILOW  = IONLES + 51
      IHIGH = IONLES + 146
      LEFT  = LEFT - 146
      IF (LEFT .LT. 0) CALL MESAGE (-8,0,NAME)
      J = ICC + 38
      DO 70 I = ILOW,IHIGH
      J = J + 1
   70 IZ(I) = IZ(J)
C
C     STORE LOAD FACTOR AND INTEGER IN LABEL PORTION OF OUTPUT
C
      IZ(IONLES+135) = ITITLE(1)
      IZ(IONLES+136) = ITITLE(2)
      IZ(IONLES+137) = ITITLE(3)
      III = PLACNT - 1
      CALL INT2AL (III,IZ(IONLES+138),ICHAR)
C
C     DEFINE DESTINATION OF OUTPUT
C
      I = ICC + IDEST
      JDEST = IZ(I)
C
C     OPEN THE ESTNL1 FILE FOR OUTPUT.
C
      FILE = ESTNL1
      CALL OPEN (*9010,ESTNL1,Z(BUFR3),OUTRW)
      CALL FNAME (ESTNL1,DUM2)
      CALL WRITE (ESTNL1,DUM2,2,EOR)
      FILE = ESTNLS
C
C     READ ELEMENT TYPE
C
   80 CALL READ (*220,*9030,ESTNLS,ELTYPE,1,NEOR,IFLAG)
C
C     FILL IN REMAINDER OF ID RECORD FOR THE ONLES FILE
C
      IZ(IONLES+3)  = ELTYPE
      IZ(IONLES+10) = NSTWDS(ELTYPE)
      IF (NSTWDS(ELTYPE) .LE. 0) CALL MESAGE (-30,91,ELTYPE)
C
C     WRITE ID RECORD FOR ONLES FILE
C
      CALL WRITE (ONLES,IZ(IONLES+1),146,EOR)
      CALL WRITE (ESTNL1,ELTYPE,1,NEOR)
C
C     READ AN ENTRY FROM THE APPENDED ESTNL FILE AND CALL THE PROPER
C     ROUTINE
C
   90 CALL READ (*9020,*210,ESTNLS,ESTBK,NWDSP2(ELTYPE),NEOR,IFLAG)
C
C               1,ROD    2,BEAM    3,TUBE   4,SHEAR   5,TWIST
      GO TO (     110,      999,      120,      999,      999,
C             6,TRIA1   7,TRBSC   8,TRPLT   9,TRMEM 10,CONROD
     1            130,      999,      999,      140,      110,
C            11,ELAS1  12,ELAS2  13,ELAS3  14,ELAS4  15,QDPLT
     2            999,      999,      999,      999,      999,
C            16,QDMEM  17,TRIA2  18,QUAD2  19,QUAD1  20,DAMP1
     3            150,      160,      170,      180,      999,
C            21,DAMP2  22,DAMP3  23,DAMP4   24,VISC  25,MASS1
     4            999,      999,      999,      999,      999,
C            26,MASS2  27,MASS3  28,MASS4  29,CONM1  30,CONM2
     5            999,      999,      999,      999,      999,
C           31,PLOTEL  32,REACT  33,QUAD3    34,BAR   35,CONE
     6            999,      999,      999,      190,      999,
C           36,TRIARG 37,TRAPRG 38,TORDRG   39,CORE?   40,CAP?
     7            999,      999,      999,      999,      999), ELTYPE
C
C     ROD, CONROD
C
  110 CALL PSROD
C
C     IF ELEMENT IS A TUBE, RESTORE SAVED EST ENTRY AND STORE UPDATED
C     STRESS VARIABLES IN PROPER SLOTS.
C
      IF (ELTYPE .NE. 3) GO TO 200
      DO 115 I = 1,16
  115 ESTBK(I) = TUBSAV(I)
      ESTBK(17) = ESTBK(18)
      ESTBK(18) = ESTBK(19)
      ESTBK(19) = ESTBK(20)
      ESTBK(20) = ESTBK(21)
      GO TO 200
C
C
C     TUBE - REARRANGE ESTBK FOR THE TUBE SO THAT IT IS IDENTICAL TO THE
C            ONE FOR THE ROD
C
C     SAVE THE EST ENTRY FOR THE TUBE EXCEPT THE 4 WORDS WHICH WILL BE
C     UPDATED BY THE THE ROD ROUTINE AND THE DISPLACEMENT VECTORS
C
  120 DO 125 I = 1,16
  125 TUBSAV(I) = ESTBK(I)
C
C     COMPUTE AREA, TORSIONAL INERTIA TERM AND STRESS COEFFICIENT
C
      D  = ESTBK(5)
      T  = ESTBK(6)
      DMT= D - T
      A  = DMT*T* PI
      FJ = .25*A*(DMT**2 + T**2)
      C  = D/2.0
C
C     MOVE THE END OF THE ESTBK ARRAY DOWN ONE SLOT SO THAT ENTRIES 7
C     THRU 32 WILL BE MOVED TO POSITIONS 8 THRU 33.
C
      M = 33
      DO 127 I = 1,26
      ESTBK(M) = ESTBK(M-1)
  127 M = M - 1
      ESTBK(5) = A
      ESTBK(6) = FJ
      ESTBK(7) = C
      GO TO 110
C
C     TRIA1
C
  130 CALL PSTRI1
      GO TO 200
C
C     TRMEM
C
  140 CALL PSTRM
      GO TO 200
C
C     QDMEM
C
  150 CALL PSQDM
      GO TO 200
C
C     TRIA2
C
  160 CALL PSTRI2
      GO TO 200
C
C     QUAD2
C
  170 CALL PSQAD2
      GO TO 200
C
C     QUAD1
C
  180 CALL PSQAD1
      GO TO 200
C
C     BAR
C
  190 CALL PSBAR
      GO TO 200
C
C     ALTER ELEMENT IDENTIFICATION FROM EXTERNAL (USER) IDENTIFICATION
C     TO INTERNAL ID., AND WRITE OUTPUT FILES.
C
  200 IY(1) = 10*IY(1) + JDEST
      CALL WRITE (ONLES,IY,NSTWDS(ELTYPE),NEOR)
      CALL WRITE (ESTNL1,ESTBK,ESTWDS(ELTYPE),NEOR)
      OSTRT(2) = OSTRT(2) + 1
      ESTT(2)  = ESTT(2)  + 1
      GO TO 90
C
C     WRITE EORS
C
  210 CALL WRITE (ONLES,0,0,EOR)
      CALL WRITE (ESTNL1,0,0,EOR)
      GO TO 80
C
C     CLOSE FILES AND WRITE TRAILERS
C
  220 CALL CLOSE (ONLES,CLSRW)
      CALL CLOSE (ESTNL1,CLSRW)
      CALL CLOSE (ESTNLS,CLSRW)
      OSTRT(1) = ONLES
      ESTT(1)  = ESTNL1
      CALL WRTTRL (OSTRT)
      CALL WRTTRL (ESTT)
      RETURN
C
C     FATAL ERRORS
C
  999 CALL MESAGE (-30,92,ELTYPE)
 9010 J = -1
      GO TO 9050
 9020 J = -2
      GO TO 9050
 9030 J = -3
      GO TO 9050
 9040 J = -5
 9050 CALL MESAGE (J,FILE,NAME)
      RETURN
      END