text:
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK                          .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c ... file visgau.f
c
c     contains documentation and code for subroutine visgau
c
      SUBROUTINE DVISEQU(NLAT,NLON,H,LEN,EYER,EYELAT,EYELON,WK,LWK,IWK,
     +                  LIWK,IERROR)
      DOUBLE PRECISION H
      DOUBLE PRECISION EYER
      DOUBLE PRECISION EYELAT
      DOUBLE PRECISION EYELON
      DOUBLE PRECISION WK
      DOUBLE PRECISION PI
      DOUBLE PRECISION DTR
      DOUBLE PRECISION XEYE
      DOUBLE PRECISION YEYE
      DOUBLE PRECISION ZEYE
c
c     subroutine visequ will display a function on the sphere
c     as a solid. ie. as a "lumpy" sphere. visequ calls subroutine
c     vsurf to produce the visible surface rendering.
c
c     requires routines visequ1 interp sptc diag stride triang vsurf
c                       vsurf1 prjct box icvmg projct
c
c     visgeo uses the ncar graphics package.
c     compile with: ncargf77 (all programs above)
c
c     execute with:  a.out
c
c     on screen display with:  ctrans -d x11 gmeta
c
c     print with:  ctrans -d ps.color gmeta > gmeta.ps
c                  lpr -p(your printer) gmeta.ps
c
c     input parameters
c
c     nlat   the number of colatitudes on the full sphere including the
c            poles. for example, nlat = 37 for a five degree grid.
c            nlat determines the grid increment in colatitude as
c            pi/(nlat-1).  if nlat is odd the equator is located at
c            grid point i=(nlat+1)/2. if nlat is even the equator is
c            located half way between points i=nlat/2 and i=nlat/2+1.
c            nlat must be at least 3. note: on the half sphere, the
c            number of grid points in the colatitudinal direction is
c            nlat/2 if nlat is even or (nlat+1)/2 if nlat is odd.
c
c     nlon   the number of distinct londitude points.  nlon determines
c            the grid increment in longitude as 2*pi/nlon. for example
c            nlon = 72 for a five degree grid. nlon must be greater
c            than or equal to 4. the efficiency of the computation is
c            improved when nlon is a product of small prime numbers.
c
c
c     h      a two dimensional array that contains the discrete
c            function to be displayed. h(i,j) is the distance from the
c            center of the sphere to the surface at the colatitude
c            point theta(i) = (i-1)*pi/(nlat-1) and longitude point
c            phi(j) = (j-1)*2*pi/nlon.
c
c     len    the first dimension of the array h as it appears in the
c            program that calls visequ.
c
c     eyer   the distance from the center of the sphere to the eye.
c
c     eyelat the colatitudinal coordinate of the eye (in degrees).
c
c     eyelon the longitudinal  coordinate of the eye (in degrees).
c
c     wk     a real work array
c
c     lwk    the dimension of the array wk as it appears in the
c            program that calls visequ. lwk must be at least
c                       46*nlat*(nlon+1).
c
c     iwk    an integer work array
c
c     liwk   the dimension of the array iwk as it appears in the
c            program that calls visequ. liwk must be at least
c                       14*nlat*(nlon+1).
c
c     ierror = 0    no error
c            = 1    the eye is positioned inside the sphere
c            = 2    lwk  is less than 46*nlat*(nlon+1)
c            = 3    liwk is less than 14*nlat*(nlon+1)
c
c
      DIMENSION H(LEN,*),WK(*)
      INTEGER IWK(*)

      N = NLAT
      M = NLON + 1
      MN = M*N
      IERROR = 2
      IF (LWK.LT.46*MN) RETURN
      IERROR = 3
      IF (LIWK.LT.14*MN) RETURN
      IERROR = 1
      DO 10 J = 1,NLON
          DO 10 I = 1,NLAT
              IF (EYER.LE.H(I,J)) RETURN
   10 CONTINUE
      IERROR = 0
      PI = 4.D0*ATAN(1.D0)
      DTR = PI/180.D0
c
c     ****     set up pointers to sub work arrays in wk and iwk
c
      NTRI = MN + MN
      NW1 = 1
      NW2 = NW1 + MN
      NCLAT = 1
      NSLAT = NCLAT + N
      NXP = 1
      NYP = NXP + MN
      NX1 = 1
      NY1 = NX1 + NTRI
      NZ1 = NY1 + NTRI
      NX2 = NZ1 + NTRI
      NY2 = NX2 + NTRI
      NZ2 = NY2 + NTRI
      NX3 = NZ2 + NTRI
      NY3 = NX3 + NTRI
      NZ3 = NY3 + NTRI
      NX = NZ3 + NTRI
      NY = NX + MN
      NZ = NY + MN
      NWRK = NX
      NITYPE = 1
      NIFLAG = NTRI + 1
      NMST = NIFLAG + MN
      NMFAC = NMST + N
c     total iwk is 7*ntri
c     total wk is 58*nlat*(nlon+1)
c     ****     mid-cell interpolation, calculation of polar values
      CALL DINTERP(H,LEN,M,N,WK(NW1),WK(NW2),IWK(NIFLAG))
c     ****     transform grid points to cartesian coordinates
      CALL DSPTC(H,LEN,M,N,WK(NCLAT),WK(NSLAT),WK(NX),WK(NY),WK(NZ))
c     ****     transform eye position to cartesian coordinates
      XEYE = EYER*SIN(DTR*EYELAT)
      YEYE = XEYE*SIN(DTR*EYELON)
      XEYE = XEYE*COS(DTR*EYELON)
      ZEYE = EYER*COS(DTR*EYELAT)
c     ****     project grid points
      CALL DPROJCT(M,N,XEYE,YEYE,ZEYE,WK(NX),WK(NY),WK(NZ),WK(NXP),
     +            WK(NYP))
c     ****     check for visibility of cell boundaries
      CALL DDIAG(M,N,WK(NXP),WK(NYP),IWK(NIFLAG))
c     ****     compute longitude stride as a function of latitude
      CALL DSTRIDE(M,N,IWK(NMST),IWK(NMFAC))
c     ****     perform triangulation
      CALL DTRIANG(M,N,WK(NX),WK(NY),WK(NZ),ITRI,WK(NX1),WK(NY1),
     +             WK(NZ1),WK(NX2),WK(NY2),WK(NZ2),WK(NX3),WK(NY3),
     +             WK(NZ3),IWK(NITYPE),IWK(NIFLAG),IWK(NMST))
c     ****     call surface plotting routine
      CALL DVSURF(XEYE,YEYE,ZEYE,ITRI,WK(NX1),WK(NY1),WK(NZ1),WK(NX2),
     +           WK(NY2),WK(NZ2),WK(NX3),WK(NY3),WK(NZ3),IWK(NITYPE),
     +           WK(NWRK),IWK(NIFLAG))
      RETURN
      END
      SUBROUTINE DINTERP(H,LEN,M,N,W1,W2,IFLAG)
      DOUBLE PRECISION H
      DOUBLE PRECISION W1
      DOUBLE PRECISION W2
      DOUBLE PRECISION STEN
c     ****     interpolates to mid points of grid cells using second
c     ****     order formula
      DIMENSION H(LEN,1),W1(N,M),W2(N+2,M+2),IFLAG(N,M),STEN(4,4)
      DATA STEN/.015625D0,2*-.078125D0,.015625D0,-.078125D0,2*.390625D0,
     +     2*-.078125D0,2*.390625D0,-.078125D0,.015625D0,2*-.078125D0,
     +     .015625D0/
c     ****     copy h to w2
      MM1 = M - 1
      DO 1 I = 1,MM1
          DO 1 J = 1,N
              W2(J,I+1) = H(J,I)
    1 CONTINUE
c     ****     add periodic points
      DO 2 J = 1,N
          W2(J,1) = W2(J,M)
          W2(J,M+1) = W2(J,2)
          W2(J,M+2) = W2(J,3)
    2 CONTINUE
      N1 = 2
      N2 = N - 2
c     ****     perform interpolation
c     ****     set w1 to zero

      DO 7 I = 1,M
          DO 7 J = 1,N
              W1(J,I) = 0.D0
    7 CONTINUE
c     ****     interpolate
      DO 8 K = 1,4
          DO 8 L = 1,4
              DO 8 I = 1,M - 1
                  DO 8 J = N1,N2
                      W1(J,I) = W1(J,I) + W2(J+L-2,I+K-1)*STEN(K,L)
    8 CONTINUE
c     ****     set up iflag array
c     ****     iflag(j,i)=0  if diagonal is (j,i) to (j+1,i+1)
c     ****     iflag(j,i)=16 if diagonal is (j+1,i), (j,i+1)
      DO 9 I = 1,M - 1
          DO 9 J = N1,N2
              IFLAG(J,I) = ICVMG(16,0,ABS(.5D0* (W2(J,I+1)+W2(J+1,
     +                     I+2))-W1(J,I))-ABS(.5D0* (W2(J,I+2)+W2(J+1,
     +                     I+1))-W1(J,I)))
    9 CONTINUE
      RETURN
      END
      SUBROUTINE DSPTC(R,LEN,M,N,CLAT,SLAT,X,Y,Z)
      DOUBLE PRECISION R
      DOUBLE PRECISION CLAT
      DOUBLE PRECISION SLAT
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      DOUBLE PRECISION Z
      DOUBLE PRECISION PI
      DOUBLE PRECISION DT
      DOUBLE PRECISION DP
      DOUBLE PRECISION CLON
      DOUBLE PRECISION SLON
c     ****     transforms from spherical to cartesian coordinates
      DIMENSION R(LEN,1),CLAT(N),SLAT(N),X(N,M),Y(N,M),Z(N,M)

      PI = 4.D0*ATAN(1.D0)
      DT = PI/ (N-1)
      DP = (PI+PI)/ (M-1)
      DO 10 J = 1,N
          CLAT(J) = COS((J-1)*DT)
          SLAT(J) = SIN((J-1)*DT)
   10 CONTINUE
      DO 20 I = 1,M - 1
          CLON = COS((I-1)*DP)
          SLON = SIN((I-1)*DP)
          DO 20 J = 1,N
              X(J,I) = R(J,I)*SLAT(J)
              Y(J,I) = X(J,I)*SLON
              X(J,I) = X(J,I)*CLON
              Z(J,I) = R(J,I)*CLAT(J)
   20 CONTINUE
      DO 30 J = 1,N
          X(J,M) = X(J,1)
          Y(J,M) = Y(J,1)
          Z(J,M) = Z(J,1)
   30 CONTINUE
      RETURN
      END
      SUBROUTINE DDIAG(M,N,XP,YP,IFLAG)
      DOUBLE PRECISION XP
      DOUBLE PRECISION YP
      DOUBLE PRECISION CP
