text:
      SUBROUTINE DCONE
C
C     DIFFERENTIAL STIFFNESS FOR THE CONICAL SHELL. FMMS-68
C
C     CALLS FROM DCONE ARE MADE TO
C           MESAGE
C           MAT
C           INVERD
C           GMMATD
C           DS1B
C
C     ECPT( 1) = ELEMENT ID                                INTEGER
C     ECPT( 2) = SIL PT A                                  INTEGER
C     ECPT( 3) = SIL PT B                                  INTEGER
C     ECPT( 4) = MATID 1                                   INTEGER
C     ECPT( 5) = TM  (MEMBRANE THICK)                      REAL
C     ECPT( 6) = MATID 2                                   INTEGER
C     ECPT( 7) = I   (MOM.OF INERTIA)                      REAL
C     ECPT( 8) = MATID 3                                   INTEGER
C     ECPT( 9) = TS  (SHEAR THICKNESS)                     REAL
C     ECPT(10) = NON-STRUCTURAL-MASS                       REAL
C     ECPT(11) = Z1                                        REAL
C     ECPT(12) = Z2                                        REAL
C     ECPT(13) = PHI  1                                    REAL
C     ECPT(14) = PHI  2                                    REAL
C     ECPT(15) = PHI  3                                    REAL
C     ECPT(16) = PHI  4                                    REAL
C     ECPT(17) = PHI  5                                    REAL
C     ECPT(18) = PHI  6                                    REAL
C     ECPT(19) = PHI  7                                    REAL
C     ECPT(20) = PHI  8                                    REAL
C     ECPT(21) = PHI  9                                    REAL
C     ECPT(22) = PHI 10                                    REAL
C     ECPT(23) = PHI 11                                    REAL
C     ECPT(24) = PHI 12                                    REAL
C     ECPT(25) = PHI 13                                    REAL
C     ECPT(26) = PHI 14                                    REAL
C     ECPT(27) = COORD. SYS. ID PT.1                       INTEGER
C     ECPT(28) = RADIUS PT. 1                              REAL
C     ECPT(29) = DISTANCE TO PT.1                          REAL
C     ECPT(30) = NULL                                      REAL
C     ECPT(31) = COORD. SYS. ID PT.2                       INTEGER
C     ECPT(32) = RADIUS PT 2                               REAL
C     ECPT(33) = DISTANCE TO PT. 2                         REAL
C     ECPT(34) = NULL                                      REAL
C     ECPT(35) = ELEMENT TEMPERATURE                       REAL
C     ECPT(36) = ELEMENT DEFORMATION                       REAL
C     ECPT(37) = ELEMENT LOADING TEMPERATURE - GRID PT A   REAL
C     ECPT(38) = ELEMENT LOADING TEMPERATURE - GRID PT B   REAL
C     ECPT(39) = DISPLACEMENT COMPONENTS AT GRID POINT A   REAL
C     ECPT(40) =                  ...                      REAL
C     ECPT(41) =                  ...                      REAL
C     ECPT(42) =                  ...                      REAL
C     ECPT(43) =                  ...                      REAL
C     ECPT(44) =                  ...                      REAL
C     ECPT(45) = DISPLACEMENT COMPONENTS AT GRID POINT B   REAL
C     ECPT(46) =                  ...                      REAL
C     ECPT(47) =                  ...                      REAL
C     ECPT(48) =                  ...                      REAL
C     ECPT(49) =                  ...                      REAL
C     ECPT(50) =                  ...                      REAL
C
      INTEGER          NECPT(100) ,NERROR(2) ,NA(10)
      DOUBLE PRECISION INT(10,4)  ,HUQ       ,U(10)     ,KQD       ,
     1                 A(5,3)     ,HYQ(10)   ,Q(8)      ,KIJ       ,
     2                 B(7,3)     ,EHT       ,FAC(10)   ,C(3,3)    ,
     3                 RA         ,E11       ,TEMP      ,ONE       ,
     4                 RB         ,E12       ,TEMP1     ,OPI       ,
     5                 ZA         ,E22       ,TEMP2     ,N2D33     ,
     6                 ZB         ,E33       ,TEMP3     ,SP2D22    ,
     7                 RASQ       ,D11       ,TEMP4     ,SP2D4     ,
     8                 RBSQ       ,D12       ,TEMP5     ,OQ        ,
     9                 PI         ,D22       ,TEMP6     ,TDIF      ,
     O                 PIOVB      ,D33       ,TEMP7     ,DEPS      ,
     1                 N          ,TS        ,NSPOPI    ,DEPP      ,
     2                 N2         ,TM        ,TWOD33    ,EPS       ,
     3                 SL         ,NSP       ,NOV4      ,EPP       ,
     4                 L2         ,NCP       ,NSPOV4    ,TE11      ,
     5                 SP         ,SP2       ,N2OV4     ,TE12      ,
     6                 CP         ,CP2       ,SD22PI    ,TE22      ,
     7                 SUM        ,SIGN      ,GSHEAR    ,TEMP48(48),
     8                 A0         ,A1        ,A2        ,A3
      DOUBLE PRECISION B0         ,B1        ,B2        ,B3        ,
     1                 C0         ,C1        ,D0        ,D1        ,
     2                 DETERM     ,CONSTD
      COMMON /DS1AAA/  NPVT       ,DUMCL(34) ,NOGO
      COMMON /DS1AET/  ECPT(100)
      COMMON /DS1ADP/  HUQ(100)   ,KQD(64)   ,KIJ(36)   ,EHT(96)   ,
     1                 E11        ,E12       ,E22       ,E33
      COMMON /MATIN /  MATID      ,INFLAG    ,ELTEMP    ,STRESS    ,
     1                 SINTH      ,COSTH
      COMMON /MATOUT/  G11        ,G12       ,G13       ,G22       ,
     1                 G23        ,G33       ,RHOY      ,ALPHA1    ,
     2                 ALPHA2     ,ALPH12    ,DUM(10)
      COMMON /CONDAD/  CONSTD(5)
      EQUIVALENCE      (G,G12)               ,(ECPT(1),NECPT(1))   ,
     1                 (ECPT(4),MATID1)      ,(ECPT(6), MATID2)    ,
     2                 (ECPT(8),MATID3)      ,(CONSTD(1),   PI)
      DATA    NA    /  6*1, 2*2, 2*4/
      DATA    FAC   /  1.0D0, 1.0D0, 2.0D0, 6.0D0, 24.0D0, 120.0D0,
     1                 720.0D0, 5040.0D0, 40320.0D0, 362880.0D0   /
      DATA    ONE   /  1.0D0  /
C
C
C     CALCULATE SHELL ORIENTATION CONSTANTS
C
      SINTH = 0.0
      COSTH = 1.0
      NINT  = NECPT(1)/1000
      N     = NECPT(1) - NINT*1000 - 1
      RA    = ECPT(28)
      ZA    = ECPT(29)
      RB    = ECPT(32)
      ZB    = ECPT(33)
      TEMP1 = RB - RA
      TEMP2 = ZB - ZA
      L2    = TEMP1**2 + TEMP2**2
      SL    = DSQRT(L2)
      IF (SL) 30,20,30
   20 NERROR(1) = NECPT(1)/1000
      NERROR(2) = N + .3D0
      CALL MESAGE (30,39,NERROR(1))
C
C     SET FLAG FOR FATAL ERROR WHILE ALLOWING ERROR MESSAGES TO
C     ACCUMULATE
C
      NOGO = 1
      RETURN
C
   30 SP = TEMP1/SL
      CP = TEMP2/SL
C
C     COMPUTE INTEGRALS I     FOR M = 0,9
C                        MN       N = 0,3
C
C     FOR EVALUATION OF INTEGRALS  A = RA,  B = SP
C
      IF (SP) 60,40,60
C
C     COMPUTE INTEGRAL FOR B = 0
C
C                            1-N
C                      PI  RA     M+1
C               I   = --------- SL    (FOR ALL M,N .GE. 0)
C                M,N    M + 1
C
C
C     M = I - 1    WHERE I IS THE DO LOOP INDEX
C     N = J - 1    WHERE J IS THE DO LOOP INDEX
C     MPLUS1 THUS EQUALS I
C
   40 DO 50 I = 1,10
      NBEGIN  = NA(I)
      DO 50 J = NBEGIN,4
   50 INT(I,J) = (PI*SL**I)/(DBLE(FLOAT(I))*RA**(J-2))
C
      GO TO 100
C
C
C     COMPUTE INTEGRALS FOR (B .NE. 0)
C
C     FIRST M = 0 CASE
C
C                             2-N     2-N
C                      PI ( RB    - RA   )
C               I     =--------------------   (N NOT EQUAL TO 2)
C                0,N       (2-N)  B
C
C
C     FOR N=2   I     = PI * (LOG RB  -  LOG RA) / B
C                0,2             E          E
C
C
   60 RASQ  = RA*RA
      RBSQ  = RB*RB
      PIOVB = PI/SP
C
      INT(1,1) = 0.5D0*PIOVB*(RBSQ - RASQ)
      INT(1,2) = PIOVB*(RB - RA)
      INT(1,3) = PIOVB*DLOG(RB/RA)
      INT(1,4) =-PIOVB*(ONE/RB - ONE/RA)
C
C
C     M = I        WHERE I IS THE DO LOOP INDEX
C     N = J - 1    WHERE J IS THE DO LOOP INDEX
C
C     WE ARE GETTING INTEGRAL(M,N)
C     M = POWER OF S
C     N = POWER OF R
C
C
C     EVALUATING AT R = RB  THEN AT R = RA
C
C                                  K   NPOW
C                  M FAC.      M   (-A) (R)
C     I  = (PI)(-----------)( SUM ------------------------) + (TERM-X)
C      MN               (M+1)   K=0  (M-K)FAC.(K)FAC.(NPOW)
C                    B        (K.NE.M-N+2)                  (K.EQ.M-N+2)
C
C
C     WHERE NPOW = M - N - K + 2
C
C
C                    M-N+2
C                (-A)     LOG(R)
C       TERM-X = --------------------
C               (M-N+2)FAC.(N-2)FAC.
C
C
C     NOTE IN DATA STATEMENT THAT 0 FACTORIAL = FAC(1)
C                                 1 FACTORIAL = FAC(2)
C                                 2 FACTORIAL = FAC(3)    ETC.
C
      DO 90 I = 1,9
      MPLUS1  = I + 1
      NBEGIN  = NA(MPLUS1)
      DO 90 J = NBEGIN,4
      SUM  = 0.0D0
      SIGN =-1.0D0
      NPOW = I - J + 3
      DO 80 KK = 1,MPLUS1
      SIGN = -SIGN
      K    = KK - 1
      IF (K .EQ. NPOW) GO TO 70
      KPOW = NPOW   - K
      IFAC = MPLUS1 - K
      TEMP = KPOW
      SUM  = SUM + SIGN*RA**K*(RB**KPOW - RA**KPOW)/
     1       (FAC(IFAC)*FAC(KK)*TEMP)
      GO TO 80
   70 SUM = SUM + SIGN*RA**NPOW*DLOG(RB/RA)/(FAC(NPOW+1)*FAC(J-2))
   80 CONTINUE
C
      INT(MPLUS1,J) = SUM*PI*FAC(MPLUS1)/SP**MPLUS1
   90 CONTINUE
  100 CONTINUE
C
C     CRANK OUT HUQ MATRIX FOR ZERO HARMONIC
C     FOR EXPLICIT FORMULATION OF HUQ, SEE MS-28, PP.15,16 AND PP.24,25.
C
      DO 105 I = 1,100
  105 HUQ(  I) = 0.0D0
      HUQ(  1) = ONE
      HUQ( 13) = ONE
      HUQ( 25) = ONE
      HUQ( 36) = ONE
      HUQ( 41) = CP/RA
      HUQ( 49) = ONE
      HUQ( 51) = ONE
      HUQ( 52) = SL
      HUQ( 63) = ONE
      HUQ( 64) = SL
      HUQ( 75) = ONE
      HUQ( 76) = SL
      HUQ( 77) = L2
      HUQ( 78) = HUQ(77)*SL
      HUQ( 86) = ONE
      HUQ( 87) = 2.0D0*SL
      HUQ( 88) = 3.0D0*HUQ(77)
      HUQ( 91) = CP/RB
      HUQ( 92) = HUQ(91)*SL
      HUQ( 99) = ONE
      HUQ(100) = SL
C
C     IF TRANSVERSE SHEAR IS ZERO
C
C     OR INERTIA           = 0.0
C     OR SHEAR MODULUS(G)  = 0.0
C     OR MATID2            = 0
C     OR MATID3            = 0
C
C     THEN (HYQ)  = (0).  THEREFORE, USE HUQ MATRIX AS IS
C
      IF (MATID2.EQ.0 .OR. MATID3.EQ.0) GO TO 130
      IF (ECPT(9).EQ.0.0 .OR. ECPT(7).EQ.0.0) GO TO 130
      INFLAG = 1
      MATID  = MATID3
      ELTEMP = ECPT(35)
      CALL MAT (ECPT(1))
      GSHEAR = G
      IF (G .EQ. 0.0) GO TO 130
      INFLAG = 2
      MATID  = MATID2
      ELTEMP = ECPT(35)
      CALL MAT (ECPT(1))
