text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 04.02.04.01/10                        *
C  *    TEST TITLE : Annotation style                      *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

      COMMON /DIALOG/ DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT,
     2                DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS
      INTEGER         DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT
      REAL            DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)

C aspect identifier
      INTEGER    PLN,      PLWSC,    PPLCI,    PMK,      PMKSC
      PARAMETER (PLN  = 0, PLWSC= 1, PPLCI= 2, PMK  = 3, PMKSC= 4)

C linetype
      INTEGER    PLSOLI,     PLDASH,     PLDOT,     PLDASD
      PARAMETER (PLSOLI = 1, PLDASH = 2, PLDOT = 3, PLDASD = 4)

C marker type
      INTEGER    PPOINT,   PPLUS,   PAST,   POMARK,   PXMARK
      PARAMETER (PPOINT=1, PPLUS=2, PAST=3, POMARK=4, PXMARK=5)

C text alignment horizontal
      INTEGER    PAHNOR,     PALEFT,     PACENT,     PARITE
      PARAMETER (PAHNOR = 0, PALEFT = 1, PACENT = 2, PARITE = 3)

C text alignment vertical
      INTEGER    PAVNOR,     PATOP,      PACAP,     PAHALF,
     1           PABASE,     PABOTT
      PARAMETER (PAVNOR = 0, PATOP  = 1, PACAP = 2, PAHALF = 3,
     1           PABASE = 4, PABOTT = 5)

C text precision
      INTEGER    PSTRP,     PCHARP,     PSTRKP
      PARAMETER (PSTRP = 0, PCHARP = 1, PSTRKP = 2)

C annotation style
      INTEGER    PUNCON,     PLDLN
      PARAMETER (PUNCON = 1, PLDLN = 2)

      INTEGER    IMPMAX, NUMLIN, PERM(20), NGSTYL, LNID, TRUANS
      PARAMETER (IMPMAX = 300)
      INTEGER    NUMAS,MINANS,MAXANS, STDDX,REGDX,IMPDX, ANSTYL
      INTEGER    ADX, THISAS, LAVRAS(20), LAVSAS(2), LAVIAS(IMPMAX)
      INTEGER    LUNSUP(10), UNSUDX, NUMLT,LNTYP,LNCOL, RNDINT
      INTEGER    PICSTR, TXCI, IX
      INTEGER    IDUM1, IDUM2

      REAL       NPCX,NPCY,NPCPWC, YINCR,YLOC, ANOFFX,ANOFFY, LWSCF
      REAL       XA(20),YA(20), XB(20),YB(20), XO(20),YO(20), NOMLW
      REAL       RDUM1, RDUM2

      CHARACTER  MSG*1600

      LOGICAL    ANSAVL(0:2)

      CALL INITGL ('04.02.04.01/10')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)
C set-up of workstation and dialogue area
      PICSTR = 101
      TXCI = 1
      CALL SETDLG (PICSTR, 801,TXCI)
      CALL WCNPC (0.0, 0.0, NPCX,NPCY, NPCPWC)

      CALL POPST (PICSTR)
C  by convention, view #1 is for picture
      CALL PSVWI (1)
C  use individual attributes
      CALL SETASF (PINDIV)
C  defaults
      CALL PSTXFN (1)
      CALL PSTXPR (PSTRKP)
      CALL PSCHXP (0.7)
      CALL PSCHSP (0.1)
C  geometric attributes for annotation text
      CALL PSATCH (0.03 * NPCPWC)
      CALL PSATAL (PARITE, PAHALF)
C  geometric attributes for text
      CALL PSCHH  (0.03)
      CALL PSTXAL (PARITE, PAHALF)
      CALL PEXST (102)
      CALL PCLST

      CALL POPST (102)

C  use <inquire annotation facilities> to determine:
C    numas  = number of available annotation styles
C    lavsas = list of available standard annotation styles
C    lavras = list of available registered (non-mandatory) annotation styles
C    lavias = list of available implementor-defined annotation styles

      CALL PQANF (SPECWT, 1, ERRIND, NUMAS, MINANS, IDUM2, RDUM1,RDUM2)
      CALL CHKINQ ('pqanf', ERRIND)
      MAXANS = MINANS
      STDDX = 0
      REGDX = 0
      IMPDX = 0
      ANSAVL(0) = .FALSE.
      ANSAVL(1) = .FALSE.
      ANSAVL(2) = .FALSE.

      DO 150 ADX = 1,NUMAS
         CALL PQANF (SPECWT,ADX, ERRIND, IDUM1, THISAS,
     1               IDUM2,RDUM1,RDUM2)
         CALL CHKINQ ('pqanf', ERRIND)
         MINANS = MIN (MINANS, THISAS)
         MAXANS = MAX (MAXANS, THISAS)
         IF  (THISAS .GT. 2) THEN
            REGDX = REGDX+1
            LAVRAS (REGDX) = THISAS
         ELSEIF (THISAS .GT. 0) THEN
            STDDX = STDDX +1
            LAVSAS (STDDX) = THISAS
            ANSAVL(THISAS) = .TRUE.
         ELSE
            IF (THISAS .EQ. 0) ANSAVL(0) = .TRUE.
            IF (IMPDX .LT. IMPMAX) THEN
               IMPDX = IMPDX+1
               LAVIAS (IMPDX) = THISAS
            ENDIF
         ENDIF
150   CONTINUE

      CALL SETMSG ('77', 'The list of available annotation styles ' //
     1             'returned by <inquire annotation facilities> '   //
     2             'should include styles number 1 and 2.')
      CALL IFPF (ANSAVL(1) .AND. ANSAVL(2))

      CALL DLINE ('Throughout this program, when an annotation '  //
     1            'text primitive is displayed, its reference '   //
     2            'point is marked with an asterisk polymarker, ' //
     3            'and its annotation point with a circle '       //
     4            'polymarker (enter anything to continue).', MSG)

      IF (.NOT. ANSAVL(PUNCON)) THEN
         CALL INMSG ('Skipping tests for annotation style #1: not ' //
     1               'available.')
         GOTO 300
      ENDIF

C  *** *** ***   style #1 (unconnected)   *** *** ***

C  set up perm to randomize order
      NUMLIN = 6
      CALL RNPERM (NUMLIN, PERM)

C  ngstyl = random integer from 1 to 3
      NGSTYL = RNDINT(1,3)

C  display interspersed 3 actual and 3 simulated examples of style #1,
C    with actual #ngstyl done incorrectly (use style #2)
      YINCR = 1.0 / (NUMLIN+1)
      YLOC  = 1 - YINCR
      CALL NUMLAB (NUMLIN, 0.2, YLOC, YINCR)
C  annotation offsets in WC
      ANOFFX = -0.2
      ANOFFY =  0.3 * YINCR
      DO 200 IX = 1,NUMLIN
C  mark reference and annotation point
         XO(IX) = 0.5
         YO(IX) = YLOC
         XA(IX) = XO(IX) - ANOFFX
         YA(IX) = YO(IX) - ANOFFY

         LNID = PERM(IX)
         IF (LNID .GT. 3) THEN
C  do simulation
            CALL PTX (XO(IX), YO(IX), 'Primitive')
         ELSE
C  do actual
            IF (LNID .EQ. NGSTYL) THEN
               ANSTYL = PLDLN
               TRUANS = IX
            ELSE
               ANSTYL = PUNCON
            ENDIF
            CALL PSANS (ANSTYL)
            CALL PATR (XA(IX),YA(IX), ANOFFX*NPCPWC, ANOFFY*NPCPWC,
     1                 'Primitive')
         ENDIF
         YLOC = YLOC-YINCR
200   CONTINUE
      CALL PSMK (PAST)
      CALL PPM  (NUMLIN, XA,YA)
      CALL PSMK (POMARK)
      CALL PPM  (NUMLIN, XO,YO)

      CALL SETMSG ('1 2 5 75 76 78 80', 'The mandatory annotation '   //
     1             'style number 1 should leave the reference point ' //
     2             'and the annotation point of the annotation text ' //
     3             'primitive unconnected.')

      CALL DCHPFV ('ANNOTATION STYLE 1: Which annotation style is ' //
     1             'different?', NUMLIN, TRUANS)
      CALL PEMST (102)

      IF (.NOT. ANSAVL(PLDLN)) GOTO 300
C  lunsup = list of (mostly) unsupported annotation styles, contains:
C    a negative unsupported type (< minimum in lavias)
C    a positive unsupported type (> maximum in lavras)
C    1,
C    2, if supported (the only non-#1 style)
C    0, if not in lavias
      LUNSUP(1) = MINANS-10
      LUNSUP(2) = MAXANS+1
      LUNSUP(3) = PUNCON
      LUNSUP(4) = PLDLN
      IF (ANSAVL(0)) THEN
         UNSUDX = 4
      ELSE
         UNSUDX = 5
         LUNSUP(5) = 0
      ENDIF

C  Display in random order annotation styles from lunsup
C  set up perm to randomize order
      CALL RNPERM (UNSUDX, PERM)

      YINCR = 1.0 / (UNSUDX+1)
      YLOC  = 1 - YINCR
      CALL NUMLAB (UNSUDX, 0.2, YLOC, YINCR)
      DO 220 IX = 1,UNSUDX
         ANSTYL = LUNSUP(PERM(IX))
         IF (ANSTYL .EQ. 2) TRUANS = IX
C  mark reference and annotation point
         XO(IX) = 0.5
         YO(IX) = YLOC
         XA(IX) = XO(IX) - ANOFFX
         YA(IX) = YO(IX) - ANOFFY

         CALL PSANS (ANSTYL)
         CALL PATR (XA(IX),YA(IX), ANOFFX*NPCPWC, ANOFFY*NPCPWC,
     1              'Primitive')
         YLOC = YLOC-YINCR
220   CONTINUE
      CALL PSMK (PAST)
      CALL PPM  (UNSUDX, XA,YA)
      CALL PSMK (POMARK)
      CALL PPM  (UNSUDX, XO,YO)

      CALL SETMSG ('5 75 76 78 79 80', 'Unavailable annotation '     //
     1             'styles should be displayed as annotation style ' //
     2             'number 1.')

      CALL DCHPFV ('UNSUPPORTED ANNOTATION STYLE: Which annotation ' //
     1             'style is different?', UNSUDX, TRUANS)
      CALL PEMST (102)

C  end_style_1:
300   CONTINUE

C  *** *** ***   style #2 (connected)   *** *** ***

      IF (.NOT. ANSAVL(PLDLN)) THEN
         CALL INMSG ('Skipping tests for annotation style #2: not ' //
     1               'available.')
         GOTO 400
      ENDIF

C  set up perm to randomize order
      NUMLIN = 6
      CALL RNPERM (NUMLIN, PERM)

C  ngstyl = random integer from 1 to 3
      NGSTYL = RNDINT(1,3)

C  display interspersed 3 actual and 3 simulated examples of style #2,
C    with actual #ngstyl done incorrectly (use style #1)
      YINCR = 1.0 / (NUMLIN+1)
      YLOC  = 1 - YINCR
      CALL NUMLAB (NUMLIN, 0.2, YLOC, YINCR)

C  set up exciting polyline attributes
      CALL PQPLF (SPECWT, 0, ERRIND, NUMLT, IDUM1,
     1            IDUM2,NOMLW,RDUM1,RDUM2, IDUM2)
      CALL CHKINQ ('pqplf', ERRIND)
      CALL PQPLF (SPECWT, RNDINT(1, ABS(NUMLT)), ERRIND, NUMLT,LNTYP,
     1            IDUM2,NOMLW,RDUM1,RDUM2, IDUM2)
      CALL CHKINQ ('pqplf', ERRIND)
      LWSCF = .02 / (NOMLW * WCPDC)
      LNCOL = 2
C  set linetype/width thru bundle, color individually
      CALL PSLN   (LNTYP+1)
      CALL PSLWSC (2*LWSCF)
      CALL PSPLCI (LNCOL)
      CALL PSPLI  (2)
      CALL PSPLR  (WKID, 2, LNTYP,LWSCF,LNCOL-1)
      CALL PSIASF (PLN,   PBUNDL)
      CALL PSIASF (PLWSC, PBUNDL)

      DO 310 IX = 1,NUMLIN
C  mark reference and annotation point
         XO(IX) = 0.5
         YO(IX) = YLOC
         XA(IX) = XO(IX) - ANOFFX
         YA(IX) = YO(IX) - ANOFFY

         LNID = PERM(IX)
         IF (LNID .GT. 3) THEN
C  do simulation
            CALL PTX (XO(IX),YO(IX), 'Primitive')
            XB(1) = XA(IX)
            YB(1) = YA(IX)
            XB(2) = XO(IX)
            YB(2) = YO(IX)
            CALL PPL (2, XB,YB)
         ELSE
C  do actual
            IF (LNID .EQ. NGSTYL) THEN
               ANSTYL = PUNCON
               TRUANS = IX
            ELSE
               ANSTYL = PLDLN
            ENDIF
            CALL PSANS (ANSTYL)
            CALL PATR (XA(IX),YA(IX), ANOFFX*NPCPWC, ANOFFY*NPCPWC,
     1                 'Primitive')
         ENDIF
         YLOC = YLOC-YINCR
310   CONTINUE
      CALL PSMK (PAST)
      CALL PPM  (NUMLIN, XA,YA)
      CALL PSMK (POMARK)
      CALL PPM  (NUMLIN, XO,YO)

      CALL SETMSG ('1 2 5 75 76 78 80', 'The mandatory annotation ' //
     1             'style number 2 should connect the reference '   //
     2             'point and annotation point with a polyline, '   //
     3             'using the current attributes.')

      CALL DCHPFV ('ANNOTATION STYLE 2: Which annotation style is ' //
     1             'different?', NUMLIN, TRUANS)
      CALL PEMST (102)

C  end_style_2:
400   CONTINUE

C  *** *** ***   registered styles   *** *** ***

C  regdx = size of lavras
      IF (REGDX .LT. 1) THEN
         CALL INMSG ('No registered annotation styles are supported.')
         GOTO 800
      ENDIF

C display and label at most 8 annotation styles in lavras
      CALL DSANST (REGDX, LAVRAS)

      CALL SETMSG ('5 75 76 78 81', 'The registered non-mandatory '   //
     1             'annotation styles (> 2) should agree with their ' //
     2             'registered description.')

      CALL DYNPF ('NON-MANDATORY REGISTERED ANNOTATION STYLES: Is '  //
     1            'each annotation style depicted according to its ' //
     2            'numeric identifier''s specification in the ISO '  //
     3            'register?', 'Y')
      CALL PEMST (102)

C  neg_type:
800   CONTINUE

C  *** *** ***   implementor-defined styles   *** *** ***

C  impdx = size of lavias
      IF (IMPDX .LT. 1) THEN
         CALL INMSG ('No implementor defined annotation styles ' //
     1               'are supported.')
         GOTO 666
      ENDIF

C display and label at most 8 annotation styles in lavias
      CALL DSANST (IMPDX, LAVIAS)

      CALL SETMSG ('5 75 76 78 82', 'Implementor-defined annotation ' //
     1             'styles (< 1) should agree with the descriptions ' //
     2             'supplied by the implementor.')

      CALL DYNPF ('IMPLEMENTOR DEFINED ANNOTATION STYLES: Is each ' //
     1            'annotation style depicted according to its '     //
     2            'numeric identifier''s specification in the '     //
     3            'implementor documentation?', 'Y')
      CALL PEMST (102)

      IF (IMPDX .LE. 8) GOTO 666

      CALL SETMSG ('82', 'All reported implementor-defined '       //
     1             'annotation styles available for a given '      //
     2             'workstation type should be documented by the ' //
     3             'implementor.')

      CALL SRTIAR (IMPDX, LAVIAS)
      WRITE (MSG, '(500I5)') (LAVIAS(IX), IX=1,IMPDX)

      CALL DYNPF ('COMPLETE DOCUMENTATION FOR IMPLEMENTOR '        //
     1            'ANNOTATION STYLES: Are all these reportedly '   //
     2            'available annotation styles documented by the ' //
     3            'implementor?' // MSG, 'Y')

C  end_annotation_style:
666   CONTINUE
C  wrap it up.
      CALL ENDIT
      END


text:
      SUBROUTINE RAMBO(N,ET,XM,P,WT)
C------------------------------------------------------
C
C                       RAMBO
C
C    RA(NDOM)  M(OMENTA)  B(EAUTIFULLY)  O(RGANIZED)
C
C    A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR
C    AUTHORS:  S.D. ELLIS,  R. KLEISS,  W.J. STIRLING
C    THIS IS VERSION 1.0 -  WRITTEN BY R. KLEISS
C
C    MODIFIED SLIGHTLY BY IAN HINCHLIFFE
C
C    N  = NUMBER OF PARTICLES (>1, IN THIS VERSION <101)
C    ET = TOTAL CENTRE-OF-MASS ENERGY
C    XM = PARTICLE MASSES ( DIM=100 )
C    P  = PARTICLE MOMENTA ( DIM=(4,100) )
C    WT = WEIGHT OF THE EVENT
C
C------------------------------------------------------
      IMPLICIT NONE
      REAL*8 Z,XM,P,Q,R,B,XM2,E,V,WT,ACC,PO2LOG,XMT,ET,RN,XMAX
      REAL*8 C,S,F,RANDOM1,RMAS,A,X,G,BQ,ACCU,WT2,WT3,WTM,P2
      REAL*8 F0,G0,X2,PI
      INTEGER*4 IWARN,IBEGIN,ITMAX,K,N,NM,ITER,I
      DIMENSION XM(100),P(4,100),Q(4,100),R(4),
     .   B(3),P2(100),XM2(100),E(100),V(100),IWARN(5),Z(100)
      DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/
      DATA PO2LOG/4.515827D-1/,PI/3.14159/
C PO2LOG IS LOG(PI/2)
C
C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT
      IF(IBEGIN.EQ.0) THEN
        IBEGIN=1
        Z(2)=PO2LOG
        DO 33 K=3,100
          Z(K)=Z(K-1)+PO2LOG-2.*LOG(FLOAT(K-2))
   33   CONTINUE
        DO 35 K=3,100
          Z(K)=(Z(K)-LOG(FLOAT(K-1)))
   35   CONTINUE
C CHECK ON THE NUMBER OF PARTICLES
      ELSE IF(N.LT.2.OR.N.GT.100) THEN
        PRINT 1001,N
        STOP
      END IF
C
C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES
      XMT=0.D0
      NM=0
      DO 17 I=1,N
        IF(XM(I).NE.0.) NM=NM+1
        XMT=XMT+ABS(XM(I))
   17 CONTINUE
      IF(XMT.GT.ET) THEN
        PRINT 1002,XMT,ET
        STOP
      END IF
C
C THE PARAMETER VALUES ARE NOW ACCEPTED
C
C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE
c      write(*,*) 'rambo--start of do 27 loop'
      DO 27 I=1,N
c        write(*,*) 'rambo--calling rn 1'
        C=2.D0*RN(1)-1.
c      write(*,*) 'rambo--rn 1', (C+1.d0)/2.d0
        S=SQRT(1.-C*C)
c        write(*,*) 'rambo--calling rn'
        F=2.*PI*RN(2)
c        write(*,*) 'rambo--rn(2)', F/(2.*PI)
        RANDOM1=RN(3)*RN(4)
        IF(RANDOM1.EQ.0) THEN
          RANDOM1=RN(3)*RN(4)
        END IF
        Q(4,I)=-LOG(RANDOM1)
        Q(3,I)=Q(4,I)*C
        Q(2,I)=Q(4,I)*S*COS(F)
        Q(1,I)=Q(4,I)*S*SIN(F)
   27 CONTINUE
C
C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION
      DO 93 I=1,4
        R(I)=0.D0
   93 CONTINUE
      DO 95 I=1,N
        DO 97 K=1,4
          R(K)=R(K)+Q(K,I)
   97   CONTINUE
   95 CONTINUE
      RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2)
      DO 98 K=1,3
        B(K)=-R(K)/RMAS
   98 CONTINUE
      G=R(4)/RMAS
      A=1.D0/(1.+G)
      X=ET/RMAS
C
C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S
      DO 43 I=1,N
        BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I)
        DO 45 K=1,3
          P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ))
   45   CONTINUE
        P(4,I)=X*(G*Q(4,I)+BQ)
   43 CONTINUE