c
c     ****     label visibility of cell sides
c
c     north side corresponds to j
c     south side corresponds to j+1
c     west  side corresponds to i
c     east  side corresponds to i+1
c
c     let iflag = b4 b3 b2 b1 b0 (in binary) then b0 through b3 are
c     either o or 1 depending on whether the east, south, north
c     or west side is either invisible or visible, respectively.
c
c     b4 is o if the diagonal is from (i,j) to (i+1,j+1) and 1 if
c     the diagonal is from (i,j+1) to (i+1,j).
c
      DIMENSION XP(N,M),YP(N,M),IFLAG(N,M)
c     ****     arithmetic statement function
      CP(J1,I1,J2,I2,J3,I3) = ((XP(J1,I1)-XP(J2,I2))*
     +                        (YP(J3,I3)-YP(J2,I2))-
     +                        (XP(J3,I3)-XP(J2,I2))*
     +                        (YP(J1,I1)-YP(J2,I2)))
      DO 100 J = 2,N - 2
          DO 100 I = 1,M - 1
              IF (IFLAG(J,I).GE.16) GO TO 20
              IF (CP(J+1,I+1,J+1,I,J,I).LE.0) GO TO 10
c     west and south are visible
              IFLAG(J,I) = IFLAG(J,I) + 10
   10         IF (CP(J,I,J,I+1,J+1,I+1).LE.0) GO TO 100
c     east and north are visible
              IFLAG(J,I) = IFLAG(J,I) + 5
              GO TO 100
   20         IF (CP(J+1,I,J,I,J,I+1).LE.0) GO TO 30
c     west and north are visible
              IFLAG(J,I) = IFLAG(J,I) + 12
   30         IF (CP(J,I+1,J+1,I+1,J+1,I).LE.0) GO TO 100
c     east and south are visible
              IFLAG(J,I) = IFLAG(J,I) + 3
  100 CONTINUE
c
c     classify the poles
c
      DO 200 I = 1,M - 1
          IFLAG(1,I) = 0
          IF (CP(2,I+1,2,I,1,I).GT.0) IFLAG(1,I) = 15
          IFLAG(N-1,I) = 0
          IF (CP(N,I,N-1,I,N-1,I+1).GT.0) IFLAG(N-1,I) = 31
  200 CONTINUE
      DO 250 J = 1,N - 1
          IFLAG(J,M) = IFLAG(J,1)
  250 CONTINUE
      RETURN
      END
      SUBROUTINE DSTRIDE(M,N,MST,MFAC)
      DOUBLE PRECISION TPHI
      DOUBLE PRECISION PI
      DOUBLE PRECISION DT
      DOUBLE PRECISION THETA
      DOUBLE PRECISION ST
      DIMENSION MFAC(*),MTRYH(3),MST(N),ICL(8)
      DATA MTRYH(1),MTRYH(2),MTRYH(3)/2,3,5/
      DATA ICL(1),ICL(2),ICL(3),ICL(4),ICL(5),ICL(6),ICL(7),ICL(8)/0,1,
     +     2,12,3,13,23,123/
c
c     find prime factors of m-1
c
      ML = M - 1
      NF = 0
      J = 0
  101 J = J + 1
      IF (J-3) 102,102,103
  102 MTRY = MTRYH(J)
      GO TO 104
  103 MTRY = MTRY + 2
  104 MQ = ML/MTRY
      MR = ML - MTRY*MQ
      IF (MR) 101,105,101
  105 NF = NF + 1
      MFAC(NF) = MTRY
      ML = MQ
      IF (ML.NE.1) GO TO 104
      IF (MFAC(NF).GT.2) GO TO 106
      NF = NF - 1
      MFAC(NF) = 4
  106 TPHI = .707D0/DBLE(M-1)
      NS2 = N/2
      MF1 = MFAC(NF)
      MST(1) = (M-1)/MF1
      PI = 4.D0*ATAN(1.D0)
      DT = PI/DBLE(N-1)
      JF = NF - 1
      DO 110 JDO = 2,NS2
          J = JDO
          THETA = (J-1)*DT
          ST = SIN(THETA)
          MF2 = MF1*MFAC(JF)
          IF (ABS(ST/MF1-TPHI).GT.ABS(ST/MF2-TPHI)) GO TO 115
          MST(J) = MST(J-1)
          GO TO 110
  115     MST(J) = (M-1)/MF2
          MF1 = MF2
          JF = JF - 1
          IF (JF.EQ.0) GO TO 120
  110 CONTINUE
  120 DO 125 JDO = J,NS2
          MST(JDO) = 1
  125 CONTINUE
      DO 130 JDO = 1,NS2
          MST(N-JDO) = MST(JDO)
  130 CONTINUE
c      write (6,135) (mst(j),j=1,n)
  135 FORMAT (' colatitude strides',/, (15i5))
c
      RETURN
      END
      SUBROUTINE DTRIANG(M,N,X,Y,Z,ITRI,X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,ITYP,
     +                  IFLAG,MST)
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      DOUBLE PRECISION Z
      DOUBLE PRECISION X1
      DOUBLE PRECISION Y1
      DOUBLE PRECISION Z1
      DOUBLE PRECISION X2
      DOUBLE PRECISION Y2
      DOUBLE PRECISION Z2
      DOUBLE PRECISION X3
      DOUBLE PRECISION Y3
      DOUBLE PRECISION Z3
c     ****     performs triangulation
      DIMENSION X(N,M),Y(N,M),Z(N,M),X1(1),Y1(1),Z1(1),X2(1),Y2(1),
     +          Z2(1),X3(1),Y3(1),Z3(1),ITYP(1),IFLAG(N,M),MST(N),ICL(8)
      DATA ICL(1),ICL(2),ICL(3),ICL(4),ICL(5),ICL(6),ICL(7),ICL(8)/0,1,
     +     2,12,3,13,23,123/

      ITRI = 0
      N1 = 2
      N2 = N - 2
      DO 100 J = N1,N2
          DO 100 I = 1,M - 1
              IF (IFLAG(J,I).GE.16) GO TO 50
              IF (MOD(IFLAG(J,I),16).LT.8) GO TO 70
              ITRI = ITRI + 1
              X1(ITRI) = X(J,I)
              Y1(ITRI) = Y(J,I)
              Z1(ITRI) = Z(J,I)
              X2(ITRI) = X(J+1,I)
              Y2(ITRI) = Y(J+1,I)
              Z2(ITRI) = Z(J+1,I)
              X3(ITRI) = X(J+1,I+1)
              Y3(ITRI) = Y(J+1,I+1)
              Z3(ITRI) = Z(J+1,I+1)
              ITYPH = 3
              IF (MOD(I-1,MST(J)).EQ.0) GO TO 60
              IF (MOD(IFLAG(J,I-1),2).EQ.0) GO TO 60
              ITYPH = ITYPH - 1
   60         IF (MOD(IFLAG(J,I),2).EQ.0) ITYPH = ITYPH + 4
              ITYP(ITRI) = ICL(ITYPH+1)
   70         IF (MOD(IFLAG(J,I),2).EQ.0) GO TO 100
              ITRI = ITRI + 1
              X1(ITRI) = X(J,I)
              Y1(ITRI) = Y(J,I)
              Z1(ITRI) = Z(J,I)
              X2(ITRI) = X(J+1,I+1)
              Y2(ITRI) = Y(J+1,I+1)
              Z2(ITRI) = Z(J+1,I+1)
              X3(ITRI) = X(J,I+1)
              Y3(ITRI) = Y(J,I+1)
              Z3(ITRI) = Z(J,I+1)
              ITYPH = 0
              IF (MOD(IFLAG(J,I),16).LT.8) ITYPH = ITYPH + 1
              IF (MOD(IFLAG(J,I+1),16).LT.8) ITYPH = ITYPH + 2
              IF (MOD(IFLAG(J-1,I),4).LT.2) ITYPH = ITYPH + 4
              ITYP(ITRI) = ICL(ITYPH+1)
              GO TO 100
   50         IF (MOD(IFLAG(J,I),16).LT.8) GO TO 20
              ITRI = ITRI + 1
              X1(ITRI) = X(J,I)
              Y1(ITRI) = Y(J,I)
              Z1(ITRI) = Z(J,I)
              X2(ITRI) = X(J+1,I)
              Y2(ITRI) = Y(J+1,I)
              Z2(ITRI) = Z(J+1,I)
              X3(ITRI) = X(J,I+1)
              Y3(ITRI) = Y(J,I+1)
              Z3(ITRI) = Z(J,I+1)
              ITYPH = 1
              IF (MOD(I-1,MST(J)).EQ.0) GO TO 10
              IF (MOD(IFLAG(J,I-1),2).EQ.0) GO TO 10
              ITYPH = 0
   10         IF (MOD(IFLAG(J,I),2).EQ.0) ITYPH = ITYPH + 2
              IF (MOD(IFLAG(J-1,I),4).LT.2) ITYPH = ITYPH + 4
              ITYP(ITRI) = ICL(ITYPH+1)
   20         IF (MOD(IFLAG(J,I),2).EQ.0) GO TO 100
              ITRI = ITRI + 1
              X1(ITRI) = X(J+1,I)
              Y1(ITRI) = Y(J+1,I)
              Z1(ITRI) = Z(J+1,I)
              X2(ITRI) = X(J+1,I+1)
              Y2(ITRI) = Y(J+1,I+1)
              Z2(ITRI) = Z(J+1,I+1)
              X3(ITRI) = X(J,I+1)
              Y3(ITRI) = Y(J,I+1)
              Z3(ITRI) = Z(J,I+1)
              ITYPH = 1
              IF (MOD(IFLAG(J,I+1),16).LT.8) ITYPH = ITYPH + 2
              IF (MOD(IFLAG(J,I),16).LT.8) ITYPH = ITYPH + 4
              ITYP(ITRI) = ICL(ITYPH+1)
  100 CONTINUE