C
C     FORM
C     (D) = I*(G)
C
      D11 = ECPT(7)*G11
      D12 = ECPT(7)*G12
      D22 = ECPT(7)*G22
      D33 = ECPT(7)*G33
C
      TS  = ECPT(9)
C
      DO 110 I = 1,10
  110 HYQ(I) = 0.0D0
      CP2 = CP*CP
      SP2 = SP*SP
      N2  = N*N
      OPI = ONE/PI
      SD22PI = SP2*D22*OPI
      OQ  = SL*TS*GSHEAR*(RA+RB)*0.5D0 + SD22PI*INT(1,3)
      OQ  = ONE/OQ
C
      HYQ(6) = OQ*INT(1,3)*SD22PI
      HYQ(7) = OQ*2.0D0*(D11 *(RA-RB) + INT(2,3)*SD22PI)
      HYQ(8) = OQ*(-D11*6.0D0*SL*RB   + 3.0D0*INT(3,3)*SD22PI)
C
      DO 120 I = 6,8
      HUQ(I+30) = HUQ(I+30) - HYQ(I)
  120 HUQ(I+80) = HUQ(I+80) - HYQ(I)
C
  130 CONTINUE
C
C     NO NEED TO COMPUTE DETERMINANT SINCE IT IS NOT USED SUBSEQUENTLY.
C
      ISING = -1
      CALL INVERD (10,HUQ(1),10,DUM,0,DETERM,ISING,TEMP48(1))
C
C     CHECK SINGULARITY
C
      IF (ISING .NE. 2) GO TO 140
      CALL MESAGE (30,40,NECPT(1))
C
C     SET FLAG FOR FATAL ERROR WHILE ALLOWING ERROR MESSAGES TO
C     ACCUMULATE
C
      NOGO = 1
      RETURN
C
  140 CONTINUE
C
C     CALCULATE GENERALIZED DISPLACEMENT VECTOR(Q)
C
C                       ***     ***
C                       *  T      *
C                       *(E )(U ) *
C                       *      A  *
C        (Q)  =  (HUQ)  *---------*
C                       *  T      *
C                       *(E )(U ) *
C                       *      B  *
C                       ***     ***
C
C           WHERE
C                                0    1    0    0    0    0
C                          T     SP   0    CP   0    0    0
C                MATRIX  (E )  = CP   0   -SP   0    0    0
C                                0    0    0    0    1    0
C                                0    0    0    SP   0    CP
C
C
      K1 = 0
      K2 = 0
  320 U(K2+1) = DBLE(ECPT(K1+40))
      U(K2+2) = DBLE(ECPT(K1+39))*SP + DBLE(ECPT(K1+41))*CP
      U(K2+3) = DBLE(ECPT(K1+39))*CP - DBLE(ECPT(K1+41))*SP
      U(K2+4) = DBLE(ECPT(K1+43))
      U(K2+5) = DBLE(ECPT(K1+42))*SP + DBLE(ECPT(K1+44))*CP
C
      IF (K1 .NE. 0) GO TO 400
      K1 = 6
      K2 = 5
      GO TO 320
C
  400  CALL GMMATD (HUQ(1),8,10,0,U(1),10,1,0,Q(1))
C
C     CALCULATE STRAIN COEFFICIENTS AND OBTAIN MATERIAL PROPERTY MATRIX
C     (E)
C
      MATID  = MATID1
      INFLAG = 2
      ELTEMP = ECPT(35)
      CALL MAT (ECPT(1))
      E11   = G11
      E12   = G12
      E22   = G22
      E33   = G33
      TDIF  = (DBLE(ECPT(38)) - DBLE(ECPT(37)))/SL
      DEPS  = DBLE(ALPHA1)*TDIF
      DEPP  = DBLE(ALPHA2)*TDIF
      EPS   = DBLE(ALPHA1)*DBLE(ECPT(37))
      EPP   = DBLE(ALPHA2)*DBLE(ECPT(37))
C
C     COMPUTE COEFFICIENTS FOR POWER SERIES OF DIFFERENTIAL STIFF. COEFF
C
      TM    = ECPT(5)
      TEMP1 = SP*Q(3) + CP*Q(5)
      TEMP2 = SP*Q(4) + CP*Q(6)
      TEMP3 = Q(4) - EPS
      TE11  = TM*E11
      TE12  = TM*E12
      TE22  = TM*E22
C
      A0    = TE12*TEMP1
      A1    = TE12*TEMP2
      A2    = TE12*CP*Q(7)
      A3    = TE12*CP*Q(8)
      B0    = TE22*TEMP1
      B1    = TE22*TEMP2
      B2    = TE22*CP*Q(7)
      B3    = TE22*CP*Q(8)
      C0    = TE11*TEMP3 - TE12*EPP
      C1    =-TE11*DEPS  - TE12*DEPP
      D0    = TE12*TEMP3 - TE22*EPP
      D1    =-TE12*DEPS  - TE22*DEPP
C
C     COMPUTE DIFFERENTIAL STIFFNESS COEFFICIENTS
C
      DO 500 I = 1,3
      IP1 = I + 1
      IP2 = I + 2
      IP3 = I + 3
      DO 500 J = I,3
      JP1 = J + 1
      A(I,J) = A0*INT(I  ,JP1) + A1*INT(IP1,JP1) + A2*INT(IP2,JP1)
     1       + A3*INT(IP3,JP1) + C0*INT(I  ,J  ) + C1*INT(IP1,J  )
      B(I,J) = B0*INT(I  ,JP1) + B1*INT(IP1,JP1) + B2*INT(IP2,JP1)
     1       + B3*INT(IP3,JP1) + D0*INT(I  ,J  ) + D1*INT(IP1,J  )
  500 C(I,J) = A(I,J) + B(I,J)
C
      J   = 1
      JP1 = 2
      DO 510 I = 2,5
      IP1 = I + 1
      IP2 = I + 2
      IP3 = I + 3
  510 A(I,J) = A0*INT(I  ,JP1) + A1*INT(IP1,JP1) + A2*INT(IP2,JP1)
     1       + A3*INT(IP3,JP1) + C0*INT(I  ,J  ) + C1*INT(IP1,J  )
C
      J   = 3
      JP1 = 4
      DO 520 I = 4,7
      IP1 = I + 1
      IP2 = I + 2
      IP3 = I + 3
  520 B(I,J) = B0*INT(I  ,JP1) + B1*INT(IP1,JP1) + B2*INT(IP2,JP1)
     1       + B3*INT(IP3,JP1) + D0*INT(I  ,J  ) + D1*INT(IP1,J  )
C
C     COMPUTE KQD
C     FOR EXPLICIT FORMULATION OF KQD, SEE MS-31, PP. 8-11
C     CASE ONE.. HARMONIC NUMBER = ZERO
C
      DO 600 I = 1,64
  600 KQD(I)  = 0.0D0
      SP2D4   = SP2*0.25D0
      KQD( 1) = CP2*B(1,3) + SP2D4*C(1,3)
      KQD( 2) = CP2*B(2,3) + 0.25D0*SP*C(1,2) + SP2D4*C(2,3)
      KQD( 9) = KQD(2)
      KQD(10) = CP2*B(3,3) + (C(1,1) + 2.0D0*SP*C(2,2)
     $        + SP2*C(3,3))*0.25D0
      KQD(46) = A(1,1)
      KQD(47) = A(2,1)*2.0D0
      KQD(48) = A(3,1)*3.0D0
      KQD(54) = KQD(47)
      KQD(55) = A(3,1)*4.0D0
      KQD(56) = A(4,1)*6.0D0
      KQD(62) = KQD(48)
      KQD(63) = KQD(56)
      KQD(64) = A(5,1)*9.0D0
C
C     CHECK HARMONIC NUMBER
C
      IF (N .EQ. 0.0D0) GO TO 800
C
C     CASE TWO.. HARMONIC NUMBER .NE. ZERO
C
      NOV4    = N*0.25D0
      NSPOV4  = NOV4*SP
      NCP     = N*CP
      N2OV4   = NOV4*N
      KQD( 3) = NSPOV4*C(1,3)
      KQD( 4) = NSPOV4*C(2,3)
      KQD( 5) = NCP*B(1,3)
      KQD( 6) = NCP*B(2,3)
      KQD( 7) = NCP*B(3,3)
      KQD( 8) = NCP*B(4,3)
      KQD(11) = NOV4*(C(1,2) + SP*C(2,3))
      KQD(12) = NOV4*(C(2,2) + SP*C(3,3))
      KQD(13) = NCP*B(2,3)
      KQD(14) = NCP*B(3,3)
      KQD(15) = NCP*B(4,3)
      KQD(16) = NCP*B(5,3)
      KQD(17) = KQD( 3)
      KQD(18) = KQD(11)
      KQD(19) = N2OV4*C(1,3)
      KQD(20) = N2OV4*C(2,3)
      KQD(25) = KQD( 4)
      KQD(26) = KQD(12)
      KQD(27) = KQD(20)
      KQD(28) = N2OV4*C(3,3)
      KQD(33) = KQD( 5)
      KQD(34) = KQD(13)
      KQD(37) = N2*B(1,3)
      KQD(38) = N2*B(2,3)
      KQD(39) = N2*B(3,3)
      KQD(40) = N2*B(4,3)
      KQD(41) = KQD( 6)
      KQD(42) = KQD(14)
      KQD(45) = KQD(38)
      KQD(46) = KQD(46) + N2*B(3,3)
      KQD(47) = KQD(47) + N2*B(4,3)
      KQD(48) = KQD(48) + N2*B(5,3)
      KQD(49) = KQD( 7)
      KQD(50) = KQD(15)
      KQD(53) = KQD(39)
      KQD(54) = KQD(47)
      KQD(55) = KQD(55) + N2*B(5,3)
      KQD(56) = KQD(56) + N2*B(6,3)
      KQD(57) = KQD( 8)
      KQD(58) = KQD(16)
      KQD(61) = KQD(40)
      KQD(62) = KQD(48)
      KQD(63) = KQD(56)
      KQD(64) = KQD(64) + N2*B(7,3)
C
C     COMPUTE HUQ FOR NTH HARMONIC
C
      DO 690 I = 1,100
  690 HUQ(  I) = 0.0D0
      HUQ(  1) = ONE
      HUQ( 13) = ONE
      HUQ( 25) = ONE
      HUQ( 36) = ONE
      HUQ( 41) = CP/RA
      HUQ( 45) = N /RA
      HUQ( 49) = ONE
      HUQ( 51) = ONE
      HUQ( 52) = SL
      HUQ( 63) = ONE
      HUQ( 64) = SL
      HUQ( 75) = ONE
      HUQ( 76) = SL
      HUQ( 77) = L2
      HUQ( 78) = HUQ(77)*SL
      HUQ( 86) = ONE
      HUQ( 87) = 2.0D0*SL
      HUQ( 88) = 3.0D0*HUQ(77)
      HUQ( 91) = CP/RB
      HUQ( 92) = HUQ(91)*SL
      HUQ( 95) = N /RB
      HUQ( 96) = HUQ(95)*SL
      HUQ( 97) = HUQ(95)*L2
      HUQ( 98) = HUQ(96)*L2
      HUQ( 99) = ONE
      HUQ(100) = SL
C
C     COMPUTE HYQ
C
      IF (MATID2.EQ.0 .OR. MATID3.EQ.0) GO TO 710
      IF (ECPT(9).EQ.0.0 .OR. ECPT(7).EQ.0.0) GO TO 710
      IF (GSHEAR .EQ. 0.0D0) GO TO 710
C
      N2D33  = N2 *D33
      SP2D22 = SP2*D22
      OQ     = SL*TS*GSHEAR*(RA+RB)*0.5D0 + INT(1,3)*(N2D33+SP2D22)*OPI
      OQ     = ONE/OQ
      NSP    = N*SP
      NCP    = N*CP
      NSPOPI = NSP*OPI
      TWOD33 = 2.0D0*D33
      TEMP1  = D12*(ONE/RB - ONE/RA)
      TEMP2  = NSPOPI*(D22 + D33)
      TEMP3  = N*NSPOPI*(TWOD33 + D22)
      TEMP4  = OQ*0.5D0*N2D33*CP*OPI
      TEMP5  = OPI*(N2*TWOD33  + SP2D22)
      TEMP6  = D12*N2*L2/RB
      TEMP7  = NSPOPI*CP*0.50D0
C
      HYQ( 1) = OQ*(TEMP1*NCP - TEMP7*INT(1,4)*(D33 + 2.0D0*D22))
      HYQ( 2) = OQ*(NCP*SL/RB*D12 - TEMP7*INT(2,4)*(3.0D0*D33 + D22) +
     1          1.5D0*NCP*OPI*INT(1,3)*D33)
      HYQ( 3) = TEMP4*INT(1,4)
      HYQ( 4) = TEMP4*INT(2,4)
      HYQ( 5) = OQ*(TEMP1*N2 - TEMP3*INT(1,4))
      HYQ( 6) = OQ*(D12*N2*SL/RB - TEMP3*INT(2,4) + TEMP5*INT(1,3))
      HYQ( 7) = OQ*(2.0D0*D11*(RA-RB) + TEMP6 + 2.0D0*INT(2,3)*TEMP5
     $        - TEMP3*INT(3,4))
      HYQ( 8) = OQ*(-D11*6.0D0*SL*RB + TEMP6*SL + 3.0D0*INT(3,3)*TEMP5
     $        - TEMP3*INT(4,4))
      HYQ( 9) =-OQ*TEMP2*INT(1,3)
      HYQ(10) = OQ*(N*SL*(D12 + D33) - TEMP2*INT(2,3))
