text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE PLTGPH(X,Y,NUM,XLAB,XUNIT,YLAB,YUNIT)
      REAL DEVCAP(23)
      REAL DEFOUT(7)
      COMMON /STATUS/DEVCAP,DEFOUT
      REAL DEVP(5)
      COMMON /DEVICE/DEVP
      REAL COLP(3)
      REAL PALETT(3,16)
      COMMON /COLOR/COLP,PALETT
      REAL TEXTP(40)
      COMMON /TEXT/TEXTP
      REAL VECTP(5)
      REAL XCUR
      REAL YCUR
      COMMON /VECTRC/VECTP,XCUR,YCUR
      INTEGER IDEX(200,2)
      INTEGER NVECT(200,2)
      REAL XSIZE(200,2)
      REAL YSIZE(200,2)
      REAL X0(2300,2)
      REAL Y0(2300,2)
      REAL X1(2300,2)
      REAL Y1(2300,2)
      COMMON /FONT/IDEX,NVECT,XSIZE,YSIZE,X0,Y0,X1,Y1
      REAL GRAPHP(100)
      COMMON /GRAPH/GRAPHP
      COMMON /MAPPAR/MAPP(11)
      REAL MAPP
      COMMON /STORAG/MEMORY(1000)
      CHARACTER*(*) XLAB,XUNIT,YLAB,YUNIT
      DIMENSION X(1),Y(1)

      IF (GRAPHP(21).EQ.1.) THEN
         CALL PLTNXY(X,Y,NUM,XLAB,XUNIT,YLAB,YUNIT)

      ELSE IF (GRAPHP(21).EQ.2.) THEN
         CALL PLTLGX(X,Y,NUM,XLAB,XUNIT,YLAB,YUNIT)

      ELSE IF (GRAPHP(21).EQ.3.) THEN
         CALL PLTLGY(X,Y,NUM,XLAB,XUNIT,YLAB,YUNIT)

      ELSE IF (GRAPHP(21).EQ.4.) THEN
         CALL PLTLXY(X,Y,NUM,XLAB,XUNIT,YLAB,YUNIT)
      END IF

      RETURN

      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
       subroutine narrow(xs,ys,xf,yf,ss,iccc)
C **      for    h-10/mt        es-1022
C...to draw an arrow at an end (or both ends) of a straight line.
C   xs,ys   x and y coordinates of the starting point of the dimension
C   xf,yf  x and y coordinates of the end point of a dimension line.
C   ss     length of the arrow
C   iccc   code showing the position to draw the arrow
C   iccc=1 to draw the arrow at the starting point side.
C   iccc=2to draw the arrow at the end point side.
C   iccc=3 to draw the arrow at the both ends of the straight line
      icod=iccc
      s=ss
        alpha=15.0*3.14159265/180.0
      if(icod.le.0) icod=3
      if(icod.gt.3) icod=3
      if(xf.eq.xs) goto 240
      theta=atan((yf-ys)/(xf-xs))
      goto 250
 240  if(yf-ys) 241,500,242
 241  theta=-3.14159265/2.0
      goto 250
 242  theta=3.14159265/2.0
  250  isw=1
       go to (301,302,310),icod
 301  isw=0
 310  check=xf-xs
      px=xs
      s=abs(s)
      py=ys
       if(check.lt.0.) s=-s
       go to 340
  302  continue
      call move(xs,ys,0)
      call move(xf,yf,1)
  320  isw=0
       check=xs-xf
      px=xf
       py=yf
      s=abs(s)
       if(check.le.0.) s=-s
 340  qx=px+s*cos(theta+alpha)
      qy=py+s*sin(theta+alpha)
      call move(px,py,0)
      call move(qx,qy,1)
      qx=px+s*cos(theta-alpha)
      qy=py+s*sin(theta-alpha)
      call move(px,py,0)
      call move(qx,qy,1)
      if(icod.eq.2)goto 342
      call move(xs,ys,0)
      call move(xf,yf,1)
      icod=2
 342  continue
       if(isw.ne.0) go to 320
 500  return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
! SPDX-Identifier: MIT
!
! Copyright 1999 J.C. Matthewman, P.J. Brown
! With contributions from W.I.F. David, J.B. Forsyth and J.H. Matthewman.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
!*****************************************************************************
!
      FUNCTION ALNINT(A,B,X,N)
!
! *** ALNINT by JCM 21 May 85 ***
!
!X
!C 9C
!H Performs linear interpolation, suitable for profile backgrounds etc.
!
!A On entry:
!A     A is a real array of dimension N which holds the arguments
!A     B is a real array of dimension N which holds the function values
!A     X is a real which holds the argument for which the function value
!A       is required
!A     N is the integer dimension of A and B
!A
!A ALINT on exit will hold the function value for argument X
!
      DIMENSION A(N), B(N)

      DO I = 2, N
        IF (A(I).GT.X) GOTO 2
      ENDDO
      I = N
    2 ALNINT = (B(I-1)*(A(I)-X)-B(I)*(A(I-1)-X))/(A(I)-A(I-1))

      END FUNCTION ALNINT
!
!*****************************************************************************
!
      SUBROUTINE C1MSCA(A,B,SCALE,NI,NJ)
!
! *** C1MSCA by PJB Apr 87 ***
!
!X
!C 12C
!H Multiplies every element of the REAL matrix A by the COMPLEX scalar SCALE.
!A On entry A is a REAL matrix of dimension NIxNJ
!A          SCALE is a COMPLEX scalar
!A On exit  B is a scaled COMPLEX matrix being SCALE times A
!N B may be the same as A
!N Does nothing if NI or NJ is zero
!
      COMPLEX B(1), SCALE

      DIMENSION A(1)
      NIJ = NI*NJ
      DO I = 1, NIJ
        B(I) = SCALE*A(I)
      ENDDO

      END SUBROUTINE C1MSCA
!
!*****************************************************************************
!
      SUBROUTINE CGMADD(A,B,C,NI,NJ)
!
! *** CGMADD by JCM 20 Nov 87 ***
!
!X
!C 12C
!H Sets COMPLEX matrix C = COMPLEX matrix A + COMPLEX matrix B.
!A On entry A and B are NI by NJ COMPLEX matrices
!A On exit C is set to the sum of A and B
!
      COMPLEX A(NI,NJ), B(NI,NJ), C(NI,NJ)

      DO I = 1, NI
        DO J = 1, NJ
          C(I,J) = A(I,J) + B(I,J)
        ENDDO
      ENDDO

      END SUBROUTINE CGMADD
!
!*****************************************************************************
!
      SUBROUTINE CGMEQ(A,B,N,M)
!
! *** CGMEQ by JCM ***
!
!X
!C 12C
!H Sets COMPLEX matrix B equal to COMPLEX matrix A.
!A On entry A is a COMPLEX matrix of DIMENSION N,M
!A On exit B, a COMPLEX matrix of the same total size, has been set equal to A
!N Would not give error even if N or M were zero.
      COMPLEX A(1), B(1)

      NE = N*M
      DO I = 1, NE
        B(I) = A(I)
      ENDDO

      END SUBROUTINE CGMEQ
!
!*****************************************************************************
!
      SUBROUTINE CGMSCA(A,B,SCALE,NI,NJ)
!
! *** CGMSCA by PJB Apr 87 ***
!
!X
!C 12C
!H Multiplies every element of a COMPLEX matrix by a COMPLEX scale.
!A On entry, A holds a COMPLEX NI by NJ matrix
!A           SCALE holds a COMPLEX scalar scaling number
!A On exit,  B holds a COMPLEX NI by NJ matrix being A times SCALE
!N A may be the same as B
!N Note the existence of CMRSCA which works with a REAL scale
!
      COMPLEX A(1), B(1), SCALE

      NIJ = NI*NJ
      DO I = 1, NIJ
        B(I) = SCALE*A(I)
      ENDDO

      END SUBROUTINE CGMSCA
!
!*****************************************************************************
!
      SUBROUTINE CGMZER(A,NI,NJ)
!
! *** CGMZER by JCM 23 Oct 87 ***
!
!X
!C 12C
!H Clears to zero a complex matrix.
!A On entry A is an NI by NJ COMPLEX matrix
!A On exit A is cleared to COMPLEX zero
!
      COMPLEX A(1)

      NIJ = NI*NJ
      DO I = 1, NIJ
        A(I) = CMPLX(0.,0.)
      ENDDO

      END SUBROUTINE CGMZER
!
!*****************************************************************************
!
      SUBROUTINE CMCONJ(A,B,I,J)
!
! *** CMCONJ by PJB Nov 89 ***
!
!X
!C 12C
!H Gives the conjugate of a COMPLEX matrix.
!A On entry A is a complex I by J matrix
!A On exit the I by J complex matrix B contains conjugate complex of A
!
      COMPLEX A(I,J), B(I,J)

      DO II = 1, I
        DO JJ = 1, J
          B(II,JJ) = CONJG(A(II,JJ))
        ENDDO
      ENDDO

      END SUBROUTINE CMCONJ
!
!*****************************************************************************
!
      SUBROUTINE CMRSCA(A,B,SCALE,NI,NJ)
!
! *** CMRSCA by JCM 23 Aug 88 ***
!
!X
!C 12C
!H Multiplies every element of a COMPLEX matrix by a REAL scale.
!A On entry, A holds a COMPLEX NI by NJ matrix
!A           SCALE holds a REAL scalar scaling number
!A On exit,  B holds a COMPLEX NI by NJ matrix being A times SCALE
!N A may be the same as B
!N Note also the existence of CGMSCA which works entirely in COMPLEX.
!
      COMPLEX A(1), B(1)

      NIJ = NI*NJ
      DO I = 1, NIJ
        B(I) = SCALE*A(I)
      ENDDO

      END SUBROUTINE CMRSCA
!
!*****************************************************************************
!
      FUNCTION DEGREE(X)
!
! *** DEGREE by JCM ***
!
!X
!C 10C
!H Converts from radians to degrees.
!A On entry X is an angle in radians
!A On exit DEGREE is the angle in degrees
!
!N The function RADIAN(X) does the degrees to radians conversion.
!N In routines where time matters, it is quicker to declare the COMMON
!N /CONSTA/ and multiply X by DEG.
      REAL            PI, RAD, DEG, TWOPI, FOURPI, PIBY2, ALOG2, SQL2X8, VALMUB
      COMMON /CONSTA/ PI, RAD, DEG, TWOPI, FOURPI, PIBY2, ALOG2, SQL2X8, VALMUB

      DEGREE = DEG*X

      END FUNCTION DEGREE