c
c     ****     triangles around north and south poles
c
      DO 200 I = 1,M - 1
          IF (MOD(IFLAG(1,I),16).LT.8) GO TO 250
          ITRI = ITRI + 1
          X1(ITRI) = X(1,I)
          Y1(ITRI) = Y(1,I)
          Z1(ITRI) = Z(1,I)
          X2(ITRI) = X(2,I)
          Y2(ITRI) = Y(2,I)
          Z2(ITRI) = Z(2,I)
          X3(ITRI) = X(2,I+1)
          Y3(ITRI) = Y(2,I+1)
          Z3(ITRI) = Z(2,I+1)
          ITYPH = 3
          IF (MOD(I-1,MST(1)).EQ.0) GO TO 260
          IF (MOD(IFLAG(1,I-1),2).EQ.0) GO TO 260
          ITYPH = ITYPH - 1
  260     IF (MOD(IFLAG(1,I+1),16).LT.8) ITYPH = ITYPH + 4
          ITYP(ITRI) = ICL(ITYPH+1)
  250     IF (MOD(IFLAG(N-1,I),16).LT.8) GO TO 200
          ITRI = ITRI + 1
          X1(ITRI) = X(N-1,I)
          Y1(ITRI) = Y(N-1,I)
          Z1(ITRI) = Z(N-1,I)
          X2(ITRI) = X(N,I)
          Y2(ITRI) = Y(N,I)
          Z2(ITRI) = Z(N,I)
          X3(ITRI) = X(N-1,I+1)
          Y3(ITRI) = Y(N-1,I+1)
          Z3(ITRI) = Z(N-1,I+1)
          ITYPH = 1
          IF (MOD(I-1,MST(N-1)).EQ.0) GO TO 210
          IF (MOD(IFLAG(N-1,I-1),2).EQ.0) GO TO 210
          ITYPH = 0
  210     IF (MOD(IFLAG(N-1,I+1),16).LT.8) ITYPH = ITYPH + 2
          IF (MOD(IFLAG(N-2,I),4).LT.2) ITYPH = ITYPH + 4
          ITYP(ITRI) = ICL(ITYPH+1)
  200 CONTINUE
      RETURN
      END
      SUBROUTINE DVSURF(XEYE,YEYE,ZEYE,NTRI,X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,
     +                 ITYPE,WORK,IWORK)
      DOUBLE PRECISION XEYE
      DOUBLE PRECISION YEYE
      DOUBLE PRECISION ZEYE
      DOUBLE PRECISION X1
      DOUBLE PRECISION Y1
      DOUBLE PRECISION Z1
      DOUBLE PRECISION X2
      DOUBLE PRECISION Y2
      DOUBLE PRECISION Z2
      DOUBLE PRECISION X3
      DOUBLE PRECISION Y3
      DOUBLE PRECISION Z3
      DOUBLE PRECISION WORK
c
c    subroutine vsurf is like subroutine hidel except the triangles
c    are categorized. vsurf is also like solid except triangles rather
c    than lines are covered.
c
c     written by paul n. swarztrauber, national center for atmospheric
c     research, p.o. box 3000, boulder, colorado, 80307
c
c    this program plots visible lines for the surface defined
c    by the input 3-d triangles with corners at (x1,y1,z1), (x2,y2,z2)
c    and (x3,y3,z3). the sides of these these triangles may or
c    may not be plotted depending on itype. if itype is 1 then the
c    side between points (x1,y1,z1) and (x2,y2,z2) is plotted if it
c    is visible. if itype is 2 then the side between (x2,y2,z2)
c    and (x3,y3,z3) is plotted. if itype is 3 then the visible portion
c    of the side between (x3,y3,z3) and (x1,y1,z1) is plotted.
c    any combination is possible by specifying itype to be one
c    of the following values: 0,1,2,3,12,13,23,123.
c
c    the length of real    array  work must be at least 14*ntri
c
c    the length of integer array iwork must be at least  6*ntri
c
c
c    the vertices of the triangles are renumbered by vsurf so that
c    their projections are orientated counterclockwise. the user need
c    only be aware that the vertices may be renumbered by vsurf.
c
      DIMENSION X1(NTRI),Y1(NTRI),Z1(NTRI),X2(NTRI),Y2(NTRI),Z2(NTRI),
     +          X3(NTRI),Y3(NTRI),Z3(NTRI),ITYPE(NTRI),WORK(14*NTRI)
      INTEGER IWORK(6*NTRI)
c
      CALL DVSURF1(XEYE,YEYE,ZEYE,NTRI,X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,
     +             ITYPE,WORK,WORK(NTRI+1),WORK(2*NTRI+1),
     +             WORK(3*NTRI+1),WORK(4*NTRI+1),WORK(5*NTRI+1),
     +             WORK(6*NTRI+1),WORK(7*NTRI+1),WORK(8*NTRI+1),
     +             WORK(9*NTRI+1),WORK(10*NTRI+1),WORK(11*NTRI+1),
     +             WORK(12*NTRI+1),WORK(13*NTRI+1),IWORK,
     +             IWORK(NTRI+1),IWORK(2*NTRI+1),IWORK(4*NTRI+1))
      RETURN
      END
      SUBROUTINE DVSURF1(XEYE,YEYE,ZEYE,NTRI,X1,Y1,Z1,X2,Y2,Z2,X3,Y3,
     +                   Z3,ITYPE,PX1,PY1,PX2,PY2,PX3,PY3,VX1,VY1,VX2,
     +                   VY2,VX3,VY3,TL,TR,KH,NEXT,ISTART,IFINAL)
      DOUBLE PRECISION XEYE
      DOUBLE PRECISION YEYE
      DOUBLE PRECISION ZEYE
      DOUBLE PRECISION X1
      DOUBLE PRECISION Y1
      DOUBLE PRECISION Z1
      DOUBLE PRECISION X2
      DOUBLE PRECISION Y2
      DOUBLE PRECISION Z2
      DOUBLE PRECISION X3
      DOUBLE PRECISION Y3
      DOUBLE PRECISION Z3
      DOUBLE PRECISION PX1
      DOUBLE PRECISION PY1
      DOUBLE PRECISION PX2
      DOUBLE PRECISION PY2
      DOUBLE PRECISION PX3
      DOUBLE PRECISION PY3
      DOUBLE PRECISION VX1
      DOUBLE PRECISION VY1
      DOUBLE PRECISION VX2
      DOUBLE PRECISION VY2
      DOUBLE PRECISION VX3
      DOUBLE PRECISION VY3
      DOUBLE PRECISION TL
      DOUBLE PRECISION TR
      DOUBLE PRECISION FNTRI
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      DOUBLE PRECISION Z
      DOUBLE PRECISION DUM1
      DOUBLE PRECISION DUM2
      DOUBLE PRECISION CPROD
      DOUBLE PRECISION PX1H
      DOUBLE PRECISION PY1H
      DOUBLE PRECISION X1HOLD
      DOUBLE PRECISION Y1HOLD
      DOUBLE PRECISION Z1HOLD
      DOUBLE PRECISION PMAX
      DOUBLE PRECISION PMIN
      DOUBLE PRECISION XMIN
      DOUBLE PRECISION XMAX
      DOUBLE PRECISION YMIN
      DOUBLE PRECISION YMAX
      DOUBLE PRECISION DMX
      DOUBLE PRECISION DMY
      DOUBLE PRECISION C
      DOUBLE PRECISION D
      DOUBLE PRECISION XMID
      DOUBLE PRECISION HDY
      DOUBLE PRECISION A
      DOUBLE PRECISION B
      DOUBLE PRECISION YMID
      DOUBLE PRECISION HDX
      DOUBLE PRECISION HGR
      DOUBLE PRECISION DXT
      DOUBLE PRECISION HR
      DOUBLE PRECISION TL1
      DOUBLE PRECISION TL2
      DOUBLE PRECISION PX4
      DOUBLE PRECISION PY4
      DOUBLE PRECISION PX5
      DOUBLE PRECISION PY5
      DOUBLE PRECISION X4
      DOUBLE PRECISION Y4
      DOUBLE PRECISION Z4
      DOUBLE PRECISION X5
      DOUBLE PRECISION Y5
      DOUBLE PRECISION Z5
      DOUBLE PRECISION X54
      DOUBLE PRECISION Y54
      DOUBLE PRECISION C17
      DOUBLE PRECISION C27
      DOUBLE PRECISION C37
      DOUBLE PRECISION C14
      DOUBLE PRECISION C25
      DOUBLE PRECISION C36
      DOUBLE PRECISION TMIN
      DOUBLE PRECISION TMAX
      DOUBLE PRECISION XPL
      DOUBLE PRECISION YPL
      DOUBLE PRECISION ZPL
      DOUBLE PRECISION XPR
      DOUBLE PRECISION YPR
      DOUBLE PRECISION ZPR
      DOUBLE PRECISION VX1T
      DOUBLE PRECISION VY1T
      DOUBLE PRECISION VZ1T
      DOUBLE PRECISION VX2T
      DOUBLE PRECISION VY2T
      DOUBLE PRECISION VZ2T
      DOUBLE PRECISION APL
      DOUBLE PRECISION BPL
      DOUBLE PRECISION CPL
      DOUBLE PRECISION DPL
      DOUBLE PRECISION VX3T
      DOUBLE PRECISION VY3T
      DOUBLE PRECISION VZ3T
      DOUBLE PRECISION DEN
      DOUBLE PRECISION TIL
      DOUBLE PRECISION TIR
      DOUBLE PRECISION TIM
      DOUBLE PRECISION THOLD
      DOUBLE PRECISION TLH
      DOUBLE PRECISION TRH
      DOUBLE PRECISION XA
      DOUBLE PRECISION YA
      DOUBLE PRECISION XB
      DOUBLE PRECISION YB