C
      DO 700 I = 1,10
      HUQ(I+30) = HUQ(I+30) - HYQ(I)
  700 HUQ(I+80) = HUQ(I+80) - HYQ(I)
C
  710 CONTINUE
C
C     AGAIN SET ISING TO -1
C
      ISING = -1
      CALL INVERD (10,HUQ(1),10,DUM,0,DETERM,ISING,TEMP48(1))
      IF (ISING .NE. 2) GO TO 720
      CALL MESAGE (30,40,NECPT(1))
C
C     SET FLAG FOR FATAL ERROR WHILE ALLOWING ERROR MESSAGES TO
C     ACCUMULATE
C
      NOGO = 1
      RETURN
C
  720 CONTINUE
C
C     COMPLETE SOLUTION BY TRANSFORMING KQD TO GLOBAL COORDINATES
C
C                     T            T
C       (K  ) = (E)(H  )(KQD)(H)(E )     FOR I = PIVOT POINT
C         IJ         I         J             J = A,B
C
C     FIRST OBTAIN PRODUCTS
C                       T
C        EHAT  =  (E)(H  )      AND STORE AT EHT(1) . . . EHT(48)
C                      A
C
C                       T
C        EHBT  =  (E)(H  )      AND STORE AT EHT(49). . . EHT(96)
C                      B
C
C              0    CP   SP   0    0
C              1    0    0    0    0
C
C              0    CP  -SP   0    0
C
C        MATRIX E =
C              0    0    0    0    SP
C
C              0    0    0    1    0
C
C              0    0    0    0    CP
C
  800 INC1 = 0
      INC2 = 0
  810 DO 820 I = 1,8
      KROW = I + INC1
      NCOL = (I-1)*10 + INC2
      EHT(KROW   ) = SP*HUQ(NCOL+2) + CP*HUQ(NCOL+3)
      EHT(KROW+ 8) =    HUQ(NCOL+1)
      EHT(KROW+16) = CP*HUQ(NCOL+2) - SP*HUQ(NCOL+3)
      EHT(KROW+24) = SP*HUQ(NCOL+5)
      EHT(KROW+32) =    HUQ(NCOL+4)
  820 EHT(KROW+40) = CP*HUQ(NCOL+5)
      IF(INC1 .GT. 0) GO TO 830
      INC1 = 48
      INC2 = 5
      GO TO 810
C
C     CHECK FOR PIVOT POINT NUMBER
C
  830 DO 840 I = 1,2
      IF (NPVT .EQ. NECPT(I+1)) GO TO 850
  840 CONTINUE
C
C     FALL THRU LOOP IMPLIES NO PIVOT POINT NUMBER
C
      CALL MESAGE (-30,34,ECPT(1))
C
  850 NPIVOT = I
      CALL GMMATD ( EHT(48*NPIVOT-47),6,8,0, KQD(1),8,8,0, TEMP48(1) )
C
C     IF N = 0 DOUBLE RESULT
C
      IF (N .NE. 0.0D0) GO TO 870
      DO 860 I = 1,48
  860 TEMP48(I) = 2.0D0*TEMP48(I)
C
  870 DO 880 J = 1,2
      CALL GMMATD (TEMP48(1),6,8,0,EHT(48*J-47),6,8,1,KIJ(1))
  880 CALL DS1B (KIJ(1),NECPT(J+1))
C
      RETURN
      END


text:
Cedars are attractive, pleasantsmelling conifers. The needles are shorter than those of Pines.  You should point them out by extending an arm, index finger directed at one and loudly exclaiming See deres a tree.

 Species

And where you can find them.

  Atlas Cedar  There are some big ones around the Memorial Union and on the north side of the Quad.
  Deodar Cedar  Also known as the California Christmas Tree. They can be grown from cuttings; intriguingly, cuttings from the trees top branches will grow tall and straight, while cuttings from the older, horizonal bottom branches will grow more outward and twisted, resulting in a drooping, gnarled, spooky tree that looks something from a Tim Burton movie.
  Incense Cedar  These are a big deal up in the mountains. The bark is reddish and a bit shreddedlooking. The tree smells nice.

 See Also
  Town Flora

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    %W% %G%
      subroutine cvsol(icv, base, xi, xo)

      implicit none

      include 'tspinc/params.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/cv.inc'

      integer icv, idat, i
      real    xi, xo, xipu, xopu, base

c     abort if icv .eq. 0
      if (icv.eq.0) then
        xo = xi
        return
      endif

c     process hardcoded curves (icv < 0)
      if (icv.lt.0) then

        idat = abs(icv)

        xipu = xi/base

        if (xipu.le.inpt(idat,1)) then
          xopu = outpt(idat,1)
        else if (xipu.ge.inpt(idat,totalpts(idat))) then
          xopu = outpt(idat,totalpts(idat))
        else
          do i = 2, totalpts(idat)
            if (xipu.lt.inpt(idat,i)) then
c             interpolate between i-1 and i
              xopu =   outpt(idat,i-1) + (xipu - inpt(idat,i-1))
     &                               *(outpt(idat,i) - outpt(idat,i-1))
     &                               /(inpt(idat,i)  - inpt(idat,i-1) )
              xo = xopu*base
              return
            endif
          enddo
        endif
        xo = xopu*base
        return
      endif

c     get index to cvin and cvout arrays
      idat = cvindx(icv)

c     compute output
      if (xi.le.cvin(idat+1)) then
        xo = cvout(idat+1)
      else if (xi.ge.cvin(idat+cvpts(icv))) then
        xo = cvout(idat+cvpts(icv))
      else
        do i = 2, cvpts(icv)
          if (xi.lt.cvin(idat+i)) then
c           interpolate between i-1 and i
            xo = cvout(idat+i-1) + (xi - cvin(idat+i-1) )
     &                            * ( cvout(idat+i) - cvout(idat+i-1) )
     &                            / ( cvin(idat +i) - cvin(idat+ i-1) )
            return
          endif
        enddo
      endif

      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:
C     LAST UPDATE 10/12/96
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE LINES(BUF,GSTOPN,VALS,NV)
      IMPLICIT NONE
C
C Purpose: Does line plots and peak fits.
C
C Calls   8: ASK    , PAPER  , PSPACE , FILNAM , FPLOT  , ERASE  ,
C            GREND  , RDCOMF
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C FIX include file:
C
      INCLUDE 'FIXPAR.COM'
      INCLUDE 'FIT.COM'
C
C Arguments:
C
      REAL BUF(NPIX*NRAST)
      LOGICAL GSTOPN
      REAL VALS(10)
      INTEGER NV
C
C Local variables:
C
      REAL LBUF(MAXDIM)
      REAL*8 ASTDER(MAXPAR)
      REAL AREA,PI,X,Y,DX,DY,VECX,VECY,DIST,SHAPE,X0,Y0,VECXN,VECYN,
     &     DISTN
      INTEGER IVAL(10),ITEM(1)
      INTEGER I,J,K,L,M,IX,IY,IRC,KUNIT,NCHAN,N,NW
      CHARACTER*10 WORD
C
C FPLOT declarations:
C
      REAL*8 A(MAXPAR)
      REAL XPOS(MAXDIM)
      INTEGER NFRAMI,IFLAG
      LOGICAL REPT,FIT,HCOPY,INIT,XAXIS,AUTO
      CHARACTER*12 CTYPE(5)
c++++++++GHOST++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      CHARACTER*80 NAMFIL
c++++++++END++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C External function:
C
      REAL PSNINT,VGTINT
      EXTERNAL PSNINT,VGTINT
c++++++++PGPLOT+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c      INTEGER PGBEG
c      EXTERNAL PGBEG
c++++++++END++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C Common block:
C
      INTEGER NPEAKS,NBAK,FEXP,NSELPK,KORIG
      LOGICAL LATCON
      DOUBLE PRECISION XORIG
      INTEGER FTYPE(20),VTABLE(MAXPAR)
      COMMON /FCOM  /NPEAKS,NBAK,FEXP,FTYPE,VTABLE,NSELPK
      COMMON /FCOM2 /XORIG,KORIG,LATCON
C
C Data:
C
      DATA  KUNIT/12/ , NFRAMI /0/
      DATA  PI /3.1415927/
      DATA  CTYPE /'Gaussian    ' , 'Lorentzian  ' , 'Pearson VII ' ,
     &             'Voigt       ' , 'Debye       '/
      DATA INIT /.FALSE./ , REPT /.FALSE./ , XAXIS /.FALSE./ ,
     &     AUTO /.FALSE./
C
C-----------------------------------------------------------------------
      IF(NLIN.EQ.0)RETURN
C
C========Initialize graphics
C
      IF(.NOT.GSTOPN)THEN
         CALL CFREE
c++++++++GHOST++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
         CALL PAPER(1)
         CALL PSPACE(0.1,0.9,0.1,0.9)
         CALL GPSTOP(0)
c++++++++PGPLOT+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c         IRC = PGBEG(0,'/XWINDOW',1,1)
c         IF(IRC.NE.1)GOTO 9999
c         CALL PGPAP(6.0,0.7)
c         CALL PGASK(.FALSE.)
c         CALL PGSVP(0.1,0.9,0.1,0.9)
c         CALL PGPAGE
c++++++++END++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
         GSTOPN = .TRUE.
      ENDIF
      IF(NV.EQ.0)THEN
         IVAL(1) = 1
         IVAL(2) = NLIN
         NV = 2
      ELSE
         DO 25 I=1,NV
            IVAL(I) = NINT(VALS(I))
 25      CONTINUE
      ENDIF
C
C========Prompt for peak fitting
C
      FIT = .FALSE.
      CALL ASK('702 Do you want to fit the data',FIT,0)
C
C========Prompt for hardcopy output
C
      HCOPY = .FALSE.
c++++++++GHOST+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      CALL ASK('Do you require hardcopy output',HCOPY,0)
      IF(HCOPY)THEN
         WRITE(6,1000)
         READ(5,'(A80)',ERR=9999,END=9999)NAMFIL
         CALL FILNAM(NAMFIL)
      ENDIF
c++++++++PGPLOT+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C========No hardcopy at present
c++++++++END++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      WRITE(6,'(/)')
C
C========Start main loop
C
      DO 60 I=1,NV,2
         IF(I+1.GT.NV)IVAL(I+1) = IVAL(I)
         DO 50 J=IVAL(I),IVAL(I+1)
C
C========Interpolate between end points
C
            VECX = XLINE(2,J) - XLINE(1,J)
            VECY = YLINE(2,J) - YLINE(1,J)
            DIST = SQRT(VECX*VECX+VECY*VECY)
            IF(INT(DIST).GT.MAXDIM)THEN
               NCHAN = MAXDIM
               DIST = FLOAT(NCHAN)
            ELSE
               NCHAN = INT(DIST)
            ENDIF
            IF(NCHAN.GT.0)THEN
               VECX = VECX/DIST
               VECY = VECY/DIST
            ELSE
               GOTO 50
            ENDIF
            IF(ABS(VECY).GT.0.0)THEN
               VECXN = 1.0
               VECYN = -VECX/VECY
               DISTN = SQRT(VECXN*VECXN+VECYN*VECYN)
               VECXN = VECXN/DISTN
               VECYN = VECYN/DISTN
            ELSE
               VECYN = 1.0
               VECXN = 0.0
            ENDIF
            DO 30 K=1,NCHAN
               X0 = XLINE(1,J) + FLOAT(K)*VECX
               Y0 = YLINE(1,J) + FLOAT(K)*VECY
               LBUF(K) = 0.0
               DO 27 L=-LWIDTH(J),LWIDTH(J)
                  X = X0 + FLOAT(L)*VECXN
                  Y = Y0 + FLOAT(L)*VECYN
                  IX = INT(X)
                  IY = INT(Y)
                  IF(IX.GE.1.AND.IX.LT.NPIX.AND.
     &               IY.GE.1.AND.IY.LT.NRAST)THEN
                     M = (IY-1)*NPIX + IX
                     DX = X - FLOAT(IX)
                     DY = Y - FLOAT(IY)
                     LBUF(K) = LBUF(K) + (1.0-DX)*(1.0-DY)*BUF(M)
     &                                 +       DX*(1.0-DY)*BUF(M+1)
     &                                 +       (1.0-DX)*DY*BUF(M+NPIX)
     &                                 +             DX*DY*BUF(M+NPIX+1)
                  ENDIF
 27            CONTINUE
 30         CONTINUE
C
C========Do plotting and fitting
C
            CALL FPLOT(LBUF,NCHAN,REPT,FIT,A,ASTDER,INIT,NFRAMI,
     &                 XPOS,XAXIS,AUTO,IFLAG)
C
C========Exit if GUI demands it (i.e. ^D received)
C
            IF(IFLAG.EQ.-2)GOTO 60