!
!*****************************************************************************
!
      FUNCTION DETER3(A)
!
! *** DETER3 by JCM 28 Jun 83 ***
!
!X
!C 12C
!H Forms the determinant of a 3 x 3 matrix.
!A On entry A is a 3x3 real matrix.
!A On exit DETER3 holds its determinant.
!
      DIMENSION A(3,3)

      D = 0.0
      J = 2
      K = 3
      DO I = 1, 3
        D = D + A(1,I)*(A(2,J)*A(3,K)-A(2,K)*A(3,J))
        J = K
        K = I
      ENDDO
      DETER3 = D

      END FUNCTION DETER3
!
!*****************************************************************************
!
      FUNCTION ERFNC(X)
!
! *** ERFNC by WIFD 22 Aug 85 ***
!
!X
!C 9C
!H Calculates the error function accurate to 3E-7, for + and - X.
!A On entry X is the argument at which the function is required.
!A On exit ERFNC holds the function.
!D See Abramovitz and Stegun p.299
!
      DATA A1, A2, A3/0.0705230784, 0.0422820123, 0.0092705272/
      DATA A4, A5, A6/0.0001520143, 0.0002765672, 0.0000430638/
!  THESE COEFFICIENTS WERE TAKEN FROM ABRAMOVITZ AND STEGUN P.299
!  AND GIVE ERFNC ACCURATE TO 3E-07 FOR BOTH +VE AND -VE X.

      Z = ABS(X)
      ZZ = Z*Z
      ZZZ = Z*Z*Z
      E = 1.0 + A1*Z + A2*ZZ + A3*ZZZ + A4*ZZ*ZZ + A5*ZZ*ZZZ + A6*ZZZ*ZZZ
      IF (E.LT.10.0) GOTO 1
      E = 0.0
      GOTO 2
    1 E = 1.0/E**16
      IF (E.LE.0.0) E = 0.0
    2 IF (X.LT.0.0) E = 2.0 - E
      ERFNC = E

      END FUNCTION ERFNC
!
!*****************************************************************************
!
      SUBROUTINE FRAC3(VEC)
!
! *** FRAC3 by JCM 11 Nov 83 ***
!
!X
!C 1C
!H Makes all 3 elements of a vector fractional
!A On entry VEC is a 1x3 real vector.
!A On exit the elements of VEC have each been put into the range 0 =< X < 1
!
      DIMENSION VEC(3)

      DO I = 1, 3
        CALL FRACT(VEC(I),A,J)
      ENDDO

      END SUBROUTINE FRAC3
!
!*****************************************************************************
!
      SUBROUTINE FRACT(X,Y,N)
!
! *** FRACT by JCM ***
!
!X
!C 9C
!H Forms the fractional part of a real number.
!A On entry X is a real number
!A On exit  X is in the range 0=< X <1
!A          Y is set so that X+Y=original X
!A          N= 0 if X was unchanged
!A           = 1 if X was >= 1
!A           =-1 if X was < 0
!
      Y = 0.
      N = 1
      IF (X.GE.1.) GOTO 1
      N = 0
      IF (X.GE.0.) RETURN
! OUT IN USUAL CASE OF X ALREADY BEING +VE FRACTION
      N = -1
    1 Y = AINT(X)
      IF (N.EQ.-1) Y = Y - 1.0
      X = X - Y

      END SUBROUTINE FRACT
!
!*****************************************************************************
!
      SUBROUTINE FT01A(IT,INV,TR,TI)
!
! *** FT01A updated by JCM FROM HARWELL ROUTINE 9 Sep 91 ***
!
!X
!C 9C
!H Modification of Harwell Fast Fourier Transform.
!
      DIMENSION TR(*), TI(*)

      REAL            PI, RAD, DEG, TWOPI, FOURPI, PIBY2, ALOG2, SQL2X8, VALMUB
      COMMON /CONSTA/ PI, RAD, DEG, TWOPI, FOURPI, PIBY2, ALOG2, SQL2X8, VALMUB

      REAL            UR,     UI
      COMMON /FFTDA / UR(15), UI(15)

    2 UM = 1.
      IF (INV.EQ.1) UM = -1.
      IO = 2
      DO I = 2, 16
        IO = IO + IO
        IF (IO-IT) 3, 4, 99
    3 ENDDO
! ERROR EXIT - IF NOT A POWER OF 2, OR TOO BIG:
   99 INV = -1
      CALL DebugErrorMessage('NOT A POWER OF 2, OR TOO BIG in FT01A')
      RETURN
    4 IO = I
      II = IO
      I1 = IT/2
      I3 = 1
   10 K = 0
      I2 = I1 + I1
   11 WR = 1.0
      WI = 0.0
      KK = K
      JO = IO
   12 IF (KK.EQ.0) GOTO 13
   14 JO = JO - 1
      KK1 = KK
      KK = KK/2
      IF (KK1.EQ.2*KK) GOTO 14
      WS = WR*UR(JO) - WI*UI(JO)
      WI = WR*UI(JO) + WI*UR(JO)
      WR = WS
      GOTO 12
   13 WI = WI*UM
      J = 0
    9 L = J*I2 + K
      L1 = L + I1
      ZR = TR(L+1) + TR(L1+1)
      ZI = TI(L+1) + TI(L1+1)
      Z = WR*(TR(L+1)-TR(L1+1)) - WI*(TI(L+1)-TI(L1+1))
      TI(L1+1) = WR*(TI(L+1)-TI(L1+1)) + WI*(TR(L+1)-TR(L1+1))
      TR(L+1) = ZR
      TR(L1+1) = Z
      TI(L+1) = ZI
      J = J + 1
      IF (J.LT.I3) GOTO 9
      K = K + 1
      IF (K.LT.I1) GOTO 11
      I3 = I3 + I3
      IO = IO - 1
      I1 = I1/2
      IF (I1.GT.0) GOTO 10
      J = 1
      UM = 1.
      IF (INV.EQ.1) UM = 1./FLOAT(IT)
    7 K = 0
      J1 = J
      DO I = 1, II
        J2 = J1/2
        K = 2*(K-J2) + J1
        J1 = J2
      ENDDO
      IF (K.GE.J) THEN
        IF (K.EQ.J) THEN
          TR(J+1) = TR(J+1)*UM
          TI(J+1) = TI(J+1)*UM
        ELSE
          ZR = TR(J+1)
          ZI = TI(J+1)
          TR(J+1) = TR(K+1)*UM
          TI(J+1) = TI(K+1)*UM
          TR(K+1) = ZR*UM
          TI(K+1) = ZI*UM
        ENDIF
      ENDIF
      J = J + 1
      IF (J.LT.IT-1) GOTO 7
      TR(1) = TR(1)*UM
      TI(1) = TI(1)*UM
      TR(IT) = TR(IT)*UM
      TI(IT) = TI(IT)*UM
  100 RETURN
      END SUBROUTINE FT01A
!
!*****************************************************************************
!
      SUBROUTINE GMADD(A,B,C,NI,NJ)
!
! *** GMADD by JCM ***
!
!X
!C 12C
!H Sets matrix C = matrix A plus matrix B.
!A On entry A and B are real matrices of dimension (NIxNJ)
!A On exit  C is a real matrix which is their sum.
!N NI and NJ must be at least 1
!
      DIMENSION A(NI,NJ), B(NI,NJ), C(NI,NJ)

      DO I = 1, NI
        DO J = 1, NJ
          C(I,J) = A(I,J) + B(I,J)
        ENDDO
      ENDDO

      END SUBROUTINE GMADD
!
!*****************************************************************************
!
      SUBROUTINE GMEQ(A,B,NI,NJ)
!
! *** GMEQ by PJB/JCM 28 Jun 83 ***
!
!X
!C 12C
!H Sets matrix B = matrix A.
!A On entry A is a real matrix of dimension NIxNJ
!A On exit  B is a real matrix equal to A
!N NI and NJ must be at least 1
!
      DIMENSION A(NI,NJ), B(NI,NJ)
      DO I = 1, NI
        DO J = 1, NJ
          B(I,J) = A(I,J)
        ENDDO
      ENDDO

      END SUBROUTINE GMEQ

!
!*****************************************************************************
!
      SUBROUTINE INGMEQ(A,B,NI,NJ)
!
! *** GMEQ by PJB/JCM 28 Jun 83 ***
!
!X
!C 12C
!H Sets matrix B = matrix A.
!A On entry A is a real matrix of dimension NIxNJ
!A On exit  B is a real matrix equal to A
!N NI and NJ must be at least 1
!
      INTEGER A,B
      DIMENSION A(NI,NJ), B(NI,NJ)
      DO I = 1, NI
        DO J = 1, NJ
          B(I,J) = A(I,J)
        ENDDO
      ENDDO

      END SUBROUTINE INGMEQ
!
!*****************************************************************************
!
      SUBROUTINE GMREV(A,B,NI,NJ)
!
! *** GMREV by PJB/JCM 28 Jun 83 ***
!
!X
!C 12C
!H Reverses the signs of the elements of an NI X NJ matrix.
!A On entry A is a real matrix of dimension NIxNJ
!A On exit  B is a real matrix holding -A
!N A and B may be the same matrix
!
      DIMENSION A(NI,NJ), B(NI,NJ)
      DO I = 1, NI
        DO J = 1, NJ
          B(I,J) = -A(I,J)
        ENDDO
      ENDDO

      END SUBROUTINE GMREV

!
      SUBROUTINE CPGMREV(A,B,NI,NJ)
!
! *** GMREV by PJB/JCM 28 Jun 83 ***
!
!X
!C 12C
!H Reverses the signs of the elements of an NI X NJ matrix.
!A On entry A is a real matrix of dimension NIxNJ
!A On exit  B is a real matrix holding -A
!N A and B may be the same matrix
!
      COMPLEX A,B
      DIMENSION A(NI,NJ), B(NI,NJ)
      DO I = 1, NI
        DO J = 1, NJ
          B(I,J) = -A(I,J)
        ENDDO
      ENDDO

      END SUBROUTINE CPGMREV
!
!*****************************************************************************
!
      LOGICAL FUNCTION GMSAME(A,B,N,TOLER)
