text:
c
c  The small main program below is an example of how to compute field
c   components with T89C.
c    See GEOPACK.DOC for an example of field line tracing.
c
      dimension parmod(10)
  1    print *, '  enter x,y,z,ps,iopt'
      read*, x,y,z,ps,iopt
      call t89c(iopt,parmod,ps,x,y,z,bx,by,bz)
      print *, bx,by,bz
      goto 1
      end

C
      SUBROUTINE T89C(IOPT,PARMOD,PS,X,Y,Z,BX,BY,BZ)
C
C
C   COMPUTES GSM COMPONENTS OF THE MAGNETIC FIELD PRODUCED BY EXTRA-
C   TERRESTRIAL CURRENT SYSTEMS IN THE GEOMAGNETOSPHERE. THE MODEL IS
C   VALID UP TO GEOCENTRIC DISTANCES OF 70 RE AND IS BASED ON THE MER-
C   GED IMP-A,C,D,E,F,G,H,I,J (1966-1974), HEOS-1 AND -2 (1969-1974),
C   AND ISEE-1 AND -2  SPACECRAFT DATA SET.
C
C   THIS IS A MODIFIED VERSION (T89c), WHICH REPLACED THE ORIGINAL ONE
C     IN 1992 AND DIFFERS FROM IT IN THE FOLLOWING:
C
C   (1)  ISEE-1,2 DATA WERE ADDED TO THE ORIGINAL IMP-HEOS DATASET
C   (2)  TWO TERMS WERE ADDED TO THE ORIGINAL TAIL FIELD MODES, ALLOWING
C          A MODULATION OF THE CURRENT BY THE GEODIPOLE TILT ANGLE
C
C
C  REFERENCE FOR THE ORIGINAL MODEL: N.A. TSYGANENKO, A MAGNETOSPHERIC MAGNETIC
C       FIELD MODEL WITH A WARPED TAIL CURRENT SHEET: PLANET.SPACE SCI., V.37,
C         PP.5-20, 1989.
C
C----INPUT PARAMETERS: IOPT - SPECIFIES THE GROUND DISTURBANCE LEVEL:
C
C   IOPT= 1       2        3        4        5        6      7
C                  CORRESPOND TO:
C    KP= 0,0+  1-,1,1+  2-,2,2+  3-,3,3+  4-,4,4+  5-,5,5+  > =6-
C
C    PS - GEODIPOLE TILT ANGLE IN RADIANS
C    X, Y, Z  - GSM COORDINATES OF THE POINT IN EARTH RADII
C
C----OUTPUT PARAMETERS: BX,BY,BZ - GSM COMPONENTS OF THE MODEL MAGNETIC
C                        FIELD IN NANOTESLAS
c
c   THE PARAMETER PARMOD(10) IS A DUMMY ARRAY.  IT IS NOT USED IN THIS
C        SUBROUTINE AND IS PROVIDED JUST FOR MAKING IT COMPATIBLE WITH THE
C           NEW VERSION (4/16/96) OF THE GEOPACK SOFTWARE.
C
C   THIS RELEASE OF T89C IS DATED  FEB 12, 1996;
C--------------------------------------------------------------------------
C
C
C              AUTHOR:     NIKOLAI A. TSYGANENKO
C                          HSTX CORP./NASA GSFC
C
       DIMENSION XI(4),F(3),DER(3,30),PARAM(30,7),A(30),PARMOD(10)
       DOUBLE PRECISION F,DER
        DATA PARAM/-116.53,-10719.,42.375,59.753,-11363.,1.7844,30.268,
     * -0.35372E-01,-0.66832E-01,0.16456E-01,-1.3024,0.16529E-02,
     * 0.20293E-02,20.289,-0.25203E-01,224.91,-9234.8,22.788,7.8813,
     * 1.8362,-0.27228,8.8184,2.8714,14.468,32.177,0.01,0.0,
     * 7.0459,4.0,20.0,-55.553,-13198.,60.647,61.072,-16064.,
     * 2.2534,34.407,-0.38887E-01,-0.94571E-01,0.27154E-01,-1.3901,
     * 0.13460E-02,0.13238E-02,23.005,-0.30565E-01,55.047,-3875.7,
     * 20.178,7.9693,1.4575,0.89471,9.4039,3.5215,14.474,36.555,
     * 0.01,0.0,7.0787,4.0,20.0,-101.34,-13480.,111.35,12.386,-24699.,
     * 2.6459,38.948,-0.34080E-01,-0.12404,0.29702E-01,-1.4052,
     * 0.12103E-02,0.16381E-02,24.49,-0.37705E-01,-298.32,4400.9,18.692,
     * 7.9064,1.3047,2.4541,9.7012,7.1624,14.288,33.822,0.01,0.0,6.7442,
     * 4.0,20.0,-181.69,-12320.,173.79,-96.664,-39051.,3.2633,44.968,
     * -0.46377E-01,-0.16686,0.048298,-1.5473,0.10277E-02,0.31632E-02,
     * 27.341,-0.50655E-01,-514.10,12482.,16.257,8.5834,1.0194,3.6148,
     * 8.6042,5.5057,13.778,32.373,0.01,0.0,7.3195,4.0,20.0,-436.54,
     * -9001.0,323.66,-410.08,-50340.,3.9932,58.524,-0.38519E-01,
     * -0.26822,0.74528E-01,-1.4268,-0.10985E-02,0.96613E-02,27.557,
     * -0.56522E-01,-867.03,20652.,14.101,8.3501,0.72996,3.8149,9.2908,
     *  6.4674,13.729,28.353,0.01,0.0,7.4237,4.0,20.0,-707.77,-4471.9,
     * 432.81,-435.51,-60400.,4.6229,68.178,-0.88245E-01,-0.21002,
     * 0.11846,-2.6711,0.22305E-02,0.10910E-01,27.547,-0.54080E-01,
     * -424.23,1100.2,13.954,7.5337,0.89714,3.7813,8.2945,5.174,14.213,
     * 25.237,0.01,0.0,7.0037,4.0,20.0,-1190.4,2749.9,742.56,-1110.3,
     * -77193.,7.6727,102.05,-0.96015E-01,-0.74507,0.11214,-1.3614,
     * 0.15157E-02,0.22283E-01,23.164,-0.74146E-01,-2219.1,48253.,
     * 12.714,7.6777,0.57138,2.9633,9.3909,9.7263,11.123,21.558,0.01,
     * 0.0,4.4518,4.0,20.0/

        DATA IOP/10/
C
         IF (IOP.NE.IOPT) THEN
C
            ID=1
            IOP=IOPT
            DO 1 I=1,30
   1        A(I)=PARAM(I,IOPT)
C
         ENDIF
C
        XI(1)=X
        XI(2)=Y
        XI(3)=Z
        XI(4)=PS
         CALL T89(ID,A,XI,F,DER)
          IF (ID.EQ.1) ID=2
        BX=F(1)
        BY=F(2)
        BZ=F(3)
        RETURN
       END
C-------------------------------------------------------------------
C
          SUBROUTINE  T89 (ID, A, XI, F, DER)
C
C        ***  N.A. Tsyganenko ***  8-10.12.1991  ***
C
C      Calculates dependent model variables and their deriva-
C  tives for given independent variables and model parame-
C  ters.  Specifies model functions with free parameters which
C  must be determined by means of least squares fits (RMS
C  minimization procedure).
C
C      Description of parameters:
C
C  ID  - number of the data point in a set (initial assignments are performed
c        only for ID=1, saving thus CPU time)
C  A   - input vector containing model parameters;
C  XI  - input vector containing independent variables;
C  F   - output double precision vector containing
C        calculated values of dependent variables;
C  DER   - output double precision vector containing
C        calculated values for derivatives of dependent
C        variables with respect to model parameters;
C
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
C      T89 represents external magnetospheric magnetic field
C  in Cartesian SOLAR MAGNETOSPHERIC coordinates (Tsyganenko N.A.,
C  Planet. Space Sci., 1989, v.37, p.5-20; the "T89 model" with the warped
c  tail current sheet) + A MODIFICATION ADDED IN APRIL 1992 (SEE BELOW)
C
C      Model formulas for the magnetic field components contain in total
c  30 free parameters (17 linear and 13 nonlinear parameters).
C      First 2 independent linear parameters A(1)-A(2) correspond to contribu-
c  tion from the tail current system, then follow A(3) and A(4) which are the
c  amplitudes of symmetric and antisymmetric terms in the contribution from
c  the closure currents; A(5) is the ring current amplitude. Then follow the
c coefficients A(6)-A(15) which define Chapman-Ferraro+Birkeland current field.
c    The coefficients c16-c19  (see Formula 20 in the original paper),
c   due to DivB=0 condition, are expressed through A(6)-A(15) and hence are not
c    independent ones.
c  A(16) AND A(17) CORRESPOND TO THE TERMS WHICH YIELD THE TILT ANGLE DEPEN-
C    DENCE OF THE TAIL CURRENT INTENSITY (ADDED ON APRIL 9, 1992)
C
C      Nonlinear parameters:
C
C    A(18) : DX - Characteristic scale of the Chapman-Ferraro field along the
c        X-axis
C    A(19) : ADR (aRC) - Characteristic radius of the ring current
c    A(20) : D0 - Basic half-thickness of the tail current sheet
C    A(21) : DD (GamRC)- defines rate of thickening of the ring current, as
c             we go from night- to dayside
C    A(22) : Rc - an analog of "hinging distance" entering formula (11)
C    A(23) : G - amplitude of tail current warping in the Y-direction
C    A(24) : aT - Characteristic radius of the tail current
c    A(25) : Dy - characteristic scale distance in the Y direction entering
c                 in W(x,y) in (13)
c    A(26) : Delta - defines the rate of thickening of the tail current sheet
c                 in the Y-direction (in T89 it was fixed at 0.01)
c    A(27) : Q - this parameter was fixed at 0 in the final version of T89;
c              initially it was introduced for making Dy to depend on X
c    A(28) : Sx (Xo) - enters in W(x,y) ; see (13)
c    A(29) : Gam (GamT) - enters in DT in (13) and defines rate of tail sheet
c              thickening on going from night to dayside; in T89 fixed at 4.0
c    A(30) : Dyc - the Dy parameter for closure current system; in T89 fixed
c               at 20.0
c  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
	 IMPLICIT  REAL * 8  (A - H, O - Z)
C
         REAL  A(1), XI(1)
C
	 DIMENSION  F(3), DER(3,30)
C
         INTEGER  ID, I, L
         DATA A02,XLW2,YN,RPI,RT/25.D0,170.D0,30.D0,0.31830989D0,30.D0/
         DATA XD,XLD2/0.D0,40.D0/
C
C   The last four quantities define variation of tail sheet thickness along X
C
      DATA SXC,XLWC2/4.D0,50.D0/
C
C   The two quantities belong to the function WC which confines tail closure
c    current in X- and Y- direction
C
      DATA DXL/20.D0/
C
C
         IF (ID.NE.1)  GOTO  3
	   DO  2  I = 1, 30
	     DO  1  L = 1, 3
  1            DER(L,I) = 0.0D0
  2        CONTINUE
C
       DYC=A(30)
       DYC2=DYC**2
       DX=A(18)
       HA02=0.5D0*A02
       RDX2M=-1.D0/DX**2
       RDX2=-RDX2M
       RDYC2=1.D0/DYC2
       HLWC2M=-0.5D0*XLWC2
       DRDYC2=-2.D0*RDYC2
       DRDYC3=2.D0*RDYC2*DSQRT(RDYC2)
       HXLW2M=-0.5D0*XLW2
       ADR=A(19)
       D0=A(20)
       DD=A(21)
       RC=A(22)
       G=A(23)
       AT=A(24)
       DT=D0
       DEL=A(26)
       P=A(25)
       Q=A(27)
       SX=A(28)
       GAM=A(29)
       HXLD2M=-0.5D0*XLD2
       ADSL=0.D0
       XGHS=0.D0
       H=0.D0
       HS=0.D0
       GAMH=0.D0
       W1=-0.5D0/DX
       DBLDEL=2.D0*DEL
       W2=W1*2.D0
       W4=-1.D0/3.D0
       W3=W4/DX
       W5=-0.5D0
       W6=-3.D0
       AK1=A(1)
       AK2=A(2)
       AK3=A(3)
       AK4=A(4)
       AK5=A(5)
       AK6=A(6)
       AK7=A(7)
       AK8=A(8)
       AK9=A(9)
       AK10=A(10)
       AK11=A(11)
       AK12=A(12)
       AK13=A(13)
       AK14=A(14)
       AK15=A(15)
        AK16=A(16)
        AK17=A(17)
       SXA=0.D0
       SYA=0.D0
       SZA=0.D0
       AK610=AK6*W1+AK10*W5
       AK711=AK7*W2-AK11
       AK812=AK8*W2+AK12*W6
       AK913=AK9*W3+AK13*W4
       RDXL=1.D0/DXL
       HRDXL=0.5D0*RDXL
       A6H=AK6*0.5D0
       A9T=AK9/3.D0
       YNP=RPI/YN*0.5D0
       YND=2.D0*YN
C
  3      CONTINUE
C
           X  = XI(1)
	   Y  = XI(2)
	   Z  = XI(3)
           TILT=XI(4)
               TLT2=TILT**2
	   SPS = DSIN(TILT)
           CPS = DSQRT (1.0D0 - SPS ** 2)
C
       X2=X*X
       Y2=Y*Y
       Z2=Z*Z
       TPS=SPS/CPS
       HTP=TPS*0.5D0
       GSP=G*SPS
       XSM=X*CPS-Z*SPS
       ZSM=X*SPS+Z*CPS
C
C   CALCULATE THE FUNCTION ZS DEFINING THE SHAPE OF THE TAIL CURRENT SHEET
C    AND ITS SPATIAL DERIVATIVES:
C
       XRC=XSM+RC
       XRC16=XRC**2+16.D0
       SXRC=DSQRT(XRC16)
       Y4=Y2*Y2
       Y410=Y4+1.D4
       SY4=SPS/Y410
       GSY4=G*SY4
       ZS1=HTP*(XRC-SXRC)
       DZSX=-ZS1/SXRC
       ZS=ZS1-GSY4*Y4
       D2ZSGY=-SY4/Y410*4.D4*Y2*Y
       DZSY=G*D2ZSGY