C
C========Interpret A and write out parameters
C
            IF(FIT.AND.IFLAG.EQ.0)THEN
               WRITE(KUNIT,1040)J,XLINE(1,J),YLINE(1,J),
     &                          XLINE(2,J),YLINE(2,J)
               WRITE(6,1040)J,XLINE(1,J),YLINE(1,J),
     &                      XLINE(2,J),YLINE(2,J)
               CALL FLUSH(6)
               M = 0
               DO 40 L=1,NPEAKS
                  IF(FTYPE(L).EQ.1.OR.FTYPE(L).EQ.2)THEN
                     A(M+3) = ABS(A(M+3))
                     IF(FTYPE(L).EQ.1)THEN
                        AREA = A(M+1)*A(M+3)*SQRT(PI)
                     ELSE
                        AREA = A(M+1)*A(M+3)*PI
                     ENDIF
                     SHAPE = 0.0
                     N = 3
                  ELSEIF(FTYPE(L).EQ.3.OR.FTYPE(L).EQ.4)THEN
                     IF(FTYPE(L).EQ.3)THEN
                        AREA = PSNINT(A(M+1),A(M+3),A(M+4))
                     ELSE
                        AREA = VGTINT(A(M+1),A(M+3),A(M+4))
                     ENDIF
                     SHAPE = A(M+4)
                     N = 4
                  ELSEIF(FTYPE(L).EQ.5)THEN
                     AREA = A(M+1)*A(M+3)*SQRT(PI)*8.0/3.0
                     SHAPE = 0.0
                     N = 3
                  ENDIF
                  WRITE(KUNIT,1050)L,(A(M+K),K=1,3),SHAPE,AREA,
     &                 CTYPE(FTYPE(L))
                  WRITE(6,1050)L,(A(M+K),K=1,3),SHAPE,AREA,
     &                 CTYPE(FTYPE(L))
                  WRITE(KUNIT,1055)(ASTDER(M+K),K=1,N)
                  WRITE(6,1055)(ASTDER(M+K),K=1,N)
                  CALL FLUSH(6)
                  M = M + N
 40            CONTINUE
C
C========Write out background polynomial coefficients
C
               WRITE(6,1060)NBAK,(K,K=0,4)
               WRITE(KUNIT,1060)NBAK,(K,K=0,4)
               WRITE(6,1070)(A(M+K),K=NBAK+1,1,-1)
               WRITE(KUNIT,1070)(A(M+K),K=NBAK+1,1,-1)
               WRITE(6,1075)(ASTDER(M+K),K=NBAK+1,1,-1)
               WRITE(KUNIT,1075)(ASTDER(M+K),K=NBAK+1,1,-1)
               IF(FEXP.GT.0)THEN
                  WRITE(6,1080)A(M+NBAK+2),A(M+NBAK+3)
                  WRITE(KUNIT,1080)A(M+NBAK+2),A(M+NBAK+3)
                  WRITE(6,1085)ASTDER(M+NBAK+2),ASTDER(M+NBAK+3)
                  WRITE(KUNIT,1085)ASTDER(M+NBAK+2),ASTDER(M+NBAK+3)
               ENDIF
               IF(LATCON)THEN
                  WRITE(6,1090)XORIG
                  WRITE(KUNIT,1090)XORIG
                  IF(KORIG.GT.0)THEN
                     WRITE(6,1095)ASTDER(KORIG)
                     WRITE(KUNIT,1095)ASTDER(KORIG)
                  ENDIF
               ENDIF
               CALL FLUSH(6)
            ENDIF
C
C========Wait until ready for next plot
C
            IF(.NOT.FIT)THEN
               WRITE(6,1100)
               CALL FLUSH(6)
               CALL RDCOMF(5,WORD,NW,VALS,NV,ITEM,1,1,10,IRC)
               IF(IRC.NE.0)THEN
                  IFLAG = -2
                  GOTO 60
               ELSEIF(NW.GT.0)THEN
                  IF(WORD(1:2).EQ.'^D')THEN
                     IFLAG = -2
                     GOTO 60
                  ENDIF
               ENDIF
            ENDIF
c++++++++GHOST++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            IF(HCOPY)CALL FILEND
c++++++++PGPLOT+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C========No hardcopy at present
c++++++++END++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 50      CONTINUE
 60   CONTINUE
 9999 RETURN
C
 1000 FORMAT(1X,'900',
     &       1X,'Enter gridfile name: ',$)
 1040 FORMAT(/1X,'200'/
     &       1X,'Line ',I3,':  (',F6.1,',',F6.1,')    to    (',
     &       F6.1,',',F6.1,')'/1X,8('-'),
     &       3X,'Height      Position     Width      Shape   '
     &       4X,'Integrated      Type')
 1050 FORMAT(/1X,'Peak ',I3,1X,5G12.5,4X,A12)
 1055 FORMAT(1X,'Std Err',2X,5G12.5)
 1060 FORMAT(/1X,'Background polynomial degree ',I3/
     &       1X,'200',6X,5(4X,'a',I1,6X))
 1070 FORMAT(1X,'Coeff',4X,5G12.5)
 1075 FORMAT(1X,'Std Err',2X,5G12.5)
 1080 FORMAT(/1X,'Exponential background'/
     &       10X,G12.5,' * exp(',G12.5,' * x)')
 1085 FORMAT(1X,'Std Err',2X,G12.5,7X,G12.5)
 1090 FORMAT(/1X,'Origin for lattice constraints = ',G12.5)
 1095 FORMAT(1X,'                       Std Err = ',G12.5)
 1100 FORMAT(1X,'806',
     &       1X,'Hit return for next frame or <ctrl-D> to exit'/
     &       1X,'200',
     &       1X,'Waiting to continue...')
 2000 FORMAT(1X,'***Numeric input expected')
      END


program main
    use FEMDomainClass

    type(FEMDomain_)::omega
    real(8),allocatable :: nod_coord(:,:)

    call ImportArray(nod_coord,FileName="test.mesh")
    call omega%Mesh%ImportNodCoord(nod_coord)
    call omega%Meshing()

end program main


! PR middle-end/27416
! { dg-do run }

  integer :: j
  j = 6
!$omp parallel num_threads (4)
  call foo (j)
!$omp end parallel
  if (j.ne.6+16) call abort
end

subroutine foo (j)
  integer :: i, j

!$omp do firstprivate (j) lastprivate (j)
  do i = 1, 16
    if (i.eq.16) j = j + i
  end do
end subroutine foo


     program demo_elips
     use M_calcomp, only : plots, plot, newpen, factor, nframe, elips
     implicit none
     character(len=:),allocatable :: lines(:)
     integer                  :: i
     real                     :: x, y
     real                     :: a
     real                     :: major_axis_length, minor_axis_length
     real                     :: start_angle, finish_angle
     integer                  :: ipen
     integer,parameter        :: END=999, MOVE=3, DRAW=2
     lines=[character(len=80) :: &
     '#---------------------------------------------------------------------#-',&
     '#  move over in the x direction the same amount you change axis length  ',&
     '#---------#---------#---------#---------#---------#---------#---------#-',&
     '|6.5      !8.0      !0.5      !0.7      ! 0.0     ! 0.0     ! 360.0   !3',&
     '|6.6      !8.0      !0.6      !0.6      ! 0.0     ! 0.0     ! 360.0   !3',&
     '|6.7      !8.0      !0.7      !0.5      ! 0.0     ! 0.0     ! 360.0   !3',&
     '|6.8      !8.0      !0.8      !0.4      ! 0.0     ! 0.0     ! 360.0   !3',&
     '|6.9      !8.0      !0.9      !0.3      ! 0.0     ! 0.0     ! 360.0   !3',&
     '|7.0      !8.0      !1.0      !0.2      ! 0.0     ! 0.0     ! 360.0   !3',&
     '#---------#---------#---------#---------#---------#---------#---------#-',&
     '#  different end angles of different signs                              ',&
     '#---------#---------#---------#---------#---------#---------#---------#-',&
     '|5.0      !8.0      !1.0      !0.2      ! 0.0     ! 0.0     !  45.0   !3',&
     '|3.0      !8.0      !1.0      !0.2      ! 0.0     ! 0.0     ! -60.0   !3',&
     '#---------#---------#---------#---------#---------#---------#---------#-',&
     '# circles                                                               ',&
     '#---------#---------#---------#---------#---------#---------#---------#-',&
     '#---------#---------#---------#---------#---------#---------#---------#-',&
     '#  end of values to call ELIPS(3f) with                                 ',&
     '#---------------------------------------------------------------------#-']
        call plots(0.0,10.0,0.0,10.0)
        do i=1,size(lines)
           write(*,'(a)')lines(i)
           if( lines(i).eq.' ' .or. index(lines(i),'#').ne.0 )cycle
           read(lines(i),9012) x,y, &
                           & major_axis_length,minor_axis_length, &
                           & a, &
                           & start_angle,finish_angle, &
                           & ipen
           9012 format(7(1X,F9.3),1X,I1)
           call elips(x,y, &
                           & major_axis_length,minor_axis_length, &
                           & a, &
                           & start_angle,finish_angle, &
                           & ipen)
        enddo
        call nframe()
        call plot(0.0,0.0,END)
     end program demo_elips


text:
program main
  integer :: i = 1
  integer :: j = 10

  do while (i < j)
     i = i + 1
  end do

end program main

max_stars_repo_name: ouankou/rose

max_issues_repo_name: ouankou/rose

max_forks_repo_name: ouankou/rose


text:
      PROGRAM MAIN

C     RESULT REPORTING
C     CHARACTER(LEN=*), PARAMETER  :: FMT1 = "(F18.15)"
      CHARACTER(LEN=*), PARAMETER  :: FMT1 = "(E12.6)"

C     ----------------------------------------------
C     Variables and data for VMF1
C     ----------------------------------------------
      DOUBLE PRECISION VMF1H, VMF1W

      CALL VMF1(0.00127683D0, 0.00060955D0, 55055D0, 0.6708665767D0,
     .  1.278564131D0, VMF1H, VMF1W)

C     REPORT RESULTS
      PRINT *, '----------------------------------------'
      PRINT *, '> VMF1 Results:'
      PRINT *, '----------------------------------------'
      WRITE(*, FMT1) DABS(VMF1H-3.424342122738070593D0)
      WRITE(*, FMT1) DABS(VMF1W-3.448299714692572238D0)

C     ALL DONE
      END


text:
Describe Users/BobbyCoyote here.

20070919 23:54:44 nbsp Welcome to the Wiki. Im sorry that your introduction to editing the wiki came as the result of such an incident, but your comment demonstrates your commitment to customer service. My only suggestion would be that it would have been easier to read if it were broken into paragraphs, but it was worth reading even without them. The {{{Footnote()}}} macro footnote(Like this.) is also useful for formatting.
I hope that this experience hasnt left you feeling too negative about the wiki. Pages like Welcome to the Wiki/Business Owner and Reviews show that there are elements of the wiki community working to find a way to value the contribution of all users while dealing with people who only leave negative comments, spam, and or incorrect information.
One of the techniques that has been tried is to move older comments off to a comment archive page. This preserves the comment, but keeps a business from having older and potentially out of date comments near the top of the comment section. There are also some editors who will push back and leave comments for those who only leave negative comments. Users/JasonAller


text:
!==============================================================================!
  subroutine Backup_Mod_Write_Bnd(comm, fh, disp, vc, var_name, array)
!------------------------------------------------------------------------------!
!   Writes a vector variable with boundary cells to backup file.               !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Comm_Type)  :: comm
  integer          :: fh, disp, vc
  character(len=*) :: var_name
  real             :: array(-comm % nb_s:-1)
!-----------------------------------[Locals]-----------------------------------!
  character(len=80) :: vn
  integer           :: vs  ! variable size
!==============================================================================!

  if(this_proc < 2) print *, '# Writing variable: ', trim(var_name)

  ! Increase variable count
  vc = vc + 1

  ! Vector with boundaries
  vn = var_name
  call Comm_Mod_Write_Text(fh, vn, disp)
  vs = comm % nb_t * SIZE_REAL
  call Comm_Mod_Write_Int (fh, vs, disp)

  call Comm_Mod_Write_Bnd_Real(comm, fh, array(-comm % nb_s:-1), disp)

  end subroutine

max_stars_repo_name: EdinSmartLab/T-Flows

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
c
c      TESTING SUITE FOR HELMHOLTZ SUBROUTINE LIBRARY:  Mar 20, 2019
c
c
c
c         Sources at S, Targets at T.
c
c         h3dformmp      forms expansion about c0
c         h3dmpeval      evaluates expansion at T
c         h3dmpmp        shifts expansion to c1
c         h3dmploc       converts to local about y1
c         h3dtaeval      evaluates expansion at T
c         h3dlocloc      shifts to local about c0
c         h3dtaeval      evaluates expansion at T
c         h3dformta      forms local expansion about c0
c                        directly from S
c         h3dtaeval      evaluates expansion at T
c
c         -------------------------------------------------
c         |       |   S   |       |       |T      |       |
c         |       |   c0  |       |       |  c3   |       |
c         |       |       |       |       |       |       |
c         -------c1-------------------------------c2------
c         |       |       |       |       |       |       |
c         |       |       |       |       |       |       |
c         |       |       |       |       |       |       |
c         -------------------------------------------------
c
      implicit real *8 (a-h,o-z)
      real *8 ztrg(3),sources(3,10)
	  real *8 c0(3),c1(3),c2(3),c3(3)
      real *8 xnodes(2000),wts(2000)