C
C CALCULATE WEIGHT AND POSSIBLE WARNINGS
      WT=PO2LOG
      IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N)
      IF(WT.LT.-180.) THEN
        IF(IWARN(1).LE.5) PRINT 1004,WT
        IWARN(1)=IWARN(1)+1
      ELSE IF(WT.GT. 174.) THEN
        IF(IWARN(2).LE.5) PRINT 1005,WT
        IWARN(2)=IWARN(2)+1
C
C RETURN FOR WEIGHTED MASSLESS MOMENTA
      ELSE IF(NM.EQ.0) THEN
        WT=DEXP(WT)
        RETURN
      END IF
C
C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X
      XMAX=SQRT(1.-(XMT/ET)**2)
      DO 63 I=1,N
        XM2(I)=XM(I)**2
        P2(I)=P(4,I)**2
   63 CONTINUE
      ITER=0
      X=XMAX
      ACCU=ET*ACC
      DO WHILE (ITER.LE.ITMAX)
        F0=-ET
        G0=0.D0
        X2=X*X
        DO 65 I=1,N
          E(I)=SQRT(XM2(I)+X2*P2(I))
          F0=F0+E(I)
          G0=G0+P2(I)/E(I)
   65   CONTINUE
        IF(ABS(F0).GT.ACCU.AND.ITER.LE.ITMAX) THEN
          ITER=ITER+1
          X=X-F0/(X*G0)
        ELSE IF(ABS(F0).LE.ACCU) THEN
          ITER=ITMAX+1
        ELSE
          PRINT 1006,ITMAX
        END IF
      END DO
      DO 67 I=1,N
        V(I)=X*P(4,I)
        DO 69 K=1,3
          P(K,I)=X*P(K,I)
   69   CONTINUE
        P(4,I)=E(I)
   67 CONTINUE
C
C CALCULATE THE MASS-EFFECT WEIGHT FACTOR
      WT2=1.D0
      WT3=0.D0
      DO 73 I=1,N
        WT2=WT2*V(I)/E(I)
        WT3=WT3+V(I)**2/E(I)
   73 CONTINUE
      WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET)
C
C RETURN FOR  WEIGHTED MASSIVE MOMENTA
      WT=WT+WTM
      IF(WT.LT.-180.) THEN
        IF(IWARN(3).LE.5) PRINT 1004,WT
        IWARN(3)=IWARN(3)+1
      ELSE IF(WT.GT. 174.) THEN
        IF(IWARN(4).LE.5) PRINT 1005,WT
        IWARN(4)=IWARN(4)+1
      END IF
      WT=DEXP(WT)
      RETURN
C
 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED')
 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',
     .  ' SMALLER THAN TOTAL ENERGY =',D15.6)
 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW')
 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY  OVERFLOW')
 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE',
     .  ' DESIRED ACCURACY =',D15.6)
C******************************************************
      END
c---------------------------------------------------------------------------
c      double precision function rn(idum)
c      integer*4 idum
c      real dum,myran_sing
c      rn = dble(myran_sing)
c      return
c      end
c-------------------------------------------------------------------------

ext: f

max_stars_repo_name: LoganAMorrison/Hazma

max_issues_repo_name: LoganAMorrison/Hazma

max_forks_repo_name: LoganAMorrison/Hazma


text:
!===========================================================================
!  Program BIGSTICK
!
!  Based upon code REDSTICK by W. E. Ormand
!  BDEFORMEDLIB1.f90
!
!  optional subroutines for dealing with "deformed" systems, usually with
!  non-rotationally invariant single particle energies.

!===============================================
!
!  master routine for deformed (non-spherically-symmetric) single particle energies
!
! CALLED BY:
!    master_readin_tbmes
!
!  SUBROUTINES CALLED:
!    readdeformedSPEs
!
subroutine deformedSPEs

use interaction_deformed
use spstate
use haiku_info
use io

implicit none
integer it
logical foundit
integer :: aerr

   allocate(speXm(-2:2, nhspsmax), stat=aerr)
   if(aerr /= 0) call memerror("deformedSPEs")
   speXm(:,:) = 0.0

!   open(unit=22,file='deformed.spe',status='old',err=101)
  write(logfile,*)' Reading in deformed single-particle energies '
  do it = 1,2
    call readdeformedSPEs(it,foundit)
    if(.not.foundit)then
      print*,' Could not find file deformed.int '
      write(logfile,*)' Could not find file deformed.int '
    end if

  end do

return

end subroutine deformedSPEs
!=======================================================
!
!  INPUT:
!   it = species
!  OUTPUT
!   foundit = logical flag to signify success
!            in opening file deformed.int
!
!  CALLED BY:
!    deformedSPEs
!
subroutine readdeformedSPEs(it,foundit)
   use verbosity
   use spstate
   use haiku_info
   use system_parameters
   use interaction_deformed

   implicit none
   integer it      ! which species
   logical foundit
   integer ttmp, mtmp,jtmp,ntmp,ltmp
   integer mtot,partot
   integer, allocatable :: mapsps(:),mapm(:),mappar(:)
   integer nspslocal
   integer i
   integer ntbme  ! # of uncoupled TBMEs in file
   integer a,b,c,d
   integer isps,jsps,ksps,lsps
   integer cpair,dpair
   real hfactor
   real vxx
   logical foundsps
   integer iref, istart
   integer cphase,dphase
   integer itbme
   real,pointer :: v(:)
   integer,pointer :: map(:)
   integer jzmin,jzmax
   integer :: aerr

   print*,' XX deformed ',it,nhspsmax
   foundit = .false.
   open(unit=23,file='deformed.int',status='old',err=101)
   foundit = .true.
!...... FIRST READ PAST AND GET MAPPING OF SINGLE PARTICLE STATES.....

   read(23,*)nspslocal
   allocate(mapsps(nspslocal),mapm(nspslocal),mappar(nspslocal), stat=aerr )
   if(aerr /= 0) then
      call memerror("readdeformedSPEs")
      stop 5
   end if
   mapsps(:) = -1
   do i = 1,nspslocal
      read(23,*)ttmp, ntmp, ltmp,jtmp,mtmp
      if(ttmp /= it)then
         mapsps(i) = -1
         cycle
      end if
      mapm(i) = mtmp
      mappar(i) = (-1)**ltmp
      foundsps = .false.
!............ SEARCH FOR SINGLE-PARTICLE STATES.....
      if(mtmp > 0)then
        do isps = 1,nhsps(it)
          if( mtmp == hspsqn(it,isps)%m .and. jtmp == hspsqn(it,isps)%j  &
          .and. ntmp == hspsqn(it,isps)%nr .and. ltmp == hspsqn(it,isps)%l )then
              foundsps = .true.
              mapsps(i) = isps+nhsps(-it)
              exit
          end if

        end do

      else
        do isps = 1,nhsps(-it)
          if( mtmp == hspsqn(-it,isps)%m .and. jtmp == hspsqn(-it,isps)%j  &
          .and. ntmp == hspsqn(-it,isps)%nr .and. ltmp == hspsqn(-it,isps)%l )then
              foundsps = .true.
              mapsps(i) = isps
              exit
          end if

        end do

      end if
      if(.not.foundsps)then   ! ERROR TRAP
         print*,' DID NOT FIND CORRESPONDING STATE '
         print*,i,ntmp,ltmp,jtmp,mtmp
         stop

      end if

   end do
!................. READ SINGLE-PARTICLE ENERGIES  (added in V7.2.8)....

   do a = 1,nspslocal
      read(23,*)vxx
      if(mapsps(a) < 0)cycle  ! wrong species of nucleon
      isps = mapsps(a)
      if(isps > nhsps(-it))then
        speXm(-it,isps-nhsps(-it) ) = vxx
      else
        speXm(it,isps) = vxx
      end if

   end do

   close(23)
   deallocate(mapsps,mapm,mappar )
   return
101 continue
   print*,' file deformed.int does not exist '
   foundit = .false.
   return
end subroutine readdeformedSPEs
!
!
!=======================================================
!
!  intermediate routine to create single-particle energies for haikus
!
!  alternate to routine makehaikuspes
!
!  CALLED BY:
!    makespeme
!
   subroutine makehaikuspes_deformed(ith,bspe)

   use bitstuff
   use diagh
   use blocks
   use haiku_info
   use spstate
   use interaction_deformed

   implicit none
   integer ith
   type (blockhspe), pointer :: bspe(:)