!
! *** GMSAME by JCM 22 Oct 86 ***
!
!X
!C 11C
!H Tells whether one vector is the same as another, to a given tolerance.
!
!A A on entry is an N-sized array, to be compared with:
!A B, also an N-sized array (A and B may of course be parts of larger arrays)
!A N on entry is the number of elements of A and B to compare
!A TOLER on entry is the number within which all the elements of A and B
!A      must agree
!A
!A GMSAME will be set .TRUE. if all elements of A and B agree within TOLER,
!A      and .FALSE. otherwise
!
      DIMENSION A(N), B(N)

      GMSAME = .FALSE.
      DO I = 1, N
        IF (ABS(A(I)-B(I)).GT.TOLER) RETURN
      ENDDO
      GMSAME = .TRUE.

      END FUNCTION GMSAME
!
!*****************************************************************************
!
      SUBROUTINE GMSCA(A,B,SCALE,NI,NJ)
!
! *** GMSCA by JCM 22 Nov 84 ***
!
!X
!C 12C
!H Multiplies every element of the matrix A by the scalar SCALE.
!A On entry A is a real matrix of dimension NIxNJ
!A          SCALE is the required multiplying scale
!A On exit  B is a real matrix whose elements are those of A times SCALE.
!N A and B may be the same matrix
!
      DIMENSION A(1), B(1)

      NIJ = NI*NJ
      DO I = 1, NIJ
        B(I) = SCALE*A(I)
      ENDDO

      END SUBROUTINE GMSCA
!
!*****************************************************************************
!
      SUBROUTINE GMSUB(A,B,C,NI,NJ)
!
! *** GMSUB by PJB/JCM 28 Jun 83 ***
!
!X
!C 12C
!H Sets matrix C = matrix A minus matrix B.
!A On entry A and B are real matrices of dimension NIxNJ
!A On exit  C is a real matrix whose elements are those of A-B
!
      DIMENSION A(1), B(1), C(1)
      NIJ = NI*NJ
      DO I = 1, NIJ
        C(I) = A(I) - B(I)
      ENDDO

      END SUBROUTINE GMSUB
!
!*****************************************************************************
!
      SUBROUTINE GMTRAN(A,B,JJ,II)
!
! *** GMTRAN by PJB ***
!
!X
!C 12C
!H Transposes a JJxII matrix A into B.
!A On entry A is a real matrix of dimension (JJxII)
!A On exit  B is a real matrix of dimension (IIxJJ) which is the transpose of A
!
      DIMENSION A(JJ,II), B(II,JJ)

      DO J = 1, JJ
        DO I = 1, II
          B(I,J) = A(J,I)
        ENDDO
      ENDDO

      END SUBROUTINE GMTRAN
!
!*****************************************************************************
!
      SUBROUTINE GMUNI(A,NI)
!
! *** GMUNI by JCM 7 Jul 83 ***
!
!X
!C 12C
!H Writes a unit matrix into the square matrix A.
!A On entry NI is the dimension of the required matrix
!A On exit  A is a square real matrix of dimension NIxNI, set to a unit matrix
!
      DIMENSION A(NI,NI)

      DO I = 1, NI
        A(I,I) = 1.
        I1 = I + 1
        IF (I1.GT.NI) RETURN
        DO J = I1, NI
          A(I,J) = 0.
          A(J,I) = 0.
        ENDDO
      ENDDO

      END SUBROUTINE GMUNI
!
!*****************************************************************************
!
      SUBROUTINE GMZER(A,NI,NJ)
!
! *** GMZER by JCM 7 Jul 83 ***
!
!X
!C 12C
!H Clears to zero the matrix A.
!A On entry NI and NJ specify A to be of dimension NIxNJ
!A On exit  A is a real matrix of dimension NIxNJ all cleared to zero.
!
      DIMENSION A(1)

      NIJ = NI*NJ
      DO I = 1, NIJ
        A(I) = 0.0
      ENDDO

      END SUBROUTINE GMZER
!
!*****************************************************************************
!
      SUBROUTINE JGMEQ(JA,JB,NI,NJ)
!
! *** JGMEQ by JCM  Jun 88 ***
!
!X
!C 12C
!H Equates an integer matrix to a given integer matrix.
!A On entry JA is an integer matrix of dimension NIxNJ
!A On exit  JB = integer matrix being a copy of JA
!
      DIMENSION JA(NI,NJ), JB(NI,NJ)

      DO I = 1, NI
        DO J = 1, NJ
          JB(I,J) = JA(I,J)
        ENDDO
      ENDDO

      END SUBROUTINE JGMEQ
!
!*****************************************************************************
!
      SUBROUTINE JGMZER(JA,NI,NJ)
!
! *** JGMZER by JCM Jun 88 ***
!
!X
!C 12C
!H Clears an integer matrix to zero.
!A On entry NI and NJ are the dimensions of the required matrix
!A On exit  the integer matrix JA, of dimensions NIxNJ is cleared to zero.
!
      INTEGER, INTENT (  OUT) :: JA(*)
      INTEGER, INTENT (IN   ) :: NI, NJ

      INTEGER I, NIJ

      NIJ = NI*NJ
      DO I = 1, NIJ
        JA(I) = 0
      ENDDO

      END SUBROUTINE JGMZER
!
!*****************************************************************************
!
      FUNCTION RADIAN(X)
!
! *** RADIAN by JCM ***
!
!X
!C 10C
!H Converts from degrees to radians.
!A On entry X is the value of an angle in degrees
!D On exit  X has been converted to radians.
!N The function DEGREE(X) does the radians to degrees conversion
!N In routines where time matters, it is quicker to declare the COMMON
!N /CONSTA/ and multiply X by RAD.
!
      REAL            PI, RAD, DEG, TWOPI, FOURPI, PIBY2, ALOG2, SQL2X8, VALMUB
      COMMON /CONSTA/ PI, RAD, DEG, TWOPI, FOURPI, PIBY2, ALOG2, SQL2X8, VALMUB

      RADIAN = RAD*X

      END FUNCTION RADIAN
!
!*****************************************************************************
!
      SUBROUTINE RCMPRD(A,B,C,I,J,K)
!
! *** RCMPRD by PJB Jun 87 ***
!
!X
!C 12C
!H Performs the multiplication of a COMPLEX by a REAL matrix.
!A On entry A is a real IxJ matrix,
!A          B is a complex JxK matrix,
!A On exit  C is returned as the complex IxK matrix A*B
!
      DIMENSION A(1)
      COMPLEX B(1), C(1)

      DO II = 1, I
        IK = II
        JK = 1
        DO KK = 1, K
          IJ = II
          C(IK) = 0.
          DO JJ = 1, J
            C(IK) = C(IK) + A(IJ)*B(JK)
            IJ = IJ + I
            JK = JK + 1
          ENDDO
          IK = IK + I
        ENDDO
      ENDDO

      END SUBROUTINE RCMPRD
!
!*****************************************************************************
!
      FUNCTION RSCALP(A,B)
!
! *** RSCALP by PJB Jun 87 ***
!
!X
!C 12C
!H Forms the scalar product of two COMPLEX vectors.
!A On entry A and B are complex vectors of dimension 3
!D The function value is the real product (A.Conjg(B) + Conjg(A).B)/2
!
      COMPLEX A(3), B(3), P

      P = CMPLX(0.,0.)
      DO I = 1, 3
        P = P + A(I)*CONJG(B(I))
      ENDDO
      RSCALP = REAL(P)

      END FUNCTION RSCALP
!
!*****************************************************************************
!
      FUNCTION SCALPR(A,B)
!
! *** SCALPR by JCM ***
!
!X
!C 12C
!H Forms the scalar product of two orthogonal vectors (or of one from
!H real space with one from reciprocal space).
!A On entry A and B are vectors of dimension 3
!D The function is set to the scalar product of A and B (not involving
!D  cell parameters) e.g. with one vector from real & one from reciprocal
!D  space
!N The function SCLPRD gives the scalar product using the cell parameters.
!
      DIMENSION A(3), B(3), C(1)

      CALL MultiplyMatrices(A,B,C,1,3,1)
      SCALPR = C(1)

      END FUNCTION SCALPR
!
!*****************************************************************************
!
      SUBROUTINE SINCOS(SN,CS,L1)
!
! *** SINCOS updated by JCM 1 Jul 86 ***
!
!X
!C 10C
!H Calculates sin from cos or vice-versa.
!A On entry SN is a sine or cosine of an angle
!A On exit  CS is the cosine or sine respectively of the same angle
!A          L1 is a character variable which is used to identify the
!A             calling routine if the absolute value of SN is greater
!A             than 1.
!N If ABS(SN)>1 the routine writes an error message and STOPS
!
      CHARACTER*(*) L1
      INTEGER         LPT, LUNI
      COMMON /IOUNIT/ LPT, LUNI

      IF (ABS(SN)-10.E-5.GT.1.) THEN
        WRITE (LPT,3000) SN, L1
 3000   FORMAT (/' Sin or Cos value',E12.4,' greater than unity',' - called from ',A6)
        CALL BMBOUT
        CALL DebugErrorMessage('Error in SINCOS() in Mvec')
        RETURN
      ENDIF
!
!>> JCC There can be an error raised here. SN can be passed
!>> As a NaN, due to a corruption in the CELL array. Please check
!>>
! NOW DETECT CASE WHERE SN IS ONLY JUST MORE THAN 1:
      IF (ABS(SN).GE.1.0) THEN
        SN = SIGN(1.0,SN)
        CS = 0.
      ELSE
! USUAL CASE:
        CS = SQRT(1.-SN*SN)
      ENDIF

      END SUBROUTINE SINCOS
!
!*****************************************************************************
!
      SUBROUTINE SPLINE(N,X,F,D)
!
! *** SPLINE from HARWELL TB04A 21 May 85 ***
!
!X
!C 9A
!H Sets up a cubic spline to fit a curve.
!A On entry N is the number of points on a curve
!A          F(I) are the N function values at the N points X(I)
!A          D(I) is set to the derivative if the spline at the point I
!A              for all N points.
!N Uses working space in COMMON /SCRAT/
!N For further details see the description of the Harwell Library
!
      DIMENSION X(N), F(N), D(N)
      INTEGER         LPT, LUNI
      COMMON /IOUNIT/ LPT, LUNI
      COMMON /SCRAT / A(2000)