c
      complex *16 pot,fld(3),opot,ofld(3)
      complex *16, allocatable :: mpole1(:,:),mpole2(:,:)
      complex *16, allocatable :: locexp1(:,:),locexp2(:,:)
      complex *16 charge(100),dipvec(3,100)
      real *8 wlege(100000),errs(2,5),err_exp(5)
      integer ipass(5)
      complex *16 zk,eye
c
      data eye/(0.0d0,1.0d0)/
c
      done=1
      pi=4.0*atan(done)
      zk = 0.2d0

      write(*,*) "=========================================="
      write(*,*) "Testing suite for helmrouts3d"

      open(unit=33,file='print_testres.txt',access='append')

c
c
c       set scale parameter neq 1 to make sure it is used correctly.
c
      bsize = 1.0d0
      rscale1 = abs(zk)*bsize/4
      rscale2 = abs(zk)*bsize/2

cc      rscale1 = 1
cc      rscale2 = 1
c
c       set center for original multipole expansion
c
      c0(1)=0.248d0
      c0(2)=0.251d0
      c0(3)=0.249d0
c
c       create two charge sources.
c
      sources(1,1)=c0(1)+0.24d0
      sources(2,1)=c0(2)+ 0.25d0
      sources(3,1)=c0(3)+ 0.25d0
      charge(1)= 1.0d0 + eye*0.2d0
      dipvec(1,1) = hkrand(0) + eye*hkrand(0)
      dipvec(2,1) = hkrand(0) + eye*hkrand(0)
      dipvec(3,1) = hkrand(0) + eye*hkrand(0)

      ns = 1
      sources(1,2)=c0(1)+0.25d0
      sources(2,2)=c0(2)-0.25d0
      sources(3,2)=c0(3)-0.25d0
      charge(2)= -1.0d0 + eye*0.1d0
      dipvec(1,2) = hkrand(0) + eye*hkrand(0)
      dipvec(2,2) = hkrand(0) + eye*hkrand(0)
      dipvec(3,2) = hkrand(0) + eye*hkrand(0)
      ns = 2
c
c       create center for shifted expansions
c       (location jiggled for good measure)
c
      c1(1)= 0.015d0
      c1(2)= 0.012d0
      c1(3)= 0.013d0
      c2(1)= 2.015d0
      c2(2)= 0.012d0
      c2(3)= 0.013d0
      c3(1)= c2(1)-0.25d0
      c3(2)= c2(2)-0.249d0
      c3(3)= c2(3)-0.251d0
c
c       create target
c
      ztrg(1)=c3(1)-0.245d0
      ztrg(2)=c3(2)-0.25d0
      ztrg(3)=c3(3)-0.25d0
      nt = 1

      ntest = 5
      do i=1,ntest
        ipass(i) = 0
      enddo
c
c       direct calculation:
c
      opot = 0
      ofld(1) = 0
      ofld(2) = 0
      ofld(3) = 0
      thresh = 1.0d-15
      nd = 1
      call h3ddirectcdg(nd,zk,sources,charge,
     1       dipvec,ns,ztrg,nt,opot,
     1       ofld,thresh)

      eps = 0.5d-12
      call h3dterms(bsize/2, zk, eps, nterms)
      call h3dterms(bsize, zk, eps, nterms2)
      call h3dterms(bsize/2, zk, eps, nterms3)

      nlege = 100
      lw7 = 100000
      call ylgndrfwini(nlege,wlege,lw7,lused7)
c
c
c       create h-expansion:
c
      ifinit = 1
      nquad = 2.0*max(nterms2,nterms3)
      call legewhts(nquad,xnodes,wts,ifinit)
      allocate(mpole1(0:nterms,-nterms:nterms))

      nd = 1
      call mpzero(nd,mpole1,nterms)
      call h3dformmpcd(nd,zk,rscale1,sources,charge,
     1       dipvec,ns,c0,nterms,mpole1,wlege,nlege)

      pot = 0
      fld(1) = 0
      fld(2) = 0
      fld(3) = 0
      nd = 1

      rconv1 = 1.0d0/sqrt(3.0d0)
      rconv2 = sqrt(3.0d0)/2.0d0
      rconv3 = 0.75d0

      call h3dmpevalg(nd,zk,rscale1,c0,mpole1,nterms,ztrg,nt,pot,
     1      fld,wlege,nlege,thresh)

      err_exp(1) = 10*max(rconv1**(nterms),eps)
      write(*,'(a,e11.4)') 'Testing formmp and mpeval, expected error='
     1   ,err_exp(1)
      call errprint(pot,opot,fld,ofld,errs(1,1))
      if(max(errs(1,1),errs(2,1)).lt.err_exp(1)) ipass(1) = 1
      write(*,*)

      if(ipass(1).ne.1) write(33,*) 'Failed formmp and mpeval test'

c
c    mpmp shift
c
      radius = sqrt(3.0d0)/2.0d0
      allocate(mpole2(0:nterms2,-nterms2:nterms2))
      call mpzero(nd,mpole2,nterms2)
      call h3dmpmp(nd,zk,rscale1,c0,mpole1,nterms,
     1       rscale2,c1,mpole2,nterms2,radius,xnodes,wts,nquad)

      pot = 0
      fld(1) = 0
      fld(2) = 0
      fld(3) = 0
      nd = 1
      call h3dmpevalg(nd,zk,rscale2,c1,mpole2,nterms,ztrg,nt,pot,
     1      fld,wlege,nlege,thresh)

      err_exp(2) = 10*max(rconv2**(nterms),eps)
      write(*,'(a,e11.4)') 'Testing mpmp, expected error='
     1   ,err_exp(2)
      call errprint(pot,opot,fld,ofld,errs(1,2))

      if(max(errs(1,2),errs(2,2)).lt.err_exp(2)) ipass(2) = 1
      write(*,*)
      if(ipass(2).ne.1) write(33,*) 'Failed mpmp test'

c
c    convert to local
c
      radius = sqrt(3.0d0)/2.0d0
      nquad = 2.2*nterms2
      call legewhts(nquad,xnodes,wts,ifinit)

      allocate(locexp2(0:nterms2,-nterms2:nterms2))
      call mpzero(nd,locexp2,nterms2)
      call h3dmploc(nd,zk,rscale2,c1,mpole2,nterms2,
     1      rscale2,c2,locexp2,nterms2,radius,xnodes,wts,nquad)

      pot = 0
      fld(1) = 0
      fld(2) = 0
      fld(3) = 0
      nd = 1
      call h3dtaevalg(nd,zk,rscale2,c2,locexp2,nterms2,ztrg,nt,
     1      pot,fld,wlege,nlege)

      err_exp(3) = 10*max(rconv3**(nterms),eps)
      write(*,'(a,e11.4)') 'Testing mploc, expected error='
     1   ,err_exp(3)
      call errprint(pot,opot,fld,ofld,errs(1,3))
      if(max(errs(1,3),errs(2,3)).lt.err_exp(3)) ipass(3) = 1
      write(*,*)
      if(ipass(3).ne.1) write(33,*) 'Failed mploc test'

c
c    shift local
c
      radius = sqrt(3.0d0)/4.0d0
      allocate(locexp1(0:nterms3,-nterms3:nterms3))
      call mpzero(nd,locexp1,nterms3)

      call h3dlocloc(nd,zk,rscale2,c2,locexp2,nterms2,
     1      rscale1,c3,locexp1,nterms3,radius,xnodes,wts,nquad)

      pot = 0
      fld(1) = 0
      fld(2) = 0
      fld(3) = 0
      nd = 1
      call h3dtaevalg(nd,zk,rscale1,c3,locexp1,nterms3,ztrg,nt,
     1      pot,fld,wlege,nlege)

      err_exp(4) = 10*max(rconv3**(nterms),eps)
      write(*,'(a,e11.4)') 'Testing locloc, expected error='
     1   ,err_exp(4)
      call errprint(pot,opot,fld,ofld,errs(1,4))
      if(max(errs(1,4),errs(2,4)).lt.err_exp(4)) ipass(4) = 1
      write(*,*)
      if(ipass(4).ne.1) write(33,*) 'Failed locloc test'

c
c
c    create local exp from sources

      call mpzero(nd,locexp1,nterms3)
      nd = 1
      call h3dformtacd(nd,zk,rscale1,sources,charge,
     1       dipvec,ns,c3,nterms3,locexp1,wlege,nlege)
      pot = 0
      fld(1) = 0
      fld(2) = 0
      fld(3) = 0
      nd = 1
      call h3dtaevalg(nd,zk,rscale1,c3,locexp1,nterms3,ztrg,nt,
     1      pot,fld,wlege,nlege)

      err_exp(5) = 10*max(rconv3**(nterms),eps)
      write(*,'(a,e11.4)') 'Testing formta and taeval, expected error='
     1   ,err_exp(5)
      call errprint(pot,opot,fld,ofld,errs(1,5))
      if(max(errs(1,5),errs(2,5)).lt.err_exp(5)) ipass(5) = 1
      write(*,*)
      if(ipass(5).ne.1) write(33,*) 'Failed formta and taeval test'

      isum = 0
      do i=1,ntest
        isum = isum+ipass(i)
      enddo

      write(*,'(a,i1,a,i1,a)') 'Successfully completed ',isum,
     1   ' out of ',ntest,' tests in helmrouts3d testing suite'
      write(33,'(a,i1,a,i1,a)') 'Successfully completed ',isum,
     1   ' out of ',ntest,' tests in helmrouts3d testing suite'
      close(33)

      stop
      end
c
c
c
c
c
C
C
      subroutine errprint(pot,opot,fld,ofld,errs)
      implicit real *8 (a-h,o-z)
      complex *16 pot,opot,fld(3),ofld(3)
      real *8 errs(2)
 1000  format(4D15.5)
      err = 0
      ddd = 0
      err = err + abs(fld(1)-ofld(1))**2
      err = err + abs(fld(2)-ofld(2))**2
      err = err + abs(fld(3)-ofld(3))**2
      ddd = ddd + abs(ofld(1))**2
      ddd = ddd + abs(ofld(2))**2
      ddd = ddd + abs(ofld(3))**2
      err = sqrt(err)
      ddd = sqrt(ddd)

      err1 = abs(pot-opot)/abs(opot)
      write(*,'(a,e11.4,a,e11.4)')
     1     'pot error=',err1,'   grad error=',err/ddd

      errs(1) = err1
      errs(2) = err/ddd

      return
      end
c


text:
!***************************************************************************!
!     George S. Fanourgakis, PNL Laboratory                                 !
!     george.fanourgakis@pnl.gov                                            !
!                                                                           !
!     Disclaimer: This code is provided "as is". No representation is made  !
!                 that is  suitable for any particular purpose              !
!                                                                           !
!     Reference: Fanourgakis, G. S. and Xantheas, S. S.                     !
!                 J. Phys. Chem. A. vol 110, page 4100-4106, year 2006      !
!                                                                           !
!***************************************************************************!
program test_ttm2f
!*** this program is an example of how to use the subroutine "ttm2f" for the
!*** energy and derivatives calculation
!use model_mod
!use ttm2f_mod
use potential_mod
use dfpmin_mod

implicit none
character(len=100) :: fconfig
integer :: iw, Nw, ia, ixyz
double precision, dimension(:), allocatable :: p, g
double precision, dimension(:,:), allocatable :: RR, dRR, dRR0
double precision :: Energy, Energy0
double precision :: delta, diff
double precision :: fret, gtol, tolx, stpmx
integer :: iter, neqs,  iprint
!***
!*** reads the filename of the configuration file
!***
debug = .false.
!print*,'Enter model (2: TTM2-F, 21:TTM21-F,   3:TTM3-F)'
!read(*,*)imodel
imodel=21
!print*,'Enter configuration filename'
IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN
  WRITE(*,*)'ERROR, FILE SHOULD BE PROVIDED ON COMMAND LINE, STOPPING'
  STOP