c
      DIMENSION X1(NTRI),Y1(NTRI),Z1(NTRI),X2(NTRI),Y2(NTRI),Z2(NTRI),
     +          X3(NTRI),Y3(NTRI),Z3(NTRI),ITYPE(NTRI),PX1(NTRI),
     +          PY1(NTRI),PX2(NTRI),PY2(NTRI),PX3(NTRI),PY3(NTRI),
     +          VX1(NTRI),VY1(NTRI),VX2(NTRI),VY2(NTRI),VX3(NTRI),
     +          VY3(NTRI),TL(NTRI),TR(NTRI),NEXT(NTRI),KH(NTRI),
     +          ISTART(2*NTRI),IFINAL(2*NTRI),LTP(3),IRD(11),IP2(11),
     +          NCT(11),NCV(11),LAST(11)
c
      DOUBLE PRECISION L2E
C*PT*WARNING* Already double-precision
      DOUBLE PRECISION LE2
C*PT*WARNING* Constant already double-precision
c
c     compute projections of 3-d points
c
      LE2 = .6931471805599453094172321d0
C*PT*WARNING* Constant already double-precision
      L2E = 1.d0/LE2
      FNTRI = NTRI
      IRMAX = .5D0*L2E*LOG(FNTRI)
      IRMAX = MIN(IRMAX,10)
      IRMP1 = IRMAX + 1
      DO 4 ICV = 1,11
          NCV(ICV) = 0
    4 CONTINUE
      NCT(1) = 0
      IP2(1) = 1
      IRD(1) = 0
      ISIZE = 4
      DO 7 IRP1 = 2,IRMP1
          IR = IRP1 - 1
          NCT(IRP1) = 0
          IP2(IRP1) = 2**IR
          IRD(IRP1) = IRD(IR) + ISIZE
          ISIZE = (IP2(IRP1)+1)**2
    7 CONTINUE
      ISXM = IRD(IRMP1) + ISIZE + 1
      DO 8 ISX = 1,ISXM
          ISTART(ISX) = 0
          IFINAL(ISX) = 0
    8 CONTINUE
      DO 6 I = 1,NTRI
          NEXT(I) = 0
    6 CONTINUE
      CALL DPRJCT(0,XEYE,YEYE,ZEYE,X,Y,Z,DUM1,DUM2)
c      write(6,127) ntri
  127 FORMAT (' ntri in hidel',i5)
      DO 86 K = 1,NTRI
          CALL DPRJCT(1,XEYE,YEYE,ZEYE,X1(K),Y1(K),Z1(K),PX1(K),PY1(K))
          CALL DPRJCT(1,XEYE,YEYE,ZEYE,X2(K),Y2(K),Z2(K),PX2(K),PY2(K))
          CALL DPRJCT(1,XEYE,YEYE,ZEYE,X3(K),Y3(K),Z3(K),PX3(K),PY3(K))
          IF (K.LT.3) THEN
c          write(6,333) xeye,yeye,zeye,x1(k),y1(k),z1(k),px1(k),py1(k)
  333         FORMAT (' xeye, etc.',8D8.1)
          END IF
   86 CONTINUE
c
c     orientate triangles counter clockwise
c
      DO 70 K = 1,NTRI
          CPROD = (PX2(K)-PX1(K))* (PY3(K)-PY1(K)) -
     +            (PY2(K)-PY1(K))* (PX3(K)-PX1(K))
c      if(cprod.eq.0.) write(6,79) k,px1(k),px2(k),px3(k),
c     -                              py1(k),py2(k),py3(k)
   79     FORMAT ('  cprod=0 at k=',i5,6D9.2)
          IF (CPROD.GE.0.D0) GO TO 70
          PX1H = PX1(K)
          PY1H = PY1(K)
          PX1(K) = PX2(K)
          PY1(K) = PY2(K)
          PX2(K) = PX1H
          PY2(K) = PY1H
          X1HOLD = X1(K)
          Y1HOLD = Y1(K)
          Z1HOLD = Z1(K)
          X1(K) = X2(K)
          Y1(K) = Y2(K)
          Z1(K) = Z2(K)
          X2(K) = X1HOLD
          Y2(K) = Y1HOLD
          Z2(K) = Z1HOLD
          ITYP = ITYPE(K)
          IF (ITYP.EQ.2) ITYPE(K) = 3
          IF (ITYP.EQ.3) ITYPE(K) = 2
          IF (ITYP.EQ.12) ITYPE(K) = 13
          IF (ITYP.EQ.13) ITYPE(K) = 12
   70 CONTINUE
c
c     set screen limits
c
      PMAX = PX1(1)
      PMIN = PX1(1)
      DO 87 K = 1,NTRI
          PMIN = DMIN1(PMIN,PX1(K),PY1(K),PX2(K),PY2(K),PX3(K),PY3(K))
          PMAX = DMAX1(PMAX,PX1(K),PY1(K),PX2(K),PY2(K),PX3(K),PY3(K))
   87 CONTINUE
      PMIN = 1.1D0*PMIN
      PMAX = 1.1D0*PMAX
      CALL SET(0.D0,1.D0,0.D0,1.D0,PMIN,PMAX,PMIN,PMAX,1)
      XMIN = DMIN1(PX1(1),PX2(1),PX3(1))
      XMAX = DMAX1(PX1(1),PX2(1),PX3(1))
      YMIN = DMIN1(PY1(1),PY2(1),PY3(1))
      YMAX = DMAX1(PY1(1),PY2(1),PY3(1))
      DO 1 I = 2,NTRI
          XMIN = DMIN1(XMIN,PX1(I),PX2(I),PX3(I))
          XMAX = DMAX1(XMAX,PX1(I),PX2(I),PX3(I))
          YMIN = DMIN1(YMIN,PY1(I),PY2(I),PY3(I))
          YMAX = DMAX1(YMAX,PY1(I),PY2(I),PY3(I))
    1 CONTINUE
      DMX = XMAX - XMIN
      DMY = YMAX - YMIN
      IF (DMX.GT.DMY) GO TO 2
      C = YMIN
      D = YMAX
      XMID = .5D0* (XMIN+XMAX)
      HDY = .5D0*DMY
      A = XMID - HDY
      B = XMID + HDY
      GO TO 3
    2 A = XMIN
      B = XMAX
      YMID = .5D0* (YMIN+YMAX)
      HDX = .5D0*DMX
      C = YMID - HDX
      D = YMID + HDX
    3 HGR = B - A
c
c     categorize triangles
c
      DO 100 I = 1,NTRI
          XMIN = DMIN1(PX1(I),PX2(I),PX3(I))
          XMAX = DMAX1(PX1(I),PX2(I),PX3(I))
          YMIN = DMIN1(PY1(I),PY2(I),PY3(I))
          YMAX = DMAX1(PY1(I),PY2(I),PY3(I))
          DXT = DMAX1(XMAX-XMIN,YMAX-YMIN)
          IF (DXT.GT.0.D0) GO TO 10
          IR = IRMAX
          GO TO 20
   10     IR = L2E*LOG(HGR/DXT)
          IR = MIN(IR,IRMAX)
   20     IRP1 = IR + 1
          NCT(IRP1) = NCT(IRP1) + 1
          HR = HGR/IP2(IRP1)
          XMID = .5D0* (XMIN+XMAX)
          ID = (XMID-A)/HR + 1.5D0
          YMID = .5D0* (YMIN+YMAX)
          JD = (YMID-C)/HR + 1.5D0
          IJD = IP2(IRP1) + 1
          ISX = ID + (JD-1)*IJD + IRD(IRP1)
          IFX = IFINAL(ISX)
          IF (IFX.GT.0) GO TO 50
          ISTART(ISX) = I
          GO TO 60
   50     NEXT(IFX) = I
   60     IFINAL(ISX) = I
  100 CONTINUE
c      write(6,106) tcat,(irp1,nct(irp1),irp1=1,irmp1)
  106 FORMAT (' time to categorize   ',D15.6,/, (' ir+1',i3,' ntri',i7))
c
c     sort triangles into boxes
c
      L = 0
      DO 30 IRP1 = 1,IRMP1
          IF (NCT(IRP1).EQ.0) GO TO 30
          IST = IRD(IRP1) + 1
          ISD = IP2(IRP1) + 1
          CALL DBOX(ISD,ISTART(IST),NEXT,L,IFINAL)
          LAST(IRP1) = L + 1
   30 CONTINUE
      DO 35 IRP1 = 1,IRMP1
          IL = IRD(IRP1) + (IP2(IRP1)+1)**2 + 1
          IF (ISTART(IL).EQ.0) ISTART(IL) = LAST(IRP1)
   35 CONTINUE
c      write(6,31) tsort,l,ntri
   31 FORMAT (' time to sort  ',D15.6,'   l',i8,'   ntri',i8)
      DO 90 K = 1,NTRI
          VX1(K) = PX2(K) - PX1(K)
          VY1(K) = PY2(K) - PY1(K)
          VX2(K) = PX3(K) - PX2(K)
          VY2(K) = PY3(K) - PY2(K)
          VX3(K) = PX1(K) - PX3(K)
          VY3(K) = PY1(K) - PY3(K)
   90 CONTINUE
      TL1 = 0.D0
      TL2 = 0.D0
      MAXS = 0
      DO 500 IR2 = 1,IRMP1
          IF (NCT(IR2).EQ.0) GO TO 500
          IST = IRD(IR2)
          ISD = IP2(IR2) + 1
          DO 490 J2 = 1,ISD
              DO 480 I2 = 1,ISD
                  IST = IST + 1
                  LS = ISTART(IST)
                  LF = ISTART(IST+1) - 1
                  IF (LF.LT.LS) GO TO 480