!
! F(I) ARE THE FUNCTION VALUES AT THE POINTS X(I) FOR I=1,N AND
! THE SPLINE DERIVATIVES D(I) ARE FOUND.  THE DIMENSION OF A MUST
! NOT BE LESS THAN 3*N.
!
      DO I = 2, N
        IF (X(I)-X(I-1)) 1, 1, 5
    1   WRITE (LPT,3000) I
 3000   FORMAT (/'  ERROR ** in SPLINE - X',I3,' out of order')
        A(1) = 1.
        GOTO 100
    5 ENDDO
      DO I = 1, N
        J = 2
        IF (I-1) 6, 10, 6
    6   J = N - 1
        IF (I.EQ.N) GOTO 10
        H1 = 1./(X(I)-X(I-1))
        H2 = 1./(X(I+1)-X(I))
        A(3*I-2) = H1
        A(3*I-1) = 2.*(H1+H2)
        A(3*I) = H2
        D(I) = 3.0*(F(I+1)*H2*H2+F(I)*(H1*H1-H2*H2)-F(I-1)*H1*H1)
        GOTO 2
   10   H1 = 1./(X(J)-X(J-1))
        H2 = 1./(X(J+1)-X(J))
        A(3*I-2) = H1*H1
        A(3*I-1) = H1*H1 - H2*H2
        A(3*I) = -H2*H2
        D(I) = 2.*(F(J)*(H2*H2*H2+H1*H1*H1)-F(J+1)*H2*H2*H2-F(J-1)*H1*H1*H1)
    2 ENDDO
      P = A(4)/A(1)
      A(5) = A(5) - P*A(2)
      A(6) = A(6) - P*A(3)
      D(2) = D(2) - P*D(1)
      DO I = 3, N
        K = 3*I - 4
        P = A(K+2)/A(K)
        A(K+3) = A(K+3) - P*A(K+1)
        D(I) = D(I) - P*D(I-1)
        IF (I.NE.N-1) GOTO 4
        P = A(K+5)/A(K)
        A(K+5) = A(K+6) - P*A(K+1)
        A(K+6) = A(K+7)
        D(N) = D(N) - P*D(N-2)
    4 ENDDO
      D(N) = D(N)/A(3*N-1)
      DO I = 3, N
        J = N + 2 - I
        D(J) = (D(J)-A(3*J)*D(J+1))/A(3*J-1)
      ENDDO
      D(1) = (D(1)-D(2)*A(2)-D(3)*A(3))/A(1)
      A(1) = 0.
  100 RETURN
      END SUBROUTINE SPLINE
!
!*****************************************************************************
!
      FUNCTION SPLINT(IX,N,U,S,D,X)
!
! *** SPLINT from HARWELL TG01B 21 May 85 ***
!
!X
!C 9C
!H Evaluates a cubic spline given spline values and first derivative
!H values at the given knots.
!
!
!A On entry IX allows the caller to take advantage of spline parameters
!A            set on a previous call in cases when an X point follows
!A            the previous X point.
!A      If IX < 0 the whole range is searched for knot interval;
!A      If IX > 0 it is assumed that X is greater than the X of the
!A      previous call and search started from there.
!A           N=The number of knots.
!A           U holds the knots.
!A           S holds the spline values.
!A           D holds the first derivative values of the spline at the
!A             knots.
!A           X=the point at which the spline value is required.
!
!P A call of SPLINE must have set up the knots in D.
!
!N The spline value is defined as zero outside the knot range, which
!N is extended by a rounding error for the purpose.
!
!     ALLOWABLE ROUNDING ERROR ON POINTS AT EXTREMES OF KNOT RANGE
!     IS 2**IEPS*MAX(\U(1)\,\U(N)\).
      DIMENSION U(N), S(N), D(N)
      DATA IFLG, IEPS/0, -19/
!       TEST WETHER POINT IN RANGE.
      IF (X.LT.U(1)) GOTO 3
      IF (X.GT.U(N)) GOTO 4
!       JUMP IF KNOT INTERVAL REQUIRES RANDOM SEARCH.
      IF (IX.LT.0 .OR. IFLG.EQ.0) GOTO 12
!       JUMP IF KNOT INTERVAL SAME AS LAST TIME.
      IF (X.LE.U(J+1)) GOTO 8
!       LOOP TILL INTERVAL FOUND.
    1 J = J + 1
   11 IF (X.GT.U(J+1)) GOTO 1
      GOTO 7
!       ESTIMATE KNOT INTERVAL BY ASSUMING EQUALLY SPACED KNOTS.
   12 J = ABS(X-U(1))/(U(N)-U(1))*(N-1) + 1
!       ENSURE CASE X=U(N) GIVES J=N-1.
      J = MIN0(J,N-1)
!       INDICATE THAT KNOT INTERVAL INSIDE RANGE HAS BEEN USED.
      IFLG = 1
!       SEARCH FOR KNOT INTERVAL CONTAINING X.
      IF (X.GE.U(J)) GOTO 11
    2 J = J - 1
      IF (X.LT.U(J)) GOTO 2
!       CALCULATE SPLINE PARAMETERS FOR JTH INTERVAL.
    7 H = U(J+1) - U(J)
      Q1 = H*D(J)
      Q2 = H*D(J+1)
      SS = S(J+1) - S(J)
      B = 3.0*SS - 2.0*Q1 - Q2
      A = Q1 + Q2 - 2.0*SS
!       CALCULATE SPLINE VALUE.
    8 Z = (X-U(J))/H
      SPLINT = ((A*Z+B)*Z+Q1)*Z + S(J)
      GOTO 100
!       TEST IF X WITHIN ROUNDING ERROR OF U(1).
    3 IF (X.LE.U(1)-2.0**IEPS*AMAX1(ABS(U(1)),ABS(U(N)))) GOTO 101
      J = 1
      GOTO 7
!       TEST IF X WITHIN ROUNDING ERROR OF U(N).
    4 IF (X.GE.U(N)+2.0**IEPS*AMAX1(ABS(U(1)),ABS(U(N)))) GOTO 101
      J = N - 1
      GOTO 7
! POINTS OUTSIDE RANGE:
  101 IFLG = 0
      SPLINT = 0.0
  100 RETURN
      END FUNCTION SPLINT
!
!*****************************************************************************
!
      SUBROUTINE TB02A(A,F,X,Y,N)
!
! *** TB02A updated by JCM from HARWELL ***
!
!X
!C 9C
!H Interpolates in non-equal interval table.
!A On entry real arrays A and F hold the arguments and function values
!A          X holds the required argument
!A          N is the dimension of A and F, and so is the number of items
!A On exit  Y holds the interpolated function value
!
      DIMENSION A(N), F(N)

      Y = 0.0
      DO L = 1, N
        PL = 1.0
        DO J = 1, N
          IF (L.EQ.J) GOTO 3
          PL = (X-A(J))*PL/(A(L)-A(J))
    3   ENDDO
        Y = Y + PL*F(L)
      ENDDO

      END SUBROUTINE TB02A
!
!*****************************************************************************
!
      SUBROUTINE TQLI(D,E,N,NP,Z,BOTH)
!
! *** TQLI by JCM from "NUMERICAL RECIPES"  Dec 89 ***
!
!X
!C 9C
!H Performs the QL algorithm for eigenvalues and vectors of a real
!H symmetric matrix previously reduced to tridiagonal form.
!A On entry D is a real vector of length NP whose first N elements are the
!A            diagonal elements of the required tridiagonal matrix;
!A          E similarly holds the sub-diagonal matrix;  E(1) is arbitrary.
!A          N is the size of all required matrices and vectors;
!A          NP is their DIMENSION in the routine, and may be larger.
!A          Z is an NPxNP array holding an NxN matrix.  If the routine
!A            TRED2 has been used to reduce the matrix to tridiagonal, Z should
!A            on entry to TQLI contain the matrix on exit from TRED2.
!A            Otherwise it should contain a unit matrix.
!A          BOTH is a logical, set TRUE if both eigenvalues and eigenvectors
!A             are wanted.
!A On exit  D holds the eigenvalues, and the Kth column of Z holds the
!A            normalised eigenvector corresponding to D(K).
!
!N Adapted from Press, Flannery, Teukolsky & Vetterling, Numerical Recipes.
!
      LOGICAL BOTH
      DIMENSION D(NP), E(NP), Z(NP,NP)

      INTEGER         IBMBER
      COMMON /CCSLER/ IBMBER

      IF (N.LE.1) RETURN
      DO I = 2, N
        E(I-1) = E(I)
      ENDDO
      E(N) = 0.
      DO L = 1, N
        ITER = 0
    1   DO M = L, N - 1
          DD = ABS(D(M)) + ABS(D(M+1))
! @@ JvdS Isn't the following test very odd? Isn't it effectively testing E(M) = 0.0 ?
          IF (ABS(E(M))+DD.EQ.DD) GOTO 2
        ENDDO
        M = N
    2   IF (M.NE.L) THEN
          IF (ITER.EQ.30) THEN
            CALL ERRMES(1,0,'too many iterations IN TQLI')
            IF (IBMBER .NE. 0) RETURN
          ENDIF
          ITER = ITER + 1
          G = (D(L+1)-D(L))/(2.*E(L))
          R = SQRT(G**2+1.)
          G = D(M) - D(L) + E(L)/(G+SIGN(R,G))
          S = 1.
          C = 1.
          P = 0.
          DO I = M - 1, L, -1
            F = S*E(I)
            B = C*E(I)
! @@ JCC It seems that the next line causes numerous crashes in the Pawley refinement.
! As such, I've added in a check on the values of G and F to prevent it happening
! This ought to be tested - why is it causing the crash: the input data must be corrupt
! somehow.
            IF (ABS(F).GT. TINY(F)  .OR. ABS(G).GT.TINY(G)) THEN
              IF (ABS(F).GE.ABS(G)) THEN
                C = G/F
                R = SQRT(C**2+1.)
                E(I+1) = F*R
                S = 1./R
                C = C*S
              ELSE
                S = F/G
                R = SQRT(S**2+1.)
                E(I+1) = G*R
                C = 1./R
                S = S*C
              ENDIF
            ELSE