ENDIF
CALL GET_COMMAND_ARGUMENT(1, fconfig)
!read(*,'(a100)')fconfig
open(10,file=fconfig,status='old')  !
call read_natoms(10, Nw)   ! read the number of water molecules (Nw)
allocate(RR(3, 3*Nw))      ! allocate the array for the coodinates  (RR)
allocate(dRR(3, 3*Nw))     ! allocate the array for the derivatives (dRR)
call readXYZ(10, Nw, RR)   ! read and sort properly the coordinates of atoms
close(10)
!***
! Now we are ready to call the potential
! Print out also the energy (in kcal/mol) and the derivatives(kcal/mol/A)
!***
!call ttm2f(Nw,RR, dRR, Energy)
call potential(Nw,RR, dRR, Energy)
write(*,'("Energy (kcal/mol) = ",f12.6)')Energy
write(*,*)
print*,'----------Derivatives (kcal/mol/A)  ---------------'
!do iw=1,3*Nw
!   write(*,'(i3,3x,3(f12.6,2x))')iw,dRR(1:3,iw)
!enddo
!the lines are stored with oxygens first and then the hydrogens
!but the hydrogens are swapped from the original positions in the file
!do iw=1,Nw
!   write(*,'(i3,3x,3(f12.6,2x))')iw,dRR(1:3,iw)
!   write(*,'(i3,3x,3(f12.6,2x))')iw,dRR(1:3,Nw+2*iw)
!   write(*,'(i3,3x,3(f12.6,2x))')iw,dRR(1:3,Nw+2*iw-1)
!enddo
!***
!***
!***   Test numerical and analytical derivatives
!***
!print*,'----------Test analytical and numerical Derivatives (kcal/mol/A)  ---------------'
!allocate(dRR0(3, 3*nw))
!delta=1.d-8
!call potential(Nw,RR, dRR0, Energy0)
!do ia=1, 3*nw
!   do ixyz=1, 3
!      RR(ixyz, ia) = RR(ixyz, ia) + delta
!      call potential(Nw,RR, dRR, Energy)
!      diff = (Energy-Energy0)/delta
!      RR(ixyz, ia) = RR(ixyz, ia) - delta
!      write(*,'(i5,2x,i2, 5x,"anal= ", f12.6,2x,"num= ", f12.6)')ia, ixyz, dRR0(ixyz, ia), diff
!   enddo
!enddo
!deallocate(dRR0)
!***
!***
!***   Minimize energy
!***
print*,'----------  GEOMETRY OPTIMIZATION -------------------'
Neqs = 3*3*Nw
allocate(p(neqs))
allocate(g(neqs))
p = reshape(RR, (/neqs/))
call min_func(neqs, p, g, fret)
print*,'fret=', fret
gtol = 1.d-8
tolx = 1.d-8
stpmx = 0.01d0
iter=1000
iprint=10

call dfpmin(neqs, p, gtol, tolx, stpmx, iter, fret, iprint, min_func)

RR = reshape(p,(/3, 3*nw/))
call potential(Nw,RR, dRR, Energy)
print*,'Final minimized Energy = ', Energy
open(15, file='minimum.xyz', status='unknown')
write(15,*)3*nw
write(15,*)Energy
do iw=1, nw
   write(15,'(a3,3(1x, f15.8))')'O', RR(1:3, iw)
   write(15,'(a3,3(1x, f15.8))')'H', RR(1:3, iw+Nw)
   write(15,'(a3,3(1x, f15.8))')'H', RR(1:3, iw+2*Nw)
enddo
close(15)

! deallocate all arrays
!***
deallocate(RR)
deallocate(dRR)

!Contains
!   subroutine min_func(n, p, g, fret)
!   implicit none
!   integer, intent(in) :: n
!   double precision, dimension(n), intent(in) :: p
!   double precision, dimension(n), intent(out) :: g
!   double precision, intent(out) :: fret
!   !...
!   integer :: Neqs, iopt, na, nw
!   double precision, dimension(:,:), allocatable :: R, dR
!
!   na = n/3
!   nw = na/3
!   allocate(R(3,na))
!   allocate(dR(3, na))
!   R =reshape(p,(/3, na/))
!
!   call potential(Nw,R, dR, Energy)
!
!   fret = Energy
!   g = reshape(dR,(/n/))
!
!   deallocate(r, dr)
!   end subroutine min_func

end program test_ttm2f

max_stars_repo_name: heindelj/pyMD

max_issues_repo_name: heindelj/pyMD

max_forks_repo_name: heindelj/pyMD


text:
Subroutine newpos(t, i)
  Implicit Double Precision (A-H, O-Z)
  Parameter (maxptn=400001)
  Common /para5/iconfg, iordsc
  Common /prec2/gx(maxptn), gy(maxptn), gz(maxptn), ft(maxptn), px(maxptn), py(maxptn), pz(maxptn), e(maxptn), xmass(maxptn), ityp(maxptn)
  Common /prec4/vx(maxptn), vy(maxptn), vz(maxptn)
  Common /prec5/eta(maxptn), rap(maxptn), tau(maxptn)
  Common /ilist5/ct(maxptn), ot(maxptn), tlarge
  Save
  dt1 = ct(i) - ft(i)
  gx(i) = gx(i) + vx(i)*dt1
  gy(i) = gy(i) + vy(i)*dt1
  gz(i) = gz(i) + vz(i)*dt1
  ft(i) = ct(i)
  If (iconfg<=3) Then
     If (ft(i)<=abs(gz(i))) Then
        eta(i) = 1000000.D0
     Else
        eta(i) = 0.5D0*log((ft(i)+gz(i))/(ft(i)-gz(i)))
     End If
     If (eta(i)<1000000.D0) Then
        tau(i) = ft(i)/cosh(eta(i))
     Else
        tau(i) = 1D-10
     End If
  End If
  Return
End Subroutine newpos

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
Previn is a junior who is currently the Vice President of Special Events of the Davis College Democrats
He is a former resident of Tercero E building

20070202 14:00:26 nbsp Yay! And btw Tercero is awesome! Users/AngelaPourtabib


text:
MODULE InputOutputUtilitiesModule

  USE KindModule, ONLY: &
    DP
  USE MeshModule, ONLY: &
    MeshType, &
    NodeCoordinate

  IMPLICIT NONE
  PRIVATE

  PUBLIC :: NodeCoordinates
  PUBLIC :: Field3D
  PUBLIC :: FromField3D
  PUBLIC :: Field4D
  PUBLIC :: FromField4D
  PUBLIC :: Opacity4D

CONTAINS

  FUNCTION NodeCoordinates( Mesh, nX, nN )

    REAL(DP) :: NodeCoordinates(nX*nN)
    TYPE(MeshType), INTENT(in) :: Mesh
    INTEGER,        INTENT(in) :: nX
    INTEGER,        INTENT(in) :: nN

    INTEGER :: i, j, iNode

    iNode = 0
    DO j = 1, nX
      DO i = 1, nN
        iNode = iNode + 1
        NodeCoordinates(iNode) &
          = NodeCoordinate( Mesh, j, i )
      END DO
    END DO

    RETURN
  END FUNCTION NodeCoordinates

  FUNCTION Field3D( F, nX, nN, nDOF, Tab )

    INTEGER,  INTENT(in) :: &
      nX(3), nN(3), nDOF, Tab(3,nDOF)
    REAL(DP), INTENT(in) :: &
      F(nDOF,nX(1),nX(2),nX(3))
    REAL(DP) :: &
      Field3D(nX(1)*nN(1),nX(2)*nN(2),nX(3)*nN(3))

    INTEGER :: iX1, iX2, iX3
    INTEGER :: iN1, iN2, iN3, iNode

    DO iX3 = 1, nX(3)
    DO iX2 = 1, nX(2)
    DO iX1 = 1, nX(1)

      DO iNode = 1, nDOF

        iN1 = Tab(1,iNode)
        iN2 = Tab(2,iNode)
        iN3 = Tab(3,iNode)

        Field3D &
          ( (iX1-1)*nN(1)+iN1, (iX2-1)*nN(2)+iN2, (iX3-1)*nN(3)+iN3 ) &
          = F(iNode,iX1,iX2,iX3)

      END DO

    END DO
    END DO
    END DO

    RETURN
  END FUNCTION Field3D

  FUNCTION FromField3D( F, nX, nN, nDOF, Tab )

    INTEGER,  INTENT(in) :: &
      nX(3), nN(3), nDOF, Tab(3,nDOF)
    REAL(DP), INTENT(in) :: &
      F(nX(1)*nN(1),nX(2)*nN(2),nX(3)*nN(3))
    REAL(DP) :: &
      FromField3D(nDOF,nX(1),nX(2),nX(3))

    INTEGER :: iX1, iX2, iX3
    INTEGER :: iN1, iN2, iN3, iNode

    DO iX3 = 1, nX(3)
    DO iX2 = 1, nX(2)
    DO iX1 = 1, nX(1)

      DO iNode = 1, nDOF

        iN1 = Tab(1,iNode)
        iN2 = Tab(2,iNode)
        iN3 = Tab(3,iNode)

        FromField3D(iNode,iX1,iX2,iX3) &
          = F( (iX1-1)*nN(1)+iN1, (iX2-1)*nN(2)+iN2, (iX3-1)*nN(3)+iN3 )

      END DO

    END DO
    END DO
    END DO

    RETURN
  END FUNCTION FromField3D

  FUNCTION Field4D( F, nX, nN, nDOF, Tab )

    INTEGER,  INTENT(in) :: &
      nX(4), nN(4), nDOF, Tab(4,nDOF)
    REAL(DP), INTENT(in) :: &
      F(nDOF,nX(1),nX(2),nX(3),nX(4))
    REAL(DP) :: &
      Field4D(nX(1)*nN(1),nX(2)*nN(2),nX(3)*nN(3),nX(4)*nN(4))

    INTEGER :: iX1, iX2, iX3, iX4
    INTEGER :: iN1, iN2, iN3, iN4, iNode

    DO iX4 = 1, nX(4)
    DO iX3 = 1, nX(3)
    DO iX2 = 1, nX(2)
    DO iX1 = 1, nX(1)

      DO iNode = 1, nDOF

        iN1 = Tab(1,iNode)
        iN2 = Tab(2,iNode)
        iN3 = Tab(3,iNode)
        iN4 = Tab(4,iNode)

        Field4D &
          ( (iX1-1)*nN(1)+iN1, (iX2-1)*nN(2)+iN2, &
            (iX3-1)*nN(3)+iN3, (iX4-1)*nN(4)+iN4 ) &
          = F(iNode,iX1,iX2,iX3,iX4)

      END DO

    END DO
    END DO
    END DO
    END DO

    RETURN
  END FUNCTION Field4D

  FUNCTION FromField4D( F, nX, nN, nDOF, Tab )

    INTEGER,  INTENT(in) :: &
      nX(4), nN(4), nDOF, Tab(4,nDOF)
    REAL(DP), INTENT(in) :: &
      F(nX(1)*nN(1),nX(2)*nN(2),nX(3)*nN(3),nX(4)*nN(4))
    REAL(DP) :: &
      FromField4D(nDOF,nX(1),nX(2),nX(3),nX(4))

    INTEGER :: iX1, iX2, iX3, iX4
    INTEGER :: iN1, iN2, iN3, iN4, iNode

    DO iX4 = 1, nX(4)
    DO iX3 = 1, nX(3)
    DO iX2 = 1, nX(2)
    DO iX1 = 1, nX(1)

      DO iNode = 1, nDOF

        iN1 = Tab(1,iNode)
        iN2 = Tab(2,iNode)
        iN3 = Tab(3,iNode)
        iN4 = Tab(4,iNode)

        FromField4D(iNode,iX1,iX2,iX3,iX4) &
          = F( (iX1-1)*nN(1)+iN1, (iX2-1)*nN(2)+iN2, &
               (iX3-1)*nN(3)+iN3, (iX4-1)*nN(4)+iN4 )

      END DO

    END DO
    END DO
    END DO
    END DO

    RETURN
  END FUNCTION FromField4D

  FUNCTION Opacity4D( O, nE, nNE, nDOFE, nX, nNX, nDOFX, TabX )

    INTEGER,  INTENT(in) :: &
      nE,    nNE,    nDOFE
    INTEGER,  INTENT(in) :: &
      nX(3), nNX(3), nDOFX, TabX(3,nDOFX)
    REAL(DP), INTENT(in) :: &
      O(nE*nNE,PRODUCT( nX*nNX ))
    REAL(DP)             :: &
      Opacity4D(nE*nNE,nX(1)*nNX(1),nX(2)*nNX(2),nX(3)*nNX(3))

    INTEGER :: iE, iX1, iX2, iX3, iOS_X
    INTEGER :: iNX1, iNX2, iNX3, iNodeX

    DO iX3 = 1, nX(3)
    DO iX2 = 1, nX(2)
    DO iX1 = 1, nX(1)

      DO iNodeX = 1, nDOFX

        iNX1 = TabX(1,iNodeX)
        iNX2 = TabX(2,iNodeX)
        iNX3 = TabX(3,iNodeX)

        iOS_X = ( (iX3-1)*nX(2)*nX(1) + (iX2-1)*nX(1) + (iX1-1) ) * nDOFX

        Opacity4D &
          (:,(iX1-1)*nNX(1)+iNX1, (iX2-1)*nNX(2)+iNX2, (iX3-1)*nNX(3)+iNX3) &
          = O(:,iOS_X+iNodeX)

      END DO

    END DO
    END DO
    END DO

    RETURN
  END FUNCTION Opacity4D

END MODULE InputOutputUtilitiesModule

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado


text: 20100308 14:10:40 nbsp Hi Umallhairnails, and Welcome to the Wiki!  As a community of individuals, its important that folks here use their real name.  Doing so allows you to establish an identity and interact with the community in a real way.  Please have a look at the welcome page (which I linked just above), as well as the Welcome To The Wiki/Business Owner Business Owner Welcome.  Then, everyone here would be much obliged if youd disable this account (via the settings link in the upper right corner) and start a new one using YourName.  Thanks so much, and again, welcome!  If you have any questions, feel free to ask.  Were here to help. Users/TomGarberson


text:
// -f ${SIM_DIR}/tests/common/vsim.f
+DPI_LIBRARY=${BUILD_DIR}/dpi/uvm_sdv_config_test
+UVM_TESTNAME=simple_dpi_config_test


text:
      subroutine desugar01(x, y, *, *)

      read *, i, n

c     forward reference to 100
      if(i.lt.0) go to 100
