text:
      SUBROUTINE MB02CY( TYPET, STRUCG, P, Q, N, K, A, LDA, B, LDB, H,
     $                   LDH, CS, LCS, DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To apply the transformations created by the SLICOT Library
C     routine MB02CX on other columns / rows of the generator,
C     contained in the arrays A and B of positive and negative
C     generators, respectively.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     TYPET   CHARACTER*1
C             Specifies the type of the generator, as follows:
C             = 'R':  A and B are additional columns of the generator;
C             = 'C':  A and B are additional rows of the generator.
C             Note:   in the sequel, the notation x / y means that
C                     x corresponds to TYPET = 'R' and y corresponds to
C                     TYPET = 'C'.
C
C     STRUCG  CHARACTER*1
C             Information about the structure of the two generators,
C             as follows:
C             = 'T':  the trailing block of the positive generator
C                     is lower / upper triangular, and the trailing
C                     block of the negative generator is zero;
C             = 'N':  no special structure to mention.
C
C     Input/Output Parameters
C
C     P       (input)  INTEGER
C             The number of rows / columns in A containing the positive
C             generators.  P >= 0.
C
C     Q       (input)  INTEGER
C             The number of rows / columns in B containing the negative
C             generators.  Q >= 0.
C
C     N       (input)  INTEGER
C             The number of columns / rows in A and B to be processed.
C             N >= 0.
C
C     K       (input)  INTEGER
C             The number of columns / rows in H.  P >= K >= 0.
C
C     A       (input/output)  DOUBLE PRECISION array, dimension
C             (LDA, N) / (LDA, P)
C             On entry, the leading P-by-N / N-by-P part of this array
C             must contain the positive part of the generator.
C             On exit, the leading P-by-N / N-by-P part of this array
C             contains the transformed positive part of the generator.
C
C     LDA     INTEGER
C             The leading dimension of the array A.
C             LDA >= MAX(1,P),    if TYPET = 'R';
C             LDA >= MAX(1,N),    if TYPET = 'C'.
C
C     B       (input/output)  DOUBLE PRECISION array, dimension
C             (LDB, N) / (LDB, Q)
C             On entry, the leading Q-by-N / N-by-Q part of this array
C             must contain the negative part of the generator.
C             On exit, the leading Q-by-N / N-by-Q part of this array
C             contains the transformed negative part of the generator.
C
C     LDB     INTEGER
C             The leading dimension of the array B.
C             LDB >= MAX(1,Q),    if TYPET = 'R';
C             LDB >= MAX(1,N),    if TYPET = 'C'.
C
C     H       (input)  DOUBLE PRECISION array, dimension
C             (LDH, K) / (LDH, Q)
C             The leading Q-by-K / K-by-Q part of this array must
C             contain part of the necessary information for the
C             Householder transformations computed by SLICOT Library
C             routine MB02CX.
C
C     LDH     INTEGER
C             The leading dimension of the array H.
C             LDH >= MAX(1,Q),    if TYPET = 'R';
C             LDH >= MAX(1,K),    if TYPET = 'C'.
C
C     CS      (input)  DOUBLE PRECISION array, dimension (LCS)
C             The leading 2*K + MIN(K,Q) part of this array must
C             contain the necessary information for modified hyperbolic
C             rotations and the scalar factors of the Householder
C             transformations computed by SLICOT Library routine MB02CX.
C
C     LCS     INTEGER
C             The length of the array CS.  LCS >= 2*K + MIN(K,Q).
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if  INFO = 0,  DWORK(1)  returns the optimal
C             value of LDWORK.
C             On exit, if  INFO = -16,  DWORK(1)  returns the minimum
C             value of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.  LDWORK >= MAX(1,N).
C             For optimum performance LDWORK should be larger.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  succesful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The Householder transformations and modified hyperbolic rotations
C     computed by SLICOT Library routine MB02CX are applied to the
C     corresponding parts of the matrices A and B.
C
C     CONTRIBUTOR
C
C     D. Kressner, Technical Univ. Chemnitz, Germany, June 2000.
C
C     REVISIONS
C
C     V. Sima, Research Institute for Informatics, Bucharest, July 2000,
C     February 2004, March 2007, Aug. 2011.
C
C     KEYWORDS
C
C     Elementary matrix operations, Householder transformation, matrix
C     operations, Toeplitz matrix.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, ONE
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      INTEGER           INFO, K, LDA, LDB, LCS, LDH, LDWORK, N, P, Q
      CHARACTER         STRUCG, TYPET
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LDA, *), B(LDB, *), CS(*), DWORK(*), H(LDH,*)
C     .. Local Scalars ..
      LOGICAL           ISLWR, ISROW
      INTEGER           I, IERR, CI, MAXWRK
      DOUBLE PRECISION  C, S, TAU
C     .. External Functions ..
      LOGICAL           LSAME
      EXTERNAL          LSAME
C     .. External Subroutines ..
      EXTERNAL          DAXPY, DLARF, DLASET, DORMLQ, DORMQR, DSCAL,
     $                  XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         MAX, MIN
C
C     .. Executable Statements ..
C
C     Decode the scalar input parameters.
C
      INFO  = 0
      ISROW = LSAME( TYPET,  'R' )
      ISLWR = LSAME( STRUCG, 'T' )
C
C     Check the scalar input parameters.
C
      IF ( .NOT.( ISROW .OR. LSAME( TYPET, 'C' ) ) ) THEN
         INFO = -1
      ELSE IF ( .NOT.( ISLWR .OR. LSAME( STRUCG, 'N' ) ) ) THEN
         INFO = -2
      ELSE IF ( P.LT.0 ) THEN
         INFO = -3
      ELSE IF ( Q.LT.0 ) THEN
         INFO = -4
      ELSE IF ( N.LT.0 ) THEN
         INFO = -5
      ELSE IF ( K.LT.0 .OR. K.GT.P ) THEN
         INFO = -6
      ELSE IF ( LDA.LT.1 .OR. ( ISROW .AND. LDA.LT.P ) .OR.
     $                   ( .NOT.ISROW .AND. LDA.LT.N ) ) THEN
         INFO = -8
      ELSE IF ( LDB.LT.1 .OR. ( ISROW .AND. LDB.LT.Q ) .OR.
     $                   ( .NOT.ISROW .AND. LDB.LT.N ) ) THEN
         INFO = -10
      ELSE IF ( LDH.LT.1 .OR. ( ISROW .AND. LDH.LT.Q ) .OR.
     $                   ( .NOT.ISROW .AND. LDH.LT.K ) ) THEN
         INFO = -12
      ELSE IF ( LCS.LT.2*K + MIN( K, Q ) ) THEN
         INFO = -14
      ELSE IF ( LDWORK.LT.MAX( 1, N ) ) THEN
         DWORK(1) = MAX( 1, N )
         INFO = -16
      END IF