c
c     define coverings
c
                  KCV = 0
                  I2M = I2 - 1
                  J2M = J2 - 1
                  DO 300 IR1 = 1,IRMP1
                      IF (NCT(IR1).EQ.0) GO TO 300
                      IF (IR1.GE.IR2) GO TO 260
                      IRDP = 2** (IR2-IR1)
                      I1S = (I2M-1)/IRDP
                      I1F = (I2M+1)/IRDP
                      IF = I2M + 1 - I1F*IRDP
                      IF (IF.GT.0) I1F = I1F + 1
                      J1S = (J2M-1)/IRDP
                      J1F = (J2M+1)/IRDP
                      JF = J2M + 1 - J1F*IRDP
                      IF (JF.GT.0) J1F = J1F + 1
                      GO TO 270
  260                 IRDP = 2** (IR1-IR2)
                      I1S = IRDP* (I2M-1)
                      I1F = IRDP* (I2M+1)
                      J1S = IRDP* (J2M-1)
                      J1F = IRDP* (J2M+1)
  270                 IJD = IP2(IR1) + 1
                      I1S = MAX(I1S+1,1)
                      I1F = MIN(I1F+1,IJD)
                      J1S = MAX(J1S+1,1)
                      J1F = MIN(J1F+1,IJD)
                      IXH = (J1S-2)*IJD + IRD(IR1)
                      IXS = I1S + IXH
                      IXF = I1F + IXH
                      DO 290 J1 = J1S,J1F
                          IXS = IXS + IJD
                          KDS = ISTART(IXS)
                          IXF = IXF + IJD
                          KDF = ISTART(IXF+1) - 1
                          IF (KDF.LT.KDS) GO TO 290
                          DO 280 KD = KDS,KDF
                              KCV = KCV + 1
                              KH(KCV) = IFINAL(KD)
  280                     CONTINUE
  290                 CONTINUE
  300             CONTINUE
                  DO 310 ICV = 1,10
                      IF (KCV.LE.NCV(ICV)) GO TO 310
                      NCV(ICV) = KCV
                      GO TO 320
  310             CONTINUE
c
c
  320             DO 470 LDO = LS,LF
                      L = IFINAL(LDO)
                      ITH = ITYPE(L)
                      IF (ITH.EQ.0) GO TO 470
                      LTP(1) = 0
                      LTP(2) = 0
                      LTP(3) = 0
                      ID1 = ITH/100
                      ITH = ITH - 100*ID1
                      ID2 = ITH/10
                      ID3 = ITH - 10*ID2
                      IF (ID1.NE.0) LTP(ID1) = 1
                      IF (ID2.NE.0) LTP(ID2) = 1
                      IF (ID3.NE.0) LTP(ID3) = 1
c     if((ith.eq.123) .or. (ith.eq.12) .or.(ith.eq.13)) ltp(1) = 1
c     if((ith.eq.123) .or. (ith.eq.23) .or.(ith.eq.12)) ltp(2) = 1
c     if((ith.eq.123) .or. (ith.eq.13) .or.(ith.eq.23)) ltp(3) = 1
                      DO 460 NS = 1,3
                          GO TO (101,102,103) NS
  101                     IF (LTP(NS).EQ.0) GO TO 460
                          PX4 = PX1(L)
                          PY4 = PY1(L)
                          PX5 = PX2(L)
                          PY5 = PY2(L)
                          X4 = X1(L)
                          Y4 = Y1(L)
                          Z4 = Z1(L)
                          X5 = X2(L)
                          Y5 = Y2(L)
                          Z5 = Z2(L)
                          GO TO 105
  102                     IF (LTP(NS).EQ.0) GO TO 460
                          PX4 = PX2(L)
                          PY4 = PY2(L)
                          PX5 = PX3(L)
                          PY5 = PY3(L)
                          X4 = X2(L)
                          Y4 = Y2(L)
                          Z4 = Z2(L)
                          X5 = X3(L)
                          Y5 = Y3(L)
                          Z5 = Z3(L)
                          GO TO 105
  103                     IF (LTP(NS).EQ.0) GO TO 460
                          PX4 = PX1(L)
                          PY4 = PY1(L)
                          PX5 = PX3(L)
                          PY5 = PY3(L)
                          X4 = X1(L)
                          Y4 = Y1(L)
                          Z4 = Z1(L)
                          X5 = X3(L)
                          Y5 = Y3(L)
                          Z5 = Z3(L)
  105                     X54 = PX5 - PX4
                          Y54 = PY5 - PY4
                          NSEG = 0
                          DO 440 KD = 1,KCV
                              K = KH(KD)
                              C17 = VX1(K)*Y54 - VY1(K)*X54
                              C27 = VX2(K)*Y54 - VY2(K)*X54
                              C37 = VX3(K)*Y54 - VY3(K)*X54
                              C14 = VY1(K)* (PX4-PX1(K)) -
     +                              VX1(K)* (PY4-PY1(K))
                              C25 = VY2(K)* (PX4-PX2(K)) -
     +                              VX2(K)* (PY4-PY2(K))
                              C36 = VY3(K)* (PX4-PX3(K)) -
     +                              VX3(K)* (PY4-PY3(K))
                              TMIN = 0.D0
                              TMAX = 1.D0
                              IF (C17) 151,152,153
  151                         TMAX = DMIN1(C14/C17,TMAX)
                              GO TO 154
  152                         IF (C14) 154,440,440
  153                         TMIN = DMAX1(C14/C17,TMIN)
  154                         IF (C27) 155,156,157
  155                         TMAX = DMIN1(C25/C27,TMAX)
                              GO TO 158
  156                         IF (C25) 158,440,440
  157                         TMIN = DMAX1(C25/C27,TMIN)
  158                         IF (C37) 159,160,161
  159                         TMAX = DMIN1(C36/C37,TMAX)
                              GO TO 162
  160                         IF (C36) 162,440,440
  161                         TMIN = DMAX1(C36/C37,TMIN)
  162                         IF (TMAX-TMIN.LT..00001D0) GO TO 440
                              XPL = X4 + TMIN* (X5-X4)
                              YPL = Y4 + TMIN* (Y5-Y4)
                              ZPL = Z4 + TMIN* (Z5-Z4)
                              XPR = X4 + TMAX* (X5-X4)
                              YPR = Y4 + TMAX* (Y5-Y4)
                              ZPR = Z4 + TMAX* (Z5-Z4)
c
c     the projections of line and plane intersect
c     now determine if plane covers line
c
                              VX1T = X2(K) - X1(K)
                              VY1T = Y2(K) - Y1(K)
                              VZ1T = Z2(K) - Z1(K)
                              VX2T = X3(K) - X1(K)
                              VY2T = Y3(K) - Y1(K)
                              VZ2T = Z3(K) - Z1(K)
                              APL = VY1T*VZ2T - VY2T*VZ1T
                              BPL = VX2T*VZ1T - VX1T*VZ2T
                              CPL = VX1T*VY2T - VX2T*VY1T
                              DPL = APL*X1(K) + BPL*Y1(K) + CPL*Z1(K)
                              VX3T = XPL - XEYE
                              VY3T = YPL - YEYE
                              VZ3T = ZPL - ZEYE
                              DEN = APL*VX3T + BPL*VY3T + CPL*VZ3T
                              TIL = 0.D0
                              IF (DEN.EQ.0.D0) GO TO 410
                              TIL = (DPL-APL*XEYE-BPL*YEYE-CPL*ZEYE)/DEN
  410                         VX3T = XPR - XEYE
                              VY3T = YPR - YEYE
                              VZ3T = ZPR - ZEYE
                              DEN = APL*VX3T + BPL*VY3T + CPL*VZ3T
                              TIR = 0.D0
                              IF (DEN.EQ.0.D0) GO TO 412
                              TIR = (DPL-APL*XEYE-BPL*YEYE-CPL*ZEYE)/DEN
  412                         IF (TIL.GE..99999D0 .AND.
     +                            TIR.GE..99999D0) GO TO 440
                              IF (TIL.LT.1.D0 .AND.
     +                            TIR.LT.1.D0) GO TO 164
                              VX3T = XPR - XPL
                              VY3T = YPR - YPL
                              VZ3T = ZPR - ZPL
                              DEN = APL*VX3T + BPL*VY3T + CPL*VZ3T
                              TIM = 0.D0
                              IF (DEN.EQ.0.D0) GO TO 414
                              TIM = (DPL-APL*XPL-BPL*YPL-CPL*ZPL)/DEN
  414                         THOLD = TMIN + TIM* (TMAX-TMIN)
                              IF (TIL.GE.1.D0) GO TO 163
                              TMAX = THOLD
                              GO TO 164
  163                         TMIN = THOLD
  164                         NSEG = NSEG + 1
                              TL(NSEG) = TMIN
                              TR(NSEG) = TMAX
  440                     CONTINUE
                          MAXS = MAX0(MAXS,NSEG)
                          IF (NSEG-1) 171,180,172
  171                     CALL LINE(PX4,PY4,PX5,PY5)
                          GO TO 460
c
c     order the segments according to left end point tl(k)
c
  172                     DO 173 K = 2,NSEG
                              DO 173 I = K,NSEG
                                  IF (TL(K-1).LE.TL(I)) GO TO 173
                                  TLH = TL(K-1)
                                  TRH = TR(K-1)
                                  TL(K-1) = TL(I)
                                  TR(K-1) = TR(I)
                                  TL(I) = TLH
                                  TR(I) = TRH
  173                     CONTINUE
c
c     eliminate segment overlap
c
                          K1 = 1
                          K2 = 1
  174                     K2 = K2 + 1
                          IF (K2.GT.NSEG) GO TO 176
                          IF (TR(K1).LT.TL(K2)) GO TO 175
                          TR(K1) = DMAX1(TR(K1),TR(K2))
                          GO TO 174
  175                     K1 = K1 + 1
                          TL(K1) = TL(K2)
                          TR(K1) = TR(K2)
                          GO TO 174
  176                     NSEG = K1
c
c     plot all segments of the line
c
  180                     DO 181 KS = 1,NSEG
                              KB = NSEG - KS + 1
                              TL(KB+1) = TR(KB)
                              TR(KB) = TL(KB)
  181                     CONTINUE
                          TL(1) = 0.D0
                          TR(NSEG+1) = 1.D0
                          NSEGP = NSEG + 1
                          DO 450 K = 1,NSEGP
                              IF (ABS(TR(K)-TL(K)).LT.
     +                            .000001D0) GO TO 450
                              XA = PX4 + TL(K)* (PX5-PX4)
                              YA = PY4 + TL(K)* (PY5-PY4)
                              XB = PX4 + TR(K)* (PX5-PX4)
                              YB = PY4 + TR(K)* (PY5-PY4)
                              CALL LINE(XA,YA,XB,YB)
  450                     CONTINUE
  460                 CONTINUE
  470             CONTINUE
  480         CONTINUE
  490     CONTINUE
  500 CONTINUE
c      write(6,903) tl1,tl2
  903 FORMAT (' time to cover',D15.6,/,' time to test ',D15.6)
c      write(6,904) maxs
  904 FORMAT (' maximum number of segments',i5)