!>> Ok: An error has occurred in the passed data
              CALL BMBOUT
              CALL DebugErrorMessage('F or G equal 0 in TQLI() in Mvec.')
              RETURN
            ENDIF
            G = D(I+1) - P
            R = (D(I)-G)*S + 2.*C*B
            P = S*R
            D(I+1) = G + P
            G = C*R - B
            IF (BOTH) THEN
              DO K = 1, N
                F = Z(K,I+1)
                Z(K,I+1) = S*Z(K,I) + C*F
                Z(K,I) = C*Z(K,I) - S*F
              ENDDO
            ENDIF
          ENDDO
          D(L) = D(L) - P
          E(L) = G
          E(M) = 0.
          GOTO 1
        ENDIF
      ENDDO

      END SUBROUTINE TQLI
!
!*****************************************************************************
!
      SUBROUTINE TRANSQ(A,ND)
!
! *** TRANSQ by JCM 12 Jul 83 ***
!
!X
!C 12C
!H Replaces a real square matrix by its transpose.
!A On entry A is an NDxND real square matrix
!A On exit  A has been replaced by its transpose
!
      DIMENSION A(ND,ND)

      DO I = 1, ND - 1
        DO J = I + 1, ND
          B = A(I,J)
          A(I,J) = A(J,I)
          A(J,I) = B
        ENDDO
      ENDDO

      END SUBROUTINE TRANSQ
!
!*****************************************************************************
!
      SUBROUTINE TRED2(A,N,NP,D,E,BOTH)
!
! *** TRED2 by JCM from "NUMERICAL RECIPES" 4 Dec 89 ***
!
!X
!C 9C
!H Performs the Householder reduction of a real symmetric matrix to
!H tridiagonal form.
!A On entry A is a real, symmetric NxN matrix held in an NPxNP array
!A          BOTH is a LOGICAL which is TRUE if both eigenvalues and eigen-
!A          vectors will subsequently be wanted.
!A On exit  D is a vector of length NP whose first N elements hold the
!A            diagonal of the tridiagonal matrix,
!A          E similarly holds the sub-diagonal elements, and E(1)=0.
!A          If BOTH, then A is replaced by the orthogonal matrix Q which
!A          effects the transformation;  if .NOT. BOTH, A holds rubbish.
!
!N For eigenvalues and eigenvectors, should be followed by a call to TQLI.
!N Adapted from Press, Flannery, Teukolsky & Vetterling, Numerical Recipes.
!
      LOGICAL BOTH
      DIMENSION A(NP,NP), D(NP), E(NP)

      IF (N.LE.1) GOTO 1
      DO I = N, 2, -1
        L = I - 1
        H = 0.
        SCALE = 0.
        IF (L.GT.1) THEN
          DO K = 1, L
            SCALE = SCALE + ABS(A(I,K))
          ENDDO
          IF (SCALE.EQ.0.) THEN
            E(I) = A(I,L)
          ELSE
            DO K = 1, L
              A(I,K) = A(I,K)/SCALE
              H = H + A(I,K)**2
            ENDDO
            F = A(I,L)
            G = -SIGN(SQRT(H),F)
            E(I) = SCALE*G
            H = H - F*G
            A(I,L) = F - G
            F = 0.
            DO J = 1, L
              IF (BOTH) A(J,I) = A(I,J)/H
              G = 0.
              DO K = 1, J
                G = G + A(J,K)*A(I,K)
              ENDDO
              IF (L.GT.J) THEN
                DO K = J + 1, L
                  G = G + A(K,J)*A(I,K)
                ENDDO
              ENDIF
              E(J) = G/H
              F = F + E(J)*A(I,J)
            ENDDO
            HH = F/(H+H)
            DO J = 1, L
              F = A(I,J)
              G = E(J) - HH*F
              E(J) = G
              DO K = 1, J
                A(J,K) = A(J,K) - F*E(K) - G*A(I,K)
              ENDDO
            ENDDO
          ENDIF
        ELSE
          E(I) = A(I,L)
        ENDIF
        D(I) = H
      ENDDO
    1 IF (BOTH) D(1) = 0.
      E(1) = 0.
      DO I = 1, N
        IF (BOTH) THEN
          L = I - 1
          IF (D(I).NE.0.) THEN
            DO J = 1, L
              G = 0.
              DO K = 1, L
                G = G + A(I,K)*A(K,J)
              ENDDO
              DO K = 1, L
                A(K,J) = A(K,J) - G*A(K,I)
              ENDDO
            ENDDO
          ENDIF
        ENDIF
        D(I) = A(I,I)
        IF (BOTH) THEN
          A(I,I) = 1.
          IF (L.GE.1) THEN
            DO J = 1, L
              A(I,J) = 0.
              A(J,I) = 0.
            ENDDO
          ENDIF
        ENDIF
      ENDDO

      END SUBROUTINE TRED2
!
!*****************************************************************************
!
      SUBROUTINE TRIAN1(A,B,C,D,I)
!
! *** TRIAN1 by JCM 25 Jan 85 ***
!
!X
!C 10C
!H Applies the cosine formula for the solution of spherical triangles.
!A The entry and exit conditions vary according to the entry value of I:
!A     I=1: sets D=angle A given sides A,B,C
!A     I=2: sets A=side A given sides B,C and D=angle A.
!A     I=3: sets D=side A given angles A,B,C
!A     I=4: sets A=angle A given angles B,C and D=side A.
!D The angles are given and returned as their cosines.
!O Writes an error message is a cosine found >1
!
      INTEGER         LPT, LUNI
      COMMON /IOUNIT/ LPT, LUNI
!
      CS = B*C
      S = SQRT((1-B*B)*(1-C*C))
      GOTO (1,2,2,1), I
    1 CS = -CS
    2 GOTO (3,4,3,4), I
    3 IF (ABS(S).GT.1.E-3) GOTO 9
      IF (ABS(A+CS).GT.1.E-3) GOTO 10
      D = 0.
      GOTO 100
    9 X = (A+CS)/S
      GOTO 5
    4 X = D*S + CS
    5 IF (ABS(X)-1..LT.10.E-3) GOTO 8
   10 WRITE (LPT,3000) A, B, C, D, I, S, CS, X
      CALL BMBOUT
      RETURN
!
    8 IF (ABS(X).GT.1.) X = SIGN(1.,X)
      GOTO (6,7,6,7), I
    6 D = X
      GOTO 100
    7 A = X
  100 RETURN
 3000 FORMAT (/' ERROR ** in spherical triangle - cosine greater',      &
     &        ' than 1.0'/' A,B,C,D=',4F10.4,'I=',I2,'S,CS,X=',3F10.4)
      END SUBROUTINE TRIAN1
!
!*****************************************************************************
!
      SUBROUTINE UNIVEC(X,U)
!
! *** UNIVEC 1 May 84 by JCM ***
!
!X
!C 12C
!H Replaces a vector by a parallel unit vector, and gives its length.
!A On entry X is a real 1x3 vector referred to the CCSL orthogonal axes
!A On exit  X is replaced by a parallel unit vector
!A          U is set to the original length of X
!D If the length would be very small, X is not changed.
!N Works in double precision
!
      DOUBLE PRECISION DU
      DIMENSION X(3)
!
      U = 0.
      DU = 0.D0
      DO I = 1, 3
        DU = DU + X(I)*X(I)
      ENDDO
      IF (DU.LT.10.D-8) GOTO 100
      U = SNGL(DSQRT(DU))
      DO I = 1, 3
        X(I) = X(I)/U
      ENDDO
  100 RETURN
      END SUBROUTINE UNIVEC
!
!*****************************************************************************
!
      SUBROUTINE VECPRD(VEC1,VEC2,PRD)
!
! *** VECPRD by JCM ***
!
!X
!C 12C
!H Calculates the vector product of two 1x3 vectors.
!A On entry VEC1 and VEC2 hold 1x3 real vectors
!A On exit PRD = the vector product of VEC1 cross VEC2
!
      DIMENSION VEC1(3), VEC2(3), PRD(3)
      J = 2
      K = 3
      DO I = 1, 3
        PRD(I) = VEC1(J)*VEC2(K) - VEC1(K)*VEC2(J)
        J = K
        K = I
      ENDDO

      END SUBROUTINE VECPRD
!
!*****************************************************************************
!

max_stars_repo_name: ccdc-opensource/dash

max_issues_repo_name: ccdc-opensource/dash

max_forks_repo_name: ccdc-opensource/dash


text:
	SUBROUTINE GQDARR  ( szdarw, szdarh, idarwd, idartp, iret )
C************************************************************************
C* GQDARR 								*
C* 									*
C* This subroutine returns the current directional arrow size, arrow 	*
C* head size, line width and arrow type.				*
C* 									*
C* GQDARR  ( SZDARW, SZDARH, IDARWD, IDARTP, IRET )			*
C*									*
C* Output parameters:							*
C* 	SZDARW		REAL	 	Directional arrow size 		*
C* 	SZDARH		REAL		Arrow head size			*
C* 	IDARWD		INTEGER		Line width			*
C*	IDARTP		INTEGER		Arrow type			*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* I. Durham/GSC	 3/98						*
C* A. Hardy/GSC          6/98		Cleaned up prolog               *
C************************************************************************
	INCLUDE		'FUNCCODE.PRM'
	INCLUDE		'ERROR.PRM'
C*
	INTEGER		isend (2), ircv (2)
C------------------------------------------------------------------------
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend (1) = 2
	isend (2) = FQDARR
C
	CALL GPUT  ( isend, 2, iret )
	IF  ( iret .ne. NORMAL )  RETURN
C
C*	Get output parameters.
C
	CALL GGET  ( iret, 1, ierr )
	IF  ( ierr .ne. NORMAL )  THEN
	    iret = ierr
	    RETURN
	END IF
C
	CALL GGETR  ( szdarw, 1, ierr )
	IF  ( ierr .eq. NORMAL ) THEN
            CALL GGETR ( szdarh, 1, ierr )
            IF  ( ierr .eq. NORMAL ) THEN
	        CALL GGET ( ircv, 2, ierr )
                IF  ( ierr .eq. NORMAL ) THEN
		    idarwd = ircv (1)
                    idartp = ircv (2)
                END IF
            END IF
        END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      module rrlw_kg04

      use shr_kind_mod, only: r8 => shr_kind_r8

!      use parkind ,only : jpim, jprb

      implicit none
      save