C
C   CALCULATE THE COMPONENTS OF THE RING CURRENT CONTRIBUTION:
C
       XSM2=XSM**2
       DSQT=DSQRT(XSM2+A02)
       FA0=0.5D0*(1.D0+XSM/DSQT)
       DDR=D0+DD*FA0
       DFA0=HA02/DSQT**3
       ZR=ZSM-ZS
       TR=DSQRT(ZR**2+DDR**2)
       RTR=1.D0/TR
       RO2=XSM2+Y2
       ADRT=ADR+TR
       ADRT2=ADRT**2
       FK=1.D0/(ADRT2+RO2)
       DSFC=DSQRT(FK)
       FC=FK**2*DSFC
       FACXY=3.0D0*ADRT*FC*RTR
       XZR=XSM*ZR
       YZR=Y*ZR
       DBXDP=FACXY*XZR
       DER(2,5)=FACXY*YZR
       XZYZ=XSM*DZSX+Y*DZSY
       FAQ=ZR*XZYZ-DDR*DD*DFA0*XSM
       DBZDP=FC*(2.D0*ADRT2-RO2)+FACXY*FAQ
       DER(1,5)=DBXDP*CPS+DBZDP*SPS
       DER(3,5)=DBZDP*CPS-DBXDP*SPS
C
C  CALCULATE THE TAIL CURRENT SHEET CONTRIBUTION:
C
       DELY2=DEL*Y2
       D=DT+DELY2
       IF (DABS(GAM).LT.1.D-6) GOTO 8
       XXD=XSM-XD
       RQD=1.D0/(XXD**2+XLD2)
       RQDS=DSQRT(RQD)
       H=0.5D0*(1.D0+XXD*RQDS)
       HS=-HXLD2M*RQD*RQDS
       GAMH=GAM*H
       D=D+GAMH
       XGHS=XSM*GAM*HS
       ADSL=-D*XGHS
   8   D2=D**2
       T=DSQRT(ZR**2+D2)
       XSMX=XSM-SX
       RDSQ2=1.D0/(XSMX**2+XLW2)
       RDSQ=DSQRT(RDSQ2)
       V=0.5D0*(1.D0-XSMX*RDSQ)
       DVX=HXLW2M*RDSQ*RDSQ2
       OM=DSQRT(DSQRT(XSM2+16.D0)-XSM)
       OMS=-OM/(OM*OM+XSM)*0.5D0
       RDY=1.D0/(P+Q*OM)
       OMSV=OMS*V
       RDY2=RDY**2
       FY=1.D0/(1.D0+Y2*RDY2)
       W=V*FY
       YFY1=2.D0*FY*Y2*RDY2
       FYPR=YFY1*RDY
       FYDY=FYPR*FY
       DWX=DVX*FY+FYDY*Q*OMSV
       YDWY=-V*YFY1*FY
       DDY=DBLDEL*Y
       ATT=AT+T
       S1=DSQRT(ATT**2+RO2)
       F5=1.D0/S1
       F7=1.D0/(S1+ATT)
       F1=F5*F7
       F3=F5**3
       F9=ATT*F3
       FS=ZR*XZYZ-D*Y*DDY+ADSL
       XDWX=XSM*DWX+YDWY
       RTT=1.D0/T
       WT=W*RTT
       BRRZ1=WT*F1
       BRRZ2=WT*F3
       DBXC1=BRRZ1*XZR
       DBXC2=BRRZ2*XZR
       DER(2,1)=BRRZ1*YZR
       DER(2,2)=BRRZ2*YZR
          DER(2,16)=DER(2,1)*TLT2
          DER(2,17)=DER(2,2)*TLT2
       WTFS=WT*FS
       DBZC1=W*F5+XDWX*F7+WTFS*F1
       DBZC2=W*F9+XDWX*F1+WTFS*F3
       DER(1,1)=DBXC1*CPS+DBZC1*SPS
       DER(1,2)=DBXC2*CPS+DBZC2*SPS
       DER(3,1)=DBZC1*CPS-DBXC1*SPS
       DER(3,2)=DBZC2*CPS-DBXC2*SPS
          DER(1,16)=DER(1,1)*TLT2
          DER(1,17)=DER(1,2)*TLT2
          DER(3,16)=DER(3,1)*TLT2
          DER(3,17)=DER(3,2)*TLT2
C
C  CALCULATE CONTRIBUTION FROM THE CLOSURE CURRENTS
C
       ZPL=Z+RT
       ZMN=Z-RT
       ROGSM2=X2+Y2
       SPL=DSQRT(ZPL**2+ROGSM2)
       SMN=DSQRT(ZMN**2+ROGSM2)
       XSXC=X-SXC
       RQC2=1.D0/(XSXC**2+XLWC2)
       RQC=DSQRT(RQC2)
       FYC=1.D0/(1.D0+Y2*RDYC2)
       WC=0.5D0*(1.D0-XSXC*RQC)*FYC
       DWCX=HLWC2M*RQC2*RQC*FYC
       DWCY=DRDYC2*WC*FYC*Y
       SZRP=1.D0/(SPL+ZPL)
       SZRM=1.D0/(SMN-ZMN)
       XYWC=X*DWCX+Y*DWCY
       WCSP=WC/SPL
       WCSM=WC/SMN
       FXYP=WCSP*SZRP
       FXYM=WCSM*SZRM
       FXPL=X*FXYP
       FXMN=-X*FXYM
       FYPL=Y*FXYP
       FYMN=-Y*FXYM
       FZPL=WCSP+XYWC*SZRP
       FZMN=WCSM+XYWC*SZRM
       DER(1,3)=FXPL+FXMN
       DER(1,4)=(FXPL-FXMN)*SPS
       DER(2,3)=FYPL+FYMN
       DER(2,4)=(FYPL-FYMN)*SPS
       DER(3,3)=FZPL+FZMN
       DER(3,4)=(FZPL-FZMN)*SPS
C
C   NOW CALCULATE CONTRIBUTION FROM CHAPMAN-FERRARO SOURCES + ALL OTHER
C
           EX=DEXP(X/DX)
           EC=EX*CPS
           ES=EX*SPS
           ECZ=EC*Z
           ESZ=ES*Z
	   ESZY2=ESZ*Y2
	   ESZZ2=ESZ*Z2
	   ECZ2=ECZ*Z
	   ESY=ES*Y
C
	   DER(1,6)=ECZ
	   DER(1,7)=ES
	   DER(1,8)=ESY*Y
	   DER(1,9)=ESZ*Z
	   DER(2,10)=ECZ*Y
	   DER(2,11)=ESY
	   DER(2,12)=ESY*Y2
	   DER(2,13)=ESY*Z2
	   DER(3,14)=EC
	   DER(3,15)=EC*Y2
	   DER(3,6)=ECZ2*W1
	   DER(3,10)=ECZ2*W5
	   DER(3,7)=ESZ*W2
	   DER(3,11)=-ESZ
	   DER(3,8)=ESZY2*W2
	   DER(3,12)=ESZY2*W6
	   DER(3,9)=ESZZ2*W3
	   DER(3,13)=ESZZ2*W4
C
C  FINALLY, CALCULATE NET EXTERNAL MAGNETIC FIELD COMPONENTS,
C    BUT FIRST OF ALL THOSE FOR C.-F. FIELD:
C
      SX1=AK6*DER(1,6)+AK7*DER(1,7)+AK8*DER(1,8)+AK9*DER(1,9)
      SY1=AK10*DER(2,10)+AK11*DER(2,11)+AK12*DER(2,12)+AK13*DER(2,13)
      SZ1=AK14*DER(3,14)+AK15*DER(3,15)+AK610*ECZ2+AK711*ESZ+AK812
     * *ESZY2+AK913*ESZZ2
       BXCL=AK3*DER(1,3)+AK4*DER(1,4)
       BYCL=AK3*DER(2,3)+AK4*DER(2,4)
       BZCL=AK3*DER(3,3)+AK4*DER(3,4)
       BXT=AK1*DER(1,1)+AK2*DER(1,2)+BXCL +AK16*DER(1,16)+AK17*DER(1,17)
       BYT=AK1*DER(2,1)+AK2*DER(2,2)+BYCL +AK16*DER(2,16)+AK17*DER(2,17)
       BZT=AK1*DER(3,1)+AK2*DER(3,2)+BZCL +AK16*DER(3,16)+AK17*DER(3,17)
       F(1)=BXT+AK5*DER(1,5)+SX1+SXA
       F(2)=BYT+AK5*DER(2,5)+SY1+SYA
       F(3)=BZT+AK5*DER(3,5)+SZ1+SZA
C
       RETURN
       END
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ext: f

max_stars_repo_name: scivision/tsyganenko

max_issues_repo_name: scivision/tsyganenko

max_forks_repo_name: johncoxon/tsyganenko


text:
program intrinsics_14
integer, parameter :: &
    k1 = selected_int_kind(5), &
    k2 = selected_int_kind(5_4), &
    k3 = selected_int_kind(5_8), &
    k4 = selected_int_kind(8), &
    k5 = selected_int_kind(10), &
    k6 = selected_real_kind(5), &
    k7 = selected_real_kind(5_4), &
    k8 = selected_real_kind(5_8), &
    k9 = selected_real_kind(8), &
    k10 = selected_char_kind("default")

integer :: i1
character(7) :: s1

s1 = "default"
i1 = 5

print*,selected_int_kind(5), k1 ! 4
print*,selected_int_kind(5_4), k2 ! 4
print*,selected_int_kind(5_8), k3 ! 4
print*,selected_int_kind(8), k4 ! 4
print*,selected_int_kind(10), k5 ! 8
print*,selected_real_kind(5), k6 ! 4
print*,selected_real_kind(5_4), k7 ! 4
print*,selected_real_kind(5_8), k8 ! 4
print*,selected_real_kind(8), k9 ! 8
print*,selected_char_kind("default"), k10 ! 1

print *, selected_int_kind(i1)
print *, selected_real_kind(i1)
print *, selected_char_kind(s1)
end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
There are plenty of talented folks out there ready to Entertainment entertain you for a fee.

Musical Entertainment
  http://www.aharp4u.com Karen Patterson  Local Music Instruction music teacher and professional harpist. Available for weddings, banquets, and special events in the Northern California area.
  KDVS Mobile DJ Unit
  Keith Richman  Music Instruction Guitar teacher also available for weddings or casuals.
  Madrigal Choir Davis High School DHS student choir, performs in full renaissance costuming
  RentaBand  Rent Davis High School DHSs band the month of December
  Users/TigWallis Tig Wallis, Singer/Songwriter
  Yolo Mambo  World / Latin quartet also available for weddings or casuals.

Sound Reinforcement
  http://daviswiki.org/Music_Scene#sound_technicians Sound companies (not really an external link) can help make your event louder, and can often help you find a DJ or band

Movies
 http://filmgeekz.com  Local Movie Website.

Childrens Entertainment
  The Ambitious Card Magical entertainment for Corporate Parties, Trade Shows, Family Parties, Auctions, Fund Raisers and Schools.
  Anitas Childrens Parties
  Once Upon a Princess  Princess and Character Parties for children 12 and under. Reserve for your holiday party today!

Adult Entertainment
  Athenas Home Novelties Goddess
  http://www.dontjustdoitloveit.com/ Passion Parties by J & C  Their parties are free. Learn the latest in our lotions, love potions, and adult novelties. The host gets free products. You get to feel, look at, and sample lots of great top of the line products and have fun.
  Passion Parties by Raechel (Independent Consultant): A wonderful Girls Night In!

 Dancing Entertainment
 Belly Dance by Nyla ~ Professional belly dance entertainment for your wedding, party or corporate event.
  Users/ZahraMorad Zahra Morad: Middle Eastern Dance Artist for family gatherings, festivals, and corporate events.

Magical Entertainment
  The Ambitious Card Magical entertainment for Corporate Parties, Trade Shows, Family Parties, Auctions, Fund Raisers and Schools.
  http://mastermind.thewallisgroup.com/ Rob Wallis, The Master Mind Fun Mindreading Show for Banquets, Fund Raisers, Parties, or Special Events.

Oddball and other options
  Catapult for hire will do parties, see Hurlers (price and time etc.. negotiable)

Could we just keep this on the Entertainment page at the bottom, with the header its got there?

 I figured the topic merrited its own page, but it certainly doesnt need to be.  Users/JanelleAlvstadMattson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
subroutine surface()

use vars
use params
use microphysics, only: micro_field, index_water_vapor
implicit none

real qvs(0:nx,1-YES3D:ny),t_s, q_s, u_h0
real taux0, tauy0, xlmo
real diag_ustar, coef, coef1
real fluxq0_coef, ustar0, u10n
real ws0
integer i,j
real(8) buffer(2), buffer1(2)

real dummy1, dummy2

! LES mode:

call t_startf ('surface')