c      write(6,250) (ncv(icv),icv=1,10)
  250 FORMAT ('  the ten largest coverings',/, (10i5))
      CALL FRAME
      END
      SUBROUTINE DPRJCT(INIT,XEYE,YEYE,ZEYE,X,Y,Z,PX,PY)
      DOUBLE PRECISION XEYE
      DOUBLE PRECISION YEYE
      DOUBLE PRECISION ZEYE
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      DOUBLE PRECISION Z
      DOUBLE PRECISION PX
      DOUBLE PRECISION PY
      DOUBLE PRECISION RADS1
      DOUBLE PRECISION RADS2
      DOUBLE PRECISION D1
      DOUBLE PRECISION D2
      DOUBLE PRECISION CX1
      DOUBLE PRECISION CY1
      DOUBLE PRECISION CX2
      DOUBLE PRECISION CY2
      DOUBLE PRECISION CZ2
      DOUBLE PRECISION CX3
      DOUBLE PRECISION CY3
      DOUBLE PRECISION CZ3
      DOUBLE PRECISION X1
      DOUBLE PRECISION Y1
      DOUBLE PRECISION Z1
      DOUBLE PRECISION RATIO
c
c     subroutine prjct projects the point x,y,z onto a plane through
c     the origin that is perpendicular to a line between the origin
c     and the eye. the projection is along the line between the eye
c     and the point x,y,z. px and py are the coordinates of the
c     projection in the plane.
c     (version 2 , 12-10-82)
c
      SAVE

      IF (INIT.NE.0) GO TO 1
      RADS1 = XEYE**2 + YEYE**2
      RADS2 = RADS1 + ZEYE**2
      D1 = SQRT(RADS1)
      D2 = SQRT(RADS2)
      CX1 = -YEYE/D1
      CY1 = XEYE/D1
      CX2 = -XEYE*ZEYE/ (D1*D2)
      CY2 = -YEYE*ZEYE/ (D1*D2)
      CZ2 = D1/D2
      CX3 = XEYE/D2
      CY3 = YEYE/D2
      CZ3 = ZEYE/D2
      RETURN
    1 X1 = CX1*X + CY1*Y
      Y1 = CX2*X + CY2*Y + CZ2*Z
      Z1 = CX3*X + CY3*Y + CZ3*Z
      RATIO = D2/ (D2-Z1)
      PX = RATIO*X1
      PY = RATIO*Y1
      RETURN
      END
      SUBROUTINE DBOX(ISD,ISTART,NEXT,L,LIST)
      DIMENSION ISTART(ISD,ISD),NEXT(1),LIST(1)

      DO 30 JD = 1,ISD
          DO 10 ID = 1,ISD
              IDX = ISTART(ID,JD)
              ISTART(ID,JD) = L + 1
              IF (IDX.EQ.0) GO TO 10
   20         L = L + 1
              LIST(L) = IDX
              IF (NEXT(IDX).EQ.0) GO TO 10
              IDX = NEXT(IDX)
              GO TO 20
   10     CONTINUE
   30 CONTINUE
      RETURN
      END
      INTEGER FUNCTION ICVMG(I1,I2,R)
      INTEGER I1,I2
      DOUBLE PRECISION R
c
c     returns i1 if i3.ge.0 and returns i2 if i3.lt.0 .
c
      ICVMG = I1
      IF (R.LT.0.D0) ICVMG = I2
      RETURN
      END
      SUBROUTINE DPROJCT(M,N,XEYE,YEYE,ZEYE,X,Y,Z,PX,PY)
      DOUBLE PRECISION XEYE
      DOUBLE PRECISION YEYE
      DOUBLE PRECISION ZEYE
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      DOUBLE PRECISION Z
      DOUBLE PRECISION PX
      DOUBLE PRECISION PY
      DOUBLE PRECISION RDUM1
      DOUBLE PRECISION RDUM2
      DOUBLE PRECISION RDUM3
      DOUBLE PRECISION RDUM4
      DOUBLE PRECISION RDUM5
c     ****     projects point (x,y,z) onto plane thru origin and perp
c     ****     to line joining origin and eye
      DIMENSION X(N,M),Y(N,M),Z(N,M),PX(N,M),PY(N,M)

      CALL DPRJCT(0,XEYE,YEYE,ZEYE,RDUM1,RDUM2,RDUM3,RDUM4,RDUM5)
      DO 100 I = 1,M
          DO 100 J = 1,N
              CALL DPRJCT(1,XEYE,YEYE,ZEYE,X(J,I),Y(J,I),Z(J,I),
     +                    PX(J,I),PY(J,I))
  100 CONTINUE
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
c { dg-do run }
      program foo
      complex(kind=8)      z0, z1, z2

      z0 = dcmplx(0.,.5)
      z1 = 1./z0
      if (z1 .ne. dcmplx(0.,-2)) call abort

      z0 = 10.*z0
      if (z0 .ne. dcmplx(0.,5.)) call abort

      z2 = cmplx(1.,2.)
      z1 = z0/z2
      if (z1 .ne. dcmplx(2.,1.)) call abort

      z1 = z0*z2
      if (z1 .ne. dcmplx(-10.,5.)) call abort
      end

ext: f

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Ralph Hexter was named Provost and Executive Vice Chancellor on November 23, 2010.  Provost Hexter will also hold an appointment as a Professor of Classics and Comparative Literature at UCD.  He came to UCD from Hampshire College in Amherst, MA, where he had served as president since 2005.  A graduate of Harvard College, Corpus Christi College at the University of Oxford and Yale University, he began his academic career in Yales classics department. He later moved to the University of Colorado at Boulder, where he directed the graduate program in comparative literature.  He also taught at UC Berkeley beginning in 1995, serving as chair of comparative literature, and later served as both dean of the Division of Arts and Humanities and executive dean of the College of Letters and Science.

Provost Hexter is one of the few openly rainbow community gay campus leaders in the U.S.; he married his longtime partner, Manfred Kollmeier, in 2007.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C  Note: currently (NPB2.3), this routine not called
C
      subroutine wtime(tim)
      real*8 tim
      dimension tarray(2)
      call etime(tarray)
      tim = tarray(1)
      return
      end

max_stars_repo_name: sindrma/mg-parallel

max_issues_repo_name: sindrma/mg-parallel

max_forks_repo_name: sindrma/mg-parallel


text:
43, spouse of one Users/JeffreyNonken JeffreyNonken, mother of two small females, cocavemate of the two cats cats Darcy and Geezer Darcy and Geezer. The Nonken Quartet have all fallen for this lovely city.

Plants, food, knives, hapkido (at the Experimental College EC instructor Master Rick Schubert), http://en.wikipedia.org/wiki/Autism autism, and a few other oddities are either hobbies, collections or sources of sanity.

Weird Things I Can Do: Im working on the call of the Western ScrubJay.  Flag me down if you see me around, and Ill teach you how.  I think there should be more people walking around the streets of Davis going SquourEEE! SquourEEE!

I also talk to Crow Migration crows.  Sometimes they answer back.  Im pretty sure I have the basics down  I started learning by listening to the crows in Orange County, CA.  It seems they have a different dialect or language than the crows in Davis, which doesnt surprise me.  The ones in OC seem to have a simpler vocabulary, fewer variety of sounds; the places I lived havent been inhabited by humans for very long, the crow communities I learned from havent been there more than 30 to 50 years; only two or three generations in Crow Years.  The Davis Crows are more complex in their communication.

You can find me most afternoons behind the meat counter at the Davis Food CoOp.  Im one of the two females there; Im the short one with the braid.

silkensteel on LiveJournal LJ.

My hair is a lot longer now.  Its grown almost six inches since I moved here.  I wonder if, despite all the nasty chemicals in Davis water, theres something that speeds up hair growth? Oh, and that picture was taken in Irvine, CA.  I need someone to take a more recent picture, one that doesnt involve beer, sausages or classical musicians.

20050707 10:29:08 nbsp I just figured out how to add the Comments thingy! Users/BlancheNonken

20050810 01:11:52 nbsp Opossums heart Blanche Users/MichelleAccurso

20050810 21:42:46 nbsp I am loving this Wiki thing. Users/BlancheNonken

20050810 21:50:47 nbsp ha ha, i really do love opossums Users/MichelleAccurso

20050811 11:42:57 nbsp Hmm, I think Lost Pets could have a part between Featured Page and Mystery Picture Users/MichelleAccurso

20050811 11:46:45 nbsp I cant figured out how to do it exactly, try asking Users/PhilipNeustrom Philip about adding it on Users/MichelleAccurso

20050811 20:35:06 nbsp oh super, Ill be working 1015 to 315! good luck getting the lost pets page up, i bet hes somewhere in the arboretum, i dont know where he could go, though my dogs travelled pretty far when she ran off, wish i had some time to head out there and search for him, poor little guy. hope youre calling and checking the pound often, I suggest physically going there, when my other dgo ran away they told me he wasnt there and I went and he was. see you tomorrow! Users/MichelleAccurso

20050812 21:09:11 nbsp i think i saw you like, walk out the door, i was bagging and wandering around and not really working, like always. poor doggie, all lost and alone, but youre right, hell seduce some woman, does he have a collar and everything? Users/MichelleAccurso

20050814 21:30:48 nbsp finally! our paths have crossed! now ive officially met the whole family. Users/MichelleAccurso

20050815 14:50:15 nbsp Did you ever find Alex? Users/SummerSong

20050815 20:15:24 nbsp id like to hear your Western ScrubJay sound Users/MichelleAccurso

20050823 18:56:15 nbsp Hey, Baskin robbins is open til like 11pm, we should all do $1 scoop night and cider! woo! Users/MichelleAccurso

20050830 18:44:53 nbsp hmmmm might be stopping by then, ill IM tarz Users/MichelleAccurso

20050914 21:35:44 nbsp ok, were in the new place  and Alex never found us.  :(  Anyway, guests with advance warning are welcome. :) Users/BlancheNonken

20051104 19:25:58 nbsp we met at the last BBQ, we shared a wurst, i will never forget. Users/AndrewBanta

20051121 21:16:53 nbsp Hey Blanche, that chickenfat piecrust thing kind of grossed me out at first, but then I thought about it and it would be fun to try. Have you got a recipe? Users/EmmaCoats