!-----------------------------------------------------------------
! rrtmg_lw ORIGINAL abs. coefficients for interval 4
! band 4:  630-700 cm-1 (low - h2o,co2; high - o3,co2)
!
! Initial version:  JJMorcrette, ECMWF, jul1998
! Revised: MJIacono, AER, jun2006
!-----------------------------------------------------------------
!
!  name     type     purpose
!  ----   : ----   : ---------------------------------------------
!fracrefao: real
!fracrefbo: real
! kao     : real
! kbo     : real
! selfrefo: real
! forrefo : real
!-----------------------------------------------------------------

      integer, parameter :: no4  = 16

      real(kind=r8) :: fracrefao(no4,9)  ,fracrefbo(no4,6)
      real(kind=r8) :: kao(9,5,13,no4)
      real(kind=r8) :: kbo(5,5,13:59,no4)
      real(kind=r8) :: selfrefo(10,no4)  ,forrefo(4,no4)

!-----------------------------------------------------------------
! rrtmg_lw COMBINED abs. coefficients for interval 4
! band 4:  630-700 cm-1 (low - h2o,co2; high - o3,co2)
!
! Initial version:  JJMorcrette, ECMWF, jul1998
! Revised: MJIacono, AER, jun2006
!-----------------------------------------------------------------
!
!  name     type     purpose
!  ----   : ----   : ---------------------------------------------
! absa    : real
! absb    : real
!fracrefa : real
!fracrefb : real
! ka      : real
! kb      : real
! selfref : real
! forref  : real
!-----------------------------------------------------------------

      integer, parameter :: ng4  = 14

      real(kind=r8) :: fracrefa(ng4,9)  ,fracrefb(ng4,6)
      real(kind=r8) :: ka(9,5,13,ng4)   ,absa(585,ng4)
      real(kind=r8) :: kb(5,5,13:59,ng4),absb(1175,ng4)
      real(kind=r8) :: selfref(10,ng4)  ,forref(4,ng4)

      equivalence (ka(1,1,1,1),absa(1,1)),(kb(1,1,13,1),absb(1,1))

      end module rrlw_kg04

max_stars_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_issues_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_forks_repo_licenses: zlib-acknowledgement FTL RSA-MD


text:
      SUBROUTINE FI7505 (N,NBITS)
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C SUBPROGRAM:    FI7505      DETERMINE NUMBER OF BITS TO CONTAIN VALUE
C   PRGMMR: CAVANAUGH        ORG: W/NMC42    DATE: 93-06-23
C
C ABSTRACT: CALCULATE NUMBER OF BITS TO CONTAIN VALUE N, WITH A
C            MAXIMUM OF 32 BITS.
C
C PROGRAM HISTORY LOG:
C   93-06-23  CAVANAUGH
C   95-10-31  IREDELL     REMOVED SAVES AND PRINTS
C
C USAGE:    CALL FI7505 (N,NBITS)
C   INPUT ARGUMENT LIST:
C     N        - INTEGER VALUE
C
C   OUTPUT ARGUMENT LIST:      (INCLUDING WORK ARRAYS)
C     NBITS    - NUMBER OF BITS TO CONTAIN N
C
C REMARKS: SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
C
C ATTRIBUTES:
C   LANGUAGE: IBM VS FORTRAN 77, CRAY CFT77 FORTRAN
C   MACHINE:  HDS, CRAY C916/256, Y-MP8/64, Y-MP EL92/256
C
C$$$
      INTEGER        N,NBITS
      INTEGER        IBITS(31)
C
      DATA           IBITS/1,3,7,15,31,63,127,255,511,1023,2047,
     *               4095,8191,16383,32767,65535,131071,262143,
     *               524287,1048575,2097151,4194303,8388607,
     *               16777215,33554431,67108863,134217727,268435455,
     *               536870911,1073741823,2147483647/
C  ----------------------------------------------------------------
C
      DO 1000 NBITS = 1, 31
          IF (N.LE.IBITS(NBITS)) THEN
              RETURN
          END IF
 1000 CONTINUE
      RETURN
      END

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
Round Table Pizza is a pizza chain located in the northeast corner of Oakshade Town Center, just off of Pole Line Road. It offers dinein, carryout, and delivery.  They offer an all you can eat pizza/saladbar lunch buffets buffet everyday 11:30 AM  1:30 PM for $6.99, and a Tues night buffets buffet from 6 PM  to 8 PM for $7.49.

A variety of coupons can be found online at http://rtpdeals.com/.

They offer free WiFi.

Once upon a time, Round Table was downtown, where Kathmandu Kitchen is now. It went out of business, remained in restaurant heaven for quite some time, and then finally reemerged in the Davis food scene when Oakshade was developed.

20050222 02:51:16 nbsp get the garlic cheese twists Users/RishiTrivedi

20051111 17:36:06 nbsp I used to work at this one, and I can say they run a very clean and professional bussiness. Far better than other round tables I have seen/worked at.  Users/AlexGourley

20060617 21:29:23 nbsp Kathmandu Kitchen used to be a Mountain Mikes not a Round Table Users/SeanSmith
  Before it was a Mountian Mikes I remember it being a Round Table.  Users/JasonAller

20060722 16:42:29 nbsp Back in the Bay Area, they have ExtraLarge, but here, large is the biggest. Its still really tasty, but when I order from a Round Table, I desire XL!  Users/SamWard

20060913 15:25:19 nbsp Maui Zaui is my favorite! Users/RachelCakes

20070215 11:41:30 nbsp Ive had troubles lately in the evening, sometimes if you get there during the 9 Oclock hour theyll refuse to serve you, claiming theyre already closed. I submitted a comment to the website and got a prompt and friendly reply from the manager stating that their policy is not to stop taking orders until theyre really closed and that doing so is actually grounds for termination. I expect things will be better given his attitude towards the problem. Great pizza, btw. Users/ThomasPomroy

20070528 11:17:59 nbsp The round table in Woodland has been advertising a lunch buffet, but the one here in Davis doesnt seem to offer that choice. Users/AlexPomeranz

20070719 18:39:38 nbsp Very good pizza over the years since childhood i have come to love this pizza , of course i like other places but round table is a good pepperoni pizza anyone want to go let me know Users/dickjones

20070803 14:39:55 nbsp Great pizza always eat there about 4 times a year i like it very much clean store good employees great pizza always the same that is what i like about this joint! Users/Brians

20080128 23:26:02 nbsp Round Table has the Big Vinnie special ($12.99best deal) which I believe is their largest pie.  Ive had the Big Vinnie Pepperoni many times.  Ask for extra sauce because it tends to be rather dry otherwise.   Users/CurlyGirl26

20080827 12:36:32 nbsp The selection for the buffet usually consists of two or three pizzas, along with a platter of fleeting garlic twists. You need to get to the food fast, though, because youll have to do battle with a thousand screaming children in order to get any food on you plate. The kitchen seems to be regularly swamped and cannot put out new pizza anywhere near as fast as people take it. Users/JoePomidor

20080921 03:24:47 nbsp everytime i go there, it feels empty... and the salad bar looks like its been open for ages. Users/2point0student

20081016 22:07:35 nbsp The combination of painintheass customers and incompetent employees is hell.  Users/RJBTDP

20081120 11:49:30 nbsp Round Table has the best pizza buffet I have ever been to. The pizza is always fresh and the salad bar is to die for. I appreciate the excellent service as most places in Davis have forgotten how to give good service. Users/jefcar

20090304 14:30:04 nbsp Definitely amazing pizza and fresh goodness. ALWAYS come to the lunch or dinner buffet or have a coupon..cause its pretty overprixed. All you can eat slad for $5.50!? what the...and 6 garlic twists for $7.10?! yeah..Hunt down coupons or hit during the times they have buffets. YUMo. I was slightly sad when I asked one of the workers about what a Big Vinnie was..and they have no idea. _. Oh well, this place is right across from me though :D  Users/MissAmyQ

20090515 20:08:03 nbsp i love round table. its my favorite pizza chain.. BUT: 1. they keep undercooking their pizzas. 2. the last time I dined in, I was standing at the register waiting to order, but the girl picked up a brook and started sweeping, watched her sweep for a while before she took my order.. so annoying. 3.Then, I watched them leave my pizza on the counter for 10 minutes before bringing it to us. The guy there was about to bring it over to us, but he changed his mind and took other peoples order instead.  Users/julpham

20090805 00:57:44 nbsp The manager who runs this Round Table is very douchey.  Users/VMCrymore

20090901 17:42:22 nbsp Great place with delicious pizza. The Garlic Cheese Twists are to die for. They have very prompt and courteous service every time I dine in.
I do wish they had a greater vegetarian selection for the buffet. They only have Olive and Mushroom pizza. If they had to pick just one veggie pizza, why Olive and Mushroom??

Note: they describe their Gourmet Veggie pizza as having lots of chopped garlic. That is an understatement. The garlic overpowers the senses and the other flavors. Whether that is a good thing or not is up to you. Users/TheShah
  If you ask them to make a Veggie pizza (not Mushroom + Olive) and theyre not already making one they will make it for you. TusharRawat TR
   I went in for the buffet today. I asked if they could make a regular veggie pizza for the buffet, but the cashier said no, Its just how they keep track of things. I was the only one there, they made 4 pizzas, only one was vegetarian (the Mushroom + Olive). Kind of a waste considering the next people who ordered a buffet didnt come in for another half an hour. Users/TheShah
     Wow, thats pretty bad. I havent been to their buffet in a while but I may see if they do that all the time if I am there again. TusharRawat TR

20100109 14:59:23 nbsp My pizza had a good amount of pepperoni and the dough was really decent. The quality of the cheese and the sauce could probably be better, but the product is worth the money. Users/EdwardRiggins

20100531 10:00:53 nbsp I come here semiregularly for lunch buffet.  I always get prompt, courteous service. Users/drpottsiv

20100807 18:43:01 nbsp Be aware that Round Table is more expensive then some of the other pizza places in town, an extra Large chicken Garlic runs about $32, and delivery times on a saterday at about 6:30 take about an hour.  Users/PatrickMeddaugh

20101111 16:46:59 nbsp The pizza is what youd expect:  average.  The interior is grocerymall corporate:  formica and industrial carpet, completely charmless if clean.  That said, the place works quite well for families with small kids and they are very good about having team parties and so forth.  And the service is generally good  nice college kids.  I hate to say it, but Ive been going here for various events for about 10 years now. Users/Flokkenfisch