!   real, pointer :: xspe(:)

!------------------------------------
   integer iblock
   integer i
   integer isps
   integer j,iword,imask
   integer iorb
   integer, pointer :: hsd(:,:)
   integer :: aerr

   allocate(bspe( hblock(ith)%nhblocks ), stat=aerr)
   if(aerr /= 0) call memerror("makehaikuspes_deformed 1")
   do iblock = 1,hblock(ith)%nhblocks
     allocate( bspe( iblock)%hspe ( hblock(ith)%list(iblock)%nhsd ), stat=aerr)
     if(aerr /= 0) call memerror("makehaikuspes_deformed 10")
     bspe(iblock)%hspe = 0.0
     do i = 1,hblock(ith)%list(iblock)%nhsd
        hsd => hblock(ith)%list(iblock)%hsd
        do isps = 1,nhsps(ith)
          iword = (isps-1)/max_bit_word+1
          j = mod(isps-1,max_bit_word)+1
          imask = ishft(1,j-1)
          if ( iand(imask,hsd(iword,i)) /=0 ) then
!              iorb = hspsqn(ith,isps)%orb
              bspe(iblock)%hspe(i) = bspe(iblock)%hspe(i) + speXm(ith,isps)
          end if

        enddo ! iop
     enddo ! i
   enddo  ! iblock
   return
   end subroutine makehaikuspes_deformed
!
!=======================================================

!
! this version is mostly obsolete;
!

subroutine deformedSPEsOLD

use interaction_deformed
use spstate
use haiku_info

implicit none
integer ith
integer isps
integer m,j,l,nr
integer i,jtmp,ltmp,mtmp,ttmp,nrtmp
real spetmp
integer :: aerr

   allocate(speXm(-2:2, nhspsmax), stat=aerr)
   if(aerr /= 0) call memerror("deformedSPEsOLD")
   speXm(:,:) = 0.0

   open(unit=22,file='deformed.spe',status='old',err=101)

   do ith = -2,2
       if(ith == 0)cycle
       do isps = 1,nhsps(ith)
          j = hspsqn(ith,isps)%j
          m = hspsqn(ith,isps)%m
          l = hspsqn(ith,isps)%l
          nr= hspsqn(ith,isps)%nr
          rewind(22)
          do i=1,1000
             read(22,*,end=102)ttmp,nrtmp, ltmp,jtmp,mtmp,spetmp
             if(ttmp /= abs(ith))cycle
             if(ltmp == l .and. jtmp == j .and. mtmp == m .and. nrtmp == nr)then
                 speXm(ith,isps) = spetmp
                 exit
             end if

          end do

       end do

   end do
   close(22)

return

101  continue
   print*,' Looking for file deformed.sps, does not exist '
   stop

102 continue
   print*,' ran out of states '
   print*,ith,l,j,m
end subroutine deformedSPEsOLD

max_stars_repo_name: noctildon/pyBigstick

max_issues_repo_name: noctildon/pyBigstick

max_forks_repo_name: noctildon/pyBigstick


!    ..............................................
!    From Hamiltonian (interaction) to HFB fileds
!    .............................................
        subroutine HFB_FIELD_COMPLEX(zro,zkapa10,zkapa01,   &
     &  zakin,zgamma,zham,zdelta10,zdelta01,NLEV)
!      ..............................................................
        USE VAPHFB_PAR
        implicit real*8 (a-h,o-y)
        implicit complex*16 (z)

        DIMENSION zro(NLEV,NLEV)
        DIMENSION zkapa10(NLEV,NLEV)
        DIMENSION zkapa01(NLEV,NLEV)
        DIMENSION zakin(NLEV,NLEV)
        DIMENSION zgamma(NLEV,NLEV)
        DIMENSION zham(NLEV,NLEV)
        DIMENSION zdelta10(NLEV,NLEV)
        DIMENSION zdelta01(NLEV,NLEV)

!     ...................... initialization
      do ia=1,NLEV
      do ic=1,NLEV
            ZGamma(ia,ic)  =zzero
            zDelta10(ia,ic)=zzero
            zDelta01(ia,ic)=zzero
            zham(ia,ic)    =zzero
      enddo
      enddo

      do iabcd=1,H%iabcd_max
          ia = H%ka(iabcd)
          ib = H%kb(iabcd)
          ic = H%kc(iabcd)
          id = H%kd(iabcd)
          ZGamma(ia,ic)=ZGamma(ia,ic) + H%ME2BM(iabcd)*ZRO(id,ib)
          ZGamma(ib,ic)=ZGamma(ib,ic) - H%ME2BM(iabcd)*ZRO(id,ia)

          ZGamma(ia,id)=ZGamma(ia,id) - H%ME2BM(iabcd)*ZRO(ic,ib)
          ZGamma(ib,id)=ZGamma(ib,id) + H%ME2BM(iabcd)*ZRO(ic,ia)

          if(ia+ib.eq.ic+id) cycle

!         .... (ab) <-> (cd)
          ZGamma(ic,ia)=ZGamma(ic,ia) + H%ME2BM(iabcd)*ZRO(ib,id)
          ZGamma(id,ia)=ZGamma(id,ia) - H%ME2BM(iabcd)*ZRO(ib,ic)

          ZGamma(ic,ib)=ZGamma(ic,ib) - H%ME2BM(iabcd)*ZRO(ia,id)
          ZGamma(id,ib)=ZGamma(id,ib) + H%ME2BM(iabcd)*ZRO(ia,ic)

!          if(abs(H%ME2BM(iabcd)).gt.1.d-6) write(201,*) ia,ib,ic,id,H%ME2BM(iabcd)
       enddo ! iabcd

!  ................... single-particle h
       do ia=1,NLEV
         do ic=1,NLEV
            zham(ia,ic)=zGamma(ia,ic)+zakin(ia,ic)
!  ..................w/o np-mixing and pairing
            if(Input%NPMix .eq. 0 .and. tnljm%t(ia) .ne. tnljm%t(ic)) then
               zham(ia,ic)  = 0.d0
               zGamma(ia,ic)= 0.d0
            endif
!  ..................w/o np-mixing and pairing
            if(Input%IPMix .eq. 0 .and. iv(tnljm%l(ia)) .ne. iv(tnljm%l(ic))) then
               zham(ia,ic)  = 0.d0
               zGamma(ia,ic)= 0.d0
            endif
!  ..................w/o k-mixing (triaxiality) and pairing
            if(Input%KMix .eq. 0 .and. tnljm%twom(ia) .ne. tnljm%twom(ic)) then
               zham(ia,ic)  = 0.d0
               zGamma(ia,ic)= 0.d0
            endif
         end do
        end do
        !Delta
      do iabcd=1,H%iabcd_max
          ia = H%ka(iabcd)
          ib = H%kb(iabcd)
          ic = H%kc(iabcd)
          id = H%kd(iabcd)
            zDelta10(ia,ib)=zDelta10(ia,ib)+0.5*H%ME2BM(iabcd)*zkapa10(ic,id)*2

            zDelta10(ib,ia)=zDelta10(ib,ia)-0.5*H%ME2BM(iabcd)*zkapa10(ic,id)*2

            zDelta01(ia,ib)=zDelta01(ia,ib)+0.5*H%ME2BM(iabcd)*zkapa01(ic,id)*2

            zDelta01(ib,ia)=zDelta01(ib,ia)-0.5*H%ME2BM(iabcd)*zkapa01(ic,id)*2

          if(ia+ib.eq.ic+id) cycle
!         .... (ab) <-> (cd)
            zDelta10(ic,id)=zDelta10(ic,id)+0.5*H%ME2BM(iabcd)*zkapa10(ia,ib)*2
            zDelta10(id,ic)=zDelta10(id,ic)-0.5*H%ME2BM(iabcd)*zkapa10(ia,ib)*2
            zDelta01(ic,id)=zDelta01(ic,id)+0.5*H%ME2BM(iabcd)*zkapa01(ia,ib)*2
            zDelta01(id,ic)=zDelta01(id,ic)-0.5*H%ME2BM(iabcd)*zkapa01(ia,ib)*2
       enddo ! iabcd
!  ..................w/o np-mixing and pairing
       do ia=1,NLEV
         do ib=1,NLEV
            if(Input%NPMix .eq. 0 .and. tnljm%t(ia) .ne. tnljm%t(ib)) then
             zDelta10(ia,ib) = 0.d0
             zDelta10(ib,ia) = 0.d0
             zDelta01(ia,ib) = 0.d0
             zDelta01(ib,ia) = 0.d0
            endif
            if(Input%IPMix .eq. 0 .and. iv(tnljm%l(ia)) .ne. iv(tnljm%l(ib))) then
             zDelta10(ia,ib) = 0.d0
             zDelta10(ib,ia) = 0.d0
             zDelta01(ia,ib) = 0.d0
             zDelta01(ib,ia) = 0.d0
            endif

            if(Input%KMix .eq. 0 .and. abs(tnljm%twom(ia)) .ne. abs(tnljm%twom(ib)) ) then
             zDelta10(ia,ib) = 0.d0
             zDelta10(ib,ia) = 0.d0
             zDelta01(ia,ib) = 0.d0
             zDelta01(ib,ia) = 0.d0
            endif

       enddo ! ib
       enddo ! ia

        end subroutine


! SPDX-FileCopyrightText: Copyright (c) 2021 Sam W
! SPDX-License-Identifier: MIT
! Copyright (c) 2021 Sam W
!
! 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.

! Entry Point
! This is similar to Pascal
program Printchara
    implicit none

    ! Declare variables
    INTEGER :: num
    INTEGER :: i
    CHARACTER(24) :: filename

    ! Ask how many 'A' should be written
    ! and the name of the file
    PRINT *, "How many 'A' should be written to a file: "
    READ(*,*) num

    PRINT *, "What is the name for the file: (Must be 24 characters): "
    READ(*,*) filename

    ! Create a file and write to it then close it
    OPEN(1, file = filename)
    do i = 0, num
        WRITE(1, "(A)", advance='no') 'A'
    end do
    CLOSE(1)

    ! If previous operation completed then press Enter to exit
    PRINT *, "Operation Complete. Press Enter to exit."
    READ(*,*)
end program Printchara


text:
      SUBROUTINE EROZ (ROTANG, XFRMT)

      REAL XFRMT (4,4), SN, CS, ROTANG

      CALL IDMAT (4, XFRMT)
      SN = SIN(ROTANG)
      CS = COS(ROTANG)
      XFRMT (1,1) =  CS
      XFRMT (2,2) =  CS
      XFRMT (1,2) = -SN
      XFRMT (2,1) =  SN

      END


text: #redirect Gospel Choir at UC Davis

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text: #redirect Acme Theatre Co.