20051211 23:58:15 nbsp Was nice seeing you at the coop, now I know to never mess with you.  Least you kick my ass. Users/AndrewBanta

20060109 21:09:59 nbsp I belive we may be spending Michelles birthday with you and Jeffery... Users/AndrewBanta

20060109 21:20:08 nbsp yes, were going to combine on birthday discounts and get cheap food :P see you thursday! Users/MichelleAccurso

20060131 00:01:55 nbsp Oh Blanche, I LOVE YOU.  So I tell you now. Users/AndrewBanta

20060218 21:47:30 nbsp Hooray!  Michelle is coming back tonight, from her hiatus in Hawaii.  We may be going to woodstocks, in which case I invite the Nonken family for her. Users/AndrewBanta

20060325 21:18:54 nbsp Andrew and Michelle need to flirt with me some more.  Bah, I feel lonely. Users/BlancheNonken

20060325 22:36:53 nbsp Ill have to swing by meat on wed. when I make my weekly treak into Davis... Users/AndrewBanta

20060617 00:24:25 nbsp Hola Blanche! I just found you on here. I will be away for the summer and I will miss you! Users/JessaRego

20070418 21:48:10 nbsp Weeeee! Hullo! Users/DavidFong


text:
!**********************************************************************
!    LICENSED MATERIALS - PROPERTY OF IBM                             *
!    "RESTRICTED MATERIALS OF IBM"                                    *
!                                                                     *
!    5765-422                                                         *
!    5765-C41                                                         *
!    (C) COPYRIGHT IBM CORP. 1995, 1997. ALL RIGHTS RESERVED.         *
!                                                                     *
!    U.S. GOVERNMENT USERS RESTRICTED RIGHTS - USE, DUPLICATION       *
!    OR DISCLOSURE RESTRICTED BY GSA ADP SCHEDULE CONTRACT WITH       *
!    IBM CORP.                                                        *
!**********************************************************************
      module init
      use pdata
      private  ! set everything private by default
      implicit none

!
!  This module initializes the blacs communication subsystem and
!  stores the row and column processor information. Most of the other
!  utility routines use the information stored by this routine.
!  while this imposes certain restraints on the user, e.g. only a
!  single blacs grid. It also simplifies the other utility routines
!  allowing a user to develop prototype message passing code more
!  quickly. The init module contains three routines as follows
!
!
!        subroutine initutils(status,nprs,npcs)
!
! This routine inializes blacs with the processor configuration
! nearest to a square grid if the optional arguments are not given.
! If nprs or npcs are provided, they are used to determine the
! processor shape.  If both are provided, their product must be equal
! to the number of processors in the current job.
!
!        subroutine arguments:
!
!          status :: integer, required, output
!                    returns initialization status
!                      0: properly initialized
!                      1: initialization failed
!
!          nprs ::   integer, optional, input
!                    number of row processor requested
!
!          npcs ::   integer, optional, input
!                    number of column processor requested
!
!      Subroutines called:
!        blacs_pinfo    from blacs library
!        blacs_gridinit from blacs library
!        blacs_gridexit from blacs library

!
!      subroutine exitutils
!
!      This subroutine exits from the current blacs grid and
!      marks the utilities as uninitialized
!

!      integer function p_context()
!
!      this function returns the blacs context for these utilities
!

      public initutils, p_context, exitutils

       contains
       subroutine initutils(status,nprs,npcs)
         integer,intent(out)             :: status
         integer, optional,intent(in) :: nprs, npcs
         integer                         :: i,j

         if(initialized) call exitutils
         call blacs_pinfo(mynode,numprocs)

         if(present(nprs) ) then
           if( present(npcs)) then
!
! both number of rows and columns are input
!
             if( nprs*npcs .gt. numprocs) then
               status = 1                      !don't have enough processors
               return
             endif
             nprows = nprs
             npcols = npcs

           else
!
! only number of rows are present
!
             if( nprs .gt. numprocs) then
               status = 1                      !don't have enough processors
               return
             endif
             nprows = nprs
             npcols = numprocs/nprs
           endif

         else
           if( present(npcs)) then
!
! only number of columns are present
!
             if( npcs .gt. numprocs) then
               status = 1                      !don't have enough processors
               return
             endif
             npcols = npcs
             nprows = numprocs/npcs
           else
!
! neither number of columns or rows present, pick the nearest to square
! grid
!
           nprows = int(sqrt(float(numprocs)))
           npcols = numprocs/nprows
           do while ( nprows*npcols .ne. numprocs)
             nprows = nprows - 1
             npcols = numprocs/nprows
           enddo
           endif

         endif

         call blacs_get(0,0,icontext)
         call blacs_gridinit(icontext,'r',nprows,npcols)
         call blacs_gridinfo(icontext,nprows,npcols,myrow,mycol)
         initialized = .true.
         status = 0
         return
       end subroutine initutils

       integer function p_context()

         p_context = -1
         if( initialized) p_context = icontext
       return
       end function p_context

       subroutine exitutils(close)
       integer, optional :: close

        if(initialized) then
         call blacs_gridexit(icontext)
         icontext = -1
         npcols = -1
         nprows = -1
         initialized = .false.
        endif
        if( present(close)) then
          call blacs_exit(close)
        endif
        return
       end subroutine exitutils

      end module init

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
      subroutine store_reorder(fl_item,
     &     label_out,label_in,
     &     iblk_out,iblk_in,
     &     tra_out,tra_in,
     &     sign_reo,occ_op0,
     &     from_to,occ_shift,nreo,nreo_i0,
     &     occ_opout,rst_opout,nj_out,
     &     occ_opin, rst_opin, nj_in,
     &     merge1,merge1inv,merge2,merge2inv,
     &     orb_info)

      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'def_contraction.h'
      include 'def_operator.h'
      include 'def_formula_item.h'
      include 'def_orbinf.h'

      integer, parameter ::
     &     ntest = 00

      type(formula_item), intent(inout) ::
     &     fl_item
      type(orbinf), intent(in) ::
     &     orb_info
      character(len=*), intent(in) ::
     &     label_out, label_in
      integer, intent(in) ::
     &     iblk_out, iblk_in,
     &     nreo, nj_out, nj_in, nreo_i0
      logical, intent(in) ::
     &     tra_out, tra_in
      integer, intent(in) ::
     &     sign_reo,
     &     from_to(2,nreo+nreo_i0),
     &     occ_shift(ngastp,2,nreo),
     &     occ_opin(ngastp,2,nj_in),
     &     rst_opin(2,orb_info%ngas,2,2,orb_info%nspin,nj_in),
     &     occ_opout(ngastp,2,nj_out),
     &     rst_opout(2,orb_info%ngas,2,2,orb_info%nspin,nj_out),
     &     occ_op0(ngastp,2,nj_in),
     &     merge1(*),merge1inv(*),
     &     merge2(*),merge2inv(*)

      integer ::
     &     ngas, nspin, lenmap

      integer, external ::
     &     len_merge_map

      if (ntest.ge.100)
     &     write(lulog,*) 'storing reordering'

      ngas = orb_info%ngas
      nspin = orb_info%nspin

      fl_item%reo%label_in  = label_in
      fl_item%reo%label_out = label_out
      fl_item%reo%nreo    = nreo
      fl_item%reo%nreo_i0   = nreo_i0
      fl_item%reo%ngas      = ngas
      fl_item%reo%nspin     = nspin
      fl_item%reo%nj_in   = nj_in
      fl_item%reo%nj_out  = nj_out
      fl_item%reo%iblk_in   = iblk_in
      fl_item%reo%iblk_out  = iblk_out
      fl_item%reo%tra_in   = tra_in
      fl_item%reo%tra_out  = tra_out
      fl_item%reo%sign      = sign_reo
      allocate(fl_item%reo%from_to(2,nreo+nreo_i0))
      fl_item%reo%from_to = from_to
      allocate(fl_item%reo%occ_opin(ngastp,2,nj_in))
      fl_item%reo%occ_opin = occ_opin
      allocate(fl_item%reo%rst_opin(2,ngas,2,2,nspin,nj_in))
      fl_item%reo%rst_opin = rst_opin
      allocate(fl_item%reo%occ_opout(ngastp,2,nj_out))
      fl_item%reo%occ_opout = occ_opout
      allocate(fl_item%reo%rst_opout(2,ngas,2,2,nspin,nj_out))
      fl_item%reo%rst_opout = rst_opout
      allocate(fl_item%reo%occ_shift(ngastp,2,nreo))
      fl_item%reo%occ_shift = occ_shift
      allocate(fl_item%reo%occ_op0(ngastp,2,nj_in))
      fl_item%reo%occ_op0 = occ_op0

      lenmap = len_merge_map(merge1,nj_in)
      allocate(fl_item%reo%merge_stp1(lenmap))
      fl_item%reo%merge_stp1 = merge1(1:lenmap)

      lenmap = len_merge_map(merge1inv,nj_in)
      allocate(fl_item%reo%merge_stp1inv(lenmap))
      fl_item%reo%merge_stp1inv = merge1inv(1:lenmap)

      lenmap = len_merge_map(merge2,nj_out)
      allocate(fl_item%reo%merge_stp2(lenmap))
      fl_item%reo%merge_stp2 = merge2(1:lenmap)

      lenmap = len_merge_map(merge2inv,nj_out)
      allocate(fl_item%reo%merge_stp2inv(lenmap))
      fl_item%reo%merge_stp2inv = merge2inv(1:lenmap)
c dbg
c        print *,'merge maps stored:'
c        print *,' 1:  ',fl_item%reo%merge_stp1
c        print *,' 1i: ',fl_item%reo%merge_stp1inv
c        print *,' 2:  ',fl_item%reo%merge_stp2
c        print *,' 2i: ',fl_item%reo%merge_stp2inv
c dbg

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


module io_read_wf_general_module

  use parallel_module, only: para, myrank_s,myrank_k,myrank_b,myrank_g
  use watch_module, only: watchb
  use wf_module, only: unk

  implicit none

  private
  public :: read_wf_general