20110727 14:25:51 nbsp recently I had a conversation with a current employee while at another pizza restaurant that we both were dining at about the way round table is run from the inside and I must say I am disturbed at what I have heard I have on some occasions had problems with my order but I now believe the management is most likely the cause my wife orders from round table 3 times a month at least but after what I have heard I dont believe we will be ordering from here anymore its is not the family community supporting place we believed it to be. Users/codybanks

20111210 20:38:06 nbsp The Ultimeat pizza is quite tasty and the employees are always very nice.  I miss the fire roasted pizzas :( Users/LoriOrf

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C       SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)
       SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, ISTAT, L)
c  ------------------------------------------------------------------------
c  Part of "poa_calfos" (STECF/IPMG)
c  Part of "NEWBILCAL" program to produce IGRF geomagnetic model (GSFC)
c
c  History:
c  --------
c  Version   Date        Author          Description
c      1.0       Jul 00  M. Rosa         POA calibration code
c      1.2       Jun 01  A. Alexov       Added IDS's common block for file
c                                        header access
c  ------------------------------------------------------------------------
c       15 Jul 00       reads from STSDAS tables
c
C ===============================================================
C       Version 1.01
C
C       Reads spherical harmonic coefficients from the specified
C       file into an array.
C
C       Input:
C           IU    - Logical unit number
C           FSPEC - File specification
C
C       Output:
C           NMAX  - Maximum degree and order of model
C           ERAD  - Earth's radius associated with the spherical
C                   harmonic coefficients, in the same units as
C                   elevation
C           GH    - Schmidt quasi-normal internal spherical
C                   harmonic coefficients
C           IER   - Error number: =  0, no error
C                                 = -2, records out of order
C                                 = FORTRAN run-time error number
C
C       A. Zunde
C       USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225
C ===============================================================
       INTEGER  I,N,M,NN,MM,NMAX,IER,IU,L
c       REAL*4   DATE,ERAD,G,H
       REAL*4   ERAD,G,H
       REAL*4   GH(*)
c       CHARACTER*10 NAME
       CHARACTER FSPEC*(*)
C
C     UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
      CHARACTER*80 CONTXT

C
C FILE I/O ACCESS MODES
C
        INTEGER   RDONLY
        PARAMETER (RDONLY = 1)

      INTEGER   TBNROW
      PARAMETER (TBNROW = 21)
      INTEGER IDIN,COLIDS(4),I,II,JJ,ISTATS(4),ISTAT,NROWS
      CHARACTER*8 COLNAM(4)
      LOGICAL NULL
      DATA COLNAM/'NDEGREE','MORDER','GCOEFF','HCOEFF'/
C
C ---------------------------------------------------------------
C       Open coefficient file. Read past first header record.
C       Read degree and order of model and Earth's radius.
C ---------------------------------------------------------------
c       IU = 81
c       OPEN (UNIT=IU,FILE='dgrf90.dat',ACCESS='SEQUENTIAL',
c     1       FORM='FORMATTED',STATUS='OLD',IOSTAT=IER,ERR=999)
cc 991   type *,'opened ',IU,FSPEC,IER
c 881   READ (IU,*) NAME
c 882   READ (IU,*,IOSTAT=IER, ERR=999) NMAX, ERAD, DATE
C ---------------------------------------------------------------
C       Read the coefficient file, arranged as follows:
C
C                                   N     M     G     H
C                                   ----------------------
C                                /   1     0    GH(1)  -
C                               /       1     1    GH(2) GH(3)
C                              /       2     0    GH(4)  -
C                             /       2     1    GH(5) GH(6)
C           NMAX*(NMAX+3)/2        /       2     2    GH(7) GH(8)
C              records              \       3     0    GH(9)  -
C                             \      .     .     .     .
C                              \       .     .     .     .
C           NMAX*(NMAX+2)          \       .     .     .     .
C           elements in GH           \  NMAX  NMAX   .     .
C
C       N and M are, respectively, the degree and order of the
C       coefficient.
C ---------------------------------------------------------------
c 883   I = 0
c       DO 2211 NN = 1, NMAX
c         DO 2233 MM = 0, NN
c           READ (IU,*,IOSTAT=IER, ERR=999) N, M, G, H
c           IF (NN .NE. N .OR. MM .NE. M) THEN
c             IER = -2
c             GOTO 999
c           ENDIF
c           I = I + 1
c           GH2(I) = G
c           IF (M .NE. 0) THEN
c             I = I + 1
c             GH2(I) = H
c           ENDIF
c2233     CONTINUE
c2211   CONTINUE

c999    CLOSE (IU)
c       RETURN
c       END