C
C     Return if there were illegal values.
C
      IF ( INFO.NE.0 ) THEN
         CALL XERBLA( 'MB02CY', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( MIN( N, K, Q ).EQ.0 ) THEN
         DWORK(1) = ONE
         RETURN
      END IF
C
C     Applying the transformations.
C
      IF ( ISROW ) THEN
C
C        The generator is row wise stored.
C
         IF ( ISLWR ) THEN
            MAXWRK = 1
C
            DO 10  I = 1, K
C
C              Apply Householder transformation avoiding touching of
C              zero blocks.
C
               CI  = N - K + I - 1
               TAU = H(1,I)
               H(1,I) = ONE
               CALL DLARF( 'Left', MIN( I, Q ), CI, H(1,I), 1, TAU, B,
     $                     LDB, DWORK )
               H(1,I) = TAU
C
C              Now apply the hyperbolic rotation under the assumption
C              that A(I, N-K+I+1:N) and B(1, N-K+I:N) are zero.
C
               C = CS(I*2-1)
               S = CS(I*2)
C
               CALL DSCAL( CI, ONE/C, A(I,1), LDA )
               CALL DAXPY( CI,  -S/C, B(1,1), LDB, A(I,1), LDA )
               CALL DSCAL( CI,     C, B(1,1), LDB )
               CALL DAXPY( CI,    -S, A(I,1), LDA, B(1,1), LDB )
C
               B(1,N-K+I) =  -S/C * A(I,N-K+I)
               A(I,N-K+I) = ONE/C * A(I,N-K+I)
C
C              All below B(1,N-K+I) should be zero.
C
               IF( Q.GT.1 )
     $            CALL DLASET( 'All', Q-1, 1, ZERO, ZERO, B(2,N-K+I),
     $                         1 )
   10       CONTINUE
C
         ELSE
C
C           Apply the QR reduction on B.
C
            CALL DORMQR( 'Left', 'Transpose', Q, N, MIN( K, Q ), H,
     $                   LDH, CS(2*K+1), B, LDB, DWORK, LDWORK, IERR )
            MAXWRK = DWORK(1)
C
            DO 20  I = 1, K
C
C              Apply Householder transformation.
C
               TAU = H(1,I)
               H(1,I) = ONE
               CALL DLARF( 'Left', MIN( I, Q ), N, H(1,I), 1, TAU, B,
     $                     LDB, DWORK )
               H(1,I) = TAU
C
C              Apply Hyperbolic Rotation.
C
               C = CS(I*2-1)
               S = CS(I*2)
C
               CALL DSCAL( N, ONE/C, A(I,1), LDA )
               CALL DAXPY( N,  -S/C, B(1,1), LDB, A(I,1), LDA )
               CALL DSCAL( N,     C, B(1,1), LDB )
               CALL DAXPY( N,    -S, A(I,1), LDA, B(1,1), LDB )
   20       CONTINUE
C
         END IF
C
      ELSE
C
C        The generator is column wise stored.
C
         IF ( ISLWR ) THEN
            MAXWRK = 1
C
            DO 30  I = 1, K
C
C              Apply Householder transformation avoiding touching zeros.
C
               CI  = N - K + I - 1
               TAU = H(I,1)
               H(I,1) = ONE
               CALL DLARF( 'Right', CI, MIN( I, Q ), H(I,1), LDH, TAU,
     $                     B, LDB, DWORK )
               H(I,1) = TAU
C
C              Apply Hyperbolic Rotation.
C
               C = CS(I*2-1)
               S = CS(I*2)
C
               CALL DSCAL( CI, ONE/C, A(1,I), 1 )
               CALL DAXPY( CI,  -S/C, B(1,1), 1, A(1,I), 1 )
               CALL DSCAL( CI,     C, B(1,1), 1 )
               CALL DAXPY( CI,    -S, A(1,I), 1, B(1,1), 1 )
C
               B(N-K+I,1) =  -S/C * A(N-K+I,I)
               A(N-K+I,I) = ONE/C * A(N-K+I,I)
C
C              All elements right behind B(N-K+I,1) should be zero.
C
               IF( Q.GT.1 )
     $            CALL DLASET( 'All', 1, Q-1, ZERO, ZERO, B(N-K+I,2),
     $                         LDB )
   30       CONTINUE
C
         ELSE
C
C           Apply the LQ reduction on B.
C
            CALL DORMLQ( 'Right', 'Transpose', N, Q, MIN( K, Q ), H,
     $                   LDH, CS(2*K+1), B, LDB, DWORK, LDWORK, IERR )
            MAXWRK = DWORK(1)
C
            DO 40  I = 1, K
C
C              Apply Householder transformation.
C
               TAU = H(I,1)
               H(I,1) = ONE
               CALL DLARF( 'Right', N, MIN( I, Q ), H(I,1), LDH, TAU, B,
     $                     LDB, DWORK )
               H(I,1) = TAU
C
C              Apply Hyperbolic Rotation.
C
               C = CS(I*2-1)
               S = CS(I*2)
C
               CALL DSCAL( N, ONE/C, A(1,I), 1 )
               CALL DAXPY( N,  -S/C, B(1,1), 1, A(1,I), 1 )
               CALL DSCAL( N,     C, B(1,1), 1 )
               CALL DAXPY( N,    -S, A(1,I), 1, B(1,1), 1 )
   40       CONTINUE
C
         END IF
C
      END IF
C
      DWORK(1) = MAX( MAXWRK, N )
C
      RETURN
C
C *** Last line of MB02CY ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
      function RD(m_N, tau, taup, term, j_chi, coupling_Array)
      implicit none
      double precision :: RD, m_N, j_chi
      double precision :: coupling_Array(14,2)
      integer :: tau,taup
      integer :: term
      integer :: c, q2
      double precision :: c5, c5p
      double precision :: c8, c8p
c
      c = 0
      q2 = 2
c
      c5 = coupling_Array(4,tau)
      c5p = coupling_Array(4,taup)
      c8 = coupling_Array(7,tau)
      c8p = coupling_Array(7,taup)
c
      if (term.eq.c) then
         RD = (j_chi*(j_chi+1))/3. * c8*c8p
      else if (term.eq.q2) then
         RD = 1/m_N**2 * (j_chi*(j_chi+1))/3. * c5*c5p
      else
         RD = 0.
      end if
c
      end function RD

ext: f


text:
      SUBROUTINE EKTRMD (NTYPE)
C
C     THIS SUBROUTINE CALCULATES THE STIFFNESS MATRIX FOR THE
C     TRIANGULAR MEMBRANE ELEMENT
C
C     DOUBLE PRECISION VERSION
C
C     ECPT LIST
C                                                      IN
C                                                      THIS
C       ECPT       DESCRIPTION                         ROUTINE   TYPE
C     ******************************************************************
C       ECPT( 1) = ELEMENT ID                          NECPT(1)  INTEGER
C       ECPT( 2) = GRID POINT A                        NGRID(1)  INTEGER
C       ECPT( 3) = GRID POINT B                        NGRID(2)  INTEGER
C       ECPT( 4) = GRID POINT C                        NGRID(3)  INTEGER
C       ECPT( 5) = THETA = ANGLE OF MATERIAL           ANGLE     REAL
C       ECPT( 6) = MATERIAL ID                         MATID     INTEGER
C       ECPT( 7) = T                                   T         REAL
C       ECPT( 8) = NON-STRUCTURAL MASS                 FMU       REAL
C       ECPT( 9) = COORD. SYSTEM ID 1                  NECPT(9)  INTEGER
C       ECPT(10) = X1                                  X1        REAL
C       ECPT(11) = Y1                                  Y1        REAL
C       ECPT(12) = Z1                                  Z1        REAL
C       ECPT(13) = COORD. SYSTEM ID 2                  NECPT(13) INTEGER
C       ECPT(14) = X2                                  X2        REAL
C       ECPT(15) = Y2                                  Y2        REAL
C       ECPT(16) = Z2                                  Z2        REAL
C       ECPT(17) = COORD. SYSTEM ID 3                  NECPT(17) INTEGER
C       ECPT(18) = X3                                  X3        REAL
C       ECPT(19) = Y3                                  Y3        REAL
C       ECPT(20) = Z3                                  Z3        REAL
C       ECPT(21) = ELEMENT TEMPERATURE                 ELTEMP    REAL
C
C     IF NTYPE = 0  COMPLETE MEMBRANE COMPUTATION IS PERFORMED
C     IF NTYPE = 1  9  3X3 MATRICES FOR THE GRID POINTS  IN ECPT
C
C
      LOGICAL          NOGO,HEAT
      INTEGER          NECPT(21)
      REAL             ECPT(21),MATBUF
      DOUBLE PRECISION G(9),C(18),TT(2),TI(9),TEMPAR(27),E,K,KOUT,KIJ,
     1                 A,PROD9,TEMP9,XSUBB,XSUBC,YSUBC,DICT5
      DOUBLE PRECISION TEMP,VOL,REELMU,FLAMDA,DELTA
      COMMON /EMGTRX/  A(225),PROD9(9),TEMP9(9),XSUBB,XSUBC,YSUBC,DICT5,
     1                 E(18),K(324),KOUT(324),KIJ(81)
      COMMON /SYSTEM/  KSYSTM(60)
      COMMON /CONDAS/  CONSTS(5)
      COMMON /EMGPRM/  DUM(19),NOGO, HEAT
      COMMON /EMGEST/  MECPT(1),NGRID(3),ANGLE,MATID1,T,FMU,DUMMY1,X1,
     1                 Y1,Z1,DUMMY2,X2,Y2,Z2,DUMMY3,X3,Y3,Z3,DUMB(80)
      COMMON /MATIN /  MATID,INFLAG,ELTEMP,STRESS,SINTH,COSTH
      COMMON /MATOUT/  G11,G12,G13,G22,G23,G33,RHO,ALPHA1,ALPHA2,ALP12,
     1                 T SUB 0,G SUB E,SIGTEN,SIGCOM,SIGSHE,
     2                 G2X211,G2X212,G2X222
      COMMON /HMTOUT/  MATBUF(4)
      EQUIVALENCE      (CONSTS(4),DEGRA),(ECPT(1),NECPT(1),MECPT(1)),
     1                 (KSYSTM(2),IOUTPT),(KSYSTM(56),IHEAT)
C
C     SET UP THE E MATRIX WHICH IS (3X2) FOR THE TRI-MEMBRANE
C
C     E(1), E(3), E(5) WILL BE THE I-VECTOR
C     E(2), E(4), E(6) WILL BE THE J-VECTOR
C     E(7), E(8), E(9) WILL BE THE K-VECTOR NOT USED IN E FOR MEMBRANE
C
C     FIRST FIND I-VECTOR = RSUBB - RSUBA  (NON-NORMALIZED)
C
      E(1) = DBLE(X2) - DBLE(X1)
      E(3) = DBLE(Y2) - DBLE(Y1)
      E(5) = DBLE(Z2) - DBLE(Z1)
C
C     NOW FIND LENGTH = X-SUB-B   COORD. IN ELEMENT SYSTEM
C
      XSUBB = DSQRT(E(1)**2 + E(3)**2 + E(5)**2)
      IF (XSUBB .LE. 1.D-6) GO TO 7770
C
C  20 NOW NORMALIZE I-VECTOR WITH X-SUB-B
C
      E(1) = E(1)/XSUBB
      E(3) = E(3)/XSUBB
      E(5) = E(5)/XSUBB
C
C     HERE WE NOW TAKE RSUBC - RSUBA AND STORE TEMPORARILY IN
C     E(2), E(4), E(6) WHICH IS WHERE THE J-VECTOR WILL FIT LATER
C
      E(2) = DBLE(X3) - DBLE(X1)
      E(4) = DBLE(Y3) - DBLE(Y1)
      E(6) = DBLE(Z3) - DBLE(Z1)
C
C     X-SUB-C  =  I . (RSUBC - RSUBA), THUS
C
      XSUBC = E(1)*E(2) + E(3)*E(4) + E(5)*E(6)
C
C     AND CROSSING THE I-VECTOR TO (RSUBC-RSUBA) GIVES THE K-VECTOR
C     (NON-NORMALIZED)
C
      E(7) = E(3)*E(6) - E(5)*E(4)
      E(8) = E(5)*E(2) - E(1)*E(6)
      E(9) = E(1)*E(4) - E(3)*E(2)
C
C
C     THE LENGTH OF THE K-VECTOR IS NOW FOUND AND EQUALS Y-SUB-C
C     COORD. IN ELEMENT SYSTEM
C
      YSUBC = DSQRT(E(7)**2 + E(8)**2 + E(9)**2)
      IF (YSUBC .LE. 1.D-6) GO TO 7780
C
C  25 NOW NORMALIZE K-VECTOR WITH YSUBC JUST FOUND
C
      E(7) = E(7)/YSUBC
      E(8) = E(8)/YSUBC
      E(9) = E(9)/YSUBC
C
C     J VECTOR = K CROSS I
C     STORE IN THE SPOT FOR J
C
      E(2) = E(5)*E(8) - E(3)*E(9)
      E(4) = E(1)*E(9) - E(5)*E(7)
      E(6) = E(3)*E(7) - E(1)*E(8)
C
C     AND JUST FOR COMPUTER EXACTNESS NORMALIZE J-VECTOR TO MAKE SURE.
C
      TEMP = DSQRT (E(2)**2 + E(4)**2 + E(6)**2)
      IF (TEMP  .LE. 0.D0) GO TO 7790
      E(2) = E(2)/TEMP
      E(4) = E(4)/TEMP
      E(6) = E(6)/TEMP
C
C     VOLUME OF ELEMENT, THETA, MU, LAMDA, AND DELTA
C
      VOL    = XSUBB*YSUBC*T/2.D0
      REELMU = 1.D0/XSUBB
      FLAMDA = 1.D0/YSUBC
      DELTA  = XSUBC/XSUBB - 1.D0
C
C     NOW FORM THE  C MATRIX   (3X6) PARTITIONED AS FOLLOWS HERE.
C         CSUBA = (3X2) STORED IN C( 1) THRU C( 6) BY ROWS
C         CSUBB = (3X2) STORED IN C( 7) THRU C(12) BY ROWS
C         CSUBC = (3X2) STORED IN C(13) THRU C(18) BY ROWS
C
      C(1)  =-REELMU
      C(2)  = 0.
      C(3)  = 0.
      C(4)  =  FLAMDA*DELTA
      C(5)  =  C(4)
      C(6)  =-REELMU
      C(7)  =  REELMU
      C(8)  = 0.
      C(9)  = 0.
      C(10) =-FLAMDA*REELMU*XSUBC
      C(11) = C(10)
      C(12) = REELMU
      C(13) = 0.
      C(14) = 0.
      C(15) = 0.
      C(16) = FLAMDA
      C(17) = FLAMDA
      C(18) = 0.
C
      IF (NTYPE .EQ. 1) GO TO 30
C
      THETA = ANGLE*DEGRA
      SINTH = SIN(THETA)
      COSTH = COS(THETA)
   30 IF (ABS(SINTH) .LT. 1.E-6) SINTH = 0.
C
C     BRANCH ON -HEAT- PROBLEM AT THIS POINT.
C
      IF (HEAT) GO TO 300
      ELTEMP = ECPT(21)
      MATID  = MATID1
      INFLAG = 2
      CALL MAT (ECPT(1))
C
C     FILL G-MATRIX WITH OUTPUT FROM MAT ROUTINE
C
      G(1) = G11
      G(2) = G12
      G(3) = G13
      G(4) = G12
      G(5) = G22
      G(6) = G23
      G(7) = G13
      G(8) = G23
      G(9) = G33
C
C     AT 50 G, E, AND C MATRICES ARE COMPLETE
C
C     AT THIS POINT THE FOLLOWING EQUATION CAN BE SOLVED FOR K-SUB-IJ
C
C                     T        T             T
C       K   = VOL . T  * E * C  * G * C  * E  * T
C        IJ          I        I        J         J
C
C     T-SUB-I WILL BE USED IN THE ABOVE ONLY IF THE PIVOT COORDINATE
C     SYSTEM ID IS NOT ZERO, OTHERWISE IT IS ASSUMED TO BE THE
C     IDENTITY MATRIX.
C
C     THE I SUBSCRIPT IMPLIES THE PIVOT POINT  1,2, OR 3 (ELEMENT SYST)
C     THE J SUBSCRIPT IMPLIES  1 THRU 3  FOR EACH CALL TO THIS ROUTINE.
C
C     DO COMPUTATIONS FOR EACH POINT IN ECPT LIST
C
      DO 60 I = 1,81
   60 KIJ(I) = 0.
      DO 200 NPVT = 1,3
      KA = 4*NPVT + 5
      NPOINT = 6*NPVT - 5
C
C                     T
C     COMPUTE   E * C   * G       AND STORE IN TEMPAR( 1 THRU 9 )
C                    I
C
      CALL GMMATD (E,  3,2,0, C(NPOINT), 3,2,1,  TEMPAR(10))
      CALL GMMATD (TEMPAR(10),3,3,0,  G, 3,3,0,  TEMPAR(1) )
C
C     NCOM WILL ALWAYS POINT TO THE COMMON 3 X 3 PRODUCT ABOVE
C     NPT1 WILL POINT TO FREE WORKING SPACE LENGTH 9
C
      NCOM = 1
      NPT1 = 10
C
C     MULTIPLY COMMON PRODUCT BY SCALER VOL
C
      DO 90 I = 1,9
   90 TEMPAR(I) = TEMPAR(I)*VOL
C
C     CHECK FOR PIVOT  CSID = 0, IF ZERO SKIP TRANSFORMATION TSUBI.
C
      IF (NECPT(KA) .EQ. 0)  GO TO 110
C
C     NOT-ZERO THUS GET TI
C
      CALL TRANSD (NECPT(KA),TI)
C
C     INTRODUCE TI INTO THE COMMON PRODUCT AND STORE AT
C     TEMPAR(10 THRU 18)
C
      CALL GMMATD (TI,3,3,1, TEMPAR(1),3,3,0, TEMPAR(10))
C
C     COMMON PRODUCT NOW STARTS AT TEMPAR(10) THUS CHANGE NCOM AND NPT1
C
      NCOM = 10
      NPT1 =  1
C
C
C  80 NOW HAVE COMMON PRODUCT STORED BEGINNING TEMPAR(NCOM),  (3X3).
C     NPT1 POINTS TO FREE WORKING SPACE LENGTH 9.
C
C     PROCEED NOW AND RUN OUT THE 3 6X6 MATRICES KIJ-SUB-1,2,3.
C
  110 NSAVE  = NPT1
      NPOINT = (NPVT-1)*27
C
C     INSERT G INTO TEMPAR
C
      DO 115 I = 1,9
115   TEMPAR(I+18) = G(I)
      DO 190 I = 1,3
      CALL GMMATD (C(6*I-5),3,2,0, E,3,2,1, TEMPAR(NSAVE))
C
C     NPT2 IS SET TO POINT TO THE BEGINNING OF THE PRODUCT  C * E * T
C                                                            J       J
C
      NPT2 = NSAVE
      NPT1 = 19
C
C     CHECK FOR ZERO CSID IN WHICH CASE TJ IS NOT NEEDED
C
      IF (NECPT(4*I +5) .EQ. 0) GO TO 120
C
C     COMMING HERE IMPLIES NEED FOR TJ
C     WILL STORE TJ IN TI
C
      CALL TRANSD (NECPT(4*I+5),TI)
      CALL GMMATD (TEMPAR(NPT2),3,3,0, TI,3,3,0, TEMPAR(19))
      NPT1 = NPT2
      NPT2 = 19
C
C  60 AT THIS POINT COMPLETE COMPUTATION FOR  K-SUB-I,J
C
  120 CALL GMMATD (TEMPAR(NCOM),3,3,0, TEMPAR(NPT2),3,3,0, TEMPAR(NPT1)
     X)
      NPT36 = NPT1 + 35
C
      DO 140 J = 1,9
      NPOINT = NPOINT + 1
      NPT2 = NPT1 + J - 1
  140 KIJ(NPOINT) = TEMPAR(NPT2)
  190 CONTINUE
  200 CONTINUE
C
      DICT5 = GSUBE
      RETURN
C
C     HEAT PROBLEM LOGIC PICKS UP HERE.  CALL HMAT FOR MATERIAL DATA.
C
  300 INFLAG = 2
      MATID  = NECPT(6)
      ELTEMP = ECPT(21)
      CALL HMAT (NECPT)
      G(1) = MATBUF(1)
      G(2) = MATBUF(2)
      G(3) = MATBUF(2)
      G(4) = MATBUF(3)
C
C     CONDENSE C MATRIX FOR HEAT PROBLEM (FORMED ABOVE)  C IS (2X3)
C
      C(2) = C( 4)
      C(3) = C( 7)
      C(4) = C(10)
      C(5) = C(13)
      C(6) = C(16)
C
C     DETERMINE THE PIVOT POINT.
C
      KQ   = 3
      KMAX = KQ*3
      DO 320 I = 1,KMAX
  320 KIJ(I) = 0.
      DO 400 NPVT = 1,3
C
C     PIVOT C MATRIX TIMES VOLUME (STORED INTO TT(1) AND TT(2).)
C
      TT(1) = VOL*C(2*NPVT-1)
      TT(2) = VOL*C(2*NPVT)
C
C     OUTPUT THE CONDUCTIVITY MATRICES
C
      NPOINT = (NPVT-1)*KQ
C
      DO 380 I = 1,3
      N2 = 2*I
      N1 = N2 - 1
      TEMPAR(1) = (G(1)*C(N1) + G(2)*C(N2))*TT(1) +
     1            (G(3)*C(N1) + G(4)*C(N2))*TT(2)
C
C     SUB-TRIANGLE (RETURN 3X3-S AS ABOVE IN STIFFNESS PORTION)
C
      KIJ(NPOINT+1) = TEMPAR(1)
      NPOINT = NPOINT + 1
  380 CONTINUE
  400 CONTINUE
      RETURN
C
C     ERROR  EXITS
C
 7770 CALL MESAGE (30,31,NECPT(1))
 7777 NOGO = .TRUE.
      RETURN
 7780 CALL MESAGE (30,32,NECPT(1))
      GO TO 7777
 7790 CALL MESAGE (30,26,NECPT(1))
      GO TO 7777
C
      END


text:
	SUBROUTINE DTICMK ( iwndw, itick, size, np, x, y, iret )
C************************************************************************
C* DTICMK								*
C*									*
C* This subroutine draws tic marks on the current graphics device.	*
C*									*
C* DTICMK ( IWNDW, ITICK, SIZE, NP, X, Y, IRET )			*
C*									*
C* Input parameters:							*
C*	IWNDW		INTEGER		Clipping window			*
C*	ITICK		INTEGER		Tick mark number		*
C*	SIZE		REAL		Tick size			*
C*	NP		INTEGER		Number of tick marks		*
C*	X (NP)		REAL		X coordinates in device units	*
C*	Y (NP)		REAL		Y coordinates in device units	*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. Vilardo/RDS	 6/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	 6/88	Cleaned up				*
C* M. desJardins/GSFC	 6/89	Added size				*
C* M. Linda/GSC		12/96	Changed X and Y to reals		*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVACT.CMN'
	INCLUDE		'DEVCHR.CMN'
C*
	REAL		x (*), y (*)
C------------------------------------------------------------------------
	iret = NORMAL
C
C*	Save line type and line width.
C
	IF  ( ( mltyp .ne. 1 ) .or. ( mlwid .ne. 1 ) ) THEN
	    jltyp = mltyp
	    jlwid = mlwid
	    CALL DSLINE ( 1, 0, 1, 0, i1, i2, i3, i4, ier )
	  ELSE
	    jltyp = 0
	    jlwid = 0
	END IF
C
C*	Set clipping window.
C
	CALL DSCLIP ( iwndw, ier )
C
C*	Compute tic mark size.
C
	ssize = size * 2 * bscalm
C
C*	Loop through tic marks.
C
	DO  i = 1, np
	    ix = NINT ( x (i) )
	    iy = NINT ( y (i) )
	    CALL ITICMK ( ix, iy, itick, ssize, ier )
	END DO
C
C*	Restore line type and line width.
C
	IF  ( ( jltyp .ne. 0 ) .or. ( jlwid .ne. 0 ) )
     +	    CALL DSLINE ( jltyp, 0, jlwid, 0, i1, i2, i3, i4, ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Always wanted to go skydiving?  Heres your chance! Skydance offers two ways to get started:

Tandem Assisted Freefall  An hour of instruction followed by a flight to between 9,000 and 14,000 feet where you and the instructor you are strapped to jump out and freefall for 3060 seconds before deploying the parachute.  Skydance also offers a tandem jump from 18,000 feet  thats 90 seconds of freefall!

Accelerated Freefall Program  A comprehensive training program starting with a full day of instruction. This is the training program that will set you on the road to making jumps by yourself and with others.

Every year they stage a jump demonstration before the July 4th Fireworks in Community Park. As they are jumping commentary is provided over the public address system. Should you have an opportunity to visit Skydance, be sure to check out their photo album, which includes some disturbing nude skydiving photos.  Spectators are always welcome; weekends are the best time to come out and watch.

For more information about the sport of skydiving, check out http://www.dropzone.com Dropzone.com.

In 2007, an episode of Mythbusters featured Skydance.

At least eight people have been killed in at least six accidents associated with SkyDance since it opened in 1987, out of an unknown total number of jumps; the last fatality was in 2013.  These fatal accidents are documented on the Tragic Events page.  Three of the fatalities were Skydance customers diving in tandem jumps with their instructors.  The fatality rate skydiving in the United States is about 1 per 60,000 jumps according to http://www.skydivingmagazine.com/faq.htm Skydiving Magazine; this is comparable to the fatality rate for driving 1,000 miles according to the http://hazmat.dot.gov/riskmgmt/riskcompare.htm Department of Transportation.

It is my understanding that Skydance runs one of the largest and best drop zones in the world.  Users/JasonAller

20060202 22:48:54 nbsp Just googled skydance skydiving death  at least three people have jumped out of a plane to their death at this place. Users/PaulThober

20060203 04:07:35 nbsp Well, thats always comforting. Users/JohnNapier

20060203 04:14:46 nbsp From what Ive been able to tell from a little internet research, both the top two comments are true. Which only goes to show that skydiving does come with some risk. That being said, it seems they also have twice yearly 30,000 foot drops. Thats just insane. Almost makes me want to spend the time and money to get my class D certificate so I could jump on one of those outings... Users/EricKlein

20060213 16:27:05 nbsp Ive jumped at dropzones throughout the US; Skydance is a great place to make your first (and subsequent) skydives.  They are one of only 5 or 6 places in the world flying the PAC750XL, a plane designed for skydiving.  Spectators are always welcome, but pets arent. Users/LisaBriggs

20060518 12:43:43 nbsp I did a tandem jump there last year. Their staff did a great job, I felt very safe and had a lot of fun. Users/GoKerby

20070731 17:36:37 nbsp I saw a tandem and some singles come down during the Plainfield volunteer firefighter pancake breakfast Users/StevenDaubert

20090419 22:52:54 nbsp I went for my first tandem jump here at 13000 feet, I felt very safe and the staff was very fun, yet professional. They are wonderful and I definitely recommend them! Users/SamanthaP

20090619 12:11:15 nbsp Can anyone add some price ranges for this? The website doesnt seem to mention it. Users/gurglemeow
  http://www.skydanceskydiving.net/tandem.php It was a few links in, in an odd location. You can click the link a the top of that page for AFF, which I assume means the program where youre not strapped to another person. users/JoePomidor

20100601 10:32:23 nbsp Here are the prices  http://www.skydance.net/skydance_skydiving_pricing_TANDEM.htm

This is a wonderful place, I am not old enough to skydive yet, but I think I will do a tandem on my 18th birthday and I already have a jump ticket! I have known Ray for a really long time and all of the people out there are super friendly.  Users/Jumpergem

20140313 10:49:44 nbsp I am very impressed by Skydance. Very professional and worth the first time tandem 13,000 foot jump (60 sec. freefall). My instructor (Greg) is exceptionally experienced and made the jump quite fun and safe. I would recommend jumping at Skydance ONLY if you really have wanted to skydive. Skydiving is not for everyone but very special for those who really have the desire. When possible, Oahu is next on my list. Users/GarySaylin

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


    MODULE random_negative_binomial_mod
! .. Implicit None Statement ..
      IMPLICIT NONE
! ..
    CONTAINS

!*********************************************************************

      FUNCTION random_negative_binomial(n,p)
! .. Use Statements ..
        USE random_standard_gamma_mod
        USE random_poisson_mod
! ..
! .. Implicit None Statement ..
        IMPLICIT NONE
! ..
! .. Function Return Value ..
        INTEGER :: random_negative_binomial
! ..
! .. Scalar Arguments ..
        REAL, INTENT (IN) :: p
        INTEGER, INTENT (IN) :: n
! ..
! .. Local Scalars ..
        REAL :: a, r, y
! ..
! .. Intrinsic Functions ..
        INTRINSIC REAL
! ..
! .. Executable Statements ..

!**********************************************************************
!     INTEGER FUNCTION IGNNBN( N, P )
!                GENerate Negative BiNomial random deviate
!                              Function
!     Generates a single random deviate from a negative binomial
!     distribution.
!                              Arguments
!     N  --> Required number of events.
!                              INTEGER N
!     JJV                      (N > 0)
!     P  --> The probability of an event during a Bernoulli trial.
!                              REAL P
!     JJV                      (0.0 < P < 1.0)
!                              Method
!     Algorithm from page 480 of
!     Devroye, Luc
!     Non-Uniform Random Variate Generation.  Springer-Verlag,
!     New York, 1986.
!**********************************************************************
!     JJV changed to call SGAMMA directly
!     REAL random_gamma
!      EXTERNAL random_gamma,random_poisson
!     Check Arguments
!     JJV changed argumnet checker to abort if N <= 0
        IF (n<=0) STOP 'N <= 0 in IGNNBN'
        IF (p<=0.0) STOP 'P <= 0.0 in IGNNBN'
        IF (p>=1.0) STOP 'P >= 1.0 in IGNNBN'

!     Generate Y, a random gamma (n,(1-p)/p) variable
!     JJV Note: the above parametrization is consistent with Devroye,
!     JJV       but gamma (p/(1-p),n) is the equivalent in our code
        r = REAL(n)
        a = p/(1.0-p)
!      y = random_gamma(a,r)
        y = random_standard_gamma(r)/a

!     Generate a random Poisson(y) variable
        random_negative_binomial = random_poisson(y)
        RETURN

      END FUNCTION random_negative_binomial

!*********************************************************************

    END MODULE random_negative_binomial_mod


text:
!     Calculation of equivalentwidth
      REAL*4 FUNCTION EW(mode,x,y,n,xb,xe)
!     mode=0 --> equivalentwidth   (ew > 0 : absorption  ew<0 : emission)
!     mode=1 --> integral
!     x,y    --> function
!     xb,xe  --> start and end value
      INTEGER*4   mode,n
      REAL*8 x(n),xb,xe, sum
      REAL*8 y(n)

      EW  = 0.D0
***      print*, x,y,n,xb,xe
      ib = 1
      ie = n
      DO i=2,n
       IF (x(i-1).LT.xb) THEN
          ib = i
       ELSE
          GOTO 10
       ENDIF
      ENDDO
 10   CONTINUE
      DO i=n-1,ib,-1
       IF (x(i+1).GT.xe) THEN
          ie = i
       ELSE
          GOTO 20
       ENDIF
      ENDDO
 20   CONTINUE

      sum = 0.D0
      IF (mode.EQ.0) THEN
         DO i=ib,ie-1
           sum = sum + DBLE(2. - ( y(i) + y(i+1) ) )  *
     *                       0.5D0 * (x(i+1) - x(i))
         END DO
         EW = sum
         RETURN
      ELSEIF (mode.EQ.1) THEN
         DO i=ib,ie-1
           sum = sum + DBLE( y(i) + y(i+1) ) *
     *                       0.5D0 * ( x(i+1) - x(i) )
         END DO
         EW = sum
         RETURN
      ELSE
         WRITE(*,*) ' EW: Unknown mode ! '
      ENDIF

      print*, 'ew ', ew

      RETURN
      END

max_stars_repo_name: andycasey/original-oracle

max_issues_repo_name: andycasey/original-oracle

max_forks_repo_name: andycasey/original-oracle


text:
!*******************************************************************************
!  File model_T.f
!  Contains module model_T
!  Defines derived-types: model

!*******************************************************************************
!  MODULE model_T
!    (Model Type Definition, for the V3FIT code)
! SECTION I.    Variable Declarations
! SECTION II.   Derived Type Declarations
! SECTION III.  Interface Blocks
! SECTION IV.   Construction Subroutine
! SECTION V.    Destruction Subroutines
! SECTION VI.   Assignment Subroutines
! SECTION VII.  Get Subroutines
! SECTION VIII. WRITE Subroutines

! SECTION XVI.  COMMENTS FOR DIFFERENT REVISIONS
!*******************************************************************************
      MODULE model_T

!*******************************************************************************
! SECTION I. VARIABLE DECLARATIONS
!*******************************************************************************

!-------------------------------------------------------------------------------
!  Type declarations - lengths of reals, integers, and complexes.
!  Frequently used mathematical constants, lots of extra precision.
!-------------------------------------------------------------------------------
      USE stel_kinds, only : rprec, iprec, cprec
      USE stel_constants, only : pi, twopi, one, zero

!-------------------------------------------------------------------------------
!  Use Statements for other structures, V3 Utilities
!-------------------------------------------------------------------------------
      USE bsc_T
      USE v3_utilities
      USE eq_T
!      USE density_T
      USE pprofile_T

!-------------------------------------------------------------------------------
!  Implicit None comes after USE statements, before other declarations
!-------------------------------------------------------------------------------
      IMPLICIT NONE

!-------------------------------------------------------------------------------
!  Make type declarations and constants Private, so there are no conflicts.
!-------------------------------------------------------------------------------
      PRIVATE rprec, iprec, cprec, pi, twopi, one, zero

!-------------------------------------------------------------------------------
!  Lengths of Character Variables
!-------------------------------------------------------------------------------
      INTEGER(iprec), PARAMETER, PRIVATE :: type_len=20

!*******************************************************************************
! SECTION II. DERIVED-TYPE DECLARATIONS
! 1)   Model
!         model
!     A model will contain an equilibrium state, and any other information
!     that is needed to compute model signals. For example, if the experiment
!     has interferometry and polarimetry diagnsotics, then there will need
!     to be a model density profile.
!*******************************************************************************
!-------------------------------------------------------------------------------
!  Declare type model
!         ---- Equilibrium State ----
!    eqstate         Type eq_state
!
!         ---- Parameterized Profiles (type pprofile)  -----
!    pp_ne          electron density
!    pp_sxrem       Soft x-ray emissivity
!    pp_te          Electron temperature
!
!         ----  Other Information ----
!  model_sxrem_type    Character variable to specify how the Soft x-ray
!                      emissivity will be calculated by the model.
!      'none'             Not part of the model.
!      'pp_sxrem'         Calculated from the parameterized profile pp_sxrem
!      'pp_ne_vmec_p'     Calculated from the pp_ne and the VMEC pressure
!
!  model_te_type    Character variable to specify how the electron temperature
!                   will be calculated by the model.
!      'none'             Not part of the model.
!      'pp_te'            Calculated from the parameterized profile pp_te
!      'pp_ne_vmec_p'     Calculated from the pp_ne and the VMEC pressure
!
!  ne_pp_unit        Scaling/conversion factor for the density. ne * ne_pp_unit
!                    is assumed to have units of m ** -3.
!
!  ne_min            Minimum electron density, m ** -3.
!
!  e_pressure_fraction    Electron pressure fraction of the total pressure
!                         Used when model_te_type = 'pp_ne_vmec_p'
!-------------------------------------------------------------------------------
      TYPE model
         TYPE (eq_state)          :: eqstate
         TYPE (pprofile)          :: pp_ne
         TYPE (pprofile)          :: pp_sxrem
         TYPE (pprofile)          :: pp_te
         CHARACTER(LEN=type_len)  :: model_sxrem_type
         CHARACTER(LEN=type_len)  :: model_te_type
         REAL(rprec)              :: ne_pp_unit
         REAL(rprec)              :: ne_min
         REAL(rprec)              :: e_pressure_fraction
      END TYPE model

!*******************************************************************************
! SECTION III. INTERFACE BLOCKS
!*******************************************************************************
!-------------------------------------------------------------------------------
!  Assignment for structures
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!  Generic construct
!-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
!  Generic destroy
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!  Generic write
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!  Interface block for testing goes here.
!-------------------------------------------------------------------------------

      CONTAINS
!*******************************************************************************
! SECTION IV. CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!  Construct a model
!-------------------------------------------------------------------------------
      SUBROUTINE model_construct(this,eqstate, pp_ne, pp_sxrem,                &
     &   pp_te, model_sxrem_type, model_te_type,ne_pp_unit,ne_min,             &
     &   e_pressure_fraction)

      IMPLICIT NONE

!  Declare Arguments
      TYPE (model), INTENT(inout)              :: this
      TYPE (eq_state), INTENT(in), OPTIONAL    :: eqstate
      TYPE (pprofile), INTENT(in), OPTIONAL    :: pp_ne
      TYPE (pprofile), INTENT(in), OPTIONAL    :: pp_sxrem
      TYPE (pprofile), INTENT(in), OPTIONAL    :: pp_te
      CHARACTER(LEN=20), INTENT(in), OPTIONAL  :: model_sxrem_type
      CHARACTER(LEN=20), INTENT(in), OPTIONAL  :: model_te_type
      REAL(rprec), INTENT(in), OPTIONAL        :: ne_pp_unit
      REAL(rprec), INTENT(in), OPTIONAL        :: ne_min
      REAL(rprec), INTENT(in), OPTIONAL        :: e_pressure_fraction

!  Declare local variables
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_construct: '

!  Start of executable code

!  Copy the OPTIONAL arguments
      IF (PRESENT(eqstate)) THEN
         this % eqstate = eqstate
      ENDIF

      IF (PRESENT(pp_ne)) THEN
         this % pp_ne = pp_ne
      ENDIF

      IF (PRESENT(pp_sxrem)) THEN
         this % pp_sxrem = pp_sxrem
      ENDIF

      IF (PRESENT(pp_te)) THEN
         this % pp_te = pp_te
      ENDIF

      IF (PRESENT(model_sxrem_type)) THEN
         this % model_sxrem_type = model_sxrem_type
      ENDIF

      IF (PRESENT(model_te_type)) THEN
         this % model_te_type = model_te_type
      ENDIF

      IF (PRESENT(ne_pp_unit)) THEN
         this % ne_pp_unit = ne_pp_unit
      ENDIF

      IF (PRESENT(ne_min)) THEN
         this % ne_min = ne_min
      ENDIF

      IF (PRESENT(e_pressure_fraction)) THEN
         this % e_pressure_fraction = e_pressure_fraction
      ENDIF

      END SUBROUTINE model_construct

!*******************************************************************************
! SECTION V. DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!  Destroy a model
!-------------------------------------------------------------------------------
      SUBROUTINE model_destroy(this)
      IMPLICIT NONE

!  Declare Arguments
      TYPE (model), INTENT(inout) :: this

!  Declare local variables
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_destroy: '

!  Start of executable code

!  As of 12-14-2004, there isn't an eq_destroy for a state.
!      CALL eq_destroy(this % eqstate)

!.........destroy the parameterized profiles
      CALL pprofile_destroy( this % pp_ne)
      CALL pprofile_destroy( this % pp_sxrem)
      CALL pprofile_destroy( this % pp_te)

      this % model_sxrem_type = 'none'
      this % model_te_type = 'none'

      this % ne_pp_unit = zero
      this % ne_min = zero
      this % e_pressure_fraction = zero

      END SUBROUTINE model_destroy

!*******************************************************************************
! SECTION VI. ASSIGNMENT SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!  Assignment for model
!-------------------------------------------------------------------------------
      SUBROUTINE model_assign(left,right)

      IMPLICIT NONE

!  Declare Arguments
      TYPE (model), INTENT (inout) :: left
      TYPE (model), INTENT (in) :: right

!  Declare local variables
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_assign: '

!  Start of executable code

      left % eqstate = right % eqstate
      left % pp_ne = right % pp_ne
      left % pp_sxrem = right % pp_sxrem
      left % pp_te = right % pp_te
      left % model_sxrem_type = right % model_sxrem_type
      left % model_te_type = right % model_te_type
      left % ne_pp_unit = right % ne_pp_unit
      left % ne_min = right % ne_min
      left % e_pressure_fraction = right % e_pressure_fraction

      END SUBROUTINE model_assign

!*******************************************************************************
! SECTION VII. GET SUBROUTINES
!*******************************************************************************
! Information needed about anything in the model should be accessed through
! these subroutines. (That's the goal. JDH 2011-07-14)
!
!  Need to get the magnetic diagnostics working through here.
!
!  For profiles, the _xcyl functions call the _s functions.
!-------------------------------------------------------------------------------
!  Plasma surface information (Pass through to eq)
!-------------------------------------------------------------------------------
      SUBROUTINE model_get_seq1_rz(this,phi,r_at_seq1,z_at_seq1)
!  Subroutine to return the R and Z positions of points on the
!  plasma surface (s = 1 surface), on a constant_phi plane.
      IMPLICIT NONE

!  Declare Arguments
      TYPE (model), INTENT (inout)             :: this
      REAL(rprec), INTENT (in)                 :: phi
      REAL(rprec), DIMENSION(:), INTENT(inout) :: r_at_seq1, z_at_seq1

!  Start of executable code. Pass through to the eq.
      CALL eq_get_seq1_rz(this % eqstate,phi,r_at_seq1,z_at_seq1)

      END SUBROUTINE model_get_seq1_rz
!-------------------------------------------------------------------------------
!  Electron Density ne
!
!     Given Cylindrical Coordinates XCYL
!-------------------------------------------------------------------------------
      FUNCTION model_get_ne_xcyl(this,xcyl)
!  FUNCTION to return the electron density, given a model and
!  Cylindrical coordinates xcyl(3) [ = R, Phi, Z]
      IMPLICIT NONE

!  Declare Arguments
      REAL(rprec) :: model_get_ne_xcyl
      TYPE (model), INTENT (inout)             :: this
      REAL (rprec), INTENT (in)                :: xcyl(3)

!  Declare local variables
      REAL(rprec) :: c_flux(3)               ! VMEC coordinates (s,u,v)
      REAL(rprec) :: s                       ! VMEC radial coordinate
      REAL(rprec) :: ne                      ! electron density
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_get_ne_s: '

!  Start of executable code.

!  Get the radial coordinate s
      c_flux = (/zero, zero, zero /)
      CALL eq_get_flux(this % eqstate,xcyl,c_flux)
      s = c_flux(1)

      ne = model_get_ne_s(this,s)
      model_get_ne_xcyl = ne
      RETURN

      END FUNCTION model_get_ne_xcyl
!-------------------------------------------------------------------------------
!  Electron Density ne
!
!     Given VMEC radial coordinate S
!-------------------------------------------------------------------------------
      FUNCTION model_get_ne_s(this,s_arg)
!  FUNCTION to return the electron density in m ** -3, given a model
!  and a value of the radial coordinate s.
!  Note the conversion factor ne_pp_unit.
      IMPLICIT NONE

!  Declare Arguments
      REAL(rprec) :: ne                      ! electron density
      TYPE (model), INTENT (inout)             :: this
      REAL (rprec), INTENT (in)                :: s_arg

!  Declare local variables
      REAL(rprec) :: model_get_ne_s
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_get_ne_s: '

!  Start of executable code.

      model_get_ne_s = MAX(this % ne_min,this % ne_pp_unit *                   &
     &   pprofile_evaluate(this % pp_ne,s_arg))
      RETURN

      END FUNCTION model_get_ne_s

!-------------------------------------------------------------------------------
!  Soft x-ray emissivity
!
!     Given Cylindrical Coordinates XCYL
!-------------------------------------------------------------------------------
      FUNCTION model_get_sxrem_xcyl(this,xcyl)
!  FUNCTION to return the soft x-ray emissivity, given a model and
!  Cylindrical coordinates xcyl(3) [ = R, Phi, Z]
      IMPLICIT NONE

!  Declare Arguments
      REAL(rprec) :: model_get_sxrem_xcyl
      TYPE (model), INTENT (inout)             :: this
      REAL (rprec), INTENT (in)                :: xcyl(3)

!  Declare local variables
      REAL(rprec) :: c_flux(3)               ! VMEC coordinates (s,u,v)
      REAL(rprec) :: s                       ! VMEC radial coordinate
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_get_sxrem_xcyl: '

!  Start of executable code.

!  Get the radial coordinate s
      c_flux = (/zero, zero, zero /)
      CALL eq_get_flux(this % eqstate,xcyl,c_flux)
      s = c_flux(1)

      model_get_sxrem_xcyl = model_get_sxrem_s(this,s)

      RETURN
      END FUNCTION model_get_sxrem_xcyl

!-------------------------------------------------------------------------------
!  Soft x-ray emissivity
!
!     Given VMEC radial coordinate S
!-------------------------------------------------------------------------------
      FUNCTION model_get_sxrem_s(this,s_arg)
!  FUNCTION to return the soft x-ray emissivity, given a model and
!  and a value of the radial coordinate s.
      IMPLICIT NONE

!  Declare Arguments
      REAL(rprec) :: model_get_sxrem_s
      TYPE (model), INTENT (inout)             :: this
      REAL (rprec), INTENT (in)                :: s_arg

!  Declare local variables
      REAL(rprec) :: ne
      REAL(rprec) :: pressure
      REAL(rprec) :: emissivity
      REAL(rprec) :: ne_norm = 1.E18         ! Normalizing constant, per Hartwell routine
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_get_ne_s: '

!  Start of executable code.
      SELECT CASE(TRIM(this % model_sxrem_type))

      CASE('pp_sxrem')
         emissivity = pprofile_evaluate(this % pp_sxrem,s_arg)

      CASE('pp_ne_vmec_p')
         IF (s_arg .le. one) THEN
            ne = model_get_ne_s(this,s_arg)
            pressure = MAX(0.,eq_get_pressure_s(this % eqstate,s_arg))
            emissivity = (ne / ne_norm) ** (3/2) * SQRT(pressure)   ! per comments
         ELSE
            emissivity = zero
         ENDIF

      CASE ('none')
         emissivity = -2.E49_rprec

      CASE DEFAULT
         WRITE(*,*) ' IN model_get_sxrem_s. model_sxrem_type is ',             &
     &      TRIM(this % model_sxrem_type)
         WRITE(*,*) '  Needs to be pp_sxrem or pp_ne_vmec_p'
         STOP

      END SELECT

      model_get_sxrem_s = emissivity

      RETURN
      END FUNCTION model_get_sxrem_s

!-------------------------------------------------------------------------------
!  Line integral for Interferometry - ne
!
!     Given Cylindrical Coordinates XCYL
!-------------------------------------------------------------------------------
      FUNCTION model_get_ip_i_xcyl(this,xcyl)
!  FUNCTION to return the interferometry integrand, given a model and
!  Cylindrical coordinates xcyl(3) [ = R, Phi, Z]
      IMPLICIT NONE

!  Declare Arguments
      REAL(rprec) :: model_get_ip_i_xcyl
      TYPE (model), INTENT (inout)             :: this
      REAL (rprec), INTENT (in)                :: xcyl(3)

!  Declare local variables
      REAL(rprec) :: c_flux(3)               ! VMEC coordinates (s,u,v)
      REAL(rprec) :: s                       ! VMEC radial coordinate
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_get_ip_i_xcyl: '

!  Start of executable code.

!  Get the radial coordinate s
      c_flux = (/zero, zero, zero /)
      CALL eq_get_flux(this % eqstate,xcyl,c_flux)
      s = c_flux(1)

      model_get_ip_i_xcyl = model_get_ip_i_s(this,s)

      RETURN
      END FUNCTION model_get_ip_i_xcyl

!-------------------------------------------------------------------------------
!  Line integral for Interferometry - ne
!
!     Given VMEC radial coordinate S
!-------------------------------------------------------------------------------
      FUNCTION model_get_ip_i_s(this,s_arg)
!  FUNCTION to return the soft x-ray emissivity, given a model and
!  and a value of the radial coordinate s.
      IMPLICIT NONE

!  Declare Arguments
      REAL(rprec) :: model_get_ip_i_s
      TYPE (model), INTENT (inout)             :: this
      REAL (rprec), INTENT (in)                :: s_arg

!  Declare local variables
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_get_ip_i_s: '

!  Start of executable code.

      model_get_ip_i_s = model_get_ne_s(this,s_arg)

      RETURN
      END FUNCTION model_get_ip_i_s
!-------------------------------------------------------------------------------
!  Electron Temperature
!
!     Given Cylindrical Coordinates XCYL
!-------------------------------------------------------------------------------
      FUNCTION model_get_te_xcyl(this,xcyl)
!  FUNCTION to return the electron temperature in eV, given a model and
!  Cylindrical coordinates xcyl(3) [ = R, Phi, Z]
      IMPLICIT NONE

!  Declare Arguments
      REAL(rprec) :: model_get_te_xcyl
      TYPE (model), INTENT (inout)             :: this
      REAL (rprec), INTENT (in)                :: xcyl(3)

!  Declare local variables
      REAL(rprec) :: c_flux(3)               ! VMEC coordinates (s,u,v)
      REAL(rprec) :: s                       ! VMEC radial coordinate
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_get_te_xcyl: '

!  Start of executable code.

!  Get the radial coordinate s
      c_flux = (/zero, zero, zero /)
      CALL eq_get_flux(this % eqstate,xcyl,c_flux)
      s = c_flux(1)

      model_get_te_xcyl = model_get_te_s(this,s)
      RETURN

      END FUNCTION model_get_te_xcyl
!-------------------------------------------------------------------------------
!  Electron Temperature
!
!     Given VMEC radial coordinate S
!-------------------------------------------------------------------------------
      FUNCTION model_get_te_s(this,s_arg)
!  FUNCTION to return the electron temperature in eV, given a model and
!  and a value of the radial coordinate s.
      IMPLICIT NONE

!  Declare Arguments
      REAL(rprec) :: model_get_te_s
      TYPE (model), INTENT (inout)             :: this
      REAL (rprec), INTENT (in)                :: s_arg

!  Declare Parameters
      REAL(rprec), PARAMETER :: eV_per_Joule = one / 1.602e-19

!  Declare local variables
      REAL(rprec) :: ne                      ! electron density (m ** -3)
      REAL(rprec) :: te                      ! electron temperature (eV)
      REAL(rprec) :: pressure                ! pressure
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_get_te_s: '

!  Start of executable code.

!  Get the radial coordinate s
      SELECT CASE(TRIM(this % model_te_type))

      CASE('pp_te')
         te = pprofile_evaluate(this % pp_te,s_arg)

      CASE('pp_ne_vmec_p')
         ne = model_get_ne_s(this,s_arg)
         pressure = MAX(0.,eq_get_pressure_s(this % eqstate,s_arg))
         te = eV_per_Joule * this % e_pressure_fraction * pressure / ne

      CASE ('none')
         te = -2.E49_rprec

      CASE DEFAULT
         WRITE(*,*) ' IN model_get_te_s. model_te_type is ',                   &
     &      TRIM(this % model_te_type)
         WRITE(*,*) '  Needs to be pp_te or pp_ne_vmec_p'
         STOP

      END SELECT

      model_get_te_s = te
      RETURN

      END FUNCTION model_get_te_s
!*******************************************************************************
! SECTION VIII.  WRITE SUBROUTINES
!*******************************************************************************
      SUBROUTINE model_write_profiles(a_model,identifier,unit)
      IMPLICIT NONE

!  Declare Arguments
      TYPE (model), INTENT (inout) :: a_model
      CHARACTER (len=*), INTENT(in), OPTIONAL :: identifier
      INTEGER, INTENT(in), OPTIONAL :: unit
!  identifier   character variable, also written out
!  unit         I/O unit number to write to

!  Declare local variables and constants
      CHARACTER(len=*), PARAMETER :: sub_name =                                &
     &  'model_write_profiles: '
      INTEGER :: iv_default = 1
      INTEGER :: iv
      INTEGER :: iou_default = 6
      INTEGER :: iou
      INTEGER :: i, nrow
      CHARACTER (len=60) :: id
      REAL(rprec) :: s, p, ne, te, sxrem

!  start of executable code
!  Check for arguments present
      IF (PRESENT(identifier)) THEN
         id = identifier
      ELSE
         id = ' '
      END IF

      IF (PRESENT(unit)) THEN
         iou = unit
      ELSE
         iou = iou_default
      END IF

!  Write Header
      WRITE(iou,900) id

      CALL pprofile_write(a_model % pp_ne,sub_name // 'ne',iou)
      CALL pprofile_write(a_model % pp_sxrem,sub_name // 'sxrem',iou)
      CALL pprofile_write(a_model % pp_te,sub_name // 'te',iou)

      WRITE(iou,1000)
      WRITE(iou,1010)
      WRITE(iou,1020) TRIM(a_model % model_te_type),                           &
     &   TRIM(a_model % model_sxrem_type)

      nrow = a_model % eqstate % fixp % ns_array(                              &
     &   a_model % eqstate % varp % ns_index) + 2
      DO i = 1,nrow + 2
         s = (i - one) / nrow
         p = eq_get_pressure_s(a_model % eqstate,s)
         ne = model_get_ne_s(a_model,s)
         te = model_get_te_s(a_model,s)
         sxrem = model_get_sxrem_s(a_model,s)
         WRITE(iou,1100) s, p, ne, te, sxrem
      END DO

      WRITE(iou,1200) id

      RETURN

900   FORMAT(/' *** Model Write Profiles. id = ',a)
1000  FORMAT(/t5,'S',t16,'PRESSURE',t30,'e DENSITY',t46,                       &
     &   'e TEMPERATURE',t63,'SXR EMISSIVITY')
1010  FORMAT(t18,'(Pascal)',t32,'(m-3)',t48,'(eV)',t65,'( - )')
1020  FORMAT(t18,'VMEC',t32,'pp_ne',t48,a12,t65,a12/80('-'))
1100  FORMAT(2x,es9.2,t16,es10.3,t30,es10.3,t46,es10.3,t63,es10.3)
1200  FORMAT(/' *** END Model Write Profiles. id = ',a)

      END SUBROUTINE model_write_profiles
!*******************************************************************************
! SECTION XVI.  COMMENTS FOR DIFFERENT REVISIONS
!*******************************************************************************
!
!  JDH 2011-07-14
!    Added model_get_seq1_rz
!
!  JDH 12-14-04.
!     First version, modified signal_T
!  JMS 7-24-07.
!     added density to the model
!
!  JDH 2011-07-14
!    Add model_get_seq1_rz
!
!  JDH 2011-07-28
!    Add model_get_sxrem_s, model_get_ne_s, model_get_pressure_s
!
!  JDH 2011-08-07
!    Correct typo in emissivity, * SQRT(pressure) instead of + SQRT(pressure)
!
!  JDH 2011-09-06
!    model_get_szr_em_s: add MAX(0.,-) for density and pressure
!
!  2011-09-08 JDH
!    Replaced model_get_sxrem_s with _xcyl. Eliminated
!    model_get_pressure_s.
!
!  2011-10-19 JDH
!    Convert from e_density to pprofile pp_ne.
!
!  2011-10-23 JDH
!    Added model_sxrem_type and pp_sxrem
!
!  2011-10-23 JDH
!    Added model_te_type and pp_te, FUNCTION model_get_te_xcyl
!
!  2011-10-25 JDH
!    Added ne_pp_unit, so that ne returned by model_get is in m ** -3.
!    Added e_pressure_fraction. Rearranged the get functions.
!    Added model_write_profiles.

      END MODULE model_T