contains

  subroutine read_wf_general( pinfo0,pinfo1,ML1,ML2,Ml3,MB0,MBZ0,MSP0,MB1,MBZ1,MSP1,type_wf )
    implicit none
    type(para),intent(in) :: pinfo0, pinfo1
    integer,intent(in) :: ML1,ML2,ML3,MB0,MBZ0,MSP0,MB1,MBZ1,MSP1
    integer,intent(in) :: type_wf
    real(8) :: ttmp(2),tt(2)
    logical :: disp
    integer :: p1,p2,p3,p4,p5,p6,irank_g,irank_g_old,g0,g1,g2,g3
    integer :: s0,s1,k0,k1,b0,b1,i1_0,i1_1,i2_0,i2_1,i3_0,i3_1
    integer :: i1_2,i1_3,i2_2,i2_3,i3_2,i3_3
    integer :: s,k,b,i1,i2,i3,irank_old,s2,s3,k2,k3,b2,b3,n,i,ng
    integer :: ig, ig_old, n11, n12
    character(13) :: filename
    character(7) :: basename
    character(5) :: crank
    integer,allocatable :: ir_new(:,:), id_new(:,:)
    integer,allocatable :: ir_old(:,:), id_old(:,:)
    real(8),allocatable :: dtmp(:)
    complex(8),allocatable :: utmp(:)
    integer,parameter :: u=3

    call write_border( 0, " read_wf_general(start)" )
    call check_disp_switch( disp, 0 )
    call watchb( ttmp, barrier='on' ); tt=0.0d0

    basename = 'wf.dat1'

    s2 = pinfo1%spin%id(myrank_s)+1
    s3 = s2 + pinfo1%spin%ir(myrank_s)-1
    k2 = pinfo1%bzsm%id(myrank_k)+1
    k3 = k2 + pinfo1%bzsm%ir(myrank_k)-1
    b2 = pinfo1%band%id(myrank_b)+1
    b3 = b2 + pinfo1%band%ir(myrank_b)-1
    g2 = pinfo1%grid%id(myrank_g)+1
    g3 = g2 + pinfo1%grid%ir(myrank_g)-1

    n = maxval( pinfo1%np(1:3) )
    allocate( ir_new(0:n-1,3) ); ir_new=0
    allocate( id_new(0:n-1,3) ); id_new=0
    do i1 = 0, ML1-1
      n = mod( i1, pinfo1%np(1) )
      ir_new(n,1) = ir_new(n,1) + 1
    end do
    do i2 = 0, ML2-1
      n = mod( i2, pinfo1%np(2) )
      ir_new(n,2) = ir_new(n,2) + 1
    end do
    do i3 = 0, ML3-1
      n = mod( i3, pinfo1%np(3) )
      ir_new(n,3) = ir_new(n,3) + 1
    end do
    do i = 1, 3
      do n = 0, pinfo1%np(i)-1
        id_new(n,i) = sum( ir_new(0:n,i) ) - ir_new(n,i)
      end do
    end do

    n = maxval( pinfo0%np(1:3) )
    allocate( ir_old(0:n-1,3) ); ir_old=0
    allocate( id_old(0:n-1,3) ); id_old=0
    do i1 = 0, ML1-1
      n = mod( i1, pinfo0%np(1) )
      ir_old(n,1) = ir_old(n,1) + 1
    end do
    do i2 = 0, ML2-1
      n = mod( i2, pinfo0%np(2) )
      ir_old(n,2) = ir_old(n,2) + 1
    end do
    do i3 = 0, ML3-1
      n = mod( i3, pinfo0%np(3) )
      ir_old(n,3) = ir_old(n,3) + 1
    end do
    do i = 1, 3
      do n = 0, pinfo0%np(i)-1
        id_old(n,i) = sum( ir_old(0:n,i) ) - ir_old(n,i)
      end do
    end do

    loop_i3: do i3 = 0, pinfo1%np(3)-1
             do i2 = 0, pinfo1%np(2)-1
             do i1 = 0, pinfo1%np(1)-1
               irank_g = i1 + i2*pinfo1%np(1) + i3*pinfo1%np(1)*pinfo1%np(2)
               if ( irank_g == myrank_g ) then
                 i1_2 = id_new(i1,1)
                 i1_3 = id_new(i1,1) + ir_new(i1,1) - 1
                 i2_2 = id_new(i2,2)
                 i2_3 = id_new(i2,2) + ir_new(i2,2) - 1
                 i3_2 = id_new(i3,3)
                 i3_3 = id_new(i3,3) + ir_new(i3,3) - 1
                 n11 = i1_3 - i1_2 + 1
                 n12 = n11*(i2_3-i2_2+1)
                 exit loop_i3
               end if
             end do
             end do
             end do loop_i3

    n=maxval(pinfo0%grid%ir)
    if ( type_wf == 1 ) then ! real-wf
      allocate( dtmp(n) ); dtmp=0.0d0
    else ! complex-wf
      allocate( utmp(n) ); utmp=(0.0d0,0.0d0)
    end if

    irank_old = -1

    do p6 = 0, pinfo0%np(6)-1
      s0 = pinfo0%spin%id(p6)+1
      s1 = s0 + pinfo0%spin%ir(p6)-1
    do p5 = 0, pinfo0%np(5)-1
      k0 = pinfo0%bzsm%id(p5)+1
      k1 = k0 + pinfo0%bzsm%ir(p5)-1
    do p4 = 0, pinfo0%np(4)-1
      b0 = pinfo0%band%id(p4)+1
      b1 = b0 + pinfo0%band%ir(p4)-1

      irank_g_old = -1

      do p3 = 0, pinfo0%np(3)-1
        i3_0 = id_old(p3,3)
        i3_1 = i3_0 + ir_old(p3,3) - 1
      do p2 = 0, pinfo0%np(2)-1
        i2_0 = id_old(p2,2)
        i2_1 = i2_0 + ir_old(p2,2) - 1
      do p1 = 0, pinfo0%np(1)-1
        i1_0 = id_old(p1,1)
        i1_1 = i1_0 + ir_old(p1,1) - 1

        irank_g_old = irank_g_old + 1
        g0 = pinfo0%grid%id(irank_g_old)+1
        g1 = g0 + pinfo0%grid%ir(irank_g_old)-1

        irank_old = irank_old + 1

        write(crank,'(i5.5)') irank_old
        filename = basename//'.'//crank

        if ( (s3<s0 .or. s1<s2) .or. (k3<k0 .or. k1<k2) .or. (b3<b0 .or. b1<b2) .or. &
              (i1_3<i1_0 .or. i1_1<i1_2) .or. (i2_3<i2_0 .or. i2_1<i2_2) .or. &
              (i3_3<i3_0 .or. i3_1<i3_2) ) cycle

        open(u,file=filename,status='old',form='unformatted')

        ng = g1 - g0 + 1
        do s = s0, s1
        do k = k0, k1
        do b = b0, b1

          if ( s<s2 .or. s3<s .or. k<k2 .or. k3<k .or. b<b2 .or. b3<b ) then
            read(u)
            cycle
          end if

          if ( type_wf == 1 ) then ! real-wf

            read(u) dtmp(1:ng)

            ! if ( s<s2 .or. s3<s .or. k<k2 .or. k3<k .or. b<b2 .or. b3<b ) cycle

            ig_old=0
            do i3 = i3_0, i3_1
            do i2 = i2_0, i2_1
            do i1 = i1_0, i1_1
              ig_old = ig_old + 1
              if ( i1<i1_2 .or. i1_3<i1 .or. i2<i2_2 .or. i2_3<i2 .or. i3<i3_2 .or. i3_3<i3 ) cycle
              ig = g2 + i1-i1_2 + n11*(i2-i2_2) + n12*(i3-i3_2)
              unk(ig,b,k,s)=dtmp(ig_old)
            end do
            end do
            end do

          else ! complex-wf

            read(u) utmp(1:ng)

            ! if ( s<s2 .or. s3<s .or. k<k2 .or. k3<k .or. b<b2 .or. b3<b ) cycle

            ig_old=0
            do i3 = i3_0, i3_1
            do i2 = i2_0, i2_1
            do i1 = i1_0, i1_1
              ig_old = ig_old + 1
              if ( i1<i1_2 .or. i1_3<i1 .or. i2<i2_2 .or. i2_3<i2 .or. i3<i3_2 .or. i3_3<i3 ) cycle
              ig = g2 + i1-i1_2 + n11*(i2-i2_2) + n12*(i3-i3_2)
              unk(ig,b,k,s)=utmp(ig_old)
            end do
            end do
            end do

          end if

        end do !b
        end do !k
        end do !s

        close(u)

      end do !p1
      end do !p2
      end do !p3

    end do !p4
    end do !p5
    end do !p6

    deallocate( id_old )
    deallocate( ir_old )
    deallocate( id_new )
    deallocate( ir_new )
    if ( allocated(dtmp) ) deallocate(dtmp)
    if ( allocated(utmp) ) deallocate(utmp)

    call watchb( ttmp, tt, barrier='on' )
    if ( disp ) then
      write(*,*) "time(read_wf_general)=",tt
    end if
    call write_border( 0, " read_wf_general(end)" )

  end subroutine read_wf_general

end module io_read_wf_general_module


text:
+incdir+../sv
../sv/spdif_tb_top.sv

max_stars_repo_name: zsassi/SPDIF-Traffic-Generator

max_issues_repo_name: zsassi/SPDIF-Traffic-Generator

max_forks_repo_name: zsassi/SPDIF-Traffic-Generator


text:
wiki:WikiPedia:Paintball Paintball has enjoyed a certain level of popularity in the greater Sacramento area throughout its history as a game/sport. Airsoft and Nerf are similar activities.  Paintball is a surprisingly fun activity for the weekend. Get a group of friends together and head on out. It can be quite pricey however.

 Where to Play

  Davis Paintball Center  all day rental package is $35, or only $15 if you have your own equipment, ladies play for free with their own equipment.  Users/AllisonEriksen recommends beginners not try to play there on Sundays, as that is when the professional teams go to practice, and they generally dont want to let inexperienced players into the games.

 Where to Get Equipment

  Big 5 Sporting Goods has a limited selection of markers (paintball guns) and facemasks.  Also their paintballs are quite a bit cheaper that what Davis Paintball Center charges on the field.
  Apparently ACE Hardware carries some paintball supplies... somewhere.  I only know this from an ad in the The California Aggie Aggie advertising a sale on paintball equipment.

 Teams
 The Paintball Club at UC Davis  is again active. All players are welcome whether inexperienced or experienced! Join us for an awesome time.