c     forward reference to 200
      if(i.lt.0) go to 200
c     forward reference to 300
      if(i.lt.0) go to 300
c     forward reference to 400
      if(i.lt.0) go to 400

c     assign statement
      assign 200 to next

c     labelled logical if with alternate return
 100  if(x.gt.0.) return 2

c     labelled logical if with computed go to
 200  if(x.eq.0.) go to (100, 200, 300), i

c     labelled logical if with assigned go to
 300  if(x.lt.0.) go to next, (100, 200, 300)

c     labelled loop
      continue

 400  do 500 i = 1, n
         print *, i
 500  continue

      end


text:
module class_Circle
  implicit none
  private

  real :: pi = 3.1415926535897931d0 ! Class-wide private constant

  type, public :: Circle
     real :: radius
   contains
     procedure, pass :: area => circle_area
     procedure, pass :: print => circle_print
  end type Circle

contains

  function circle_area(this) result(area)
    type(Circle), intent(in) :: this
    real :: area
    area = pi * this%radius**2
  end function circle_area

  subroutine circle_print(this)
    type(Circle), intent(in) :: this
    real :: area
    area = this%area()  ! Call the type-bound function
    print *, 'Circle: r = ', this%radius, ' area = ', area
  end subroutine circle_print

end module class_Circle

program circle_test
  use class_Circle
  implicit none

  type(Circle) :: c     ! Declare a variable of type Circle.
  c = Circle(1.5)       ! Use the implicit constructor, radius = 1.5.
  call c%print          ! Call the type-bound subroutine
end program circle_test

max_stars_repo_name: jrblevin/scicomp

max_issues_repo_name: jrblevin/scicomp

max_forks_repo_name: jrblevin/scicomp


text:
! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! 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 and documentation.

! 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.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE OUTPUT4_MATRIX_MSGS ( OUNT )

! Writes messages to F06 file regarding matrices that will be output in OUTOUT4 format

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG
      USE IOUNT1, ONLY                :  ERR, F06, LEN_INPUT_FNAME, MOU4, OU4, OU4FIL, OU4_MSG, OU4STAT
      USE OUTPUT4_MATRICES, ONLY      :  ACT_OU4_MYSTRAN_NAMES, ACT_OU4_MYSTRAN_NAMES, ACT_OU4_OUTPUT_NAMES,                       &
                                         ALLOW_OU4_MYSTRAN_NAMES, ALLOW_OU4_OUTPUT_NAMES, NUM_OU4_REQUESTS, NUM_OU4_VALID_NAMES,   &
                                         OU4_FILE_UNITS
      USE PARAMS, ONLY                :  SUPINFO

      USE OUTPUT4_MATRIX_MSGS_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN(ALLOW_OU4_MYSTRAN_NAMES))                                                                                &
                                      :: MYSTRAN_NAMES_I(NUM_OU4_VALID_NAMES)
                                                                     ! OUTPUT4 matrix names written to one OU4 unit

      CHARACTER(LEN(ALLOW_OU4_MYSTRAN_NAMES))                                                                                &
                                      :: OUTPUT_NAMES_I(NUM_OU4_VALID_NAMES)
                                                                     ! OUTPUT4 matrix names written to one OU4 unit

      INTEGER(LONG), INTENT(IN)       :: OUNT(2)             ! File units to write messages to. Input to subr UNFORMATTED_OPEN.
      INTEGER(LONG)                   :: I,J                         ! DO loop indices
      INTEGER(LONG)                   :: NOU4_UNITS                  ! Number of OU4 units for files requested
      INTEGER(LONG)                   :: NOU4_FILES                  ! Number of matrices to write to a specific OPi file

! **********************************************************************************************************************************
      DO I=1,MOU4                                          ! Initially set all to close status of 'DELETE'
         OU4STAT(I) = 'DELETE'                             ! (should have been done in module IOUNT1 but do it here to make sure)
      ENDDO

      NOU4_UNITS = 0                                       ! Reset to "KEEP' close status of units where OUTPUT4 files written
      DO I=1,NUM_OU4_REQUESTS
         DO J=1,MOU4
            IF (OU4_FILE_UNITS(I) == OU4(J)) THEN
               OU4STAT(J) = 'KEEP'
               EXIT
            ENDIF
         ENDDO
      ENDDO

      NOU4_UNITS = 0                                       ! Count the number of OUTPUT4 files that need to be kept
      DO I=1,MOU4
         IF (OU4STAT(I) == 'KEEP') THEN
            NOU4_UNITS = NOU4_UNITS + 1
         ENDIF
      ENDDO

      IF (NOU4_UNITS > 0) THEN
         WRITE(F06,*)
         WRITE(ERR,288) NUM_OU4_REQUESTS, NOU4_UNITS
            WRITE(F06,288) NUM_OU4_REQUESTS, NOU4_UNITS
         DO I=1,MOU4                                       ! Open and then close as "KEEP' the requested OUTPUT4 files
            NOU4_FILES = 0
            IF (OU4STAT(I) == 'KEEP') THEN
               WRITE(F06,289) OU4(I), OU4FIL(I)(1:LEN_INPUT_FNAME+3)
               DO J=1,NUM_OU4_REQUESTS
                  IF (OU4_FILE_UNITS(J) == OU4(I)) THEN
                     NOU4_FILES = NOU4_FILES + 1
                     MYSTRAN_NAMES_I(NOU4_FILES) = ACT_OU4_MYSTRAN_NAMES(J)
                     OUTPUT_NAMES_I(NOU4_FILES)  = ACT_OU4_OUTPUT_NAMES(J)
                  ENDIF
               ENDDO
               DO J=1,NOU4_FILES
                  IF (OUTPUT_NAMES_I(J) == MYSTRAN_NAMES_I(J)) THEN
                     WRITE(F06,290) J, OUTPUT_NAMES_I(J)
                  ELSE
                     WRITE(F06,291) J, OUTPUT_NAMES_I(J), MYSTRAN_NAMES_I(J)
                  ENDIF
               ENDDO
               WRITE(F06,*)
               CALL FILE_OPEN  ( OU4(I), OU4FIL(I), OUNT, 'REPLACE', OU4_MSG(I), 'NEITHER', 'UNFORMATTED', 'WRITE', 'REWIND',      &
                                 'Y', 'N', 'Y')
               CALL FILE_CLOSE ( OU4(I), OU4FIL(I), 'KEEP', 'Y' )
            ENDIF
         ENDDO
         WRITE(F06,*)
      ENDIF

! **********************************************************************************************************************************
  288 FORMAT(' *INFORMATION: THE FOLLOWING ',I3,' MATRICES HAVE BEEN REQUESTED TO BE WRITTEN TO ',I2,' OUTPUT4 FILES IN THE ORDER',&
                           ' LISTED BELOW:',/)

  289 FORMAT(14X,' OUTPUT4 file on unit ',I3,' has been created as: ',A,' and will contain the matrices:')

  290 FORMAT(23X,'(',I2,') ',A)

  291 FORMAT(23X,'(',I2,') ',A,3X,': this is MYSTRAN matrix ',A)

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

      END SUBROUTINE OUTPUT4_MATRIX_MSGS

max_stars_repo_name: luzpaz/MYSTRAN

max_issues_repo_name: luzpaz/MYSTRAN

max_forks_repo_name: luzpaz/MYSTRAN


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

      subroutine somsed ( pmsa   , fl     , ipoint , increm , noseg  ,
     &                    noflux , iexpnt , iknmrk , noq1   , noq2   ,
     &                    noq3   , noq4   )
!>\file
!>       Total of all sedimenting substances

!
!     Description of the module :
!
! Name    T   L I/O   Description                                    Units
! ----    --- -  -    -------------------                            -----
! DMCFy   R*4 1 I  conversion factor for gX->dry matter subst y   [gDM/gX]
! FLXy    R*4 1 I  sedimentation flux substance y                [gX/m3/d]
! TDMSED  R*4 1 O  total dry matter sedimentation flux          [gDM/m2/d]
! TIMSED  R*4 1 O  total inorganic mattter sedimentation flux   [gDM/m2/d]

!     Logical Units : -
!     Modules called : -
!     Name     Type   Library

!     ------   -----  ------------

      IMPLICIT NONE

      REAL     PMSA  ( * ) , FL    (*)
      INTEGER  IPOINT( 40 ) , INCREM(40) , NOSEG , NOFLUX,
     +         IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4

      INTEGER  IP(40)
      REAL     FLX1  , FLX2  , FLX3  , FLX1S2  , FLX2S2  , FLX3S2  ,
     J         FLPOC , FLPOM , FLALGC,
     J         FLALGM, DMCF1 , DMCF2 , DMCF3 , TIMSED, TDMSED, POCSED,
     J         C1    , C2    , C3    , V1    , V2    , V3    , CTOT,
     J         FLPOC1, FLPOC2, FLPOC3, FLPOC4, DMPOC1, DMPOC2, DMPOC3,
     J         DMPOC4, POC1  , POC2  , POC3  , POC4  , CPTOT ,
     J         CP1   , VP1   ,
     J         CP2   , VP2   ,
     J         CP3   , VP3   ,
     J         CP4   , VP4
      INTEGER  IFLUX , ISEG  , IKMRK1, IKMRK2, IQ    , IVAN  , INAAR
      INTEGER  IKMRKN, IKMRKV

      IP = IPOINT
!
      IFLUX = 0
      DO 9000 ISEG = 1 , NOSEG
!!    CALL DHKMRK(1,IKNMRK(ISEG),IKMRK1)
!!    IF (IKMRK1.EQ.1) THEN
      IF (BTEST(IKNMRK(ISEG),0)) THEN
      CALL DHKMRK(2,IKNMRK(ISEG),IKMRK2)
      IF ((IKMRK2.EQ.0).OR.(IKMRK2.EQ.3)) THEN
!

      FLX1    = PMSA(IP(1 ))
      FLX2    = PMSA(IP(2 ))
      FLX3    = PMSA(IP(3 ))
      FLX1S2  = PMSA(IP(4 ))
      FLX2S2  = PMSA(IP(5 ))
      FLX3S2  = PMSA(IP(6 ))
      FLPOC1  = PMSA(IP(7 ))
      FLPOC2  = PMSA(IP(8 ))
      FLPOC3  = PMSA(IP(9 ))
      FLPOC4  = PMSA(IP(10))
      FLALGC  = PMSA(IP(11))
      FLALGM  = PMSA(IP(12))
      DMCF1   = PMSA(IP(13))
      DMCF2   = PMSA(IP(14))
      DMCF3   = PMSA(IP(15))
      DMPOC1  = PMSA(IP(16))
      DMPOC2  = PMSA(IP(17))
      DMPOC3  = PMSA(IP(18))
      DMPOC4  = PMSA(IP(19))

!*******************************************************************************
!**** Calculations connected to the sedimentation
!***********************************************************************

!    Calculate som sedimentation of dry matter
      TIMSED = (FLX1 + FLX1S2) * DMCF1 +
     &         (FLX2 + FLX2S2) * DMCF2 +
     &         (FLX3 + FLX3S2) * DMCF3
      FLPOC  = FLPOC1 + FLPOC2 + FLPOC3 + FLPOC4
      FLPOM  = FLPOC1*DMPOC1 + FLPOC2*DMPOC2 + FLPOC3*DMPOC3
     J                                       + FLPOC4*DMPOC4

      TDMSED = TIMSED + FLPOM + FLALGM

      POCSED = FLPOC + FLALGC

      PMSA (IP(34)) = TDMSED
      PMSA (IP(35)) = TIMSED
      PMSA (IP(36)) = POCSED
      PMSA (IP(37)) = FLPOC
      PMSA (IP(38)) = FLPOM

      ENDIF
      ENDIF
      IFLUX = IFLUX + NOFLUX
      IP    = IP    + INCREM
!
 9000 CONTINUE
!

!.....Exchangeloop over de horizontale richting
      IP = IPOINT
      DO 8000 IQ=1,NOQ1+NOQ2
         PMSA(IP(39)) = 0.0
         PMSA(IP(40)) = 0.0
         IP = IP + INCREM
 8000 CONTINUE

!.....Exchangeloop over de verticale richting
      DO 7000 IQ = NOQ1+NOQ2+1 , NOQ1+NOQ2+NOQ3

         PMSA(IP(39)) = 0.0
         PMSA(IP(40)) = 0.0
         IVAN  = IEXPNT(1,IQ)
         INAAR = IEXPNT(2,IQ)

!        Zoek eerste kenmerk van- en naar-segmenten

         IF ( IVAN .GT. 0 .AND. INAAR .GT. 0 ) THEN
         CALL DHKMRK(1,IKNMRK(IVAN ),IKMRKV)
         CALL DHKMRK(1,IKNMRK(INAAR),IKMRKN)
         IF (IKMRKV.EQ.1.AND.IKMRKN.EQ.1) THEN