C
C Open table
C
        CALL UTTOPN(FSPEC,RDONLY,IDIN,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='ERROR opening dgrf table '//FSPEC
                GO TO 998
        ENDIF
C
C get NMAX and ERAD header keywords
C
        CALL UTHGTI(IDIN,'NMAX',NMAX,ISTAT)
        CALL UTPGTI(IDIN,TBNROW,NROWS,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='ERROR reading NMAX key in dgrf table '//FSPEC
                GO TO 999
        ENDIF
        CALL UTHGTR(IDIN,'ERAD',ERAD,ISTAT)
        CALL UTPGTI(IDIN,TBNROW,NROWS,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='ERROR reading ERAD key in dgrf table '//FSPEC
                GO TO 999
        ENDIF
C
C get number of rows
C
        CALL UTPGTI(IDIN,TBNROW,NROWS,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='ERROR reading dgrf table '//FSPEC
                GO TO 999
        ENDIF
C
C Get column ids.
        CALL UTCFND(IDIN,COLNAM,4,COLIDS,ISTAT)
        IF(ISTAT.NE.0)THEN
                CONTXT='ERROR locating needed columns in dgrf table '//
     *                  FSPEC
                GO TO 999
        ENDIF

c get the data
        I = 0
        II = 1
        DO 2214 NN = 1, NMAX
           DO 2235 MM = 0, NN
              CALL UTRGTI(IDIN,COLIDS(1),1,II,N,NULL,ISTATS(1))
              CALL UTRGTI(IDIN,COLIDS(2),1,II,M,NULL,ISTATS(2))
              CALL UTRGTR(IDIN,COLIDS(3),1,II,G,NULL,ISTATS(3))
              CALL UTRGTR(IDIN,COLIDS(4),1,II,H,NULL,ISTATS(4))
              DO 30 JJ=1,4
                 IF(ISTATS(JJ).NE.0)THEN
                    CONTXT='ERROR reading dgrf table '//FSPEC
                    GO TO 999
                 ENDIF
30            CONTINUE
              IF (NN .NE. N .OR. MM .NE. M) THEN
                IER = -2
                GOTO 999
              ENDIF
              I = I + 1
              GH(I) = G
              IF (M .NE. 0) THEN
                I = I + 1
                GH(I) = H
              ENDIF
              II = II + 1
 2235      CONTINUE
 2214   CONTINUE

C close the table
        CALL UTTCLO(IDIN,ISTAT)
        ISTAT=0
        GO TO 1000
C
999     CALL UTTCLO(IDIN,ISTAT)
998     CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
        ISTAT=1
1000    RETURN
        END

max_stars_repo_name: iraf-community/stsdas


text:
      DOUBLE PRECISION FUNCTION alnrel(a)
C-----------------------------------------------------------------------
C            EVALUATION OF THE FUNCTION LN(1 + A)
C-----------------------------------------------------------------------
C     .. Scalar Arguments ..
      DOUBLE PRECISION a
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION p1,p2,p3,q1,q2,q3,t,t2,w,x
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,dble,dlog
C     ..
C     .. Data statements ..
      DATA p1/-.129418923021993D+01/,p2/.405303492862024D+00/,
     +     p3/-.178874546012214D-01/
      DATA q1/-.162752256355323D+01/,q2/.747811014037616D+00/,
     +     q3/-.845104217945565D-01/
C     ..
C     .. Executable Statements ..
C--------------------------
      IF (abs(a).GT.0.375D0) GO TO 10
      t = a/ (a+2.0D0)
      t2 = t*t
      w = (((p3*t2+p2)*t2+p1)*t2+1.0D0)/ (((q3*t2+q2)*t2+q1)*t2+1.0D0)
      alnrel = 2.0D0*t*w
      RETURN
C
   10 x = 1.D0 + dble(a)
      alnrel = dlog(x)
      RETURN

      END

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy


subroutine antest_maps0_glb(mype,theta0f,z0f,theta2f,z2f,theta3f,z3f)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    antest_maps0_glb
!   prgmmr:
!
! abstract: this routine creates output maps of background error correlations
!
! program history log:
!   2009-09-21  lueken - added subprogram doc block
!   2010-03-30  zhu    - use nvars from control_vectors
!   2010-12-05  pondeca - change variable names: "tv" is now "t" , and "st" is "sf"
!   2013-10-19  todling - metguess now holds background
!   2016-04-19  pondeca - comment out call to "ansmoothrf(howrk)". must update to use
!                         bundle type as argument
!
!   input argument list:
!    mype
!    theta0f
!    z0f
!    theta2f
!    z2f
!    theta3f
!    z3f
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind,r_single
  use anberror, only: kvar_start,kvar_end,var_names,indices,pf2aP1,pf2aP2,pf2aP3
  use gridmod, only: nsig,nsig1o,nlon,nlat,istart,jstart,lat2,lon2
  use constants, only: zero_single,zero,one,rd_over_cp,r100
  use mpimod, only: ierror,mpi_real4,mpi_real8,mpi_sum,mpi_comm_world
  use guess_grids, only: ntguessig,ges_prsl
  use patch2grid_mod, only: vpatch2grid
  use control_vectors, only: nvars
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use mpeu_util, only: die
  implicit none

  integer(i_kind),intent(in   ) :: mype
  real(r_single) ,intent(in   ) :: theta0f(pf2aP1%nlatf,pf2aP1%nlonf,nsig1o)
  real(r_single) ,intent(in   ) :: z0f(pf2aP1%nlatf,pf2aP1%nlonf,nsig1o)
  real(r_single) ,intent(in   ) :: theta2f(pf2aP2%nlatf,pf2aP2%nlonf,nsig1o)
  real(r_single) ,intent(in   ) :: z2f(pf2aP2%nlatf,pf2aP2%nlonf,nsig1o)
  real(r_single) ,intent(in   ) :: theta3f(pf2aP3%nlatf,pf2aP3%nlonf,nsig1o)
  real(r_single) ,intent(in   ) :: z3f(pf2aP3%nlatf,pf2aP3%nlonf,nsig1o)

  character(len=*),parameter::myname='antest_maps0_glb'
  real(r_kind),dimension(nlat,nlon,nsig1o):: hwork
  real(r_kind) tempf(nlat,nlon), &
               tempc  (pf2aP1%nlatf,pf2aP1%nlonf), &
               tempcp2(pf2aP2%nlatf ,pf2aP2%nlonf ), &
               tempcp3(pf2aP3%nlatf ,pf2aP3%nlonf )
  real(r_single) outwork(nlon,nlat),outwork0(nlon,nlat)
  character(80) ref_plotcor
  character(80) var_plotcor
  integer(i_kind) i_plotcor,j_plotcor,k_plotcor

  real(r_kind) h00,h000
  integer(i_kind) lunin,i,j,k,ivar,iglob,jglob,ivar_plot,k_plot,ielm
  integer(i_kind) it,mm1,ier,istatus

  real(r_kind),dimension(:,:  ),pointer:: ges_z_it=>NULL()
  real(r_kind),dimension(:,:,:),pointer:: ges_tv_it=>NULL()

!*********************************************************************
!          variable names expected for var_plotcor are
!
!    sf  -- stream function
!    vp  -- velocity potential
!    ps  -- surface pressure
!    t   -- virtual temperature
!    q   -- specific humidity
!    oz  -- ozone
!    sst -- sea surface temperature
!    stl -- skin temp over land
!    sti -- skin temp over ice
!    cw  -- cloud water
!*********************************************************************
! Make choices here!
! i_plotcor=500
! j_plotcor=500
! k_plotcor=25
! var_plotcor='sf'
!Note: Must call this subroutine from anprewgt_reg.f90
!Make sure statement has been uncommented!
! End of choice section
!*********************************************************************

  do ielm=1,5
     select case (ielm)
        case(1)   ; var_plotcor='sf'
        case(2)   ; var_plotcor='vp'
        case(3)   ; var_plotcor='t'
        case(4)   ; var_plotcor='q'
        case(5)   ; var_plotcor='ps'
     end select

     ref_plotcor=var_plotcor
     it=ntguessig

     ier=0
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'z' ,ges_z_it,   istatus)
     ier=ier+istatus
     call gsi_bundlegetpointer (gsi_metguess_bundle(it),'tv',ges_tv_it,  istatus)
     ier=ier+istatus
     if(ier/=0) call die(myname,'missing fields, ier= ', ier)

     lunin=91
     if(mype==0) then
        open(lunin,file="cormaps_"//trim(var_plotcor),form='unformatted')
        rewind lunin
     end if

     ivar_plot=0
     do ivar=1,nvars
        if(trim(var_names(ivar))==trim(var_plotcor)) then
           ivar_plot=ivar
           exit
        end if
     end do

     if(ivar_plot==0) then
        write(6,*)' in antest_maps0, variable ',trim(var_plotcor),'  not found.  program stops'
        call mpi_finalize(ierror)
        stop
     end if

     hwork=zero
     do k_plotcor=10,50,20
        if(var_plotcor=='ps') then
           k_plot=kvar_start(ivar_plot)
        else
           k_plot=kvar_start(ivar_plot)+k_plotcor-1
        end if

        if(k_plot>=indices%kps.and.k_plot<=indices%kpe) then
           if( var_plotcor /= 'hires' ) then

              if(k_plotcor<50.or.var_plotcor=='ps') then
                 do i_plotcor=nlat/8,nlat-nlat/8+1,nlat/8
                    do j_plotcor= nlon/4,nlon,nlon/2
                       hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one
                    end do
                 end do
                 do i_plotcor=nlat/6,nlat-nlat/6+1,nlat/6
                    do j_plotcor=   1,nlon,nlon/2
                       hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one
                    end do
                 end do
              else
                 do i_plotcor=nlat/4,nlat-nlat/4+1,nlat/4
                    do j_plotcor= 1,nlon,nlon/4
                       hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one
                    end do
                 end do
              end if

              i_plotcor=1
              j_plotcor=1
              hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one

              i_plotcor=nlat
              j_plotcor=1
              hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one

           else

              do j_plotcor= 1,nlon,nlon/16
                 do i_plotcor=(nlat+2)/8,nlat-(nlat+2)/8-1,(nlat+2)/16
                    hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one
                 end do
              end do

              i_plotcor=1
              j_plotcor=1
              hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one

              i_plotcor=nlat
              j_plotcor=1
              hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one

              i_plotcor=(nlat+2)/32
              do j_plotcor= 1, nlon,nlon/4
                 hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one
              end do
              i_plotcor=nlat-(nlat+2)/32
              do j_plotcor= 1, nlon,nlon/4
                 hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one
              end do

              i_plotcor=(nlat+2)/16
              do j_plotcor= 1, nlon,nlon/8
                 hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one
              end do
              i_plotcor=nlat-(nlat+2)/16
              do j_plotcor= 1, nlon,nlon/8
                 hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)=one
              end do

           end if
        end if
        if(var_plotcor=='ps') exit
     end do

!    call ansmoothrf(hwork) !Must replace hwork with bundle / MPondeca 19Apr2016

     if(mype==0) write(lunin) ref_plotcor,var_plotcor,j_plotcor,i_plotcor,k_plotcor, &
                  nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1

     if(mype==0) write(6,*) ' refvar= ',trim(ref_plotcor),' corvar= ',trim(var_plotcor), &
            '  i,j,k_plotcor =', j_plotcor,i_plotcor,k_plotcor, ' nlon,nlat,nsig=', &
                  nlon,nlat,kvar_end(ivar_plot)-kvar_start(ivar_plot)+1

     !---------------------in case we haven't normalized, divide by value of correlation point

     h00=zero
     if(k_plot>=indices%kps.and.k_plot<=indices%kpe) h00=hwork(i_plotcor,j_plotcor,k_plot-indices%kps+1)
     call mpi_allreduce(h00,h000,1,mpi_real8,mpi_sum,mpi_comm_world,ierror)
     hwork=hwork/h000

!     output original pot temp  (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     do k=1,nsig
        outwork=zero_single
        do j=2,lon2-1
           jglob=jstart(mm1)-2+j
           do i=2,lat2-1
              iglob=istart(mm1)-2+i
              outwork(jglob,iglob)=ges_tv_it(i,j,k)/(ges_prsl(i,j,k ,it)/r100)**rd_over_cp
           end do
        end do
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

!             output "smoothed pot temp"

     do k=kvar_start(ivar_plot),kvar_end(ivar_plot)
        outwork=zero_single
        if(k>=indices%kps.and.k<=indices%kpe) then
           do j=1,pf2aP1%nlonf
              do i=1,pf2aP1%nlatf
                 tempc(i,j)=theta0f(i,j,k-indices%kps+1)
              end do
           end do
           do j=1,pf2aP2%nlonf
              do i=1,pf2aP2%nlatf
                 tempcp2(i,j)=theta2f(i,j,k-indices%kps+1)
              end do
           end do
           do j=1,pf2aP3%nlonf
              do i=1,pf2aP3%nlatf
                 tempcp3(i,j)=theta3f(i,j,k-indices%kps+1)
              end do
           end do
           call vpatch2grid(tempf,tempc,tempcp2,tempcp3)
           do j=1,nlon
              do i=1,nlat
                 outwork(j,i)=tempf(i,j)
              end do
           end do
        end if
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

     do k=kvar_start(ivar_plot),kvar_end(ivar_plot)
        outwork=zero_single
        if(k>=indices%kps.and.k<=indices%kpe) then
           do j=1,nlon
              do i=1,nlat
                 outwork(j,i)=hwork(i,j,k-indices%kps+1)
              end do
           end do
        end if

!             very slow way to move field from local processor to processor 0

        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

!     output original terrain (slow way to get full 2d field)  -- this is reference field

     mm1=mype+1
     outwork=zero_single
     do j=2,lon2-1
        jglob=jstart(mm1)-2+j
        do i=2,lat2-1
           iglob=istart(mm1)-2+i
           outwork(jglob,iglob)=ges_z_it(i,j)
        end do
     end do
     call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
     if(mype==0) write(lunin) outwork0

!             output "smoothed terrain"

     PRINT*,'IN ANPREWGT,KPS,KPE=',indices%KPS,indices%KPE
     do k=1, 1  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!kvar_start(ivar_plot),kvar_end(ivar_plot)
        outwork=zero_single
        if(k>=indices%kps.and.k<=indices%kpe) then
           do j=1,pf2aP1%nlonf
              do i=1,pf2aP1%nlatf
                 tempc(i,j)=z0f(i,j,k)!  theta0f(i,j,k-indices%kps+1)
              end do
           end do
           do j=1,pf2aP2%nlonf
              do i=1,pf2aP2%nlatf
                 tempcp2(i,j)=z2f(i,j,k)!  theta0f(i,j,k-indices%kps+1)
              end do
           end do
           do j=1,pf2aP3%nlonf
              do i=1,pf2aP3%nlatf
                 tempcp3(i,j)=z3f(i,j,k)!  theta0f(i,j,k-indices%kps+1)
              end do
           end do
           call vpatch2grid(tempf,tempc,tempcp2,tempcp3)
           do j=1,nlon
              do i=1,nlat
                 outwork(j,i)=tempf(i,j)
              end do
           end do
        end if
        call mpi_reduce(outwork,outwork0,nlon*nlat,mpi_real4,mpi_sum,0,mpi_comm_world,ierror)
        if(mype==0) write(lunin) outwork0
     end do

     close(lunin)
!    if(mype>-1000) then
!       call mpi_finalize(i)
!       stop
!    end if

  end do

end subroutine antest_maps0_glb

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


text:
Bower Bird is closed. There is a childrens clothing store in its place now.

Nothing that you need, but everything you want.

Bower Birds are collectors, so if you have kindred tendencies this gift shop may be for you. They have antiques, French soaps and lotions, unique jewelry, ribbon, and various trinkets for sale.

Come visit Deborahs resident Zebra Finches!

Walked by the store today, it was cleared out and empty.