if(.not.SFC_FLX_FXD) then

  if(sstxy(1,1).le.-100.) then
     print*,'surface: sst is undefined. Quitting...'
     call task_abort()
  end if

  if(OCEAN) then

    if(LES) then

       call oceflx(pres(1),u0(1)+ug, v0(1)+vg,t0(1)-gamaz(1), q0(1),t0(1),z(1),&
                          sstxy(1,1)+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
                          fluxq0_coef, ustar0, u10n)
       if(SFC_TAU_FXD) then
         u_h0 = max(1.,sqrt((u0(1)+ug)**2+(v0(1)+vg)**2))
         taux0 = -(u0(1)+ug)/u_h0*tau0*rhow(1)
         tauy0 = -(v0(1)+vg)/u_h0*tau0*rhow(1)
       else
         tau0=sqrt( taux0**2 +  tauy0**2)/rhow(1)
       end if

       if(doFixedWindSpeedForSurfaceFluxes) then
         call oceflx(pres(1),WindSpeedForFluxes,0.,t0(1)-gamaz(1), q0(1),t0(1),z(1),&
                          sstxy(1,1)+t00, fluxt0, fluxq0, dummy1, dummy2, q_s, &
                          fluxq0_coef, ustar0, u10n)
       end if

       fluxbt(:,:) = fluxt0
       fluxbq(:,:) = fluxq0
       fluxbu(:,:) = taux0/rhow(1)
       fluxbv(:,:) = tauy0/rhow(1)
       ! extra parameters for surface fluxes of water isotopologues
       fluxbq_coef(:,:) = fluxq0_coef
       qsat_surf(:,:) = q_s
       ustar(:,:) = ustar0
       u10arr(:,:) = u10n !bloss(2018-02): For aerosol surface flux in M2005_PA

    end if ! LES

    if(CEM) then

       qvs(0:nx,1-YES3D:ny) = micro_field(0:nx,1-YES3D:ny,1,index_water_vapor)

       do j=1,ny
         do i=1,nx

           if(doFixedWindSpeedForSurfaceFluxes) then
             ! fixed wind speed
             call oceflx(pres(1),WindSpeedForFluxes, 0., &
                  t(i,j,1)-gamaz(1),qv(i,j,1),t(i,j,1),z(1), &
                  sstxy(i,j)+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
                  fluxq0_coef, ustar0, u10n)
           else
             ! normal call
             call oceflx(pres(1),0.5*(u(i+1,j,1)+u(i,j,1))+ug, &
                  0.5*(v(i,j+YES3D,1)+v(i,j,1))+vg, &
                  t(i,j,1)-gamaz(1),qv(i,j,1),t(i,j,1),z(1), &
                  sstxy(i,j)+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
                  fluxq0_coef, ustar0, u10n)
           end if

           fluxbt(i,j) = fluxt0
           fluxbq(i,j) = fluxq0
           ! extra parameters for surface fluxes of water isotopologues
           fluxbq_coef(i,j) = fluxq0_coef
           qsat_surf(i,j) = q_s
           ustar(i,j) = ustar0
           u10arr(i,j) = u10n

           call oceflx(pres(1),u(i,j,1)+ug, &
                       0.25*(v(i-1,j+YES3D,1)+v(i-1,j,1)+v(i,j+YES3D,1)+v(i,j,1))+vg, &
                       0.5*(t(i-1,j,1)+t(i,j,1))-gamaz(1),0.5*(qvs(i-1,j)+qvs(i,j)), &
                       0.5*(t(i-1,j,1)+t(i,j,1)),z(1), &
                       0.5*(sstxy(i-1,j)+sstxy(i,j))+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
                       fluxq0_coef, ustar0, u10n)
           if(SFC_TAU_FXD) then
             u_h0 = max(1.,sqrt((u(i,j,1)+ug)**2+ &
                     (0.25*(v(i-1,j+YES3D,1)+v(i-1,j,1)+v(i,j+YES3D,1)+v(i,j,1))+vg)**2))
             taux0 = -(u(i,j,1)+ug)/u_h0*tau0*rhow(1)
           end if
           fluxbu(i,j) = taux0/rhow(1)

           call oceflx(pres(1),0.25*(u(i+1,j-YES3D,1)+u(i,j-YES3D,1)+u(i+1,j,1)+u(i,j,1))+ug, &
                       v(i,j,1)+vg, &
                       0.5*(t(i,j-YES3D,1)+t(i,j,1))-gamaz(1),0.5*(qvs(i,j-YES3D)+qvs(i,j)), &
                       0.5*(t(i,j-YES3D,1)+t(i,j,1)),z(1), &
                       0.5*(sstxy(i,j-YES3D)+sstxy(i,j))+t00, fluxt0, fluxq0, taux0, tauy0, q_s, &
                       fluxq0_coef, ustar0, u10n)
           if(SFC_TAU_FXD) then
             u_h0 = max(1.,sqrt( &
                       (0.25*(u(i+1,j-YES3D,1)+u(i,j-YES3D,1)+u(i+1,j,1)+u(i,j,1))+ug)**2+ &
                       (v(i,j,1)+vg)**2))
             tauy0 = -(v(i,j,1)+vg)/u_h0*tau0*rhow(1)
           end if
           fluxbv(i,j) = tauy0/rhow(1)

         end do
       end do

    end if ! CEM

  end if ! OCEAN

  if(LAND) then

       if(doFixedWindSpeedForSurfaceFluxes) then
         write(*,*) 'Fixed Wind Speed not yet implemented for land fluxes'
         STOP 'in surface.f90'
       end if

            if(LES) then

               coef = (1000./pres0)**(rgas/cp)
               coef1 = (1000./pres(1))**(rgas/cp)
               t_s = (sstxy(1,1)+t00)*coef
               q_s = soil_wetness*qsatw(sstxy(1,1)+t00,pres(1))
               call landflx(pres(1),(t0(1)-gamaz(1))*coef1, t_s,     &
                      q0(1), q_s, u0(1)+ug, v0(1)+vg, z(1), z0,      &
                      fluxt0, fluxq0, taux0, tauy0, xlmo)
               if(SFC_TAU_FXD) then
                 u_h0 = max(1.,sqrt((u0(1)+ug)**2+(v0(1)+vg)**2))
                 taux0 = -(u0(1)+ug)/u_h0*tau0*rhow(1)
                 tauy0 = -(v0(1)+vg)/u_h0*tau0*rhow(1)
               else
                 tau0=sqrt( taux0**2 +  tauy0**2)/rhow(1)
               end if

               fluxbt(:,:) = fluxt0
               fluxbq(:,:) = fluxq0
               fluxbu(:,:) = taux0/rhow(1)
               fluxbv(:,:) = tauy0/rhow(1)

            end if ! LES

            if(CEM) then

              coef = (1000./pres0)**(rgas/cp)
              coef1 = (1000./pres(1))**(rgas/cp)
              qvs(0:nx,1-YES3D:ny) = micro_field(0:nx,1-YES3D:ny,1,index_water_vapor)

              do j=1,ny
               do i=1,nx

               t_s = (sstxy(i,j)+t00)*coef
               q_s = soil_wetness*qsatw(sstxy(i,j)+t00,pres(1))
               call landflx(pres(1),(t(i,j,1)-gamaz(1))*coef1, t_s,   &
                      qv(i,j,1), q_s, 0.5*(u(i+1,j,1)+u(i,j,1))+ug,     &
                        0.5*(v(i,j+YES3D,1)+v(i,j,1))+vg, z(1), z0,        &
                      fluxt0, fluxq0, taux0, tauy0, xlmo)
               fluxbt(i,j) = fluxt0
               fluxbq(i,j) = fluxq0

               t_s = (0.5*(sstxy(i-1,j)+sstxy(i,j))+t00)*coef
               q_s = soil_wetness*qsatw(0.5*(sstxy(i-1,j)+sstxy(i,j))+t00,pres(1))
               call landflx(pres(1),(0.5*(t(i-1,j,1)+t(i,j,1))-gamaz(1))*coef1, t_s,   &
                      0.5*(qvs(i-1,j)+qvs(i,j)), q_s, u(i,j,1)+ug,     &
                        0.25*(v(i-1,j+YES3D,1)+v(i-1,j,1)+v(i,j+YES3D,1)+v(i,j,1))+vg, &
                       z(1), z0, fluxt0, fluxq0, taux0, tauy0, xlmo)
               if(SFC_TAU_FXD) then
                   u_h0 = max(1.,sqrt((u(i,j,1)+ug)**2+ &
                        (0.25*(v(i-1,j+YES3D,1)+v(i-1,j,1)+v(i,j+YES3D,1)+v(i,j,1))+vg)**2))
                   taux0 = -(u(i,j,1)+ug)/u_h0*tau0*rhow(1)
               end if
               fluxbu(i,j) = taux0/rhow(1)

               t_s = (0.5*(sstxy(i,j-YES3D)+sstxy(i,j))+t00)*coef
               q_s = soil_wetness*qsatw(0.5*(sstxy(i,j-YES3D)+sstxy(i,j))+t00,pres(1))
               call landflx(pres(1),(0.5*(t(i,j-YES3D,1)+t(i,j,1))-gamaz(1))*coef1, t_s,   &
                      0.5*(qvs(i,j-YES3D)+qvs(i,j)), q_s,  &
                      0.25*(u(i+1,j-YES3D,1)+u(i,j-YES3D,1)+u(i+1,j,1)+u(i,j,1))+ug,     &
                      v(i,j,1)+vg, &
                      z(1), z0, fluxt0, fluxq0, taux0, tauy0, xlmo)
               if(SFC_TAU_FXD) then
                  u_h0 = max(1.,sqrt( &
                       (0.25*(u(i+1,j-YES3D,1)+u(i,j-YES3D,1)+u(i+1,j,1)+u(i,j,1))+ug)**2+ &
                       (v(i,j,1)+vg)**2))
                  tauy0 = -(v(i,j,1)+vg)/u_h0*tau0*rhow(1)
               end if
               fluxbv(i,j) = tauy0/rhow(1)

               end do
              end do

            end if ! CEM

  end if ! LAND

end if! .not.SFC_FLX_FXD

if(SFC_FLX_FXD) then

  u_h0 = max(1.,sqrt((u0(1)+ug)**2+(v0(1)+vg)**2))

  if(.not.SFC_TAU_FXD) then
    if(OCEAN) z0 = 0.0001  ! for LAND z0 should be set in namelist (default z0=0.035)

    tau0 = diag_ustar(z(1),  &
                bet(1)*(fluxt0+epsv*(t0(1)-gamaz(1))*fluxq0),u_h0,z0)**2

  end if ! .not.SFC_TAU_FXD

  if(LES) then
    taux0 = -(u0(1)+ug)/u_h0*tau0
    tauy0 = -(v0(1)+vg)/u_h0*tau0
    fluxbu(:,:) = taux0
    fluxbv(:,:) = tauy0
  else
    fluxbu(:,:) = -(u(1:nx,1:ny,1)+ug)/u_h0*tau0
    fluxbv(:,:) = -(v(1:nx,1:ny,1)+vg)/u_h0*tau0
  end if

  fluxbt(:,:) = fluxt0
  fluxbq(:,:) = fluxq0

end if ! SFC_FLX_FXD

!
! Homogenize the surface scalar fluxes if needed for sensitivity studies
!
   if(dosfchomo) then

	fluxt0 = 0.
	fluxq0 = 0.
	do j=1,ny
         do i=1,nx
	   fluxt0 = fluxt0 + fluxbt(i,j)
	   fluxq0 = fluxq0 + fluxbq(i,j)
         end do
        end do
	fluxt0 = fluxt0 / float(nx*ny)
	fluxq0 = fluxq0 / float(nx*ny)
        if(dompi) then
            buffer(1) = fluxt0
            buffer(2) = fluxq0
            call task_sum_real8(buffer,buffer1,2)
	    fluxt0 = buffer1(1) /float(nsubdomains)
	    fluxq0 = buffer1(2) /float(nsubdomains)
        end if ! dompi
	fluxbt(:,:) = fluxt0
	fluxbq(:,:) = fluxq0

   end if

shf_xy(:,:) = shf_xy(:,:) + fluxbt(:,:) * dtfactor
lhf_xy(:,:) = lhf_xy(:,:) + fluxbq(:,:) * dtfactor

call t_stopf ('surface')

end

! ----------------------------------------------------------------------
!
! DISCLAIMER : this code appears to be correct but has not been
!              very thouroughly tested. If you do notice any
!              anomalous behaviour then please contact Andy and/or
!              Bjorn
!
! Function diag_ustar:  returns value of ustar using the below
! similarity functions and a specified buoyancy flux (bflx) given in
! kinematic units
!
! phi_m (zeta > 0) =  (1 + am * zeta)
! phi_m (zeta < 0) =  (1 - bm * zeta)^(-1/4)
!
! where zeta = z/lmo and lmo = (theta_rev/g*vonk) * (ustar^2/tstar)
!
! Ref: Businger, 1973, Turbulent Transfer in the Atmospheric Surface
! Layer, in Workshop on Micormeteorology, pages 67-100.
!
! Code writen March, 1999 by Bjorn Stevens
!
! Code corrected 8th June 1999 (obukhov length was wrong way up,
! so now used as reciprocal of obukhov length)

      real function diag_ustar(z,bflx,wnd,z0)

      implicit none
      real, parameter      :: vonk =  0.4   ! von Karmans constant
      real, parameter      :: g    = 9.81   ! gravitational acceleration
      real, parameter      :: am   =  4.8   !   "          "         "
      real, parameter      :: bm   = 19.3   !   "          "         "
      real, parameter      :: eps  = 1.e-10 ! non-zero, small number

      real, intent (in)    :: z             ! height where u locates
      real, intent (in)    :: bflx          ! surface buoyancy flux (m^2/s^3)
      real, intent (in)    :: wnd           ! wind speed at z
      real, intent (in)    :: z0            ! momentum roughness height

      integer :: iterate
      real    :: lnz, klnz, c1, x, psi1, zeta, rlmo, ustar

      lnz   = log(z/z0)
      klnz  = vonk/lnz
      c1    = 3.14159/2. - 3.*log(2.)

      ustar =  wnd*klnz
      if (bflx /= 0.0) then
        do iterate=1,4
          rlmo   = -bflx * vonk/(ustar**3 + eps)   !reciprocal of
                                                   !obukhov length
          zeta  = z*rlmo
          if (zeta > 0.) then
            ustar =  vonk*wnd  /(lnz + am*zeta)
          else
            x     = sqrt( sqrt( 1.0 - bm*zeta ) )
            psi1  = 2.*log(1.0+x) + log(1.0+x*x) - 2.*atan(x) + c1
            ustar = wnd*vonk/(lnz - psi1)
          end if
        end do
      end if

      diag_ustar = ustar

      return
      end function diag_ustar