text:
C    @(#)solton.f	20.13 5/27/99
      subroutine solton

      include 'ipfinc/parametr.inc'

      include 'ipfinc/addtbx.inc'
      include 'ipfinc/alpha.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/chgprm.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/epridc.inc'
      include 'ipfinc/gendrp.inc'
      include 'ipfinc/ikk.inc'
      include 'ipfinc/intbus.inc'
      include 'ipfinc/itrhis.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/ltcsln.inc'
      include 'ipfinc/prt.inc'
      include 'ipfinc/slnopt.inc'
      include 'ipfinc/tbx.inc'
      include 'ipfinc/tbxsrt.inc'
      include 'ipfinc/update.inc'
      include 'ipfinc/usranl.inc'
      include 'ipfinc/ycomp.inc'
      include 'ipfinc/comdrx.inc'
      include 'ipfinc/xblim.inc'
      include 'ipfinc/miscfile.inc'

      common /drptmp/ drptmp

      common /ycpsln/ dytot, ditot, kownty, ycpscr(6,MXYCMP)

      common /ownflg/ ownflg
      logical ownflg

      common /is_batch / is_batch

      common /sol_temp/ idchgx, kswsln, itlow

      integer fichsv, crtsv, status, chkerr, error, open_file, ckzoneup
      logical batch, chkbch, epdc, opened

      external nrpqv

      character prmfil*60, tmpfil*60, tempc*60, bigbuf*512, ljstfy*512

C     Initialize solution arrays.

      do i = 1, MAXLTC
         ltcsln(i) = 0
      enddo
      do i = 1, MAXTBX
         tbxsrt(i) = 0
         tbxrat(i) = 1.0
         do j = 1, 6
            tbxslp(j,i) = 0.0
         enddo
      enddo
      do i = 1, kxtot
         xblim(1,i) = -9999.0
         xblim(2,i) = -9999.0
         xvlim(1,i) = 0.0
         xvlim(2,i) = 0.0
         xvrlim(1,i) = 0.0
         xvrlim(2,i) = 0.0
      enddo
      do i = 1, MXYCMP
         do j = 1, 6
            ycpscr(j,i) = 0.0
         enddo
      enddo
      if (nepbus .gt. 0) then
         epdc = .true.
      else
         epdc = .false.
      endif
      lprtsv=lprtsw
      fichsv=fichsw
      crtsv=crtsw
      crtsw=0
      koptsw = 0
      lskp = 1
      ownflg = .false.             ! Set flag to force updating OWNLOS
      update(1) = 1                ! Set flag to force updating ZSUM
      ickzone = ckzoneup(status)   ! Determine whether it is necessary
C                                  ! to update ZSUM each solution
      call presln

      if ( chkbch (batch)) then
         crtsw=0
      else
         crtsw=crtsv
      endif

      if (iopton(20) .eq. 1) then
         crtsw=crtsv
         write(outbuf, 59)
   59    format('0     Base case solution voltages used by request    ')
         call prtout (1)
         call nrbase
         crtsw=0
         go to 200
      else if (iopton(20) .eq. 2) then
         crtsw=crtsv
         write(outbuf, 60)
   60    format('0     Alternate solution voltages obtained by ',
     &          '/LOAD_SOLUTION command   ')
         call prtout (1)
         call nrbase
         crtsw=0
         go to 200
      else if (iopton(20) .eq. 3) then
         crtsw=crtsv
         write(outbuf, 61)
   61    format('0     Alternate starting voltages obtained by ',
     &          '/LOAD_SOLUTION command   ')
         call prtout (1)
      endif

      lprtsv=lprtsw
      fichsv=fichsw
      lprtsw=1
      if (kspare(16) .ge. 0) fichsw=1
C
C     Check override of default decoupled solution iteration.
C
      if (iopton(1) .eq. 2 .and. numdrp .gt. 0) then
         iopton(1) = 6
         write (errbuf(1), 62) iopton(1)
   62    format('/GEN_DROP has forced ', i2, ' DECOUPLED iterations.')
         call prterx ('W',1)
      endif

      call sumsol
C
C     Open scratch file
C
      inquire (unit = lunscr, opened = opened)
      if (.not. opened) then
        status = open_file (lunscr, ' ', 'U', 'W', iostat)
      endif

      if (epdc) then
C
C        Initialization of variables.
C
         call epdcin
C
C        Calculation of per unit constants and the per unitizing of
C        data.
C
         call lfdcpu
C
C        Determination of VECBUS.
C
         call solvec (jdw1, jdw2, jdw3, jdw4)
C
C        Topology checking.
C
         call dctopo
C
C        Generation of nodal conductance matrix.
C
         call bldgdc
C
C        Printing of DC input data.
C
C        CALL PRTDCI
C
      end if

      call opsln1

      kerr=chkerr('F')
      if (kerr .gt. 0) then
         write (errbuf(1), 80)
   80    format(' Case aborted by fatal errors.')
         if (is_batch .eq. 0) then
            call prterx ('E',1)
            return
         else
            call prterx ('F',1)
            call erexit()
         endif
      endif

      lprtsw=lprtsv
      fichsw=fichsv
      crtsw=1
      call slnhdr
      drptmp = drptot
      if (iopton(1).gt.0) then
         if (numdrp .gt. 0) then
            call comdrp ('DC Iteration', 1, num1, drptmp, num2, drop2,
     1                    status, nrpqv, itrx)
         endif
C
C        Set switches and counters to store the decoupled solution
C        history data.
C
         itrtyp = 1
         call savsln
         call dcstrt  (kerr)
C
C        Reset switches and counters to store the Newton_Raphson
C        solution history data.
C
         itrhis(1) = ittot - 1
         itrhis(2) = 0
         itrtyp = 2
         if (kerr.ne.0) call retsln
      else
C
C        Otherwise, obtain X-bus sensitivities.
C
         if (kxtot .gt. 0) call getxsn(error)
      endif
      if (iopton(2).gt.0) then
         if (numdrp .gt. 0) then
            call comdrp ('NX Iteration', 1, num1, drptmp, num2, drop2,
     1                    status, nrpqv, itrx)
         endif

         call savsln
         call cursol (kerr)
         if (kerr .ne. 0) call retsln
      endif
      call nrcnvt
      lprtsv=lprtsw
      fichsv=fichsw
      lprtsw=1
      if (kspare(16) .ge. 0) fichsw=1

      ittot=0
      kswsln = 1
      idchgx = 0

      if (numdrp .gt. 0) then
         drop2 = 1.0e10
         itx = 1
         do while (itx .le. 4 .and. drop2 .ge. 10.0)
            call comdrp ('NR Iteration', 1, num1, drptmp, num2, drop2,
     1                    status, nrpqv, itrx)
            if (itx .gt. 1) ittot = 3
            call nrsoln (ittot, iopton(3), idchgx, kswsln, itlow)
            if (itrhis(2) .eq. 0) itrhis(2) = itrhis(1) + ittot
            itx = itx + 1
         enddo
      else
         call nrsoln (ittot, iopton(3), idchgx, kswsln, itlow)
         if (itrhis(2) .eq. 0) itrhis(2) = itrhis(1) + ittot
      endif
C
C     Perform parameter perturbation summary.
C
      if (index (buf,'CHANGE_PAR') .ne. 0 .or.
     1    index (buf,'CHANGEPAR') .ne. 0) then
         if (lskp .eq. 3) then
            write (errbuf(1), 83)
   83       format ('The following / CHANGE_PARAMETERS ',
     &              'commands aborted by failed solution.')
            call prterx ('W', 1)
C
C           Restore solution
C
            call retsln
            error = 1
C
C           Skip current and subsequent / CHANGE_PARAMETER records.
C
            do while (buf(1:1) .eq. '/' .and.
     &               (index (buf,'CHANGE_PAR') .ne. 0 .or.
     &                index (buf,'CHANGEPAR') .ne. 0))

               write (errbuf(1), 84) buf(1:80)
   84          format (' COMMAND IGNORED (', a, ')')
               call prterx ('W', 1)
               card = buf(1:1)

               read (inp, 105, end=85) buf
               card = buf(1:1)
               last = lastch (buf)
               do while (buf(last:last) .eq. '-')
                  write (errbuf(1), 103) buf(1:80)
                  call prterx ('W', 1)
                  read (inp, 105, end=85) buf
                  card = buf(1:1)
                  last = lastch (buf)
               enddo
C
C              Examine next input record
C
               do while (card .eq. '.')
                  read (inp, 105, end=85) buf
                  card = buf(1:1)
               enddo
            enddo
            go to 86

   85       buf = '(STOP) SOLTON'
            card = buf(1:1)

   86       call ctlpow
         else

            bigbuf = buf
            call space (1)
            write (outbuf, 87) buf(1:80)
   87       format (' Solution Command (', a, ')')
            last = lastch (bigbuf)
            do while (bigbuf(last:last) .eq. '-')
               read (inp, fmt='(a)', end=88) buf
               write (outbuf, 87) buf(1:80)
               call prtout (1)
               if (buf(1:1) .ne. '.')
     &           bigbuf(last:) = ' ' // ljstfy (buf)
               last = lastch (bigbuf)
               buf = bigbuf
            enddo
   88       continue
C
C           Save initial solution in case we get a failed solution.
C
            call savsln
C
C           Open file for parameters
C
            inquire (unit=inp, name=prmfil)
            i = index (prmfil, ']')
            j = index (prmfil(i+1:), '.') + i
            if (j .eq. 0) then
               j = lastch (prmfil) + 1
               prmfil(j:j) = '.'
            endif
            tempc = prmfil
            prmfil = tempc(1:j) // 'QVPT'
            status = open_file (23, prmfil, 'F', 'W', iostat)
            numprm = 0
C
C           Open OUTPUT files.  If none are explicitly defined,
C           create default names.
C
            do iuser = 1, numusr
               if (usrfil(iuser) .eq. ' ') then
                  write (usrfil(iuser), 92) prmfil(1:j),
     1               'USER_REPORT_', iuser
   92             format (a, a, i2.2)
               endif
               lun(iuser) = 26
               close(unit=lun(iuser),err=93)
   93          open (unit = lun(iuser),
     1               file = usrfil(iuser),
     2               status = 'NEW',
     3               form = 'FORMATTED',
     5               iostat = ios,
     6               err = 94)
               go to 97

   94          last = lastch(usrfil(iuser))
               write (errbuf(1), 95) ios, usrfil(iuser)(1:last),
     1            lun(iuser)
   95          format (' Error ', i3, ' opening file ', a,
     1             ' on logical unit ', i2)
               call prterx ('W', 1)

   97          close(unit=lun(iuser))
               numupx(iuser) = 0
            enddo
C
C           Write user-specified report if requested.
C
            if (numusr .gt. 0) then
               call postsl
               update(1) = ickzone
               call updzon()
               call prtusr (lun, numupx)
            endif

  100       call getprm (bigbuf, error)
C
C           We start at iteration 3 because bus switching is permitted.
C           The bus history logic has a maximum iteration limit of 50,
C           which could easily be violated with several parameter studie
C
            ittot = 3
            call nrsoln (ittot, iopton(3), idchgx, kswsln, itlow)
            if (itrhis(2) .eq. 0) itrhis(2) = itrhis(1) + ittot
            if (lskp .ne. 3 .and. numdrp .gt. 0) then
               itrtot = 1
               if (numdrp .gt. 0) then
                  drop2 = 1.0e10
                  itx = 1
                  do while (itx .le. 4 .and. drop2 .ge. 10.0)
                     call comdrp ('NR Iteration', 1, num1, drptmp,
     1                             num2, drop2, status, nrpqv, itrx)
                     ittot = 3
                     call nrsoln (ittot, iopton(3), idchgx, kswsln,
     1                            itlow)
                     if (itrhis(2) .eq. 0) then
                        itrhis(2) = itrhis(1) + ittot
                     endif
                     itx = itx + 1
                  enddo
               endif
            else
               ittot = 3
               call nrsoln (ittot, iopton(3), idchgx, kswsln, itlow)
               if (itrhis(2) .eq. 0) itrhis(2) = itrhis(1) + ittot
            endif
            kerr=chkerr('F')
            if (lskp .eq. 3 .or. kerr .gt. 0) then
               write (errbuf(1), 101)
  101          format ('The following / CHANGE_PARAMETERS ',
     &                 'commands aborted by failed solution.')
               call prterx ('W', 1)
C
C              Restore solution
C
               call retsln
               error = 1
C
C              Skip current and subsequent / CHANGE_PARAMETER records.
C
               do while (card .eq. '/' .and.
     &                  (index (buf,'CHANGE_PAR') .ne. 0 .or.
     &                  index (buf,'CHANGEPAR') .ne. 0))
                  write (errbuf(1), 103) buf(1:80)
  103             format (' COMMAND IGNORED (', a, ')')
                  call prterx ('W', 1)

                  read (inp, 105, end=107) buf
  105             format (a)
                  card = buf(1:1)
                  last = lastch (buf)
                  do while (buf(last:last) .eq. '-')
                     write (errbuf(1), 103) buf(1:80)
                     call prterx ('W', 1)
                     read (inp, 105, end=107) buf
                     card = buf(1:1)
                     last = lastch (buf)
                  enddo
C
C                 Examine next input record
C
                  do while (card .eq. '.')
                     read (inp, 105, end=107) buf
                     card = buf(1:1)
                  enddo
               enddo
               go to 108

  107          buf = '(STOP) SOLTON'
               card = buf(1:1)

  108          call ctlpow
            else
               call prtprm (xpoint, ypoint, error)
               call savsln
C
C              Write user-specified report if requested.
C
               if (numusr .gt. 0) then
                  call postsl
                  update(1) = ickzone
                  call updzon()
                  call prtusr (lun, numupx)
               endif
            endif

            if (index (buf,'CHANGE_PAR') .ne. 0 .or.
     1         index (buf,'CHANGEPAR') .ne. 0) then

               bigbuf = buf
               call space (1)
               write (outbuf, 87) buf(1:80)
               last = lastch (bigbuf)
               do while (bigbuf(last:last) .eq. '-')
                  read (inp, fmt='(a)', end=109) buf
                  write (outbuf, 87) buf(1:80)
                  call prtout (1)
                  if (buf(1:1) .ne. '.')
     &              bigbuf(last:) = ' ' // ljstfy (buf)
                  last = lastch (bigbuf)
                  buf = bigbuf
               enddo
  109          continue
               go to 100
            else
C
C              Close CHANGE_PARAMETER file
C
               if (numprm .gt. 0) then
                  inquire (unit=23, name=prmfil)
                  l = lastch (prmfil)
                  write (outbuf, 110) numprm, prmfil(1:l)
  110             format ('0 ', i4, ' points written to file ',
     1                a)
                  call prtout (1)
               endif
               close (unit=23, err=112)
  112          continue
               prmfil = ' '
               numprm = 0
C
C              Close USER_REPORT file
C
               do iuser = 1, numusr
                  if (numupx(iuser) .gt. 0) then
                     inquire (file=usrfil(iuser), name=tmpfil)
                     l = lastch (tmpfil)
                     write (outbuf, 120) numupx(iuser), tmpfil(1:l)
  120                format ('0 ', i4, ' records written to file ',
     1                   a)
                     call prtout (1)
                  endif
                  usrfil(iuser) = ' '
               enddo
            endif
         endif
      endif
C
C     Print out solution summary if failed solution.
C
      if (lskp .eq. 3 .or. iopton(22) .ne. 0) then
         call itrsum
      endif

      lprtsw=lprtsv
      fichsw=fichsv
C
C     Perform post-solution update:
C
C     1. Update P- and Q-injections for slack buses, PV buses.
C     2. Update d-c values.
C     3. Update Taps for LTC transformers, Xij and Bis for variable
C        var-compensation.
C
      call postsl
C
C     Summarize Generator Dropping
C
      crtsw=0
      if (numdrp .gt. 0) then
         call sumdrp
      endif

      if (index (buf,'BUSSEN') .ne. 0 .or.
     1    index (buf,'BUS_SEN') .ne. 0 .or.
     2    index (buf,'LINESE') .ne. 0 .or.
     3    index (buf,'LINE_SE') .ne. 0 .or.
     4    index (buf,'LOSSSE') .ne. 0 .or.
     5    index (buf,'LOSS_SE') .ne. 0 .or.
     6    index (buf,'TRANSF') .ne. 0) then
         call opsln3
      endif

      call nrrstr

C     Restore temporary "BF" NTYP created in DCSTRT.

      do i = 1, ntbxad
         jt = ntotb - ntbxad + i
         kt = tbx(2,jt)
         ntypu(kt) = ltbxad(i)
      enddo
      ntotb = ntotb - ntbxad

  200 crtsw=crtsv
      ownflg = .false. ! Flag OWNLOS as being obsolete
      update(1) = 1    ! Flag ZSUM as being obsolete
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      SUBROUTINE IFP3
C
C     DATA PROCESSING AND GENERATION OF THE AXIS-SYMETRIC-CONICAL SHELL
C
C          CARDS          TYPE       REC.ID-BIT CARDS-FILE,  CARDS-FILE
C     ===  =======      ===========  ========== ===========  ==========
C      1   AXIC     --  AX.SY.SHELL     515- 5
C      2   CCONEAX  --  AX.SY.SHELL    8515-85  CCONE-GEOM2,
C      3   FORCEAX  --  AX.SY.SHELL    2115-21  FORCE-GEOM3,
C      4   FORCE    --  STANDARD       4201-42  FORCE-GEOM3,
C      5   GRAV     --  STANDARD       4401-44   GRAV-GEOM3,
C      6   LOAD     --  STANDARD       4551-61   LOAD-GEOM3,
C      7   MOMAX    --  AX.SY.SHELL    3815-38  MOMNT-GEOM3,
C      8   MOMENT   --  STANDARD       4801-48  MOMNT-GEOM3,
C      9   MPCADD   --  STANDARD       4891-60 MPCADD-GEOM4,
C     10   MPCAX    --  AX.SY.SHELL    4015-40    MPC-GEOM4,
C     11   OMITAX   --  AX.SY.SHELL    4315-43   OMIT-GEOM4,
C     12   POINTAX  --  AX.SY.SHELL    4915-49    MPC-GEOM4, GRID-GEOM1
C     13+  RFORCE   --  STANDARD       5509-55 RFORCE-GEOM3,
C     14   RINGAX   --  AX.SY.SHELL    5615-56    SPC-GEOM4, GRID-GEOM1
C     15   SECTAX   --  AX.SY.SHELL    6315-63    MPC-GEOM4, GRID-GEOM1
C     16   SEQGP    --  STANDARD       5301-53  SEQGP-GEOM1,
C     17   SPCADD   --  STANDARD       5491-59 SPCADD-GEOM4,
C     18   SPCAX    --  AX.SY.SHELL    6215-62    SPC-GEOM4,
C     19   SUPAX    --  AX.SY.SHELL    6415-64 SUPORT-GEOM4,
C     20   TEMPAX   --  AX.SY.SHELL    6815-68   TEMP-GEOM3,
C     21   TEMPD    --  STANDARD       5641-65  TEMPD-GEOM3,
C     22   CTRIAAX  --  AX.TR.CR       7012-70  CTRIA-GEOM2
C     23   CTRAPAX  --  AX.TRA.CR      7042-74  CTRAP-GEOM2
C
      IMPLICIT INTEGER (A-Z)
      EXTERNAL        LSHIFT    ,RSHIFT    ,ANDF      ,ORF       ,
     1                COMPLF
      LOGICAL         SECD      ,NOGO      ,RECOFF    ,PIEZ
      REAL            NPHI      ,NPHI1     ,NISQ      ,NI        ,
     1                A1        ,A2        ,A3        ,A4        ,
     2                ANGLE     ,RADDEG    ,PI        ,DIFPHI    ,
     3                RZ        ,T1        ,T2        ,COEF      ,
     4                CONSTS    ,SUM       ,TWOPI
      DIMENSION       GEOM(4)   ,Z(8)      ,NUM(11)   ,INUM(11)  ,
     1                MSG1(2)   ,MSG2(2)
      CHARACTER       UFM*23    ,UWM*25    ,UIM*29    ,SFM*25
      COMMON /XMSSG / UFM       ,UWM       ,UIM       ,SFM
      COMMON /SYSTEM/ IBUFSZ    ,NOUT      ,NOFLAG    ,DUMDUM(8) ,
     1                NLINES    ,DUM1(26)  ,NBPC      ,NBPW      ,
     2                DUM37(37) ,IPIEZ
      COMMON /TWO   / TWO(32)
      COMMON /CONDAS/ CONSTS(5)
      COMMON /IFP3LV/            RECID(3)  ,RECID1(3) ,RECIDX(3) ,
     1                IEND      ,REC(3)    ,REC1(3)   ,TRAIL(7)  ,
     2                IT        ,AXTRL(7)  ,OPENFL(6) ,N         ,
     3                A1        ,CSID      ,NI        ,NISQ      ,
     4                A2        ,IBUFF1    ,IBUFF2    ,IBUFF3    ,
     5                A3        ,BUFF      ,NOGO      ,OP        ,
     6                A4        ,IHEADR    ,IBITR     ,IFILE     ,
     7                NOREG     ,LAST      ,IERRTN    ,ICONT     ,
     8                NOAXIC    ,RINGID    ,OUTBUF    ,VEOR      ,
     9                ISTART    ,IRETRN    ,FLAG      ,IAMT      ,
     O                SUM       ,IBIT      ,SETID     ,SORC      ,
     1                IBEGIN    ,MPCON     ,NWORDS    ,NNN       ,
     2                ANGLE     ,K3OR6     ,NPHI1     ,ZPT       ,
     3                NMOVE     ,CSSET     ,NOPONT    ,NON       ,
     4                IPHI      ,RECOFF    ,NPHI      ,N3OR5     ,
     5                ION       ,NPLUS1    ,NOSECT    ,COEF      ,
     6                IPT       ,COMPON    ,ICORE     ,ISCRAT    ,
     7                ICORE1    ,NCARDS    ,I1        ,IAT       ,
     8                I2        ,T1        ,T2        ,NFILE     ,
     9                NADD      ,NCARD
      COMMON /IFP3CM/ FILE(6)   ,INAME(12) ,CDTYPE(50),AXIC1(3)  ,
     1                CCONEX(3) ,FORCEX(3) ,FORCE(3)  ,GRAV(3)   ,
     2                LOAD(3)   ,MOMAX(3)  ,MOMENT(3) ,MPCADD(3) ,
     3                MPCAX(3)  ,OMITAX(3) ,POINTX(3) ,PRESAX(3) ,
     4                RINGAX(3) ,SECTAX(3) ,SEQGP(3)  ,SPCAX(3)  ,
     5                SUPAX(3)  ,TEMPAX(3) ,TEMPD(3)  ,PLOAD(3)  ,
     6                MPC(3)    ,SPC(3)    ,GRID(3)   ,SUPORT(3) ,
     7                NEG111(3) ,T65535(3) ,TEMP(3)   ,OMIT(3)   ,
     8                SPCADD(3) ,ONE       ,ZERO      ,IHEADB(96),
     9                CTRIAA(3) ,CTRAPA(3) ,ICONSO    ,RFORCE(3)
      COMMON /OUTPUT/ DUMMY(96) ,IHEAD(96)
      COMMON /ZZZZZZ/ RZ(1)
      EQUIVALENCE     (CONSTS(1),PI     )  ,(CONSTS(2),TWOPI  )  ,
     1                (CONSTS(4),RADDEG )  ,(Z(1)     ,RZ(1)  )  ,
     2                (GEOM(1)  ,FILE(1))  ,(SCRTCH   ,FILE(5))  ,
     3                (AXIC     ,FILE(6))  ,(NUM(11)  ,B      )  ,
     4                (NOEOR    ,INPRWD    ,ZERO)                ,
     5                (EOR      ,CLORWD    ,OUTRWD    ,ONE    )
      DATA    INUM  / 1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H /
      DATA    IFIAT / 4HFIAT/, IFIST  /  4HFIST/,I5,I6 / 5,6     /
      DATA    MSG1  / 4HIFP3 , 4HBEGN /, MSG2  / 4HIFP3, 4HEND   /
C
      CALL CONMSG (MSG1,2,0)
C
C     RIGHT-JUSTIFY INUM AND CALL IT NUM
C
      DO 1 I = 1,11
    1 NUM(I) = RSHIFT(INUM(I),NBPW-NBPC)
C
C     INITIAL CHECK TO MAKE SURE TRAILER BITS ARE ALL OFF FOR GEOM1,
C     GEOM2, GEOM3, GEOM4.
C
      DO 10 I = 1,96
   10 IHEAD(I) = IHEADB(I)
C
      IF (NOFLAG) 30,20,30
   20 NOGO = .FALSE.
      GO TO 40
   30 NOGO = .TRUE.
C
   40 OPENFL(1) = 0
      OPENFL(2) = 0
      OPENFL(3) = 0
      OPENFL(4) = 0
      OPENFL(5) = 0
      OPENFL(6) = 0
      DO 110 I = 1,4
      TRAIL(1) = GEOM(I)
      CALL RDTRL (TRAIL(1))
      IF (TRAIL(1)) 50,50,70
   50 CALL PAGE2 (3)
      IMSG = 1061
      WRITE  (NOUT,55) SFM,IMSG
   55 FORMAT (A25,I5)
      WRITE  (NOUT,60) GEOM(I),INAME(2*I-1),INAME(2*I),IFIAT
   60 FORMAT (5X,11HFILE NUMBER,I4,3H ( ,2A4,12H) IS NOT IN ,A4,1H.)
      NOGO = .TRUE.
      GO TO 110
C
   70 DO 100 J = 2,7
      IF (TRAIL(J)) 80,100,80
   80 CALL PAGE2 (3)
      IMSG = 1062
      WRITE  (NOUT,55) SFM,IMSG
      WRITE  (NOUT,90) GEOM(I),INAME(2*I-1),INAME(2*I)
   90 FORMAT (5X,'FILE NUMBER',I4,3H ( ,2A4,') HAS TRAILER BIT ON.  ',
     1       'FILE SHOULD BE CLEAN AT ENTRY TO IFP3.')
      NOGO = .TRUE.
      GO TO 110
  100 CONTINUE
  110 CONTINUE
C
C     PROCEED TO SETUP CORE AND OPEN AXIC FILE
C     ICORE1 WILL ALWAYS EQUAL THE GROSS OPEN CORE TO IFP3 AT START
C
      ICORE1 = KORSZ(Z)
      IBUFF1 = ICORE1 - IBUFSZ - 2
      IBUFF2 = IBUFF1 - IBUFSZ
      IBUFF3 = IBUFF2 - IBUFSZ
      ICORE  = IBUFF3 - 1
      ICRQ   = 100 - ICORE
      IF (ICORE .LT. 100) GO TO 1310
C
C     OPEN  AXIC FILE
C
      CALL PRELOC (*1330,Z(IBUFF1),AXIC)
      OPENFL(6) = 1
      AXTRL(1)  = AXIC
      CALL RDTRL (AXTRL(1))
C
C     READ AXIC CARD
C
      CALL LOCATE (*130,Z(IBUFF1),AXIC1(1),FLAG)
      CALL READ (*1600,*130,AXIC,Z(1),2,EOR,FLAG)
      N    = Z(1)
      CSID = Z(2)
      NNN  = N
      NCARD= 1
      ASSIGN 140 TO IERRTN
      GO TO 1420
C
C     MISSING REQUIRED AXIC CARD
C
  130 ASSIGN 140 TO IERRTN
      NNN   = 0
      NCARD = 1
      GO TO 1510
  140 N = NNN
      NPLUS1 = N + 1
C
C
C     GEOM2  PROCESSING
C     =================
C
C     OPEN GEOM2
C
      IFILE= GEOM(2)
      I    = 2
      OP   = OUTRWD
      BUFF = IBUFF2
      ASSIGN 150 TO IRETRN
      GO TO 1350
C
C     CCONEAX CARDS
C
  150 REC(1) = CCONEX(1)
      REC(2) = CCONEX(2)
      REC(3) = CCONEX(3)
      NCARD  = 2
C
C     IF THERE IS NO CCONEAX CARD, THEN GO TO 1750 AND LOOK FOR
C     CTRAPAX OR CTRIAAX CARDS
C
      ICONB  = 0
      ICONSO = 0
      CALL LOCATE (*1750,Z(IBUFF1),REC(1),FLAG)
C
C     INPUT IS IN 4-WORD CARDS
C     OUTPUT IS N+1 4-WORD CARDS FOR EACH CARD INPUT
C
C     RECORD HEADER FOR CCONES
C
      ASSIGN 160 TO IHEADR
      GO TO 1470
C
  160 CALL READ (*1600,*220,AXIC,Z(1),4,NOEOR,IAMT)
C
C     CHECK RING ID-S FOR SIZE
C
      NNN = Z(3)
      ASSIGN 170 TO IERRTN
      GO TO 1440
  170 NNN = Z(4)
      ASSIGN 180 TO IERRTN
      GO TO 1440
C
C     CHECK CCONEAX ID FOR 1-9999 ALLOWABLE RANGE
C
  180 IF (Z(1).GT.0 .AND. Z(1).LT.10000) GO TO 200
      CALL PAGE2 (3)
      IMSG = 361
      WRITE  (NOUT,185) UFM,IMSG
  185 FORMAT (A23,I4)
      WRITE  (NOUT,190) Z(1)
  190 FORMAT (5X,'CCONEAX ID =',I10,'.  OUT OF 1 TO 9999 PERMISSIBLE ',
     1       'RANGE')
      NOGO = .TRUE.
C
  200 Z(1) = Z(1)*1000
      DO 210 I = 1,NPLUS1
      Z(1) = Z(1) + 1
      Z(3) = Z(3) + 1000000
      Z(4) = Z(4) + 1000000
      IF (NOGO) GO TO 210
      CALL WRITE (GEOM(2),Z(1),4,NOEOR)
  210 CONTINUE
      GO TO 160
C
C     OUT OF CCONEAX CARDS
C
  220 IF (IAMT) 230,240,230
C
C     GO TO 356 FOR RECORD ERROR
C
  230 ASSIGN 260 TO IERRTN
      GO TO 1490
C
C     WRITE EOR AND PUT BITS IN TRAILER
C
  240 ASSIGN 250 TO IRETRN
      ICONSO = 1
      GO TO 1270
  250 ICONB = 1
      GO TO 1750
C
C     CLOSE GEOM2
C
  260 I = 2
      ASSIGN 270 TO IRETRN
      GO TO 1380
C
C     GEOM3 PROCESSING
C     ================
C
C     OPEN GEOM3
C
  270 IFILE= GEOM(3)
      I    = 3
      OP   = OUTRWD
      BUFF = IBUFF2
      ASSIGN 280 TO IRETRN
      GO TO 1350
C
C     FORCE, FORCEAX, MOMNT, AND MOMNTAX CARDS
C
  280 RECID(1) = FORCE(1)
      RECID(2) = FORCE(2)
      RECID(3) = FORCE(3)
      RECIDX(1)= FORCEX(1)
      RECIDX(2)= FORCEX(2)
      RECIDX(3)= FORCEX(3)
      NCARD    = 3
      ASSIGN 620 TO ICONT
C
C     SET NOREG = 0 OR 1, DEPENDING ON PRESSENCE OF RECID
C     SET NOAXIC= 0 OR 1, DEPENDING ON PRESSENCE OF RECIDX
C
  290 IBIT = RECIDX(2)
      ASSIGN 300 TO IBITR
      GO TO 1460
  300 NOAXIC = NON
      IBIT = RECID(2)
      ASSIGN 310 TO IBITR
      GO TO 1460
  310 NOREG = NON
C
      REC(1) = RECID(1)
      REC(2) = RECID(2)
      REC(3) = RECID(3)
C
      IF (NOAXIC) 340,320,340
  320 IF (NOREG ) 330,610,330
C
C     TRANSFER FORCE OR MOMENT RECORD DIRECTLY.
C     THERE ARE NO FORCEAX OR MOMAX CARDS RESPECTIVELY.
C
  330 ASSIGN 610 TO IRETRN
      GO TO 1230
C
C     AT 410 READ IN ALL FORCEAX OR MOMNTAX CARDS AND PUT OUT ON GEOM(3)
C     IF NOREG=0,AND ON SCRTCH IF NOREG NON-ZERO.FIRST WRITE 3-WORD-
C     REC ID ON GEOM3.
C
  340 ASSIGN 350 TO IHEADR
      GO TO 1470
C
C     OPEN SCRATCH IF NEEDED
C
  350 IF (NOREG) 360,370,360
  360 I    = 5
      OP   = OUTRWD
      BUFF = IBUFF3
      ASSIGN 370 TO IRETRN
      GO TO 1350
  370 CALL LOCATE (*1530,Z(IBUFF1),RECIDX(1),FLAG)
  380 CALL READ (*1600,*440,AXIC,Z(1),8,NOEOR,IAMT)
C
C     CHECK RING ID
C
      ASSIGN 390 TO IERRTN
      NNN = Z(2)
      GO TO 1440
C
C     CHECK HARMONIC NUMBER AND FOR A SEQUENCE OF HARMONICS
C
  390 IF (Z(4) .EQ. 0) GO TO 396
      II   = 1
      NH1  = 0
      NH2  = 0
      SECD = .TRUE.
      WORD = 4
      DO 391 IJ = 1,2
      DO 392 IX = 1,4
      CHR = RSHIFT(LSHIFT(Z(WORD),NBPC*IABS(IX-4)),NBPW-NBPC)
      IF (CHR .EQ. B) GO TO 392
      DO 393 I = 1,10
      K = I-1
      IF (NUM(I) .EQ. CHR) GO TO 394
  393 CONTINUE
      SECD = .FALSE.
      II   = 1
      GO TO 392
  394 IF (SECD) GO TO 395
      NH1 = NH1 + II*K
      II  = II*10
      GO TO 392
  395 NH2 = NH2 + II*K
      II  = II*10
  392 CONTINUE
  391 WORD = WORD -1
      IF (NH1 .LE. NH2) GO TO 398
      WORD = NH1
      NH1  = NH2
      NH2  = WORD
  398 NNN  = NH1
      ASSIGN 397 TO IERRTN
      GO TO 1420
  396 NH1  = Z(3)
      NH2  = Z(3)
  397 NNN  = NH2
      ASSIGN 400 TO IERRTN
      GO TO 1420
  400 Z(4) = Z(5)
      Z(5) = Z(6)
      Z(6) = Z(7)
      Z(7) = Z(8)
      NH1  = NH1 + 1
      NH2  = NH2 + 1
      SUM  = Z(2)
      MUS  = Z(2)
      DO 430 I = NH1,NH2
      Z(2) = MUS + I*1000000
      Z(3) = 0
C
C     OUTPUT TO GEOM(3) IF NOREG = 0
C     OUTPUT TO SCRTCH  IF NOREG = NON-ZERO
C
      IF (NOGO ) GO TO 380
      IF (NOREG) 420,410,420
  410 NFILE = GEOM(3)
      GO TO 430
  420 NFILE = SCRTCH
  430 CALL WRITE (NFILE,Z(1),7,NOEOR)
      GO TO 380
C
C     OUT OF CARDS
C
  440 IF (IAMT) 450,460,450
C
C     CHECK FOR RECORD INCONSISTANCY ERROR.
C
  450 REC(1) = RECIDX(1)
      REC(2) = RECIDX(2)
      REC(3) = RECIDX(3)
      ASSIGN 460 TO IERRTN
      GO TO 1490
C
  460 IF (NOREG) 470,590,470
C
C     CLOSE THE SCRTCH FILE AND THEN MERGE SCRTCH WITH AXIC
C     ON TO GEOM3
C
  470 I = 5
      ASSIGN 480 TO IRETRN
      GO TO 1380
C
C     OPEN SCRTCH FILE FOR INPUT AND LOCATE FORCE OR MOMENT CARDS ON
C     AXIC FILE.
C
  480 ASSIGN 490 TO IRETRN
      OP = INPRWD
      GO TO 1350
  490 CALL LOCATE (*1560,Z(IBUFF1),RECID(1),FLAG)
      IF (NOGO) GO TO 610
C
      CALL READ (*1600,*600,AXIC,Z(1),7,NOEOR,IAMT)
      CALL READ (*1610,*1610,SCRTCH,Z(8),7,NOEOR,IAMT)
C
  500 IF (Z(1) .LE. Z(8)) GO TO 510
C
      NFILE  = SCRTCH
      OUTBUF = 8
      GO TO 520
C
  510 NFILE  = AXIC
      OUTBUF = 1
C
  520 IF (NOGO) GO TO 610
      CALL WRITE (GEOM(3),Z(OUTBUF),7,NOEOR)
      CALL READ (*1620,*540,NFILE,Z(OUTBUF),7,NOEOR,IAMT)
      GO TO 500
C
C     OK ALL WORDS PROCESSED FOR FILE-NFILE
C
  540 IF (NFILE .EQ. AXIC) GO TO 550
      NFILE  = AXIC
      OUTBUF = 1
      GO TO 560
  550 NFILE  = SCRTCH
      OUTBUF = 8
  560 IF (NOGO) GO TO 610
      CALL WRITE (GEOM(3),Z(OUTBUF),7,NOEOR)
      CALL READ (*1620,*580,NFILE,Z(OUTBUF),7,NOEOR,IAMT)
      GO TO 560
C
C     CLOSE SCRTCH, WRITE EOR, AND PUT BITS IN TRAILER.
C
  580 I = 5
      ASSIGN 590 TO IRETRN
      GO TO 1380
  590 ASSIGN 610 TO IRETRN
      GO TO 1270
C
C     RECORD LENGTH ERROR
C
  600 REC(1) = RECID(1)
      REC(2) = RECID(2)
      REC(3) = RECID(3)
      ASSIGN 610 TO IERRTN
      GO TO 1490
C
  610 GO TO ICONT, (620,650)
C
C     GRAV CARD
C
  620 REC(1) = GRAV(1)
      REC(2) = GRAV(2)
      REC(3) = GRAV(3)
      ASSIGN 630 TO IRETRN
      GO TO 1230
C
C     LOAD CARD
C
  630 REC(1) = LOAD(1)
      REC(2) = LOAD(2)
      REC(3) = LOAD(3)
      ASSIGN 640 TO IRETRN
      GO TO 1230
C
C     MOMENT AND MOMAX CARDS
C
  640 RECID(1)  = MOMENT(1)
      RECID(2)  = MOMENT(2)
      RECID(3)  = MOMENT(3)
      RECIDX(1) = MOMAX(1)
      RECIDX(2) = MOMAX(2)
      RECIDX(3) = MOMAX(3)
      NCARD = 7
      ASSIGN 650 TO ICONT
      GO TO 290
C
C     PRESAX CARD
C
  650 CALL LOCATE (*722,Z(IBUFF1),PRESAX(1),FLAG)
C
C     RECORD HEADER FOR PRESAX CARDS IS FORMED HERE
C
      REC(1) = PRESAX(1)
      REC(2) = PRESAX(2)
      REC(3) = PRESAX(3)
      NCARD  = 13
      ASSIGN 660 TO IHEADR
      GO TO 1470
C
  660 CALL READ (*1600,*700,AXIC,Z(1),6,NOEOR,IAMT)
C
C     CREATE N+1 CARDS OF SAME LENGTH AS INPUT CARD.
C
C     CHECK RING ID-S IN FIELDS 3 AND 4 FOR PROPER SIZE.
C
C     CHECK FOR PIEZOELECTRIC
C
      PIEZ = .FALSE.
      IF (IPIEZ.EQ.1 .AND. Z(3).LT.0) PIEZ = .TRUE.
      IF (.NOT. PIEZ) GO TO 661
      Z(3) = -Z(3)
  661 CONTINUE
      NNN = Z(3)
      ASSIGN 670 TO IERRTN
      GO TO 1440
  670 NNN = Z(4)
      ASSIGN 680 TO IERRTN
      GO TO 1440
C
  680 DIFPHI = ABS(RZ(I6) - RZ(I5))
      DO 690 I = 1,NPLUS1
      Z(7) = I - 1
      Z(3) = Z(3) + 1000000
      IF (PIEZ) Z(3) = -Z(3)
      Z(4) = Z(4) + 1000000
      IF (NOGO) GO TO 690
      IF (DIFPHI .EQ. 0.0) GO TO 690
      IF (I.GT.1 .AND. ABS(DIFPHI-360.).LT.1.E-6) GO TO 690
      CALL WRITE (GEOM(3),Z(1),7,NOEOR)
      IF (PIEZ) Z(3) = -Z(3)
  690 CONTINUE
      GO TO 660
C
C     OUT OF PRESAX CARDS
C
  700 IF (IAMT) 710,720,710
C
C     CHECK FOR RECORD INCONSISTANCY ERROR.
C
  710 ASSIGN 722 TO IERRTN
      REC(1) = PRESAX(1)
      REC(2) = PRESAX(2)
      REC(3) = PRESAX(3)
      GO TO 1490
C
C     WRITE EOR AND PUT BITS IN TRAILER
C
  720 ASSIGN 722 TO IRETRN
      GO TO 1270
C
C     RFORCE CARD
C
 722  CALL LOCATE (*730,Z(IBUFF1),RFORCE(1),FLAG)
      REC(1) = RFORCE(1)
      REC(2) = RFORCE(2)
      REC(3) = RFORCE(3)
      NCARD  = 24
      ASSIGN 723 TO IHEADR
      GO TO 1470
C
C     PROCESS RFORCE DATA
C
  723 CALL READ (*1600,*725,AXIC,Z(1),7,NOEOR,IAMT)
      IF (Z(2).EQ.0 .AND. Z(3).EQ.0 .AND. Z(5).EQ.0 .AND. Z(6).EQ.0)
     1    GO TO 7240
      WRITE  (NOUT,724) UFM,Z(1)
  724 FORMAT (A23,' 336, RFORCE DATA IN SET NO.',I8,
     1       ' CONTAINS ILLEGAL DIRECTION FOR AXISYMMETRIC PROBLEM')
      NOGO = .TRUE.
      GO TO 723
 7240 Z(2) = 0
      Z(3) = 0
      Z(5) = 0
      Z(6) = Z(7)
      Z(7) = 0
      CALL WRITE (GEOM(3),Z(1),7,NOEOR)
      GO TO 723
C
C     END OF RFORCE CARDS
C
  725 IF (IAMT) 726,727,726
C
C     RECORD INCONSISTENCY ERROR
C
  726 ASSIGN 730 TO IERRTN
      REC(1) = RFORCE(1)
      REC(2) = RFORCE(2)
      REC(3) = RFORCE(3)
      GO TO 1490
C
C     WRITE EOR AND BITS IN TRAILER
C
  727 ASSIGN 730 TO IRETRN
      GO TO 1270
C
C     TEMPD CARD
C
  730 REC(1) = TEMPD(1)
      REC(2) = TEMPD(2)
      REC(3) = TEMPD(3)
      ASSIGN 740 TO IRETRN
      IF (NOGO) GO TO 740
      CALL LOCATE (*740,Z(IBUFF1),REC(1),FLAG)
      CALL WRITE (IFILE,REC(1),3,NOEOR)
      VEOR = 0
  735 CALL READ (*1600,*738,AXIC,Z(1),ICORE,NOEOR,IAMT)
      IAMT = ICORE
  736 DO 737 I = 1,IAMT,2
  737 Z(I) = Z(I) + 100000000
      CALL WRITE (IFILE,Z(1),IAMT,0)
      DO 739 I = 1,IAMT,2
  739 Z(I) = Z(I) + 100000000
      CALL WRITE (IFILE,Z(1),IAMT,VEOR)
      IF (VEOR) 1290,735,1290
  738 VEOR = 1
      GO TO 736
C
C     TEMPAX CARD
C
  740 CALL LOCATE (*1210,Z(IBUFF1),TEMPAX(1),FLAG)
C
C     RECORD HEADER ON GEOM3 FOR TEMP CARDS
C
      REC(1) = TEMP(1)
      REC(2) = TEMP(2)
      REC(3) = TEMP(3)
      NCARD  = 20
      ASSIGN 750 TO IHEADR
      GO TO 1470
C
C     AT 604(?) SET UP SCRATCH FILE.
C
  750 I    = 5
      BUFF = IBUFF3
      OP   = OUTRWD
      ASSIGN 760 TO IRETRN
      GO TO 1350
C
C     PICK UP FIRST TEMPAX CARD = 4 WORDS.
C
  760 LAST = 0
      CALL READ (*1600,*1200,AXIC,Z(1),4,NOEOR,IAMT)
  770 K = 0
      SETID = Z(1)
      RINGID= Z(2)
C
C     CHECK RING ID FOR PROPER RANGE OF VALUE
C
      NNN = RINGID
      ASSIGN 780 TO IERRTN
      GO TO 1440
C
  780 IAT = 3
  790 K   = K + 1
      IAT = IAT + 2
      ICRQ= IAT + 3 - ICORE
      IF (ICORE .LT. IAT+3) GO TO 1310
C
C     ALL TEMPAX CARDS HAVING SAME SET AND RING ID MUST BE ABLE TO
C     HAVE 2 WORDS EACH FIT IN CORE.
C
      Z(IAT  ) = Z(3)
      Z(IAT+1) = Z(4)
C
      CALL READ (*1600,*1130,AXIC,Z(1),4,NOEOR,IAMT)
C
C     DOES THIS CARD HAVE SAME SET AND RING ID AS LAST IN CURRENT SERIES
C
      IF (Z(1) .NE. SETID ) GO TO 800
      IF (Z(2) .NE. RINGID) GO TO 800
      GO TO 790
C
C     WE HAVE A  K X 2  ARRAY OF  PHI-S  AND T-S.
C
C     CONVERT ALL  PHIS SUCH THAT (0.LE. PHI .LT.TWOPI)
C
  800 IEND   = IAT + 1
      IBEGIN = 5
C
      DO 840 I = IBEGIN,IEND,2
      ANGLE = RZ(I)
      IF (ANGLE) 810,840,830
  810 IF (ANGLE) 820,840,840
  820 ANGLE = ANGLE + 360.0
      GO TO 810
C
  830 IF (ANGLE .LT. 360.0) GO TO 840
      ANGLE = ANGLE - 360.0
      GO TO 830
C
  840 RZ(I) = ANGLE*RADDEG
C
C     SIMPLE SORT FOR THE K X 2  MATRIX.
C     SORT IS PERFORMED ON COLUMN 1 ONLY
C
      IF (K .EQ. 1) GO TO 950
      ISTART = IBEGIN + 2
      DO 900 I = ISTART,IEND,2
      IAT = I - 2
      IF (RZ(I) .GE. RZ(IAT)) GO TO 900
C
C     ROW NOT HIGH ENOUGH.  MOVE IT UP.
C
  850 IAT = IAT - 2
      IF (IAT-IBEGIN) 870,870,860
  860 IF (RZ(I) .LT. RZ(IAT)) GO TO 850
      IAT = IAT + 2
      GO TO 880
  870 IAT = IBEGIN
C
C     THE ELEMENTS (I) AND (I+1) WILL BE MOVED UP TO POSITIONS (IAT) AND
C     (IAT+1) AND ELEMENTS (IAT) THRU (I-1) WILL BE  MOVED DOWN 1 ROW.
C
C     FIRST SAVE THE ROW BEING MOVED UP
C
  880 RZ(IEND+1) = RZ(I)
      RZ(IEND+2) = RZ(I+1)
      NMOVE = I - IAT
      IAT   = I + 2
      DO 890 J = 1,NMOVE
      IAT = IAT - 1
  890 RZ(IAT) = RZ(IAT-2)
C
C     REPLACE SAVE ROW IN NEW SLOT
C
      RZ(IAT-2) = RZ(IEND+1)
      RZ(IAT-1) = RZ(IEND+2)
C
  900 CONTINUE
C
C     CHECK FOR ANY DUPLICATE ANGLES AND REMOVE THEM...
C
      IBEGIN = IBEGIN + 2
  910 DO 920 I = IBEGIN,IEND,2
      IF (Z(I) .EQ. Z(I-2)) GO TO 930
  920 CONTINUE
      GO TO 950
C
C     DUPLICATE, SHRINK LIST UP OVER IT.
C
  930 IEND = IEND - 2
      K=K-1
      DO 940 J = I,IEND,2
      Z(J  ) = Z(J+2)
  940 Z(J+1) = Z(J+3)
      IBEGIN = I
      IF (IBEGIN - IEND) 910,950,950
C
C     SET UP K + 1  CARD
C
  950 RZ(IEND+1) = RZ(I5) + TWOPI
      RZ(IEND+2) = RZ(I6)
C
C     THERE ARE K CARDS NOW WITH SETID, AND RINGID, NOT INCLUDING THE
C     K + 1ST CARD
C
C     N+1 TEMP CARDS FOR S SET (PUT ON GEOM3)
C     N+1 TEMP CARDS FOR C SET (PUT ON SCRTCH FOR NOW)
C
C     NOTE FMMS-52  (10/04/67) PAGE -9- FOR FOLLOWING...
C
      CSSET = 1
      SETID = SETID + 100000000
C
C     CSSET = 0 FOR C-SET  AND NON-ZERO FOR S-SET.
C
      IBEGIN = K + K + 7
      ICRQ   = IBEGIN + 2 - ICORE
      IF ((IBEGIN+2) .GT. ICORE) GO TO 1310
C
  960 NADD = 0
      Z(IBEGIN) = SETID
      DO 1100 I = 1,NPLUS1
      NADD = NADD + 1000000
C
C     NI IS REAL
C
      NI   = I - 1
      NISQ = (I-1)**2
      Z(IBEGIN+1) = RINGID + NADD
      IPHI = 3
      IT   = 4
      SUM  = 0.0E0
      IF (NI   ) 1010,970,1010
  970 IF (CSSET) 1000,980,1000
  980 DO 990 IK = 1,K
      IPHI = IPHI + 2
      IT   = IT   + 2
  990 SUM  = SUM  + (RZ(IT)+RZ(IT+2))*(RZ(IPHI+2)-RZ(IPHI))
 1000 RZ(IBEGIN+2) = 0.25*SUM/PI
      GO TO 1060
C
C     NON-ZERO NI
C
 1010 IF (K .EQ. 1) GO TO 1050
      DO 1040 IK = 1,K
      IPHI = IPHI + 2
      IT   = IT   + 2
      NPHI = NI*RZ(IPHI  )
      NPHI1= NI*RZ(IPHI+2)
C
      IF (CSSET) 1030,1020,1030
C
C     C-SET
C
 1020 A1 =  SIN(NPHI1)
      A2 = -SIN(NPHI )
      A3 =  COS(NPHI1)
      A4 = -COS(NPHI )
      GO TO 1040
C
C     S-SET
C
 1030 A1 = -COS(NPHI1)
      A2 =  COS(NPHI )
      A3 =  SIN(NPHI1)
      A4 = -SIN(NPHI )
C
C
 1040 SUM = SUM + (((RZ(IT)*RZ(IPHI+2) - RZ(IT+2)*RZ(IPHI))*
     1      (A1 + A2)/NI) + ((RZ(IT+2) - RZ(IT))*
     2      (A3 + A4 + NPHI1*A1 + NPHI*A2)/NISQ))/
     3      (RZ(IPHI+2) - RZ(IPHI))
C
 1050 RZ(IBEGIN+2) = SUM/PI
C
 1060 IF (NOGO ) GO TO 1105
      IF (CSSET) 1070,1080,1070
 1070 NFILE = GEOM(3)
      GO TO 1090
 1080 NFILE = SCRTCH
 1090 CALL WRITE (NFILE,Z(IBEGIN),3,NOEOR)
 1100 CONTINUE
 1105 IF (CSSET) 1110,1120,1110
 1110 CSSET = 0
      SETID = SETID + 100000000
      GO TO 960
C
C     THIS SERIES OF TEMPAX CARDS COMPLETE GO FOR MORE IF LAST = 0
C
 1120 IF (LAST) 1140,770,1140
 1130 LAST = 1
      GO TO 800
C
C     ALL TEMPAX CARDS COMPLETE. CLOSE SCRATCH, OPEN SCRATCH
C     AND COPY SCRATCH TO GEOM3.
C
 1140 IF (NOGO) GO TO 1210
      CALL WRITE (SCRTCH,Z(1),0,EOR)
      CALL CLOSE (SCRTCH,CLORWD)
      CALL OPEN  (*1640,SCRTCH,Z(IBUFF3),INPRWD)
C
      VEOR = 0
 1150 CALL READ (*1610,*1170,SCRTCH,Z(1),ICORE,NOEOR,IAMT)
      IAMT = ICORE
 1160 CALL WRITE (GEOM(3),Z(1),IAMT,VEOR)
      IF (VEOR) 1180,1150,1180
 1170 VEOR = 1
      GO TO 1160
C
C     ALL  TEMPAX  CARDS  PROCESSED.
C
 1180 CALL CLOSE (SCRTCH,CLORWD)
C
C     PUT BITS IN TRAILER FOR TEMP CARDS WRITTEN
C
      REC(1) = TEMP(1)
      REC(2) = TEMP(2)
      REC(3) = TEMP(3)
      ASSIGN 1210 TO IRETRN
      GO TO 1290
C
C     RECORD LENGTH ERROR
C
 1200 REC(1) = TEMPAX(1)
      REC(2) = TEMPAX(2)
      REC(3) = TEMPAX(3)
      ASSIGN 1210 TO IERRTN
      GO TO 1490
C
C     CLOSE GEOM3
C
 1210 I = 3
      ASSIGN 1220 TO IRETRN
      GO TO 1380
C
C     CTRIAAX CARD
C
 1700 REC(1) = CTRIAA (1)
      REC(2) = CTRIAA (2)
      REC(3) = CTRIAA (3)
      NCARD  = 43
      CALL LOCATE (*1800,Z(IBUFF1),REC(1),FLAG)
C
C     RECORD HEADER FOR CTRIAAX
C
      ASSIGN 1710 TO IHEADR
      ICONB  = 2
      ICONSO = 1
      GO TO 1470
 1710 CALL READ (*1600,*1770,AXIC,Z(1),6,NOEOR,IAMT)
      Z(1) = Z(1)*1000
      DO 1720 I = 1,NPLUS1
      Z(1) = Z(1) + 1
      Z(3) = Z(3) + 1000000
      Z(4) = Z(4) + 1000000
      Z(5) = Z(5) + 1000000
      IF (NOGO) GO TO 1720
      CALL WRITE (GEOM(2),Z(1),6,NOEOR)
 1720 CONTINUE
      GO TO 1710
C
C     OUT OF CTRIAAX CARD
C
 1770 IF (IAMT) 1730,1740,1730
 1730 ASSIGN 260 TO IERRTN
      GO TO 1490
C
C     PUT BITS IN TRILER
C
 1740 ASSIGN 260 TO IRETRN
      GO TO 1270
 1800 IF (ICONSO .EQ. 1) GO TO 1740
      ASSIGN 260 TO IERRTN
C
C     MISSING REQUIRED CCONEAX OR CTRIAAX OR CTRAPAX CARD
C
      CALL PAGE2 (3)
      IMSG = 362
      WRITE (NOUT,185) UFM,IMSG
      WRITE (NOUT,1910) CDTYPE(3),CDTYPE(4),CDTYPE(43),CDTYPE(44),
     1                  CDTYPE(45),CDTYPE(46)
 1910 FORMAT (5X,'MINIMUM PROBLEM REQUIRES ',2A4,2H,  ,2A4,4H OR ,2A4,
     1       ' CARD.  NONE FOUND')
      NOGO = .TRUE.
      GO TO IERRTN, (260,240)
C
C     CTRAPAX CARD
C     ============
C
 1750 REC(1) = CTRAPA (1)
      REC(2) = CTRAPA (2)
      REC(3) = CTRAPA (3)
      CALL LOCATE (*1700,Z(IBUFF1),REC(1),FLAG)
      ICONB  = 1
C
C     RECORD HEADER FOR CTRAPAX
C
      ASSIGN 1751 TO IHEADR
      ICONSO = 1
      GO TO 1470
 1751 CALL READ (*1600,*1753,AXIC,Z(1),7,NOEOR,IAMT)
      Z(1) = Z(1)*1000
      DO 1752 I = 1,NPLUS1
      Z(1) = Z(1) + 1
      Z(3) = Z(3) + 1000000
      Z(4) = Z(4) + 1000000
      Z(5) = Z(5) + 1000000
      Z(6) = Z(6) + 1000000
      IF (NOGO) GO TO 1752
      CALL WRITE (GEOM(2),Z(1),7,NOEOR)
 1752 CONTINUE
      GO TO 1751
C
C     OUT OF CTRAPAX CARD
C
 1753 IF (IAMT) 1754,1755,1754
 1754 ASSIGN 260 TO IERRTN
      GO TO 1490
C
C     PUT BITS IN TRILER
C
 1755 ASSIGN 260 TO IRETRN
      IF (NOGO) GO TO 1300
      CALL WRITE (IFILE,Z(1),IAMT,EOR)
      I1 = (REC(2)-1)/16 + 2
      I2 = REC(2) - (I1-2)*16 + 16
      TRAIL (I1) = ORF(TRAIL(I1),TWO(I2))
      GO TO 1700
C
C     GEOM4 AND GEOM1 PROCESSING IS PERFORMED IN IFP3B ROUTINE
C                                                =====
C
 1220 CALL IFP3B
      GO TO 1570
C
C     UTILITY SECTION FOR IFP3
C     AXIS-SYMETRIC-CONICAL-SHELL DATA GENERATOR.
C     ==========================================
C
C     COMMON CODE FOR TRANSFER OF RECORD FROM AXIC FILE TO SOME
C     OTHER FILE
C
 1230 CALL LOCATE (*1300,Z(IBUFF1),REC(1),FLAG)
      IF (NOGO) GO TO 1300
      CALL WRITE (IFILE,REC(1),3,NOEOR)
 1260 CALL READ (*1600,*1270,AXIC,Z(1),ICORE,NOEOR,IAMT)
      IAMT = ICORE
      CALL WRITE (IFILE,Z(1),IAMT,NOEOR)
      GO TO 1260
 1270 IF (NOGO) GO TO 1300
      IF (IFILE .EQ. GEOM(3)) GO TO 1280
      IF (IFILE.EQ.GEOM(2) .AND. ICONB.EQ.1) GO TO 1300
 1280 CALL WRITE (IFILE,Z(1),IAMT,EOR)
C
C     PUT BITS IN TRAILER
C
 1290 I1 = (REC(2)-1)/16 + 2
      I2 =  REC(2) - (I1-2)*16 + 16
      TRAIL(I1) = ORF(TRAIL(I1),TWO(I2))
C
 1300 GO TO IRETRN, (250,260,610,630,640,722,730,740,1210)
C
C     OUT OF CORE
C
 1310 CALL PAGE2 (4)
      IMSG = 363
      WRITE  (NOUT, 185) IMSG
      WRITE  (NOUT,1320) ICRQ
 1320 FORMAT (5X,'INSUFFICIENT CORE TO PROCESS AXIC DATA IN SUBROUTINE',
     1       'IFP3', /5X,'ADDITIONAL CORE NEEDED =',I8,' WORDS.')
      NOGO = .TRUE.
C
C     GO TO FATAL ERROR RETURN
C
      GO TO 1570
C
C     AXIC FILE NOT IN FIST
C
 1330 CALL PAGE2 (3)
      IMSG = 1061
      WRITE (NOUT,55) SFM,IMSG
      WRITE (NOUT,60) AXIC,INAME(11),INAME(12),IFIST
      NOGO = .TRUE.
C
C     GO TO FATAL ERROR RETURN
C
      GO TO 1570
C
C     OPEN A FILE AND GET THE TRAILER
C
 1350 IF (NOGO) GO TO 1360
      CALL OPEN (*1370,FILE(I),Z(BUFF),OP)
      OPENFL(I) = 1
      IF (I .GT. 4) GO TO 1360
C
C     WRITE THE HEADER RECORD
C
      CALL WRITE (FILE(I),INAME(2*I-1),2,EOR)
      TRAIL(1) = FILE(I)
      CALL RDTRL (TRAIL(1))
C
 1360 GO TO IRETRN, (150,280,370,760,490)
C
 1370 CALL PAGE2 (3)
      IMSG = 1061
      WRITE (NOUT,55) SFM,IMSG
      WRITE (NOUT,60) FILE(I),INAME(2*I-1),INAME(2*I),IFIST
      NOGO = .TRUE.
      GO TO 1570
C
C     CLOSE A FILE
C
 1380 IF (OPENFL(I)) 1390,1410,1390
 1390 IF (I.GT. 4) GO TO 1400
      CALL WRITE (FILE(I),T65535(1),3,EOR)
 1400 CALL CLOSE (FILE(I),CLORWD)
      OPENFL(I) = 0
      IF (I .GT. 4) GO TO 1410
      CALL WRTTRL (TRAIL(1))
 1410 GO TO IRETRN, (270,590,1220,480)
C
C     HARMONIC NUMBER ... ON CARD TYPE ...... IS OUT OF RANGE 0 TO 998
C
 1420 IF (NNN.LT.999 .AND. NNN.GE.0 .AND. NNN.LE.N)
     1    GO TO IERRTN, (140,400,397)
      CALL PAGE2 (3)
      IMSG = 364
      WRITE  (NOUT,185 ) UFM,IMSG
      WRITE  (NOUT,1430) NNN,CDTYPE(2*NCARD-1),CDTYPE(2*NCARD),N
 1430 FORMAT (5X,'HARMONIC NUMBER ',I6,4H ON ,2A4,' CARD OUT OF 0 TO ',
     1        I4,' ALLOWABLE RANGE.')
      NOGO = .TRUE.
      GO TO IERRTN, (140,400,397)
C
C     RING ID OUT OF PERMISSABLE RANGE OF 1 TO 999999
C
 1440 IF (NNN.GT.0 .AND. NNN.LE.999999)
     1    GO TO IERRTN, (170,180,390,670,680,780)
      CALL PAGE2 (3)
      IMSG = 365
      WRITE  (NOUT,185 ) UFM,IMSG
      WRITE  (NOUT,1450) NNN,CDTYPE(2*NCARD-1),CDTYPE(2*NCARD)
 1450 FORMAT (5X,'RING ID',I10,4H ON ,2A4,' CARD OUT OF 1 TO 999999 ',
     1       'ALLOWABLE RANGE')
      NOGO = .TRUE.
      GO TO IERRTN, (170,180,390,670,680,780)
C
C     CHECK BIT-IBIT IN TRAILER AND RETURN NON = ZERO OR NON-ZERO...
C
 1460 I1 = (IBIT-1)/16 + 2
      I2 = IBIT - (I1-2)*16 + 16
      NON = ANDF(AXTRL(I1),TWO(I2))
      GO TO IBITR, (300,310)
C
C     WRITE 3 WORD RECORD HEADER
C
 1470 IF (NOGO) GO TO 1480
      CALL WRITE (IFILE,REC(1),3,NOEOR)
 1480 GO TO IHEADR, (160,350,660,723,750,1710,1751)
C
C     END-OF-RECORD ON AXIC FILE
C
 1490 CALL PAGE2 (3)
      IMSG = 1063
      WRITE  (NOUT,55) SFM,IMSG
      WRITE  (NOUT,1500) CDTYPE(2*NCARD-1),CDTYPE(2*NCARD)
 1500 FORMAT (5X,'EOR ON AXIC FILE WHILE READING ',2A4,'CARD RECORDS.')
      NOGO = .TRUE.
      GO TO IERRTN, (260,460,610,722,730,1210)
C
C     MISSING REQUIRED CARD
C
 1510 CALL PAGE2 (3)
      IMSG = 362
      WRITE  (NOUT,185 ) UFM,IMSG
      WRITE  (NOUT,1520) CDTYPE(2*NCARD-1),CDTYPE(2*NCARD)
 1520 FORMAT (5X,'MINIMUM PROBLEM REQUIRES ',2A4,' CARD.  NONE FOUND.')
      NOGO = .TRUE.
      GO TO IERRTN, (260,140)
C
C     AXIC TRAILER BIT ON BUT CAN NOT LOCATE RECORD
C
 1530 CALL PAGE2 (3)
      IMSG = 1064
      WRITE  (NOUT,55) SFM,IMSG
      WRITE  (NOUT,1540) CDTYPE(2*NCARD-1),CDTYPE(2*NCARD)
 1540 FORMAT (5X,2A4,' CARD COULD NOT BE LOCATED ON AXIC FILE AS ',
     1       'EXPECTED')
 1550 NOGO = .TRUE.
      GO TO 610
 1560 CALL PAGE2 (2)
      WRITE (NOUT,1540) RECID(1),RECID(2),RECID(3)
      GO TO 1550
C
C     CLOSE ANY OPEN FILES AND RETURN
C
 1570 DO 1590 I = 1,6
      IF (OPENFL(I)) 1580,1590,1580
 1580 CALL CLOSE (FILE(I),CLORWD)
      OPENFL(I) = 0
 1590 CONTINUE
      IF (NOGO) NOFLAG = 32767
      CALL CONMSG (MSG2,2,0)
      RETURN
C
C     EOF ENCOUNTERED READING AXIC FILE.
C
 1600 NFILE = AXIC
      IN  = 11
      IN1 = 12
      GO TO 1620
 1610 NFILE = SCRTCH
      IN  = 9
      IN1 = 10
 1620 CALL PAGE2 (3)
      IMSG = 3002
      WRITE  (NOUT,55) SFM,IMSG
      WRITE  (NOUT,1630) INAME(IN),INAME(IN1),NFILE
 1630 FORMAT (5X,'EOF ENCOUNTERED WHILE READING DATA SET ',2A4,' (FILE',
     1        I4,') IN SUBROUTINE IFP3')
      NOGO = .TRUE.
      GO TO 1570
C
 1640 I = 5
      GO TO 1370
      END


text:
c  File: DRDTGFI1.[F|FOR] CONTAINS DRDTGFI1 AND DTGFCN.
      program DRDTGFI1
c>> 2001-07-16 DRDTGFI1 Krogh Added exponent 0 to some constants.
c>> 2001-05-22 DRDTGFI1 Krogh Minor change for making .f90 version.
c>> 1997-07-01 DRDTGFI1 Krogh Reversed subscripts in B (CLL suggestion)
c>> 1997-06-19 DRDTGFI1 Krogh Minor changes to get proper C conversion.
c>> 1997-06-18 DRDTGFI1 CLL
c>> 1996-03-04 DRDTGFI1 CLL
c>> 1996-02-02 DRDTGFI1 CLL
c>> 1995-10-31 DRDTGFI1 CLL
c   Demo driver for DTGFI, DTGGRD, DTGPD, etc.
c     ------------------------------------------------------------------
c--D replaces "?": DR?TGFI1, ?TGFI, ?TGGRD, ?TGPD, ?TGFCN
c     ------------------------------------------------------------------
      integer mb, mp, meval, mt
      parameter(mb= 16, mp= 28, mt= 336)
      parameter(meval = 5)
      integer Bdry(4,mb), i, ieval, INFO(3), IP(mp), kf
      integer mode, ncont, np, nt, TRIANG(mt)
      double precision DZ(2,mp), dzout(2)
      double precision dztemp(2), DZTRUE(2)
      double precision q(2)
      double precision savwrk(28)
      double precision W(mp), X(mp)
      double precision Y(mp)
      double precision Z(mp), zout, ztrue
      logical  WANTDZ
      character title*75
      data X /
     *-0.76059D0,-0.02286D0,-0.44790D0, 0.15068D0,-0.87287D0,-0.23390D0,
     * 0.06093D0,-0.84142D0,-0.69173D0,-0.56613D0,-0.42243D0,-0.17249D0,
     *-0.06484D0, 0.48286D0,-0.88784D0, 0.10277D0, 0.69087D0,-0.84292D0,
     * 0.08784D0,-0.95068D0, 0.02496D0, 0.94973D0, 0.04588D0,-0.51667D0,
     *-0.77561D0, 0.90000D0,-0.70830D0,-0.40500D0 /
      data Y /
     *-0.31421D0, 0.75657D0, 0.14321D0, 0.42353D0,-0.62983D0, 0.12326D0,
     * 0.34054D0,-0.55144D0, 0.34158D0,-0.67143D0, 0.43087D0, 0.96081D0,
     *-0.68130D0, 0.12095D0,-0.10663D0, 0.29219D0, 0.31028D0, 0.12934D0,
     * 0.10709D0,-0.42307D0, 0.49895D0, 0.68597D0,-0.78215D0,-0.12362D0,
     *-0.88827D0,-0.60000D0,-0.88620D0, 0.08600D0 /
      data np / mp /
      data wantdz / .true. /
c     ------------------------------------------------------------------
      write(*,'(a/a)') ' Program DRDTGFI1.  Demo driver for DTGFI,',
     *                 '                    DTGGRD, DTGPD, etc.'
      call DTGGRD(X,Y,NP,IP,W,TRIANG,MT,Bdry,MB,NT, INFO)
      if (INFO(1) .NE. 0) then
         write(*,'(a,i5)')
     *      ' ERROR return FROM DTGGRD. INFO(1) =', INFO(1)
         stop
      endif
c
c     do 200 ncont = 0,1
      do 200 ncont = 1,1
c     do 100 KF= 0,4
      do 100 KF= 2,4,2
         write(*,'(//a,i1,a)') ' New Case.  Interpolation with C',
     *      ncont,' continuity.'
         do 20 I=1,NP
               call DTGFCN(KF,X(I),Y(I),title,Z(I),DZtemp(1),DZtemp(2))
   20    continue
         write(*,'(/a)') title
         call DTGPD(X,Y,Z,DZ,NP, TRIANG, NT, IP)
         savwrk(1) = 0.0d0
         write(*,'(/3x,a//a,a/16x,a/16x,a/)')
     *       ' VALUES AND PARTIAL DERIVS INTERPOLATED ALONG A PATH.',
     *       '    I    X     Y',
     *       '    Z_INTERP     Z_TRUE     Z_ERR',
     *       '  DZ1_INTERP   DZ1_TRUE   DZ1_ERR',
     *       '  DZ2_INTERP   DZ2_TRUE   DZ2_ERR'
         do 30 ieval = 0, meval-1
            q(1) = -1.0d0 + ieval*2.0d0/(meval-1)
            q(2) =  0.8d0 * q(1)
            call DTGFI( X, Y, Z, DZ, TRIANG, NT, Bdry, Mb,
     *                 ncont, Q, zout, WANTDZ, DZout, MODE, SAVWRK)
            if(mode .ge. 0) then
               call DTGFCN(KF,  q(1), q(2),
     *             title, ztrue, dztrue(1), dztrue(2))
               write(*, '(1X,I4,2F6.2,2F11.6,E10.2)') ieval,q,
     *              zout, ztrue, zout - ztrue
               write(*,'(17x,2F11.6,E10.2)') dzout(1), dztrue(1),
     *              dzout(1) - dztrue(1)
               write(*,'(17x,2F11.6,E10.2)') dzout(2), dztrue(2),
     *              dzout(2) - dztrue(2)
            else
               write(*,'(1x,i4,2f6.2,a)')
     *            ieval, q, '   Error.'
            endif
   30    continue
  100 continue
  200 continue
      end
c     ==================================================================
      subroutine DTGFCN  (KF,X,Y, TITLE, Z,ZX,ZY)
c>> 1995-09-26 DRDTGFI1 CLL Editing for inclusion into MATH77.
c     C.L.Lawson, JPL, 1976 Dec 10.  Edited comments 1979 Mar 3.
c     This subr evaluates a function and its first partial derivs as
c     selected by KF.  KF can be from 0 to 4.
c     Input is KF, X, and Y.  Output is TITLE, Z, ZX, and ZY.
c     ------------------------------------------------------------------
      integer KF
      double precision X,Y,Z,ZX,ZY
      character title*75
      character title1(0:4)*56, title2(0:4)*19
      integer I
      data (title1(i),title2(i),i=0,4) /
     * '  CONSTANT FUNCTION     Z = 2', ' ',
     * '  LINEAR FUNCTION       Z = ( 1 + 2*X + 3*Y ) / 6', ' ',
     * '  QUADRATIC FUNCTION    Z = ( -1 + 2*X - 3*Y + 4*X**2 - ',
     *                             'X*Y + 9*Y**2 ) / 10',
     * '  CUBIC FUNCTION    Z = ( 9*X**3 - 2*(X**2)*Y + 3*X*Y**2',
     *                             ' - 4 * Y**3 ) / 10',
     * '  EXPONENTIAL FUNCTION   Z = EXP( -2 * (X**2 + Y**2) )', ' '/
c     ------------------------------------------------------------------
      title = title1(kf) // title2(kf)
       if(kf .eq. 0) then
c                                  KF=0  CONSTANT FCN.
         Z=2.0d0
         ZX=0.0d0
         ZY=0.0d0
       elseif(kf .eq. 1) then
c                                  KF=1  LINEAR FCN.
         Z=(1.0d0 + 2.0d0*X + 3.0d0*Y) / 6.0d0
         ZX=2.0d0/6.0d0
         ZY=3.0d0/6.0d0
       elseif(kf .eq. 2) then
c                                  KF=2  QUADRATIC FCN.
         Z=(-1.0d0 + 2.0d0*X - 3.0d0*Y +4.0d0*X**2 -X*Y + 9.0d0*Y**2)*
     *       0.1d0
         ZX =(2.0d0 + 8.0d0*X -Y)* 0.1d0
         ZY =(-3.0d0 -X + 18.0d0*Y)* 0.1d0
       elseif(kf .eq. 3) then
c                                  KF=3  CUBIC FCN.
         Z=(9.0d0*X**3 - 2.0d0*(X**2)*Y + 3.0d0*X*Y**2 - 4.0d0*Y**3)*
     *       0.1d0
         ZX=(27.0d0*X**2 - 4.0d0*X*Y + 3.0d0*Y**2)* 0.1d0
         ZY=(-2.0d0*X**2 + 6.0d0*X*Y - 12.0d0*Y**2)* 0.1d0
       elseif(kf .eq. 4) then
c                                  KF=4  EXPONENTIAL FCN.
         Z=  EXP(-(X**2 + Y**2) * 2.0d0)
c
c              NOTE THAT THE INFLECTION POINT OF THIS FCN IN ANY
c              RADIAL DIRECTION FROM THE ORIGIN IS AT R = .5
c
         ZX= -4.0d0*X*Z
         ZY= -4.0d0*Y*Z
       endif
      return
      end