!            Water-water uitwisseling

             C1 = PMSA(IPOINT(20)+(IVAN-1)*INCREM(20))
             C2 = PMSA(IPOINT(21)+(IVAN-1)*INCREM(21))
             C3 = PMSA(IPOINT(22)+(IVAN-1)*INCREM(22))
             CP1 = PMSA(IPOINT(23)+(IVAN-1)*INCREM(23))
             CP2 = PMSA(IPOINT(24)+(IVAN-1)*INCREM(24))
             CP3 = PMSA(IPOINT(25)+(IVAN-1)*INCREM(25))
             CP4 = PMSA(IPOINT(26)+(IVAN-1)*INCREM(26))
             V1 = PMSA(IP(27))
             V2 = PMSA(IP(28))
             V3 = PMSA(IP(29))
             VP1 = PMSA(IP(30))
             VP2 = PMSA(IP(31))
             VP3 = PMSA(IP(32))
             VP4 = PMSA(IP(33))
             CTOT = C1 + C2 + C3
             CPTOT = CP1 + CP2 + CP3 + CP4
             IF ( CTOT .GT. 0.0 )
     J       PMSA(IP(39)) = ( C1*V1+C2*V2+C3*V3 ) / CTOT
             IF ( CPTOT .GT. 0.0 )
     J       PMSA(IP(40)) = ( CP1*VP1+CP2*VP2+CP3*VP3+CP4*VP4 ) / CPTOT
         ENDIF
         ENDIF

         IP = IP + INCREM

 7000 CONTINUE

      RETURN
!
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
c
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     *                                                               *
c     *                  copyright (c) 1998 by UCAR                   *
c     *                                                               *
c     *       University Corporation for Atmospheric Research         *
c     *                                                               *
c     *                      all rights reserved                      *
c     *                                                               *
c     *                      SPHEREPACK version 3.2                   *
c     *                                                               *
c     *       A Package of Fortran77 Subroutines and Programs         *
c     *                                                               *
c     *              for Modeling Geophysical Processes               *
c     *                                                               *
c     *                             by                                *
c     *                                                               *
c     *                  John Adams and Paul Swarztrauber             *
c     *                                                               *
c     *                             of                                *
c     *                                                               *
c     *         the National Center for Atmospheric Research          *
c     *                                                               *
c     *                Boulder, Colorado  (80307)  U.S.A.             *
c     *                                                               *
c     *                   which is sponsored by                       *
c     *                                                               *
c     *              the National Science Foundation                  *
c     *                                                               *
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c
c
c     This version of gaqd implements the method presented in:
c     P. N. swarztrauber, Computing the points and weights for
c     Gauss-Legendre quadrature, SIAM J. Sci. Comput.,
c     24(2002) pp. 945-954.
c
c     The w and lwork arrays are dummy and included only to
c     permit a simple pluggable exchange with the
c     old gaqd in previous versions of SPHEREPACK
c
c
c
      subroutine gaqd(nlat,theta,wts,w,lwork,ierror)
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c     gauss points and weights are computed using the fourier-newton
c     described in "on computing the points and weights for
c     gauss-legendre quadrature", paul n. swarztrauber, siam journal
c     on scientific computing that has been accepted for publication.
c     This routine is faster and more accurate than older program
c     with the same name.
c
c     subroutine gaqd computes the nlat gaussian colatitudes and weights
c     in double precision. the colatitudes are in radians and lie in the
c     in the interval (0,pi).
c
c     input parameters
c
c     nlat    the number of gaussian colatitudes in the interval (0,pi)
c             (between the two poles).  nlat must be greater than zero.
c
c     w       unused double precision variable that permits a simple
c             exchange with the old routine with the same name
c             in spherepack.
c
c     lwork   unused variable that permits a simple exchange with the
c             old routine with the same name in spherepack.
c
c     output parameters
c
c     theta   a double precision array with length nlat
c             containing the gaussian colatitudes in
c             increasing radians on the interval (0,pi).
c
c     wts     a double precision array with lenght nlat
c             containing the gaussian weights.
c
c     ierror = 0 no errors
c            = 1 if nlat.le.0
c
c  *****************************************************************
c
      double precision theta(nlat),wts(nlat),w,
     1 x,pi,pis2,dtheta,dthalf,cmax,zprev,zlast,zero,
     2 zhold,pb,dpb,dcor,sum,cz
c
c     check work space length
c
      ierror = 1
      if (nlat.le.0) return
      ierror = 0
c
c     compute weights and points analytically when nlat=1,2
c
      if (nlat.eq.1) then
      theta(1) = dacos(0.0d0)
      wts(1) = 2.0d0
      return
      end if
      if (nlat.eq.2) then
      x = dsqrt(1.0d0/3.0d0)
      theta(1) = dacos(x)
      theta(2) = dacos(-x)
      wts(1) = 1.0d0
      wts(2) = 1.0d0
      return
      end if
      eps = sqrt(dzeps(1.0d0))
      eps = eps*sqrt(eps)
      pis2 = 2.0d0*datan(1.0d0)
      pi = pis2+pis2
      mnlat = mod(nlat,2)
      ns2 = nlat/2
      nhalf = (nlat+1)/2
      idx = ns2+2
c
      call cpdp (nlat,cz,theta(ns2+1),wts(ns2+1))
c
      dtheta = pis2/nhalf
      dthalf = dtheta/2.0d0
      cmax = .2d0*dtheta
c
c     estimate first point next to theta = pi/2
c
      if(mnlat.ne.0) then
      zero = pis2-dtheta
      zprev = pis2
      nix = nhalf-1
      else
      zero = pis2-dthalf
      nix = nhalf
      end if
 9    it = 0
 10   it = it+1
      zlast = zero
c
c     newton iterations
c
      call tpdp (nlat,zero,cz,theta(ns2+1),wts(ns2+1),pb,dpb)
      dcor = pb/dpb
      sgnd = 1.0
      if(dcor .ne. 0.0d0) sgnd = dcor/dabs(dcor)
      dcor = sgnd*min(dabs(dcor),cmax)
      zero = zero-dcor
      if(dabs(zero-zlast).gt.eps*dabs(zero)) go to 10
      theta(nix) = zero
      zhold = zero
c      wts(nix) = (nlat+nlat+1)/(dpb*dpb)
c
c     yakimiw's formula permits using old pb and dpb
c
      wts(nix) = (nlat+nlat+1)/(dpb+pb*dcos(zlast)/dsin(zlast))**2
      nix = nix-1
      if(nix.eq.0) go to 30
      if(nix.eq.nhalf-1)  zero = 3.0*zero-pi
      if(nix.lt.nhalf-1)  zero = zero+zero-zprev
      zprev = zhold
      go to 9
c
c     extend points and weights via symmetries
c
 30   if(mnlat.ne.0) then
      theta(nhalf) = pis2
      call tpdp (nlat,pis2,cz,theta(ns2+1),wts(ns2+1),pb,dpb)
      wts(nhalf) = (nlat+nlat+1)/(dpb*dpb)
      end if
      do i=1,ns2
      wts(nlat-i+1) = wts(i)
      theta(nlat-i+1) = pi-theta(i)
      end do
      sum = 0.0d0
      do i=1,nlat
      sum = sum+wts(i)
      end do
      do i=1,nlat
      wts(i) = 2.0d0*wts(i)/sum
      end do
      return
      end
      subroutine cpdp(n,cz,cp,dcp)
c
c     computes the fourier coefficients of the legendre
c     polynomial p_n^0 and its derivative.
c     n is the degree and n/2 or (n+1)/2
c     coefficients are returned in cp depending on whether
c     n is even or odd. The same number of coefficients
c     are returned in dcp. For n even the constant
c     coefficient is returned in cz.
c
      double precision cp(n/2+1),dcp(n/2+1),
     1 t1,t2,t3,t4,cz
      ncp = (n+1)/2
      t1 = -1.0d0
      t2 = n+1.0d0
      t3 = 0.0d0
      t4 = n+n+1.0d0
      if(mod(n,2).eq.0) then
      cp(ncp) = 1.0d0
      do j = ncp,2,-1
      t1 = t1+2.0d0
      t2 = t2-1.0d0
      t3 = t3+1.0d0
      t4 = t4-2.0d0
      cp(j-1) = (t1*t2)/(t3*t4)*cp(j)
      end do
      t1 = t1+2.0d0
      t2 = t2-1.0d0
      t3 = t3+1.0d0
      t4 = t4-2.0d0
      cz = (t1*t2)/(t3*t4)*cp(1)
      do j=1,ncp
      dcp(j) = (j+j)*cp(j)
      end do
      else
      cp(ncp) = 1.0d0
      do j = ncp-1,1,-1
      t1 = t1+2.0d0
      t2 = t2-1.0d0
      t3 = t3+1.0d0
      t4 = t4-2.0d0
      cp(j) = (t1*t2)/(t3*t4)*cp(j+1)
      end do
      do j=1,ncp
      dcp(j) = (j+j-1)*cp(j)
      end do
      end if
      return
      end
      subroutine tpdp (n,theta,cz,cp,dcp,pb,dpb)
c
c     computes pn(theta) and its derivative dpb(theta) with
c     respect to theta
c
      double precision cp(n/2+1),dcp(n/2+1),cz,
     1  pb,dpb,fn,theta,cdt,sdt,cth,sth,chh
c
      fn = n
      cdt = dcos(theta+theta)
      sdt = dsin(theta+theta)
      if(mod(n,2) .eq.0) then
c
c     n even
c
      kdo = n/2
      pb = .5d0*cz
      dpb = 0.0d0
      if(n .gt. 0) then
      cth = cdt
      sth = sdt
      do 170 k=1,kdo
c      pb = pb+cp(k)*cos(2*k*theta)
      pb = pb+cp(k)*cth
c      dpb = dpb-(k+k)*cp(k)*sin(2*k*theta)
      dpb = dpb-dcp(k)*sth
      chh = cdt*cth-sdt*sth
      sth = sdt*cth+cdt*sth
      cth = chh
  170 continue
      end if
      else
c
c     n odd
c
      kdo = (n+1)/2
      pb = 0.0d0
      dpb = 0.0d0
      cth = dcos(theta)
      sth = dsin(theta)
      do 190 k=1,kdo
c      pb = pb+cp(k)*cos((2*k-1)*theta)
      pb = pb+cp(k)*cth
c      dpb = dpb-(k+k-1)*cp(k)*sin((2*k-1)*theta)
      dpb = dpb-dcp(k)*sth
      chh = cdt*cth-sdt*sth
      sth = sdt*cth+cdt*sth
      cth = chh
  190 continue
      end if
      return
      end
      real function dzeps (x)
      double precision  x
c
c     estimate unit roundoff in quantities of size x.
c
      double precision a,b,c,eps
c
c     this program should function properly on all systems
c     satisfying the following two assumptions,
c        1.  the base used in representing floating point
c            numbers is not a power of three.
c        2.  the quantity  a  in statement 10 is represented to
c            the accuracy used in floating point variables
c            that are stored in memory.
c     the statement number 10 and the go to 10 are intended to
c     force optimizing compilers to generate code satisfying
c     assumption 2.
c     under these assumptions, it should be true that,
c            a  is not exactly equal to four-thirds,
c            b  has a zero for its last bit or digit,
c            c  is not exactly equal to one,
c            eps  measures the separation of 1.0 from
c                 the next larger floating point number.
c     the developers of eispack would appreciate being informed
c     about any systems where these assumptions do not hold.
c
c     this version dated 4/6/83.
c
      a = 4.0d0/3.0d0
   10 b = a - 1.0d0
      c = b + b + b
      eps = abs(c-1.0d0)
      if (eps .eq. 0.0d0) go to 10
      dzeps = eps*dabs(x)
      return
      end


text:
      function RMP2(m_N, tau, taup, term, j_chi, coupling_Array)
      implicit none
      double precision :: RMP2, m_N, j_chi
      double precision :: coupling_Array(14,2)
      integer :: tau,taup
      integer :: term
      integer :: c, q2
      double precision :: c3, c1p
      double precision :: c12, c15, c11p
c
      c = 0
      q2 = 2
c
      c3 = coupling_Array(2,tau)
      c1p = coupling_Array(1,taup)
      c12 = coupling_Array(11,tau)
      c15 = coupling_Array(14,tau)
      c11p = coupling_Array(10,taup)
c
      if (term.eq.c) then
         RMP2 = c3*c1p + (j_chi*(j_chi+1))/3. * c12*c11p
      else if (term.eq.q2) then
         RMP2 = 1/m_N**2 * (j_chi*(j_chi+1))/3. * (-c15*c11p)
      else
         RMP2 = 0.
      end if
c
      end function RMP2


text:
module tscalar_module
    use, intrinsic :: iso_fortran_env, only : int64
    use :: data_storage_module, only : data_storage

    ! All public objects in tensor are public here
    use :: tensor_module

    implicit none
    public

    type, extends(tensor) :: tscalar
    end type tscalar

    interface tscalar
        module procedure constructor_empty
        module procedure constructor
    end interface tscalar
contains
    function constructor_empty() result(this)
        type(tscalar) :: this

        call this%clear()
        this%rank = 0
    end function constructor_empty

    function constructor(storage, datatype) result(this)
        class(data_storage), intent(in) :: storage
        integer, intent(in) :: datatype
        type(tscalar) :: this

        integer(int64), dimension(0) :: empty

        this = tscalar()
        call this%setup(storage, datatype, empty)
    end function constructor
end module tscalar_module

max_stars_repo_name: cheshyre/ntcl-tensor

max_issues_repo_name: cheshyre/ntcl-tensor

max_forks_repo_name: cheshyre/ntcl-tensor