! ----------------------------------------------------------------------

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
C     This File is Automatically generated by ALOHA
C     The process calculated in this file is:
C     P(-1,3)*Gamma(-1,-2,1)*Gamma(3,2,-2) - P(-1,3)*Gamma(-1,2,-3)*Gam
C     ma(3,-3,-2)*ProjM(-2,1) - P(-1,3)*Gamma(-1,2,-3)*Gamma(3,-3,-2)*P
C     rojP(-2,1)
C
      SUBROUTINE FFVS19L2_1(P2, V3, S4, COUP, M1, W1, P1, COEFF)
      IMPLICIT NONE
      COMPLEX*16 CI
      PARAMETER (CI=(0D0,1D0))
      COMPLEX*16 V3(*)
      REAL*8 M1
      INCLUDE 'coef_specs.inc'
      COMPLEX*16 COEFF(MAXLWFSIZE,0:VERTEXMAXCOEFS-1,MAXLWFSIZE)
      COMPLEX*16 P2(0:3)
      REAL*8 W1
      COMPLEX*16 P3(0:3)
      COMPLEX*16 P1(0:3)
      COMPLEX*16 COUP
      COMPLEX*16 S4(*)
      P3(0) = V3(1)
      P3(1) = V3(2)
      P3(2) = V3(3)
      P3(3) = V3(4)
      P1(0) = +P2(0)+V3(1)+S4(1)
      P1(1) = +P2(1)+V3(2)+S4(2)
      P1(2) = +P2(2)+V3(3)+S4(3)
      P1(3) = +P2(3)+V3(4)+S4(4)
      COEFF(1,0,1)= COUP*2D0 * M1*S4(5)*(V3(6)*P3(2)-CI*(V3(8)*P3(0))
     $ +CI*(V3(5)*P3(3))-V3(7)*P3(1))
      COEFF(2,0,1)= COUP*2D0 * M1*S4(5)*(P3(0)*(-1D0)*(V3(7)+CI*(V3(6))
     $ )+(P3(1)*(-CI*(V3(8))+CI*(V3(5)))+(P3(2)*(V3(5)-V3(8))+P3(3)
     $ *(V3(7)+CI*(V3(6))))))
      COEFF(3,0,1)= COUP*2D0 * S4(5)*(P1(1)*(P3(0)*(V3(7)+CI*(V3(6)))
     $ +(P3(1)*(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(8)-V3(5))-P3(3)
     $ *(V3(7)+CI*(V3(6))))))+(P1(2)*(P3(0)*(+CI*(V3(7))-V3(6))+(P3(1)
     $ *(V3(5)-V3(8))+(P3(2)*(-CI*(V3(5))+CI*(V3(8)))+P3(3)*(V3(6)-CI
     $ *(V3(7))))))+(P1(0)*(V3(6)*P3(2)-CI*(V3(8)*P3(0))+CI*(V3(5)
     $ *P3(3))-V3(7)*P3(1))+P1(3)*(V3(7)*P3(1)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(6)*P3(2)))))
      COEFF(4,0,1)= COUP*2D0 * S4(5)*(P1(0)*(P3(0)*(-1D0)*(V3(7)+CI
     $ *(V3(6)))+(P3(1)*(-CI*(V3(8))+CI*(V3(5)))+(P3(2)*(V3(5)-V3(8))
     $ +P3(3)*(V3(7)+CI*(V3(6))))))+(P1(3)*(P3(0)*(-1D0)*(V3(7)+CI
     $ *(V3(6)))+(P3(1)*(-CI*(V3(8))+CI*(V3(5)))+(P3(2)*(V3(5)-V3(8))
     $ +P3(3)*(V3(7)+CI*(V3(6))))))+(P1(1)*(V3(7)*P3(1)-CI*(V3(5)*P3(3)
     $ )+CI*(V3(8)*P3(0))-V3(6)*P3(2))+P1(2)*(V3(8)*P3(0)-CI*(V3(7)
     $ *P3(1))+CI*(V3(6)*P3(2))-V3(5)*P3(3)))))
      COEFF(1,1,1)= 0D0
      COEFF(2,1,1)= 0D0
      COEFF(3,1,1)= COUP*2D0 * S4(5)*(V3(6)*P3(2)-CI*(V3(8)*P3(0))+CI
     $ *(V3(5)*P3(3))-V3(7)*P3(1))
      COEFF(4,1,1)= COUP*2D0 * S4(5)*(P3(0)*(-1D0)*(V3(7)+CI*(V3(6)))
     $ +(P3(1)*(-CI*(V3(8))+CI*(V3(5)))+(P3(2)*(V3(5)-V3(8))+P3(3)
     $ *(V3(7)+CI*(V3(6))))))
      COEFF(1,2,1)= 0D0
      COEFF(2,2,1)= 0D0
      COEFF(3,2,1)= COUP*2D0 * S4(5)*(P3(0)*(V3(7)+CI*(V3(6)))+(P3(1)
     $ *(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(8)-V3(5))-P3(3)*(V3(7)+CI
     $ *(V3(6))))))
      COEFF(4,2,1)= COUP*2D0 * S4(5)*(V3(7)*P3(1)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(6)*P3(2))
      COEFF(1,3,1)= 0D0
      COEFF(2,3,1)= 0D0
      COEFF(3,3,1)= COUP*2D0 * S4(5)*(P3(0)*(+CI*(V3(7))-V3(6))+(P3(1)
     $ *(V3(5)-V3(8))+(P3(2)*(-CI*(V3(5))+CI*(V3(8)))+P3(3)*(V3(6)-CI
     $ *(V3(7))))))
      COEFF(4,3,1)= COUP*2D0 * S4(5)*(V3(8)*P3(0)-CI*(V3(7)*P3(1))+CI
     $ *(V3(6)*P3(2))-V3(5)*P3(3))
      COEFF(1,4,1)= 0D0
      COEFF(2,4,1)= 0D0
      COEFF(3,4,1)= COUP*2D0 * S4(5)*(V3(7)*P3(1)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(6)*P3(2))
      COEFF(4,4,1)= COUP*2D0 * S4(5)*(P3(0)*(-1D0)*(V3(7)+CI*(V3(6)))
     $ +(P3(1)*(-CI*(V3(8))+CI*(V3(5)))+(P3(2)*(V3(5)-V3(8))+P3(3)
     $ *(V3(7)+CI*(V3(6))))))
      COEFF(1,0,2)= COUP*2D0 * M1*S4(5)*(P3(0)*(V3(7)-CI*(V3(6)))
     $ +(P3(1)*(+CI*(V3(5)+V3(8)))+(P3(2)*(-1D0)*(V3(5)+V3(8))+P3(3)
     $ *(V3(7)-CI*(V3(6))))))
      COEFF(2,0,2)= COUP*2D0 * M1*S4(5)*(V3(7)*P3(1)-CI*(V3(5)*P3(3))
     $ +CI*(V3(8)*P3(0))-V3(6)*P3(2))
      COEFF(3,0,2)= COUP*2D0 * S4(5)*(P1(0)*(P3(0)*(V3(7)-CI*(V3(6)))
     $ +(P3(1)*(+CI*(V3(5)+V3(8)))+(P3(2)*(-1D0)*(V3(5)+V3(8))+P3(3)
     $ *(V3(7)-CI*(V3(6))))))+(P1(3)*(P3(0)*(+CI*(V3(6))-V3(7))+(P3(1)
     $ *(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(V3(5)+V3(8))+P3(3)*(+CI
     $ *(V3(6))-V3(7)))))+(P1(1)*(V3(6)*P3(2)-CI*(V3(8)*P3(0))+CI
     $ *(V3(5)*P3(3))-V3(7)*P3(1))+P1(2)*(V3(8)*P3(0)-CI*(V3(7)*P3(1))
     $ +CI*(V3(6)*P3(2))-V3(5)*P3(3)))))
      COEFF(4,0,2)= COUP*2D0 * S4(5)*(P1(1)*(P3(0)*(+CI*(V3(6))-V3(7))
     $ +(P3(1)*(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(V3(5)+V3(8))+P3(3)*(
     $ +CI*(V3(6))-V3(7)))))+(P1(2)*(P3(0)*(V3(6)+CI*(V3(7)))+(P3(1)*(
     $ -1D0)*(V3(5)+V3(8))+(P3(2)*(-1D0)*(+CI*(V3(5)+V3(8)))+P3(3)
     $ *(V3(6)+CI*(V3(7))))))+(P1(0)*(V3(7)*P3(1)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(6)*P3(2))+P1(3)*(V3(7)*P3(1)-CI*(V3(5)*P3(3))
     $ +CI*(V3(8)*P3(0))-V3(6)*P3(2)))))
      COEFF(1,1,2)= 0D0
      COEFF(2,1,2)= 0D0
      COEFF(3,1,2)= COUP*2D0 * S4(5)*(P3(0)*(V3(7)-CI*(V3(6)))+(P3(1)
     $ *(+CI*(V3(5)+V3(8)))+(P3(2)*(-1D0)*(V3(5)+V3(8))+P3(3)*(V3(7)
     $ -CI*(V3(6))))))
      COEFF(4,1,2)= COUP*2D0 * S4(5)*(V3(7)*P3(1)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(6)*P3(2))
      COEFF(1,2,2)= 0D0
      COEFF(2,2,2)= 0D0
      COEFF(3,2,2)= COUP*2D0 * S4(5)*(V3(6)*P3(2)-CI*(V3(8)*P3(0))+CI
     $ *(V3(5)*P3(3))-V3(7)*P3(1))
      COEFF(4,2,2)= COUP*2D0 * S4(5)*(P3(0)*(+CI*(V3(6))-V3(7))+(P3(1)
     $ *(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(V3(5)+V3(8))+P3(3)*(+CI
     $ *(V3(6))-V3(7)))))
      COEFF(1,3,2)= 0D0
      COEFF(2,3,2)= 0D0
      COEFF(3,3,2)= COUP*2D0 * S4(5)*(V3(8)*P3(0)-CI*(V3(7)*P3(1))+CI
     $ *(V3(6)*P3(2))-V3(5)*P3(3))
      COEFF(4,3,2)= COUP*2D0 * S4(5)*(P3(0)*(V3(6)+CI*(V3(7)))+(P3(1)
     $ *(-1D0)*(V3(5)+V3(8))+(P3(2)*(-1D0)*(+CI*(V3(5)+V3(8)))+P3(3)
     $ *(V3(6)+CI*(V3(7))))))
      COEFF(1,4,2)= 0D0
      COEFF(2,4,2)= 0D0
      COEFF(3,4,2)= COUP*2D0 * S4(5)*(P3(0)*(+CI*(V3(6))-V3(7))+(P3(1)
     $ *(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(V3(5)+V3(8))+P3(3)*(+CI
     $ *(V3(6))-V3(7)))))
      COEFF(4,4,2)= COUP*2D0 * S4(5)*(V3(7)*P3(1)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(6)*P3(2))
      COEFF(1,0,3)= COUP*2D0 * S4(5)*(P1(1)*(P3(0)*(V3(7)+CI*(V3(6)))
     $ +(P3(1)*(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(-1D0)*(V3(5)+V3(8))
     $ +P3(3)*(V3(7)+CI*(V3(6))))))+(P1(2)*(P3(0)*(+CI*(V3(7))-V3(6))
     $ +(P3(1)*(V3(5)+V3(8))+(P3(2)*(-1D0)*(+CI*(V3(5)+V3(8)))+P3(3)*(
     $ +CI*(V3(7))-V3(6)))))+(P1(0)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(7)*P3(1))+P1(3)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))
     $ +CI*(V3(8)*P3(0))-V3(7)*P3(1)))))
      COEFF(2,0,3)= COUP*2D0 * S4(5)*(P1(0)*(P3(0)*(V3(7)+CI*(V3(6)))
     $ +(P3(1)*(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(-1D0)*(V3(5)+V3(8))
     $ +P3(3)*(V3(7)+CI*(V3(6))))))+(P1(3)*(P3(0)*(-1D0)*(V3(7)+CI
     $ *(V3(6)))+(P3(1)*(+CI*(V3(5)+V3(8)))+(P3(2)*(V3(5)+V3(8))-P3(3)
     $ *(V3(7)+CI*(V3(6))))))+(P1(1)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(7)*P3(1))+P1(2)*(V3(8)*P3(0)-CI*(V3(6)*P3(2))
     $ +CI*(V3(7)*P3(1))-V3(5)*P3(3)))))
      COEFF(3,0,3)= COUP*2D0 * M1*S4(5)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))
     $ +CI*(V3(8)*P3(0))-V3(7)*P3(1))
      COEFF(4,0,3)= COUP*2D0 * M1*S4(5)*(P3(0)*(V3(7)+CI*(V3(6)))
     $ +(P3(1)*(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(-1D0)*(V3(5)+V3(8))
     $ +P3(3)*(V3(7)+CI*(V3(6))))))
      COEFF(1,1,3)= COUP*2D0 * S4(5)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(7)*P3(1))
      COEFF(2,1,3)= COUP*2D0 * S4(5)*(P3(0)*(V3(7)+CI*(V3(6)))+(P3(1)
     $ *(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(-1D0)*(V3(5)+V3(8))+P3(3)
     $ *(V3(7)+CI*(V3(6))))))
      COEFF(3,1,3)= 0D0
      COEFF(4,1,3)= 0D0
      COEFF(1,2,3)= COUP*2D0 * S4(5)*(P3(0)*(V3(7)+CI*(V3(6)))+(P3(1)
     $ *(-1D0)*(+CI*(V3(5)+V3(8)))+(P3(2)*(-1D0)*(V3(5)+V3(8))+P3(3)
     $ *(V3(7)+CI*(V3(6))))))
      COEFF(2,2,3)= COUP*2D0 * S4(5)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(7)*P3(1))
      COEFF(3,2,3)= 0D0
      COEFF(4,2,3)= 0D0
      COEFF(1,3,3)= COUP*2D0 * S4(5)*(P3(0)*(+CI*(V3(7))-V3(6))+(P3(1)
     $ *(V3(5)+V3(8))+(P3(2)*(-1D0)*(+CI*(V3(5)+V3(8)))+P3(3)*(+CI
     $ *(V3(7))-V3(6)))))
      COEFF(2,3,3)= COUP*2D0 * S4(5)*(V3(8)*P3(0)-CI*(V3(6)*P3(2))+CI
     $ *(V3(7)*P3(1))-V3(5)*P3(3))
      COEFF(3,3,3)= 0D0
      COEFF(4,3,3)= 0D0
      COEFF(1,4,3)= COUP*2D0 * S4(5)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(7)*P3(1))
      COEFF(2,4,3)= COUP*2D0 * S4(5)*(P3(0)*(-1D0)*(V3(7)+CI*(V3(6)))
     $ +(P3(1)*(+CI*(V3(5)+V3(8)))+(P3(2)*(V3(5)+V3(8))-P3(3)*(V3(7)
     $ +CI*(V3(6))))))
      COEFF(3,4,3)= 0D0
      COEFF(4,4,3)= 0D0
      COEFF(1,0,4)= COUP*2D0 * S4(5)*(P1(0)*(P3(0)*(+CI*(V3(6))-V3(7))
     $ +(P3(1)*(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(5)-V3(8))+P3(3)
     $ *(V3(7)-CI*(V3(6))))))+(P1(3)*(P3(0)*(+CI*(V3(6))-V3(7))+(P3(1)
     $ *(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(5)-V3(8))+P3(3)*(V3(7)-CI
     $ *(V3(6))))))+(P1(1)*(V3(7)*P3(1)-CI*(V3(8)*P3(0))+CI*(V3(5)
     $ *P3(3))-V3(6)*P3(2))+P1(2)*(V3(8)*P3(0)-CI*(V3(6)*P3(2))+CI
     $ *(V3(7)*P3(1))-V3(5)*P3(3)))))
      COEFF(2,0,4)= COUP*2D0 * S4(5)*(P1(1)*(P3(0)*(+CI*(V3(6))-V3(7))
     $ +(P3(1)*(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(5)-V3(8))+P3(3)
     $ *(V3(7)-CI*(V3(6))))))+(P1(2)*(P3(0)*(V3(6)+CI*(V3(7)))+(P3(1)
     $ *(V3(8)-V3(5))+(P3(2)*(-CI*(V3(5))+CI*(V3(8)))-P3(3)*(V3(6)+CI
     $ *(V3(7))))))+(P1(0)*(V3(7)*P3(1)-CI*(V3(8)*P3(0))+CI*(V3(5)
     $ *P3(3))-V3(6)*P3(2))+P1(3)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(7)*P3(1)))))
      COEFF(3,0,4)= COUP*2D0 * M1*S4(5)*(P3(0)*(+CI*(V3(6))-V3(7))
     $ +(P3(1)*(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(5)-V3(8))+P3(3)
     $ *(V3(7)-CI*(V3(6))))))
      COEFF(4,0,4)= COUP*2D0 * M1*S4(5)*(V3(7)*P3(1)-CI*(V3(8)*P3(0))
     $ +CI*(V3(5)*P3(3))-V3(6)*P3(2))
      COEFF(1,1,4)= COUP*2D0 * S4(5)*(P3(0)*(+CI*(V3(6))-V3(7))+(P3(1)
     $ *(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(5)-V3(8))+P3(3)*(V3(7)-CI
     $ *(V3(6))))))
      COEFF(2,1,4)= COUP*2D0 * S4(5)*(V3(7)*P3(1)-CI*(V3(8)*P3(0))+CI
     $ *(V3(5)*P3(3))-V3(6)*P3(2))
      COEFF(3,1,4)= 0D0
      COEFF(4,1,4)= 0D0
      COEFF(1,2,4)= COUP*2D0 * S4(5)*(V3(7)*P3(1)-CI*(V3(8)*P3(0))+CI
     $ *(V3(5)*P3(3))-V3(6)*P3(2))
      COEFF(2,2,4)= COUP*2D0 * S4(5)*(P3(0)*(+CI*(V3(6))-V3(7))+(P3(1)
     $ *(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(5)-V3(8))+P3(3)*(V3(7)-CI
     $ *(V3(6))))))
      COEFF(3,2,4)= 0D0
      COEFF(4,2,4)= 0D0
      COEFF(1,3,4)= COUP*2D0 * S4(5)*(V3(8)*P3(0)-CI*(V3(6)*P3(2))+CI
     $ *(V3(7)*P3(1))-V3(5)*P3(3))
      COEFF(2,3,4)= COUP*2D0 * S4(5)*(P3(0)*(V3(6)+CI*(V3(7)))+(P3(1)
     $ *(V3(8)-V3(5))+(P3(2)*(-CI*(V3(5))+CI*(V3(8)))-P3(3)*(V3(6)+CI
     $ *(V3(7))))))
      COEFF(3,3,4)= 0D0
      COEFF(4,3,4)= 0D0
      COEFF(1,4,4)= COUP*2D0 * S4(5)*(P3(0)*(+CI*(V3(6))-V3(7))+(P3(1)
     $ *(-CI*(V3(5))+CI*(V3(8)))+(P3(2)*(V3(5)-V3(8))+P3(3)*(V3(7)-CI
     $ *(V3(6))))))
      COEFF(2,4,4)= COUP*2D0 * S4(5)*(V3(6)*P3(2)-CI*(V3(5)*P3(3))+CI
     $ *(V3(8)*P3(0))-V3(7)*P3(1))
      COEFF(3,4,4)= 0D0
      COEFF(4,4,4)= 0D0
      END

max_stars_repo_name: JaySandesara/madminer

max_issues_repo_name: JaySandesara/madminer

max_forks_repo_name: JaySandesara/madminer


text:
      SUBROUTINE DINT1 (A,B,ANSWER,WORK,IOPT)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 1996-03-31 DINT1  Krogh  Removed unused variable in common.
c>> 1995-11-20 DINT1  Krogh  Converted from SFTRAN to Fortran 77.
c>> 1994-10-19 DINT1  Krogh  Changes to use M77CON
c>> 1994-07-07 DINT1  Snyder set up for CHGTYP.
C>> 1993-05-18 DINT1  Krogh -- Changed "END" to "END PROGRAM"
C>> 1991-09-20 DINT1  Krogh converted '(1)' dimensioning to '(*)'.
C>> 1987-11-19 DINT1  Snyder  Initial code.
C
c--D replaces "?": ?INT1, ?INTA, ?intc, ?intec, ?INTF, ?INTOP, ?INTNC
C     ******************************************************************
C
C     THIS SUBROUTINE ATTEMPTS TO CALCULATE THE INTEGRAL OF A FUNCTION
C     OVER THE INTERVAL A TO B TO WITHIN A TOLERANCE LEVEL SPECIFIED
C     BY THE USER.  THE FUNCTION IS PROVIDED BY THE USER VIA A
C     SUBPROGRAM REFERENCED BY CALL DINTF (F,X,IFLAG), OR BY REVERSE
C     COMMUNICATION.  IN THE CASE OF ONE-DIMENSIONAL QUADRATURE, IFLAG
C     IS ALWAYS ZERO.  ALL ABSCISSAE ARE WITHIN THE INTERVAL OF
C     INTEGRATION.
C
C     THE RESULT IS OBTAINED USING QUADRATURE FORMULAE DUE TO
C     T. N. L. PATTERSON, MATHEMATICS OF COMPUTATION, VOLUME 22,
C     PAGES 847-856, 1968.
C
C     *****    WARNING   ***********************************************
C
C     THE RELIABILITY AND EFFICIENCY OF THIS PROGRAM ARE STRONGLY
C     INFLUENCED BY DISCONTINUITIES IN THE FUNCTION OR IT'S
C     DERIVATIVES, INTEGRABLE SINGULARITIES ON THE INTERVAL OF
C     INTEGRATION, AND NON-INTEGRABLE SINGULARITIES NEAR THIS INTERVAL
C     (INCLUDING COMPLEX VALUES).  THE EFFICIENCY AND RELIABILITY
C     OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY
C     SUBDIVIDING THE INTERVAL OF INTEGRATION AT THE DISCONTINUITY,
C     SINGULARITY OR CLOSEST POINT TO A POLE, AND SUMMING THE ANSWERS.
C     A CHANGE OF VARIABLE TO ELIMINATE OR REDUCE THE STRENGTH OF THE
C     SINGULARITY WILL SIGNIFICANTLY IMPROVE PERFORMANCE.
C
C     *****    FORMAL ARGUMENTS    *************************************
C
C A       LOWER LIMIT OF INTERVAL OF INTEGRATION.
      DOUBLE PRECISION A
C B       UPPER LIMIT OF INTERVAL OF INTEGRATION.
      DOUBLE PRECISION B
C ANSWER  IS THE ESTIMATE OF THE INTEGRAL.  WHEN REVERSE COMMUNICATION
C         IS SPECIFIED IT IS USED TO PASS FUNCTION VALUES FROM
C         THE USER PROGRAM TO DINTA.
      DOUBLE PRECISION ANSWER
C WORK    IS A WORKING STORAGE VECTOR USED BY OPTIONS (SEE IOPT).
C         WHEN THE INTEGRATION IS COMPLETE, WORK(1) CONTAINS THE
C         ESTIMATED ABSOLUTE ERROR.  WHEN REVERSE COMMUNICATION
C         IS SPECIFIED, WORK(1) IS USED TO PASS ABSCISSAE FROM
C         DINTA TO THE USER PROGRAM.
      DOUBLE PRECISION WORK(*)
C IOPT    IS A VECTOR OF INTEGERS USED TO RETURN THE STATUS AND TO
C         SPECIFY OPTIONS.
      INTEGER IOPT(*)
C         IOPT(1) IS USED TO RETURN A STATUS INDICATOR TO THE USER.
C         VALUES OF THE STATUS INDICATOR ARE
C         -1 - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE RELATIVE
C              ERROR CRITERIA ARE SATISFIED.
C         -2 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR THE
C              RELATIVE ERROR CRITERIA ARE SATISFIED, BUT THE RELATIVE
C              TO THE OBTAINABLE ERROR CRITERION IS SATISFIED.
C         -3 - NORMAL TERMINATION, BUT NONE OF THE ERROR CRITERIA ARE
C              SATISFIED.
C         +4 - BAD IOPT VALUE.
C         +5 - TOO MANY FUNCTION VALUES NEEDED.
C         +6 - APPARENT NON-INTEGRABLE SINGULARITY NEAR ANSWER.
C         ENTRIES IN IOPT STARTING AT IOPT(2) ARE DESCRIBED BELOW.
C         IOPT(I)    ENTRY IN IOPT(I) MEANS
C           0        NO MORE OPTIONS, BEGIN INTEGRATION.
C           1        IOPT(I+1) CONTAINS THE UNIT NUMBER FOR OUTPUT.
C           2        IOPT(I+1) IS DIAGNOSTIC PRINT LEVEL
C                    0 - NO PRINTING
C                    1 - MINIMUM PRINTING - ERROR MESSAGES (DEFAULT)
C                    2 - PANEL BOUNDARIES AND ANSWERS
C                    3 - ERROR ESTIMATES FOR EACH QUADRATURE FORMULA
C                    4 - DETAILED OUTPUT (DIFFERENCE LINES, ETC).
C           3        WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE TOLERANCE,
C                    WORK(IOPT(I+1)+1) CONTAINS TOLERANCE RELATIVE TO
C                    THE LOCALLY OBTAINABLE PRECISION, AND
C                    WORK(IOPT(I+1)+2) CONTAINS TOLERANCE RELATIVE TO
C                    THE VALUE OF THE INTEGRAL.  THE TOLERANCE RELATIVE
C                    TO THE LOCALLY OBTAINABLE PRECISION IS SPECIFIED AS
C                    THE FRACTION OF LOCALLY OBTAINABLE DIGITS THAT ARE
C                    PERMITTED TO BE WRONG, AND IS INTERNALLY BOUNDED
C                    BETWEEN 0.0 AND 1.0.  IF ANY OF THE ERROR CONTROL
C                    CRITERIA ARE SATISFIED, THE ANSWER IS ACCEPTED.  IF
C                    THIS OPTION IS NOT USED, THE ABSOLUTE AND RELATIVE
C                    TOLERANCES ARE SET TO ZERO, AND THE TOLERANCE
C                    RELATIVE TO THE LOCALLY OBTAINABLE PRECISION IS SET
C                    TO 0.25.
C           4        WORK(IOPT(I+1)) CONTAINS ABSOLUTE ERROR COMMITTED
C                    COMPUTING F.  WORK(IOPT(I+1)) MAY BE CHANGED DURING
C                    THE INTEGRATION.
C           5        WORK(IOPT(I+1)) CONTAINS RELATIVE ERROR EXPECTED TO
C                    BE COMMITTED COMPUTING F.  CHANGES TO THIS VALUE
C                    DURING INTEGRATION WILL NOT BE DETECTED.
C           6        USE REVERSE COMMUNICATION -
C                    CALL DINT1 (A,B,ANSWER,WORK,IOPT)
C                  1 CALL DINTA (ANSWER,WORK,IFLAG)
C                    IF (IFLAG.NE.0) GO TO 2
C                    ANSWER=F(WORK(1))
C                    GO TO 1
C                  2 CONTINUE
C                    IFLAG AT THIS POINT CONTAINS THE SAME VALUE AS
C                    THE STATUS FLAG WOULD (IF SELECTED) FOR THE
C                    FORWARD COMMUNICATION METHOD.
C           7        SET MINIMUM INDEX OF QUADRATURE FORMULA TO USE
C                    BEFORE SUBDIVISION TO IOPT(I+1).
C           8        NO EFFECT.  MAY BE USED TO PASS INFORMATION INTO
C                    DINTF.  INCREMENT I BY 2.
C           9        IOPT(I+1) IS THE MAXIMUM NUMBER OF FUNCTION VALUES
C                    TO USE TO INTEGRATE THE ENTIRE PROBLEM.  IF
C                    IOPT(I+1) .LE. 0, THE NUMBER OF FUNCTION VALUES
C                    IS NOT CONTROLLED.
C          10        IOPT(I+1) IS USED TO RETURN THE NUMBER OF FUNCTION
C                    VALUES USED TO INTEGRATE THE ENTIRE PROBLEM.
C          11        WORK(IABS(IOPT(I+1))) IS THE LOCATION OF A
C                    SINGULARITY OR DISCONTINUITY.  IF THE LOCATION IS
C                    INSIDE THE INTERVAL, THE INTERVAL WILL BE
C                    IMMEDIATELY SUBDIVIDED.  IF IOPT(I+1).GT.0, A
C                    T**2 SUBSTITUTION BASED AT WORK(...) WILL BE USED.
C                    IF IOPT(I+1) .LT. 0 A T**4 SUBSTITUTION WILL BE
C                    USED.  IF IOPT(I+1) .EQ. 0, NO SUBSTITUTION WILL BE
C                    USED.
C          12        WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE ERROR IN THE
C                    LOWER LIMIT OF THE CURRENT DIMENSION.  THE ERROR IN
C                    THE UPPER LIMIT IS IN WORK(IOPT(I+1)+1).
C          13        IS USED ONLY FOR MULTIDIMENSIONAL PROBLEMS.
C
C     ALL OPTIONS ARE SET TO THEIR NOMINAL VALUES BEFORE THE OPTION
C     VECTOR IS PROCESSED.
C
C     *****    EXTERNAL REFERENCES     *********************************
C
C DINTA   TO PERFORM THE INTEGRATION.
C DINTOP  TO SET OPTIONS.
C
C     *****    COMMON STORAGE     **************************************
C
C                    Common blocks are: DINTNC, DINTC, and DINTEC
C     COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR
C     EACH DIMENSION OF A MULTIPLE QUADRATURE.  COMMON /DINTC/
C     CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE
C     QUADRATURE.  THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE
C     ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE
C     IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND
C     SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL.  A PAD OF LOGICAL
C     VARIABLES IS INCLUDED AT THE END OF /DINTC/.  THE DIMENSION OF
C     THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END
C     OF THE COMMON BLOCK ARE ALTERED.
C
C     DECLARATIONS OF COMMON /DINTNC/ VARIABLES.
C
      DOUBLE PRECISION AINIT, BINIT, FNCVAL, S, TP
      DOUBLE PRECISION FER, FER1, RELOBT, TPS, XJ, XJP
      INTEGER     FEA,       FEA1,      INC,       INC2,      IPRINT,
     1 ISTOP(2,2),JPRINT,    KDIM,      KK,        KMAXF,     NDIM,
     2 NFINDX,    NFMAX,     NFMAXM,    RELTOL,    REVERM,    REVERS,
     3 WHEREM
      LOGICAL NEEDH
C
C     DECLARATIONS OF COMMON /DINTC/ VARIABLES.
C
c--D Next line special: S => D, X => Q, D => D, P => D
      DOUBLE PRECISION ACUM, PACUM, RESULT(2)
C     139 $.TYPE.$ VARIABLES
      DOUBLE PRECISION
     1 AACUM,     ABSCIS,    DELMIN,    DELTA,     DIFF,      DISCX(2),
     2 END(2),    ERRINA,    ERRINB,    FAT(2),    FSAVE,
     3 FUNCT(24), F1,        F2,        LOCAL(4),  PAACUM,    PF1,
     4 PF2,       PHISUM,    PHTSUM,    PX,        SPACE(6),
     5 STEP(2),   START(2),  SUM,       T,         TA,        TASAVE,
     6 TB,        TEND,      WORRY(2),  X,         X1,
     7 X2,        XT(17),    FT(17),    PHI(34)
c Note XT, FT, and PHI above are last, because they must be in adjacent
c locations in DINTC.
C     30 $DSTYP$ VARIABLES
      DOUBLE PRECISION
     1 ABSDIF,    COUNT,     EDUE2A,    EDUE2B,    EP,        EPNOIZ,
     2 EPS,       EPSMAX,    EPSMIN,    EPSO,      EPSR,      EPSS,
     3 ERR,       ERRAT(2),  ERRC,      ERRF,      ERRI,      ERRT(2),
     4 ESOLD,     EXTRA,     PEPSMN,    RE,        RELEPS,    REP,
     5 REPROD,    RNDC,      TLEN,      XJUMP
C     29 INTEGER VARIABLES
      INTEGER     DISCF,     DISCHK,    ENDPTS,    I,         INEW,
     1 IOLD,      IP,        IXKDIM,    J,         J1,        J1OLD,
     2 J2,        J2OLD,     K,         KAIMT,     KMAX,      KMIN,
     3 L,         LENDT,     NFEVAL,    NFJUMP,    NSUB,      NSUBSV,
     4 NXKDIM,    PART,      SEARCH,    TALOC,     WHERE,     WHERE2
C     11 TO 18 LOGICALS (7 ARE PADDING).
      LOGICAL     DID1,      FAIL,      FATS(2),   FSAVED,    HAVDIF,
     1 IEND,      INIT,      ROUNDF,    XCDOBT(2), PAD(7)
C
C     THE COMMON BLOCKS.
C
      COMMON /DINTNC/
c        1       2       3     4        5       6       7        8
     W AINIT,  BINIT,  FNCVAL, S,      TP,     FER,    FER1,   RELOBT,
c       9      10       11      12      13       1       2        3
     X TPS,    XJ,     XJP,    FEA,    FEA1,   KDIM,    INC,    INC2,
c     4 (2,2)    8       9     10       11      12       13      14
     Y ISTOP,  JPRINT, IPRINT, KK,     KMAXF,  NDIM,   NFINDX, NFMAX,
c        15     16       17      18      19      20
     Z NFMAXM, RELTOL, REVERM, REVERS, WHEREM, NEEDH
      COMMON /DINTC/
     1 ACUM,   PACUM,  RESULT
      COMMON /DINTC/
c        1     2 (4)     6      7        8       9      10     11 (2)
     1 AACUM,  LOCAL,  ABSCIS, TA,     DELTA,  DELMIN, DIFF,   DISCX,
c     13 (2)     15      16    17 (2)   19     20 (24) 44
     2 END,    ERRINA, ERRINB, FAT,    FSAVE,  FUNCT,  F2,
c       45      46     47       48      49     50      51 (6)
     3 PAACUM, PF1,    PF2,    PHISUM, PHTSUM, PX,     SPACE,
c      57 (2)  59 (2)   61     62        63    64       65
     4 STEP,   START,  SUM,    T,      TASAVE, TB,     TEND,
c      66 (2)  68      69      70      71       72
     5 WORRY,  X1,     X2,     X,      F1,     COUNT,
c      73 (17) 90 (17) 107 (34)
     6 XT,     FT,     PHI
      COMMON /DINTC/
c       141     142    143     144      145     146
     1 ABSDIF, EDUE2A, EDUE2B, EP,     EPNOIZ, EPSMAX,
c       147     148     149    150 (2)  152     153
     2 EPSO,   EPSR,   EPSS,   ERRAT,  ERRC,   ERRF,
c     154 (2)   156     157     158     159    160
     3 ERRT,   ESOLD,  EXTRA,  PEPSMN, RELEPS, REP,
c       161     162     163
     4 RNDC,   TLEN,   XJUMP,
c       164    165      166    167    168       169
     5 ERRI,   ERR,    EPSMIN, EPS,    RE,     REPROD
      COMMON /DINTC/
c       170     171     172
     1 DISCF,  DISCHK, ENDPTS, INEW,   IOLD,   IP,     IXKDIM,
     2 J,      J1,     J1OLD,  J2,     J2OLD,  KMAX,   KMIN,
     3 L,      LENDT,  NFJUMP, NSUBSV, NXKDIM, TALOC,  WHERE2,
c      1       2          3      4       5         6      7       8
     4 I,      K,      KAIMT,  NSUB,   PART,   SEARCH, WHERE, NFEVAL
      COMMON /DINTC/
     1 DID1,   FAIL,   FATS,   FSAVED, HAVDIF, IEND,   INIT,   ROUNDF,
     2 XCDOBT, PAD
      SAVE /DINTNC/, /DINTC/
C
C     THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT.  ALL ARE SET
C     IN DINTOP.  THE MEANING ATTACHED TO THESE VARIABLES CAN BE
C     FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP.
      DOUBLE PRECISION
     1  EMEPS,  EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF,
     2  ESMALL, ENZER,  EDELM1, ENINF
      COMMON /DINTEC/
     1  EMEPS,  EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF,
     2  ESMALL, ENZER,  EDELM1, ENINF
      SAVE /DINTEC/
C
C     *****    PROCEDURES     ******************************************
C
      WHERE=0
      NDIM=1
      KDIM=1
      NFEVAL=0
      CALL DINTOP (IOPT,WORK)
      IF (IOPT(1).EQ.0) THEN
C
C        CALL DINTA TO DO THE INTEGRATION.
C
         AINIT=A
         BINIT=B
         IF (REVERS.EQ.0) CALL DINTA (ANSWER,WORK,IOPT)
      END IF
      RETURN
C
      END


text:
!##############################################################################
!# Tutorial 003e: Random number generator
!##############################################################################

module tutorial003e

  ! Include basic Feat-2 modules
  use fsystem
  use genoutput
  use random

  implicit none
  private

  public :: start_tutorial003e

contains

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

  subroutine start_tutorial003e

    ! Declare some variables
    type(t_random) :: rrng
    integer :: i
    integer(I32), dimension(4) :: Ival
    real(DP), dimension(4) :: Dval

    ! Print a message
    call output_lbrk()
    call output_separator (OU_SEP_STAR)
    call output_line ("This is FEAT-2. Tutorial 003e")
    call output_separator (OU_SEP_MINUS)

    ! =================================
    ! Generate random numbers
    ! =================================

    ! Initialise the pseudo-random generator.
    call rng_init(rrng, 7635_I32)

    ! Generate some integer numbers.
    do i=1,4
      call rng_get(rrng, Ival(i))
    end do

    ! Generate some double numbers.
    do i=1,4
      call rng_get(rrng, Dval(i))
    end do

    ! =================================
    ! Output
    ! =================================

    do i=1,4
#ifdef USE_LARGEINT
      call output_line (" " // trim(sys_siL(int(Ival(i),I64),10)),bnolinebreak=(i .ne. 4))
#else
      call output_line (" " // trim(sys_siL(Ival(i),10)), bnolinebreak=(i .ne. 4))
#endif
    end do

    do i=1,4
      call output_line (" " // sys_sdL(Dval(i),2), bnolinebreak=(i .ne. 4))
    end do

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      subroutine testTangent(matData,F,stre,cc,finite)

      IMPLICIT NONE

      integer          finite

      double precision matData(*), F(3,*), stre(*), cc(6,*)

      logical spectralDecomposition

!      call unit_s_p(stre)
!
!      do i=1, 6
!        c(i) = dlamb2(j,i)
!      enddo
!
!      call mult_4_ss_p(cc,stre,c)
!
!      call set_s_f(stre,stre,lamb2(j))

c      call sq_s(stre,cc,b)
c      !call set_4_f(cc,cc,1.d0)

      return

      end

ext: for

max_stars_repo_name: chennachaos/mpap

max_issues_repo_name: chennachaos/mpap

max_forks_repo_name: chennachaos/mpap


text:
! Bug report from Rice: 19-f77-scale.f
! The identifier 'scale' has no predefined meaning in Fortran 77.
! ROSE always treats it as the F95 intrinsic function named 'scale'.

      scale = 0.0
      end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


!                Parallel Sparse BLAS   GPU plugin
!      (C) Copyright 2013
!
!                         Salvatore Filippone
!                         Alessandro Fanfarillo
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!

subroutine psb_s_ell_cssm(alpha,a,x,beta,y,info,trans)

  use psb_base_mod
  use psb_s_ell_mat_mod, psb_protect_name => psb_s_ell_cssm
  implicit none
  class(psb_s_ell_sparse_mat), intent(in) :: a
  real(psb_spk_), intent(in)          :: alpha, beta, x(:,:)
  real(psb_spk_), intent(inout)       :: y(:,:)
  integer(psb_ipk_), intent(out)       :: info
  character, optional, intent(in)      :: trans

  character :: trans_
  integer(psb_ipk_)   :: i,j,k,m,n, nnz, ir, jc, nxy
  real(psb_spk_), allocatable :: tmp(:,:), acc(:)
  logical            :: tra, ctra
  Integer(Psb_ipk_)  :: err_act
  character(len=20)  :: name='s_ell_cssm'
  logical, parameter :: debug=.false.

  info = psb_success_
  call psb_erractionsave(err_act)
  if (present(trans)) then
    trans_ = trans
  else
    trans_ = 'N'
  end if
  if (.not.a%is_asb()) then
    info = psb_err_invalid_mat_state_
    call psb_errpush(info,name)
    goto 9999
  endif

  if (a%is_dev()) call a%sync()
  tra  = (psb_toupper(trans_) == 'T')
  ctra = (psb_toupper(trans_) == 'C')
  m = a%get_nrows()

  if (.not. (a%is_triangle().and.a%is_sorted())) then
    info = psb_err_invalid_mat_state_
    call psb_errpush(info,name)
    goto 9999
  end if

  if (size(x,1)<m) then
    info = 36
    call psb_errpush(info,name,i_err=(/3*ione,m/))
    goto 9999
  end if

  if (size(y,1)<m) then
    info = 36
    call psb_errpush(info,name,i_err=(/5*ione,m/))
    goto 9999
  end if

  nxy = min(size(x,2),size(y,2))

  if (alpha == szero) then
    if (beta == szero) then
      do i = 1, m
        y(i,:) = szero
      enddo
    else
      do  i = 1, m
        y(i,:) = beta*y(i,:)
      end do
    endif
    return
  end if

  if (beta == szero) then
    allocate(acc(nxy), stat=info)
    if (info /= psb_success_) then
      info = psb_err_alloc_dealloc_
      call psb_errpush(info,name)
      goto 9999
    endif

    call inner_ellsm(tra,ctra,a%is_lower(),a%is_unit(),a%get_nrows(),nxy,&
         & size(a%ja,2,kind=psb_ipk_),a%irn,a%idiag,&
         & a%ja,size(a%ja,1,kind=psb_ipk_),a%val,size(a%val,1,kind=psb_ipk_),&
         & x,size(x,1,kind=psb_ipk_),y,size(y,1,kind=psb_ipk_),acc,info)

    if (info /= 0) then
      info = psb_err_invalid_mat_state_
      call psb_errpush(info,name)
      goto 9999
    endif

    if (alpha == sone) then
      ! do nothing
    else if (alpha == -sone) then
      do  i = 1, m
        y(i,:) = -y(i,:)
      end do
    else
      do  i = 1, m
        y(i,:) = alpha*y(i,:)
      end do
    end if
  else
    allocate(tmp(m,nxy),acc(nxy), stat=info)
    if (info /= psb_success_) then
      info = psb_err_alloc_dealloc_
      call psb_errpush(info,name)
      goto 9999
    endif

    call inner_ellsm(tra,ctra,a%is_lower(),a%is_unit(),a%get_nrows(),nxy,&
         & size(a%ja,2,kind=psb_ipk_),a%irn,a%idiag,&
         & a%ja,size(a%ja,1,kind=psb_ipk_),a%val,size(a%val,1,kind=psb_ipk_),&
         & x,size(x,1,kind=psb_ipk_),tmp,size(tmp,1,kind=psb_ipk_),acc,info)

    if (info == 0) &
         & call psb_geaxpby(m,nxy,alpha,tmp,beta,y,info)

    if (info /= 0) then
      info = psb_err_invalid_mat_state_
      call psb_errpush(info,name)
      goto 9999
    endif

  end if

  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)
  return

contains

  subroutine inner_ellsm(tra,ctra,lower,unit,n,nc,nxy,irn,idiag,ja,ldj,val,ldv,&
       & x,ldx,y,ldy,acc,info)
    implicit none
    logical, intent(in)                  :: tra,ctra,lower,unit
    integer(psb_ipk_), intent(in)        :: n,nc,ldj,ldv,nxy,ldx,ldy
    integer(psb_ipk_), intent(in)        :: irn(*),idiag(*), ja(ldj,*)
    real(psb_spk_), intent(in)          :: val(ldv,*)
    real(psb_spk_), intent(in)          :: x(ldx,nxy)
    real(psb_spk_), intent(out)         :: y(ldy,nxy), acc(nxy)
    integer(psb_ipk_), intent(out)       :: info

    integer(psb_ipk_) :: i,j,k,m, ir, jc

    !
    ! The only error condition here is if
    ! the matrix is non-unit and some idiag value is illegal.
    !
    info = 0

    if (.not.(tra.or.ctra)) then

      if (lower) then

        if (unit) then
          do i=1, n
            acc = szero
            do j=1,irn(i)
              acc = acc + val(i,j)*y(ja(i,j),:)
            end do
            y(i,:) = x(i,:) - acc
          end do
        else if (.not.unit) then
          do i=1, n
            acc = szero
            do j=1,idiag(i)-1
              acc = acc + val(i,j)*y(ja(i,j),:)
            end do
            if (idiag(i) <= 0) then
              info = -1
              return
            endif
            y(i,:) = (x(i,:) - acc)/val(i,idiag(i))
          end do
        end if

      else if (.not.lower) then

        if (unit) then

          do i=n, 1, -1
            acc = szero
            do j=1,irn(i)
              acc = acc + val(i,j)*y(ja(i,j),:)
            end do
            y(i,:) = x(i,:) - acc
          end do

        else if (.not.unit) then

          do i=n, 1, -1
            acc = szero
            do j=idiag(i)+1, irn(i)
              acc = acc + val(i,j)*y(ja(i,j),:)
            end do
            if (idiag(i) <= 0) then
              info = -1
              return
            endif
            y(i,:) = (x(i,:) - acc)/val(i,idiag(i))
          end do

        end if

      end if

    else if (tra) then

      do i=1, n
        y(i,:) = x(i,:)
      end do

      if (lower) then

        if (unit) then

          do i=n, 1, -1
            acc = y(i,:)
            do j=1,irn(i)
              jc    = ja(i,j)
              y(jc,:) = y(jc,:) - val(i,j)*acc
            end do

          end do

        else if (.not.unit) then

          do i=n, 1, -1
            if (idiag(i) <= 0) then
              info = -1
              return
            endif
            y(i,:) = y(i,:)/val(i,idiag(i))
            acc  = y(i,:)
            do j=1,idiag(i)
              jc    = ja(i,j)
              y(jc,:) = y(jc,:) - val(i,j)*acc
            end do
          end do

        end if

      else if (.not.lower) then

        if (unit) then

          do i=1, n
            acc  = y(i,:)
            do j=1, irn(i)
              jc    = ja(i,j)
              y(jc,:) = y(jc,:) - val(i,j)*acc
            end do
          end do

        else if (.not.unit) then

          do i=1, n
            if (idiag(i) <= 0) then
              info = -1
              return
            endif
            y(i,:) = y(i,:)/val(i,idiag(i))
            acc  = y(i,:)
            do j=idiag(i)+1, irn(i)
              jc    = ja(i,j)
              y(jc,:) = y(jc,:) - val(i,j)*acc
            end do
          end do

        end if

      end if

    else if (ctra) then

      do i=1, n
        y(i,:) = x(i,:)
      end do

      if (lower) then

        if (unit) then

          do i=n, 1, -1
            acc = y(i,:)
            do j=1,irn(i)
              jc    = ja(i,j)
              y(jc,:) = y(jc,:) - (val(i,j))*acc
            end do

          end do

        else if (.not.unit) then

          do i=n, 1, -1
            if (idiag(i) <= 0) then
              info = -1
              return
            endif
            y(i,:) = y(i,:)/(val(i,idiag(i)))
            acc  = y(i,:)
            do j=1,idiag(i)
              jc    = ja(i,j)
              y(jc,:) = y(jc,:) - (val(i,j))*acc
            end do
          end do

        end if

      else if (.not.lower) then

        if (unit) then

          do i=1, n
            acc  = y(i,:)
            do j=1, irn(i)
              jc    = ja(i,j)
              y(jc,:) = y(jc,:) - (val(i,j))*acc
            end do
          end do

        else if (.not.unit) then

          do i=1, n
            if (idiag(i) <= 0) then
              info = -1
              return
            endif
            y(i,:) = y(i,:)/(val(i,idiag(i)))
            acc  = y(i,:)
            do j=idiag(i)+1, irn(i)
              jc    = ja(i,j)
              y(jc,:) = y(jc,:) - (val(i,j))*acc
            end do
          end do

        end if

      end if
    end if
  end subroutine inner_ellsm
end subroutine psb_s_ell_cssm


text:
C Utility routines for image manipulation
C
C

       integer function img_index( defn, i, j, k, t, v )
C      -------------------------------------------------
C
C Return index to data specified by defn for the selected point
C
C Given:
C   definition record
       integer    defn(*)
C   vector of specified point
       integer    i, j, k, t, v
C
C-
       include 'ic.inc'
       integer  indx

       indx = (t-defn(14))*(defn(3)*defn(4)*defn(5)*defn(7)) +
     *        (k-defn(12))*(defn(3)*defn(4)*defn(7)) +
     *        (j-defn(10))*(defn(3)*defn(7)) +
     *        (i-defn(8))*(defn(7)) +
     *        v
       if (indx.gt.0 .and. indx.le.defn(2)) then
          img_index = indx
       else
          img_index = -1
       endif
       end
C
C
C
       real*4 function img_interp( defn, d, x, y, z, t, v )
C      ----------------------------------------------------
C
C Return index to data specified by defn for the selected point
C
C Given:
C   definition record
       integer    defn(*)
C   input data
       real*4     d(*)
C   vector of specified point
       real*4     x, y, z, t
       integer    v
C
C-
       include 'ic.inc'

       integer  img_index
       integer  i1, j1, k1, l1
       integer  i2, j2, k2, l2
       real*4   dx, dy, dz, dt, val
       real*4   rx1, ry1, rz1, rt1
       real*4   rx2, ry2, rz2, rt2

C find bounding values
       i1 = x
       j1 = y
       k1 = z
       l1 = t
       i2 = min( i1+1, defn(9) )
       j2 = min( j1+1, defn(11) )
       k2 = min( k1+1, defn(13) )
       l2 = min( l1+1, defn(15) )
       if (i1.lt.defn(8)  .or. i1.gt.defn(9)  .or.
     *     j1.lt.defn(10) .or. j1.gt.defn(11) .or.
     *     k1.lt.defn(12) .or. k1.gt.defn(13) .or.
     *     l1.lt.defn(14) .or. l1.gt.defn(15)     ) then
         image_defn(19) = defn(19)
         img_interp = blank
         return
       endif
       dx = x - i1
       dy = y - j1
       dz = z - k1
       dt = t - l1
       rx2 = d(img_index(defn,i2,j1,k1,l1,v))
       rx1 = d(img_index(defn,i1,j1,k1,l1,v))
       ry1 = rx1 + dx*(rx2-rx1)
       if (defn(1).eq.1) then
         val = ry1
       else
         rx2 = d(img_index(defn,i2,j2,k1,l1,v))
         rx1 = d(img_index(defn,i1,j2,k1,l1,v))
         ry2 = rx1 + dx*(rx2-rx1)
         rz1 = ry1 + dy*(ry2-ry1)
         if (defn(1).eq.2) then
            val = rz1
         else
           rx2 = d(img_index(defn,i2,j1,k2,l1,v))
           rx1 = d(img_index(defn,i1,j1,k2,l1,v))
           ry1 = rx1 + dx*(rx2-rx1)
           rx2 = d(img_index(defn,i2,j2,k2,l1,v))
           rx1 = d(img_index(defn,i1,j2,k2,l1,v))
           ry2 = rx1 + dx*(rx2-rx1)
           rz2 = ry1 + dy*(ry2-ry1)
           rt1 = rz1 + dz*(rz2-rz1)
           if (defn(1).eq.3) then
              val = rt1
           else
              rx2 = d(img_index(defn,i2,j1,k1,l2,v))
              rx1 = d(img_index(defn,i1,j1,k1,l2,v))
              ry1 = rx1 + dx*(rx2-rx1)
              rx2 = d(img_index(defn,i2,j2,k1,l2,v))
              rx1 = d(img_index(defn,i1,j2,k1,l2,v))
              ry2 = rx1 + dx*(rx2-rx1)
              rz1 = ry1 + dy*(ry2-ry1)
              rx2 = d(img_index(defn,i2,j1,k2,l2,v))
              rx1 = d(img_index(defn,i1,j1,k2,l2,v))
              ry1 = rx1 + dx*(rx2-rx1)
              rx2 = d(img_index(defn,i2,j2,k2,l2,v))
              rx1 = d(img_index(defn,i1,j2,k2,l2,v))
              ry2 = rx1 + dx*(rx2-rx1)
              rz2 = ry1 + dy*(ry2-ry1)
              rt2 = rz1 + dz*(rz2-rz1)
              val = rt1 + dt*(rt2-rt1)
           endif
         endif
       endif
       img_interp = val
       end
C
C
       subroutine img_getvec1( defn, d, v, vec )
C      ----------------------------------------
C
C Return a 1D vector of values
C
C Given:
C   definition record
       integer    defn(*)
C   data
       real*4     d(*)
C   vector pointer
       integer    v
C
C Returned:
C   vector of data values
       real*4     vec(*)
C
C-
       include 'ic.inc'
       integer  img_index
       integer  i, n, m

       do m=1,len_defn
         image_defn(m) = defn(m)
       enddo
       if (v.ge.v1 .and. v.le.v2) then
         m = 0
         do i=x1,x2
            n = img_index(defn,i,defn(10),defn(12),defn(14),v)
            m = m + 1
            vec(m) = d(n)
         enddo
       endif
       end
C
C
       subroutine img_getslice( defn, d, defn1, nd, vd, s )
C      ---------------------------------------------------
C
C Return an interpolated vector as specified by defn1
C
C Gievn:
C   definition record of image
       integer    defn(*)
C   input data
       real*4     d(*)
C   definition record for slice
       integer    defn1(*)
C   length of slice
       integer    nd
C Returned:
C   output data slice
       real*4     vd(*)
C
C Updated:
C   status
       integer    s
C
C-

       include 'ic.inc'
       real*4   x, y, z, t
       real*4   img_interp
       integer  n, m

       do m=1,len_defn
         image_defn(m) = defn1(m)
       enddo
       do n=1,nd
          x = float(x1) +  float(n-1)*float(x2-x1)/float(nd-1)
          y = float(y1) +  float(n-1)*float(y2-y1)/float(nd-1)
          z = float(z1) +  float(n-1)*float(z2-z1)/float(nd-1)
          t = float(t1) +  float(n-1)*float(t2-t1)/float(nd-1)
          vd(2*n-1) = sqrt((x-x1)**2+(y-y1)**2+(z-z1)**2+(t-t1)**2)
          vd(2*n) = img_interp(defn, d, x, y, z, t, v1)
       enddo
       end
C
C
       subroutine img_getvalue( defn, d, vecf, v )
C      -------------------------------------------
C
C Return an interpolated value from image data
C
C Given:
C   definition record
       integer    defn(*)
C   data
       real*4     d(*)
C   vector of position in data
       integer    vecf(*)
C
C Returned:
C   data value
       real*4     v
C
C-
       include 'ic.inc'
       real*4   img_interp
       integer  m

       do m=1,len_vecf
         image_vecf(m) = vecf(m)
       enddo
       v = img_interp( defn, d, xp, yp, zp, tp, vp )
       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
C    @(#)mvbude.f	20.3 2/13/96
        subroutine mvbude (ibus, ic)
        integer ibus, ic
c
c       Delete bus record(ibus) and store in array delete(ndelete)
c
c       Input parameters:
c
c          ibus : current bus index
c          ic   : index to chgcrd()
c
        include 'ipfinc/parametr.inc'

        include 'ipfinc/blank.inc'
        include 'ipfinc/bus.inc'
        include 'ipfinc/changr.inc'
        include 'ipfinc/delete.inc'
        include 'ipfinc/prt.inc'

        integer   status, rename_bus
        character record * 120, bus_name * 8

        call bcdbus (ibus, record)  ! obtain bus record
        call stor_del (record, ic)  ! store this record in delete()
c
c       Remove bus from bus hash table by renaming
c
        bus_name = srtlst
        bus_base = 9999.0
        status = rename_bus (ibus, bus_name, bus_base)
c
c       Flag ibus as deleted
c
        do i = 1, 11
           kbsdta(i,ibus) = 0
        enddo
        capcor(1,ibus) = 0.0
        capcor(2,ibus) = -9.0e10

        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
module parff_parsed_character_m
    use parff_parsed_value_m, only: parsed_value_t

    implicit none
    private
    public :: parsed_character_t

    type, extends(parsed_value_t) :: parsed_character_t
        character(len=1) :: value_
    end type
end module

max_stars_repo_name: everythingfunctional/parff

max_issues_repo_name: everythingfunctional/parff

max_forks_repo_name: everythingfunctional/parff


text:
#include "../defines.inc"
subroutine BDsim(R,U,TIME,TTOT,DT,BROWN, &
    inTON,IDUM,PARA,SIMTYPE,COLLISION_TIME,COL_DIST, &
    COL_TYPE)
!
!     External subroutine to perform a Brownian dynamics simulation.
!
!     Andrew Spakowitz
!     Written 11-11-13
      use polydispersity, only : length_of_chain
      use mt19937, only : rnorm
      use params, only : dp
      implicit none
      real(dp) R(3,WLC_P__NT)  ! Bead positions
      real(dp) U(3,WLC_P__NT)  ! Tangent vectors
      real(dp) TIME     ! Time of BD simulation
      real(dp) TTOT     ! Final time of BD simulation

!     Variables in the simulation

      real(dp) RS(3,WLC_P__NT) ! R during the step
      real(dp) US(3,WLC_P__NT) ! R during the step
      real(dp) DT,DT0       ! Time step size
      integer RK                ! Runge-Kutta index
      real(dp) DRDT(3,WLC_P__NT,4) ! Position rate of change
      real(dp) DUDT(3,WLC_P__NT,4) ! Position rate of change
      integer I,J,IB            ! Index Holders
      real(dp) doTU

!     Variables for use in the force calculations

      real(dp) FELAS(3,WLC_P__NT) ! Elastic force
      real(dp) FPONP(3,WLC_P__NT) ! self-int force
      real(dp) TELAS(3,WLC_P__NT) ! Elastic force
      real(dp) TPONP(3,WLC_P__NT) ! self-int force
      real(dp) FORCE    ! External force

!     Variables in the simulation

      real(dp) EB,EPAR,EPERP
      real(dp) GAM,ETA
      real(dp) XIR,XIU
      real(dp) LBOX     ! Box edge length
      real(dp) LHC      ! Length of HC int
      real(dp) VHC      ! HC strength
      real(dp) PARA(10)
      integer SIMTYPE           ! Simulation method (WLC = 1,SSWLC = 2,GC = 3)

!     Variables used for the Brownian forces

      real(dp) FRAND(3,WLC_P__NT) ! Random force
      real(dp) TRAND(3,WLC_P__NT) ! Random force
      real(dp) MAGR,MAGU ! Mag of Brownian forces
      integer BROWN             ! Logic for BD forces
      integer inTON             ! Include polymer interactions
      integer IDUM              ! Seed for the generator

!     Variables for the timestep switch

      integer SWDT

!     Variable to hold time of first collisions between each bead
      real(dp) COLLISION_TIME(WLC_P__NT,WLC_P__NT)
      real(dp) COL_DIST ! l1 dist to trigger collision
      integer COL_TYPE ! algorithm to use for collision detection

!     Load the input parameters

      EB = PARA(1)
      EPAR = PARA(2)
      EPERP = PARA(3)
      GAM = PARA(4)
      ETA = PARA(5)
      XIR = PARA(6)
      XIU = PARA(7)
      LBOX = PARA(8)
      LHC = PARA(9)
      VHC = PARA(10)

      MAGR = sqrt(XIR*2.0/DT)
      MAGU = sqrt(XIU*2.0/DT)
      DT0 = DT
      SWDT = 0

!     Setup the geometric parameters and initialize random forces

      IB = 1
      do 10 I = 1,WLC_P__NP
         do 20 J = 1,length_of_chain(I)
            RS(1,IB) = R(1,IB)
            RS(2,IB) = R(2,IB)
            RS(3,IB) = R(3,IB)
            US(1,IB) = U(1,IB)
            US(2,IB) = U(2,IB)
            US(3,IB) = U(3,IB)
            FELAS(1,IB) = 0.
            FELAS(2,IB) = 0.
            FELAS(3,IB) = 0.
            FRAND(1,IB) = 0.
            FRAND(2,IB) = 0.
            FRAND(3,IB) = 0.
            FPONP(1,IB) = 0.
            FPONP(2,IB) = 0.
            FPONP(3,IB) = 0.
            TELAS(1,IB) = 0.
            TELAS(2,IB) = 0.
            TELAS(3,IB) = 0.
            TRAND(1,IB) = 0.
            TRAND(2,IB) = 0.
            TRAND(3,IB) = 0.
            TPONP(1,IB) = 0.
            TPONP(2,IB) = 0.
            TPONP(3,IB) = 0.
            IB = IB + 1
 20      continue
 10   continue

!     Begin the time integration

      do while (TIME < TTOT)

         call CHECK_COLLISIONS(R, WLC_P__NT, COLLISION_TIME, COL_DIST, TIME, COL_TYPE)

!     Calculate the random forces and torques for use in this
!     timestep calculation if BROWN = 1

         RK = 1
         do while (RK <= 4)

 130        continue

            if (BROWN == 1.AND.RK == 1) then
               IB = 1
               do 30 I = 1,WLC_P__NP
                  do 40 J = 1,length_of_chain(I)
                     FRAND(1,IB) = MAGR*rnorm()
                     FRAND(2,IB) = MAGR*rnorm()
                     FRAND(3,IB) = MAGR*rnorm()
                     TRAND(1,IB) = MAGU*rnorm()
                     TRAND(2,IB) = MAGU*rnorm()
                     TRAND(3,IB) = MAGU*rnorm()
                     IB = IB + 1
 40               continue
 30            continue
            endif

!     Calculate the four Runge-Kutta derivatives

!     Calculate the elastic forces (same as free chain)

            call force_elas(FELAS,TELAS,R,U,EB,EPAR,EPERP,GAM,ETA,SIMTYPE)

!     Calculate the self forces

            if (inTON == 1) then
               call force_ponp(FPONP,R,LHC,VHC,LBOX,GAM,DT,XIR,SWDT)

!     If timestep is switch, reset coords and redo step

               if (SWDT == 1) then
                  print*, "Time-step switch", DT,RK,TIME
                  SWDT = 0
                  MAGR = sqrt(XIR*2.0/DT)
                  MAGU = sqrt(XIU*2.0/DT)
                  IB = 1
                  do 60 I = 1,WLC_P__NP
                     do 65 J = 1,length_of_chain(I)
                        R(1,IB) = RS(1,IB)
                        R(2,IB) = RS(2,IB)
                        R(3,IB) = RS(3,IB)
                        U(1,IB) = US(1,IB)
                        U(2,IB) = US(2,IB)
                        U(3,IB) = US(3,IB)
                        IB = IB + 1
 65                  continue
 60               continue
                  RK = 1
                  goto 130
               endif
            endif

!     Calculate the change in the position vector

            IB = 1
            do 70 I = 1,WLC_P__NP
               do 80 J = 1,length_of_chain(I)
                  DRDT(1,IB,RK) = (FELAS(1,IB) + FPONP(1,IB))/XIR
                  DRDT(2,IB,RK) = (FELAS(2,IB) + FPONP(2,IB))/XIR
                  DRDT(3,IB,RK) = (FELAS(3,IB) + FPONP(3,IB))/XIR
                  DUDT(1,IB,RK) = (TELAS(1,IB) + TPONP(1,IB))/XIU
                  DUDT(2,IB,RK) = (TELAS(2,IB) + TPONP(2,IB))/XIU
                  DUDT(3,IB,RK) = (TELAS(3,IB) + TPONP(3,IB))/XIU

                  if (BROWN == 0) then
                     doTU = DUDT(1,IB,RK)*U(1,IB) + DUDT(2,IB,RK)*U(2,IB) + DUDT(3,IB,RK)*U(3,IB)
                     DUDT(1,IB,RK) = DUDT(1,IB,RK)-doTU*U(1,IB)
                     DUDT(2,IB,RK) = DUDT(2,IB,RK)-doTU*U(2,IB)
                     DUDT(3,IB,RK) = DUDT(3,IB,RK)-doTU*U(3,IB)
                  endif
                  IB = IB + 1
 80            continue
 70         continue

            if (BROWN == 1) then
               IB = 1
               do 90 I = 1,WLC_P__NP
                  do 100 J = 1,length_of_chain(I)
                     DRDT(1,IB,RK) = DRDT(1,IB,RK) + FRAND(1,IB)/XIR
                     DRDT(2,IB,RK) = DRDT(2,IB,RK) + FRAND(2,IB)/XIR
                     DRDT(3,IB,RK) = DRDT(3,IB,RK) + FRAND(3,IB)/XIR
                     DUDT(1,IB,RK) = DUDT(1,IB,RK) + TRAND(1,IB)/XIU
                     DUDT(2,IB,RK) = DUDT(2,IB,RK) + TRAND(2,IB)/XIU
                     DUDT(3,IB,RK) = DUDT(3,IB,RK) + TRAND(3,IB)/XIU

                     doTU = DUDT(1,IB,RK)*U(1,IB) + DUDT(2,IB,RK)*U(2,IB) + DUDT(3,IB,RK)*U(3,IB)
                     DUDT(1,IB,RK) = DUDT(1,IB,RK)-doTU*U(1,IB)
                     DUDT(2,IB,RK) = DUDT(2,IB,RK)-doTU*U(2,IB)
                     DUDT(3,IB,RK) = DUDT(3,IB,RK)-doTU*U(3,IB)

                     IB = IB + 1
 100               continue
 90            continue
            endif

!     If SIMTYPE = 1 (WLC), calculate the constraint forces

            if (SIMTYPE == 1) then
               call concalc(R,DRDT,XIR,GAM,DT,RK,BROWN)
            endif

!     Step forward using the RK algorithm

            call RKstep(RS,R,US,U,DRDT,DUDT,RK,DT)

            RK = RK + 1

         ENDdo

         TIME = TIME + DT

!     Swap old variables for new ones

         DT = DT0
         MAGR = sqrt(XIR*2.0/DT)
         MAGU = sqrt(XIU*2.0/DT)

         IB = 1
         do 110 I = 1,WLC_P__NP
            do 120 J = 1,length_of_chain(I)
               RS(1,IB) = R(1,IB)
               RS(2,IB) = R(2,IB)
               RS(3,IB) = R(3,IB)
               US(1,IB) = U(1,IB)
               US(2,IB) = U(2,IB)
               US(3,IB) = U(3,IB)
               IB = IB + 1
 120        continue
 110     continue

      ENDdo

      RETURN
END

max_stars_repo_name: SpakowitzLab/BasicWLC

max_issues_repo_name: riscalab/wlcsim

max_forks_repo_name: riscalab/wlcsim


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 	RADEX
!	Original code by Van der Tak et al. 2007
!	All publications using this code should reference their release paper:
!	A&A 468, 627 (2007)
!
!   This amended version by Jon Holdship has been updated to Modern Fortran
!   with a view to removing common blocks and making an F2PY module for python.
!
!	This code has been tested against the original but the authors do not
!	guarantee there are no errors. We advise users to check results they intend
!	to publish with another version of RADEX or another code.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM Radex
USE IO
USE Solver
USE Background
IMPLICIT NONE
    !Main program: controls program flow and drives subroutines

    INTEGER :: niter,nlines   ! iteration counter
    INTEGER :: imore=1   ! are we running again?
    LOGICAL :: conv    ! are we converged?

    ! Begin executable statements
    write(*,*)
    write(*,*)'   Welcome to Radex, Moden Fortran Edition'
    write(*,*)

    method=1

    DO WHILE (imore.eq.1)
        !     Get input parameters
        IF (DEBUG) write(*,*) 'calling getinputs'

        CALL getinputs

        !     Read data file
        IF (DEBUG) write(*,*) 'calling readdata'
        CALL ReadData

        write(*,*) "Beginning Calculation"

        !     Calculate background radiation field
        IF (DEBUG) write(*,*) 'calling backrad'
        CALL backrad

        niter = 0
        conv  = .false.

        !Set up rate matrix, splitting it in radiative and collisional parts
        !Invert rate matrix to get `thin' starting condition
        IF (DEBUG) write(*,*) 'calling matrix'
        CALL matrix(niter,conv)

        !Start iterating
        DO niter=1,maxiter
            !Invert rate matrix using escape probability for line photons
            CALL matrix(niter,conv)
            IF (conv) THEN
                write(*,*) 'Finished in ',niter,' iterations.'
                EXIT
            END IF
        END DO

        IF (.NOT. conv) write(*,*) '   Warning: Calculation did not converge in ',maxiter&
            &,' iterations.'

        !Prepare output by calculating final quantities
        IF (DEBUG) write(*,*) 'calculating output summary variables'
        CALL CalcOutputArrays(nlines)

        !Write output
        IF (DEBUG) write(*,*) 'calling output'
        CALL output(niter)

        !     See if user wants more, else call it a day
        51   format(A,$)
        write(*,51) '  Another calculation [0/1] ? '
        read(*,*) imore
        write(13,52) imore
        52   format(i2)
    END DO
    write(*,*) '   Have a nice day.'
    !     Done! Now close log file ...
    close(13)
    !     ...and output file.
    close(8)
END PROGRAM Radex

max_stars_repo_name: uclchem/SpectralRadex

max_issues_repo_name: jonholdship/SpectralRadex

max_forks_repo_name: jonholdship/SpectralRadex


text:
C     The 12/15/2000 version was changed on 2/22/2002.
C     The old version is saved as OPENBIG_V01
C     The changes are:
C     1.  No more Lahey code
C     2.  SHARE='DENYNONE' is added, which is the default for MS Fortran but not for new Compaq
C     3. 	IOSTAT=IERROR is added to give more info when error occurs
C     4.  Routine will retry if user requests.

      SUBROUTINE OPENBIG(IUNIT,FILENAME,S)

      CHARACTER*(*) FILENAME,S

    5 CONTINUE

      OPEN(IUNIT,FILE=FILENAME,STATUS=S,BLOCKSIZE=32000,ACCESS='SEQUENTIAL',FORM='BINARY',
     &     SHARE='DENYNONE',IOSTAT=IERROR,ERR=10)
      RETURN
C
   10 CONTINUE
      WRITE(*,1001) FILENAME,IUNIT,S,IERROR
 1001 FORMAT('ERROR IN OPENBIG ',A50,1X,I3,1X,A10,I7)
      WRITE(*,*) 'ENTER 0 TO TRY AGAIN, ENTER 1 TO STOP'
	READ(*,*) IGO
	IF(IGO.EQ.1) STOP 'PROGRAM STOPPED BY USER'
	GOTO 5
      END


text:
MODULE wlDependentVariablesModule

  USE wlKindModule, ONLY: dp

  implicit none
  PRIVATE

  TYPE, PUBLIC :: ValueTypeDV
    REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: Values
  END TYPE

  TYPE, PUBLIC :: DVIDType
    INTEGER :: iPressure
    INTEGER :: iEntropyPerBaryon
    INTEGER :: iInternalEnergyDensity
    INTEGER :: iElectronChemicalPotential
    INTEGER :: iProtonChemicalPotential
    INTEGER :: iNeutronChemicalPotential
    INTEGER :: iProtonMassFraction
    INTEGER :: iNeutronMassFraction
    INTEGER :: iAlphaMassFraction
    INTEGER :: iHeavyMassFraction
    INTEGER :: iHeavyChargeNumber
    INTEGER :: iHeavyMassNumber
    INTEGER :: iHeavyBindingEnergy
    INTEGER :: iThermalEnergy
    INTEGER :: iGamma1
  END TYPE

  TYPE, PUBLIC :: DependentVariablesType
    INTEGER :: nVariables
    INTEGER, DIMENSION(3) :: nPoints
    CHARACTER(LEN=32), DIMENSION(:), ALLOCATABLE :: Names
    CHARACTER(LEN=32), DIMENSION(:), ALLOCATABLE :: Units
    REAL(dp), DIMENSION(:), ALLOCATABLE :: Offsets
    INTEGER, DIMENSION(:,:,:), ALLOCATABLE   :: Repaired
    REAL(dp), DIMENSION(:), ALLOCATABLE :: minValues
    REAL(dp), DIMENSION(:), ALLOCATABLE :: maxValues
    TYPE(ValueTypeDV), DIMENSION(:), ALLOCATABLE :: Variables
    TYPE(DVIDType) :: Indices
  END TYPE

  PUBLIC AllocateDependentVariables
  PUBLIC DeAllocateDependentVariables

CONTAINS

  SUBROUTINE AllocateDependentVariables( DV, nPoints, nVariables )

    TYPE(DependentVariablesType)      :: DV
    INTEGER,               INTENT(in) :: nVariables
    INTEGER, DIMENSION(3), INTENT(in) :: nPoints

    INTEGER :: i

    DV % nPoints    = nPoints
    DV % nVariables = nVariables

    ALLOCATE( DV % Names( nVariables ) )
    ALLOCATE( DV % Units( nVariables ) )
    ALLOCATE( DV % Offsets( nVariables ) )
    ALLOCATE( DV % Variables( nVariables ) )
    ALLOCATE( DV % minValues( nVariables ) )
    ALLOCATE( DV % maxValues( nVariables ) )

    ALLOCATE( DV % Repaired(1:nPoints(1), 1:nPoints(2), 1:nPoints(3)) )

    DO i = 1, nVariables
      ALLOCATE &
        ( DV % Variables(i) &
             % Values(1:nPoints(1), 1:nPoints(2), 1:nPoints(3)) )
    END DO

  END SUBROUTINE AllocateDependentVariables

  SUBROUTINE DeAllocateDependentVariables( DV )

    TYPE(DependentVariablesType) :: DV

    INTEGER :: i

    DO i = 1, DV % nVariables
      DEALLOCATE( DV % Variables(i) % Values )
    END DO

    DEALLOCATE( DV % Repaired )

    DEALLOCATE( DV % Variables )
    DEALLOCATE( DV % Offsets )
    DEALLOCATE( DV % Units )
    DEALLOCATE( DV % Names )
    DEALLOCATE( DV % minValues )
    DEALLOCATE( DV % maxValues )

  END SUBROUTINE DeAllocateDependentVariables

END MODULE wlDependentVariablesModule

max_stars_repo_name: srichers/weaklib