text:
C      FGHEVEN
C      F. Gogtas, G.G. Balint-Kurti and C.C. Marston,
C      QCPE Program No. 647, (1993).
C
C     ***************************************************************
C     This program solves one dimensional Schrodinger equation for
C     bound state eigenvalues and eigenfunctions corresponding to a
C     potential V(x).
C
C     The Fourier Grid Hamiltonian method is used in the program.
C     This method is fully described in:
C     C.C. Marston and G.G.Balint-Kurti, J.Chem. Phys.,91,3571(1989).
C     and
C     G.G. Balint-Kurti, R.N. Dixon and C.C. Marston, Internat. Rev.
C     Phys. Chem.,11, 317 (1992).
C
C     The program uses an even number of grid points. The Hamiltonian
C     matrix elements are calculated using analytic formula described
C     in the second of the above references.
C
C     The analytical Hamiltonian expression given in the reference
C     contains a small error.  The formula should read:
C
C                H(i,j) = {(h**2)/(4*m*(L**2)} *
C                          [(N-1)(N-2)/6 + (N/2)] + V(Xi),   if i=j
C
C                H(i,j) = {[(-1)**(i-j)] / m } *
C                          { h/[2*L*sin(pi*(i-j)/N)]}**2 ,     if i#j
C
C     The eigenvalues of the Hamiltonian matrix which lie below the
C     asymptotic (large x) value of V(x) correspond to the bound state
C     energies. The corresponding eigenvectors are the representation
C     of the bound state wavefunctions on a regular one dimensional grid.
C
C     The user must supply a self contained subroutine VSUB(X,V)
C     which is called with a value of X and returns the value of the
C     potential in V.
C
C     As supplied the program uses a Morse potential with parameters chosen
C     to represent the HF molecule.
C
C     The Program should work on any computer
C     Programming language used                    : Fortran77
C     Memory required to execute with typical data : 1 Mbyte
C
C     To run the program on a SUN workstation (or any Unix based system)
C     First type:
C        f77 -o fgheven fgheven.f
C     to create an executable program (fgheven).
C     Then type:
C        fgheven > fgheven.out &
C     to run the program.
C
C     The dimensions of the arrays generally depend on the number of grid
C     points used.  This number NX is set in a PARAMETER statement.
C     The following arrays represent  :
C
C        XA  : X-coordinates i.e. position on a grid.
C
C        AR  : The Hamiltonian Matrix
C
C        WCH : The eigenvalues for respective energy levels.
C
C        ZR  : The eigenvectors (X,Y); where
C                                      X : Wavefunction.
C                                      Y : Energy level.
C    The folowing data constans represent :
C        ZMA, ZMB : Mass of atoms.
C        R0       : Equilibrium bond separation
C        NPRIN    : 0 or 1 for eigenvalue or eigenvector respectively
C        NWRIT    : Number of eigenvales/eigenvectors to be printed
C        RMIN     : Starting point of grid
C        RMAX     : End point of grid
C        ZL       : Grid length
C        DX       : Grid spacings
C
C       All quantities in au (unless otherwise stated)
C     ****************************************************************
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (NX=1800)
C
C.....Declare arrays.
      DIMENSION XA(NX),AR(NX,NX),WCH(NX),ZR(NX,NX)
      DIMENSION FV1(NX),FV2(NX)
      DIMENSION POT(NX)
C
C.....Variable data input.
C     DATA ZMA/1836.9822D0/
C     DATA ZMB/34629.61319D0/
C     DATA R0/1.7329D0/
C     DATA DHART/0.2250073497D0/
C     DATA BETA/1.1741D0/
      DATA PI/3.141592653589793D0/
      DATA NPRIN/0/
      DATA NWRIT/20/
      INTEGER :: L
      LOGICAL :: MORLJ

C     NAMELIST DATA
      NAMELIST /TTDATA/ ZMA, ZMB, R0, BETA, DHART, NPRIN, NWRIT, L,
     &RMAX, RMIN, MORLJ, C6, C8, C10, C12
      READ(5,NML=TTDATA)

C
C....Test that NX is even
      ITEST=MOD(NX,2)
      IF(ITEST.NE.0) THEN
        WRITE(6,*)' **** NX MUST BE EVEN-FATAL ERROR  ****'
        STOP
      END IF
C
C.....Set up grid
      WRITE(6,*)'Grid paremeters:'
      WRITE(6,*)'   Number of grid poins    = ',NX
C      RMIN=R0*0.2D0
C      RMAX=R0*3.0D0
      ZL=(RMAX-RMIN)
      WRITE(6,*)'   Grid length             = ',ZL
      DX=ZL/DFLOAT(NX)
      WRITE(6,*)'   Grid spacings           = ',DX
C.....Print mass of atoms A and B
      WRITE(6,*)'Mass of atoms:'
      WRITE(6,*)'   Atomic mass of atom A      = ',ZMA
      WRITE(6,*)'   Atomic mass of atom B      = ',ZMB
      WRITE(6,*)'   R0                         = ',R0
      WRITE(6,*)'   BETA                       = ',BETA
      WRITE(6,*)'   DHART                      = ',DHART
C.....Compute reduced mass.
      ZMU=ZMA*ZMB/(ZMA+ZMB)
      WRITE(6,*)'   Reduced mass of AB         = ',ZMU
C
C.....Compute contants:
      PSQ=PI*PI
      CONST1=PSQ/(ZMU*(ZL**2))
      NFACT1=(NX-1)*(NX-2)
      NFACT2=(NX-2)/2
      CONST2=CONST1*(DFLOAT(NFACT1)/6.0D0 +1.0D0+DFLOAT(NFACT2))
      DARG=PI/DFLOAT(NX)
      VMIN=0.0D10
C
C.....Now compute Hamiltonian matrix:
      X=RMIN
      DO 003 I=1,NX
         XA(I)=X
       DO 004 J=1,I
          IJD=(I-J)
        IF(IJD.EQ.0)THEN
          AR(I,J)=CONST2
        ELSE
          RATIO=1.0D0/DSIN(DARG*DFLOAT(IJD))
          AR(I,J)=((-1)**IJD)*CONST1*(RATIO**2)
        END IF
 004  CONTINUE
C.....Find the potential value at x
      IF(MORLJ) THEN
       CALL VSUB_M(X,VV,R0,BETA,DHART,L,ZMU)
      ELSE
       CALL VSUB_LJ(X,VV,R0,C6,C8,C10,C12,L,ZMU)
      ENDIF
      POT(I) = VV
      IF(VMIN .gt. VV) VMIN = VV
C.....Add the potential value when the kronecker delta function
C.....equals to one, i.e. when I and J are equal
        AR(I,I)=AR(I,I)+VV
        X=X+DX
 003  CONTINUE
C.....Now fill out Hamiltonian matrix.
      DO 006 I=1,NX
       DO 007 J=1,I
        AR(J,I)=AR(I,J)
007    CONTINUE
006   CONTINUE
C.....Now call eigenvalue solver.
      call flush(6)
      CALL RS(NX,NX,AR,WCH,NPRIN,ZR,FV1,FV2,IERR)
C.....Print out eigenvalues and eigenvectors.
      WRITE(6,12)
12    FORMAT(/)
      WRITE(6,*)' The first',NWRIT,' energy levels for AB molecule '
      DO 009 I=1,NWRIT
       WRITE(6,*)' ENERGY LEVEL NO ',I,' EIGENVALUE= ',WCH(I)
       IF (WCH(I) .LT. 0.0D0 .AND. WCH(I) .GT. VMIN) THEN
        WRITE(10,*) I-1.0, WCH(I)
       ENDIF
 009  CONTINUE
      WRITE(6,12)
      WRITE(6,*)' THE CORRESPONDING EIGENFUNCTIONS ARE:'
      DO 010 I=1,NWRIT
       WRITE(6,*)' ENERGY LEVEL NO ',I,' EIGENVALUE= ',WCH(I)
       IF (NPRIN.EQ.1) THEN
       DO 011 J=1,NX
         WRITE(6,*)' R = ',XA(J),' WAVEFUNCTION=',ZR(J,I)
 011   CONTINUE
       ENDIF
 010  CONTINUE
      OPEN(unit=99,file="pot.dat")
      DO 013 I=1,NX
         WRITE(99,*)XA(I),POT(I)
 013  CONTINUE
      CLOSE(99)
      STOP
      END
C
C
      SUBROUTINE VSUB_M(X,V,R0,BETA,DHART,L,MU)
C     **************************************************************
C.....SUBROUTINE TO COMPUTE POTENTIAL ENERGY AT SEPERATION X (IN
C.....ATOMIC UNITS). THE VALUE OF THE POTENTIAL IS RETURNED IN V,
C.....AGAIN IN ATOMIC UNITS.THE POTENTIAL ENERGY IS ASSUMED TO
C..... BE A MORSE CURVE:V(X) = DHART * (EXP(-BETA*(X-R0)) - 1)**2
C..... THE FOLLOWING DATA CONSTANTS REPRESENT :
C.....    DHART : THE DISSOCIATION ENERGY FOR THE GIVEN MOLECULE
C.....    BETA  : THE EXPONENTIAL PARAMETER
C.....    R0    : THE BOND LENGTH AT EQUILIBRIUM FOR GIVEN MOLECULE
C     **************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION :: MU
      INTEGER :: L
C      DATA DHART/0.2250073497D0/
C      DATA BETA/1.1741D0/
C      DATA R0/1.7329D0/
C.....Set morse curve parameters.
      ARG=X-R0
      ARG=-BETA*ARG
      EF=DEXP(ARG)
      EF=1.0D0-EF
      EF=EF*EF
      V=DHART*EF-DHART
      V = V + DBLE(L)*(DBLE(L)+1.0)/(2.0*MU*X**2)
      RETURN
      END
CC
      SUBROUTINE VSUB_LJ(X,V,R0,C6,C8,C10,C12,L,MU)
C     **************************************************************
C.....SUBROUTINE TO COMPUTE POTENTIAL ENERGY AT SEPERATION X (IN
C.....ATOMIC UNITS). THE VALUE OF THE POTENTIAL IS RETURNED IN V,
C.....AGAIN IN ATOMIC UNITS.THE POTENTIAL ENERGY IS ASSUMED TO
C..... BE A MORSE CURVE:V(X) = DHART * (EXP(-BETA*(X-R0)) - 1)**2
C..... THE FOLLOWING DATA CONSTANTS REPRESENT :
C.....    DHART : THE DISSOCIATION ENERGY FOR THE GIVEN MOLECULE
C.....    BETA  : THE EXPONENTIAL PARAMETER
C.....    R0    : THE BOND LENGTH AT EQUILIBRIUM FOR GIVEN MOLECULE
C     **************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION :: MU
      INTEGER :: L
C      DATA DHART/0.2250073497D0/
C      DATA BETA/1.1741D0/
C      DATA R0/1.7329D0/
C.....Set morse curve parameters.
      V =C12/X**12-C10/X**6-C8/X**6-C6/X**6
      V = V + DBLE(L)*(DBLE(L)+1.0)/(2.0*MU*X**2)
      RETURN
      END
C
C
      SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
C     ****************************************************************
C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
C     OF A REAL SYMMETRIC MATRIX.
C
C     ON INPUT :
C
C        N,M  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C        DIMENSION STATEMENT,
C        N  IS THE ORDER OF THE MATRIX  A,
C        A  CONTAINS THE REAL SYMMETRIC MATRIX,
C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
C        ONLY EIGENVALUES ARE DESIRED,  OTHERWISE IT IS SET TO
C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
C
C     ON OUTPUT :
C
C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER,
C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO,
C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN
C        ERROR COMPLETION CODE DESCRIBED IN SECTION 2B OF THE
C        DOCUMENTATION.  THE NORMAL COMPLETION CODE IS ZERO,
C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C    *******************************************************************
      IF (N .LE. NM) GO TO 10
      IERR = 10 * N
      GO TO 50
C
   10 IF (MATZ .NE. 0) GO TO 20
C     ********** FIND EIGENVALUES ONLY **********
      CALL  TRED1(NM,N,A,W,FV1,FV2)
      CALL  TQLRAT(N,W,FV2,IERR)
      GO TO 50
C     ********** FIND BOTH EIGENVALUES AND EIGENVECTORS **********
   20 CALL  TRED2(NM,N,A,W,FV1,Z)
      CALL  TQL2(NM,N,W,FV1,Z,IERR)
   50 RETURN
C     ********** LAST CARD OF RS **********
      END
C
C
      SUBROUTINE TRED1(NM,N,A,D,E,E2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NM,N),D(N),E(N),E2(N)
C     ***************************************************************
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX
C     TO A SYMMETRIC TRIDIAGONAL MATRIX USING
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT :
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX,
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C     ON OUTPUT :
C
C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
C          FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
C          TRIANGLE.  THE FULL UPPER TRIANGLE OF A IS UNALTERED,
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO,
C
C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     ****************************************************************
C
      DO 100 I = 1, N
  100 D(I) = A(I,I)
C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
      DO 300 II = 1, N
         I = N + 1 - II
         L = I - 1
         H = 0.0D0
         SCALE = 0.0D0
         IF (L .LT. 1) GO TO 130
C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
         DO 120 K = 1, L
  120    SCALE = SCALE + DABS(A(I,K))
C
         IF (SCALE .NE. 0.0D0) GO TO 140
  130    E(I) = 0.0D0
         E2(I) = 0.0D0
         GO TO 290
C
  140    DO 150 K = 1, L
            A(I,K) = A(I,K) / SCALE
            H = H + A(I,K) * A(I,K)
  150    CONTINUE
C
         E2(I) = SCALE * SCALE * H
         F = A(I,L)
         G = -DSIGN(DSQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         A(I,L) = F - G
         IF (L .EQ. 1) GO TO 270
         F = 0.0D0
C
        DO 240 J = 1, L
            G = 0.0D0
C     ********** FORM ELEMENT OF A*U **********
            DO 180 K = 1, J
  180       G = G + A(J,K) * A(I,K)
C
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
  200       G = G + A(K,J) * A(I,K)
C     ********** FORM ELEMENT OF P **********
  220       E(J) = G / H
            F = F + E(J) * A(I,J)
  240    CONTINUE
C
         H = F / (H + H)
C     ********** FORM REDUCED A **********
         DO 260 J = 1, L
            F = A(I,J)
            G = E(J) - H * F
            E(J) = G
C
            DO 260 K = 1, J
               A(J,K) = A(J,K) - F * E(K) - G * A(I,K)
  260    CONTINUE
C
  270    DO 280 K = 1, L
  280    A(I,K) = SCALE * A(I,K)
C
  290    H = D(I)
         D(I) = A(I,I)
         A(I,I) = H
  300 CONTINUE
C
      RETURN
C     ********** LAST CARD OF TRED1 **********
      END
C
C
      SUBROUTINE TRED2(NM,N,A,D,E,Z)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(NM,N),D(N),E(N),Z(NM,N)
C     ****************************************************************
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A
C     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT :
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX,
C
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C     ON OUTPUT :
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO,
C
C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
C          PRODUCED IN THE REDUCTION,
C
C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     *************************************************************
      DO 100 I = 1, N
         DO 100 J = 1, I
            Z(I,J) = A(I,J)
  100 CONTINUE
C
      IF (N .EQ. 1) GO TO 320
C     ********** FOR I=N STEP -1 UNTIL 2 DO -- **********
      DO 300 II = 2, N
         I = N + 2 - II
         L = I - 1
         H = 0.0D0
         SCALE = 0.0D0
         IF (L .LT. 2) GO TO 130
C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
         DO 120 K = 1, L
  120    SCALE = SCALE + DABS(Z(I,K))
C
         IF (SCALE .NE. 0.0D0) GO TO 140
  130    E(I) = Z(I,L)
         GO TO 290
C
  140    DO 150 K = 1, L
            Z(I,K) = Z(I,K) / SCALE
            H = H + Z(I,K) * Z(I,K)
  150    CONTINUE
C
         F = Z(I,L)
         G = -DSIGN(DSQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         Z(I,L) = F - G
         F = 0.0D0
C
         DO 240 J = 1, L
            Z(J,I) = Z(I,J) / H
            G = 0.0D0
C     ********** FORM ELEMENT OF A*U **********
            DO 180 K = 1, J
  180       G = G + Z(J,K) * Z(I,K)
C
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
  200       G = G + Z(K,J) * Z(I,K)
C     ********** FORM ELEMENT OF P **********
  220       E(J) = G / H
            F = F + E(J) * Z(I,J)
  240    CONTINUE
C
         HH = F / (H + H)
C     ********** FORM REDUCED A **********
         DO 260 J = 1, L
            F = Z(I,J)
            G = E(J) - HH * F
            E(J) = G
C
            DO 260 K = 1, J
               Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
  260    CONTINUE
C
  290    D(I) = H
  300 CONTINUE
C
  320 D(1) = 0.0D0
      E(1) = 0.0D0
C     ********** ACCUMULATION OF TRANSFORMATION MATRICES **********
      DO 500 I = 1, N
         L = I - 1
         IF (D(I) .EQ. 0.0D0) GO TO 380
C
         DO 360 J = 1, L
            G = 0.0D0
C
            DO 340 K = 1, L
  340       G = G + Z(I,K) * Z(K,J)
C
            DO 360 K = 1, L
               Z(K,J) = Z(K,J) - G * Z(K,I)
  360    CONTINUE
C
  380    D(I) = Z(I,I)
         Z(I,I) = 1.0D0
         IF (L .LT. 1) GO TO 500
C
         DO 400 J = 1, L
            Z(I,J) = 0.0D0
            Z(J,I) = 0.0D0
  400    CONTINUE
C
  500 CONTINUE
C
      RETURN
C     ********** LAST CARD OF TRED2 **********
      END
C
C
      SUBROUTINE TQLRAT(N,D,E2,IERR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(N),E2(N)
      REAL*8 MACHEP
C     *****************************************************************
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
C     ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
C
C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
C     TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD.
C
C     ON INPUT :
C
C        N IS THE ORDER OF THE MATRIX,
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX'
C
C        E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE
C          INPUT MATRIX IN ITS LAST N-1 POSITIONS.  E2(1) IS ARBITRARY.
C
C      ON OUTPUT :
C
C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
C          THE SMALLEST EIGENVALUES,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     **************************************************************
C
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
C
      MACHEP = 2.D0**(-26)
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E2(I-1) = E2(I)
C
      F = 0.0D0
      B = 0.0D0
      C = 0.0D0
      E2(N) = 0.0D0
C
      DO 290 L = 1, N
         J = 0
         H = MACHEP * (DABS(D(L)) + DSQRT(E2(L)))
         IF (B .GT. H) GO TO 105
         B = H
         C = B * B
C     ********** LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT **********
  105    DO 110 M = L, N
            IF (E2(M) .LE. C) GO TO 120
C     ********** E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP **********
  110    CONTINUE
         WRITE(6,*)' **** FATAL ERROR IN TQLRAT **** '
         WRITE(6,*)' **** FALLEN THROUGH BOTTOM OF LOOP 110 *** '
         STOP
C
  120    IF (M .EQ. L) GO TO 210
  130    IF (J .EQ. 30) GO TO 1000
         J = J + 1
C     ********** FORM SHIFT **********
         L1 = L + 1
         S = DSQRT(E2(L))
         G = D(L)
         P = (D(L1) - G) / (2.0D0 * S)
         R = DSQRT(P*P+1.0D0)
         D(L) = S / (P + DSIGN(R,P))
         H = G - D(L)
C
         DO 140 I = L1, N
  140    D(I) = D(I) - H
C
         F = F + H
C     ********** RATIONAL QL TRANSFORMATION **********
         G = D(M)
         IF (G .EQ. 0.0D0) G = B
         H = G
         S = 0.0D0
         MML = M - L
C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
         DO 200 II = 1, MML
            I = M - II
            P = G * H
            R = P + E2(I)
            E2(I+1) = S * R
            S = E2(I) / R
            D(I+1) = H + S * (H + D(I))
            G = D(I) - E2(I) / G
            IF (G .EQ. 0.0D0) G = B
            H = G * P / R
  200    CONTINUE
C
         E2(L) = S * G
         D(L) = H
C     ********** GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST **********
         IF (H .EQ. 0.0D0) GO TO 210
         IF (DABS(E2(L)) .LE.DABS(C/H)) GO TO 210
         E2(L) = H * E2(L)
         IF (E2(L) .NE. 0.0D0) GO TO 130
  210    P = D(L) + F
C     ********** ORDER EIGENVALUES **********
         IF (L .EQ. 1) GO TO 250
C     ********** FOR I=L STEP -1 UNTIL 2 DO -- *********
         DO 230 II = 2, L
            I = L + 2 - II
            IF (P .GE. D(I-1)) GO TO 270
            D(I) = D(I-1)
  230    CONTINUE
C
  250    I = 1
  270    D(I) = P
  290 CONTINUE
C
      GO TO 1001
C     ********** SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS **********
 1000 IERR = L
 1001 RETURN
C     ********** LAST CARD OF TQLRAT **********
      END
C
C
      SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(N),E(N),Z(NM,N)
      REAL*8 MACHEP
C     **************************************************************
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
C     WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
C
C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
C     FULL MATRIX TO TRIDIAGONAL FORM.
C
C     ON INPUT :
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX,
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
C          THE IDENTITY MATRIX.
C
C      ON OUTPUT :
C
C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
C          UNORDERED FOR INDICES 1,2,...,IERR-1,
C
C        E HAS BEEN DESTROYED,
C
C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
C          EIGENVALUES,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ***************************************************************
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.

      MACHEP = 2.D0**(-26)
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      F = 0.0D0
      B = 0.0D0
      E(N) = 0.0D0
C
      DO 240 L = 1, N
         J = 0
         H = MACHEP * (DABS(D(L)) + DABS(E(L)))
         IF (B .LT. H) B = H
C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
         DO 110 M = L, N
            IF (DABS(E(M)) .LE. B) GO TO 120
C     ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP **********
  110    CONTINUE
  120    IF (M .EQ. L) GO TO 220
  130    IF (J .EQ. 30) GO TO 1000
         J = J + 1
C     ********** FORM SHIFT **********
         L1 = L + 1
         G = D(L)
         P = (D(L1) - G) / (2.0D0 * E(L))
         R = DSQRT(P*P+1.0D0)
         D(L) = E(L) / (P + DSIGN(R,P))
         H = G - D(L)
         DO 140 I = L1, N
  140    D(I) = D(I) - H
C
         F = F + H
C     ********** QL TRANSFORMATION **********
         P = D(M)
         C = 1.0D0
         S = 0.0D0
         MML = M - L
C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
         DO 200 II = 1, MML
            I = M - II
            G = C * E(I)
            H = C * P
            IF (DABS(P) .LT. DABS(E(I))) GO TO 150
            C = E(I) / P
            R = DSQRT(C*C+1.0D0)
            E(I+1) = S * P * R
            S = C / R
            C = 1.0D0 / R
            GO TO 160
  150       C = P / E(I)
            R = DSQRT(C*C+1.0D0)
            E(I+1) = S * E(I) * R
            S = 1.0D0 / R
            C = C * S
  160       P = C * D(I) - S * G
            D(I+1) = H + S * (C * G + S * D(I))
C     ********** FORM VECTOR **********
            DO 180 K = 1, N
               H = Z(K,I+1)
               Z(K,I+1) = S * Z(K,I) + C * H
               Z(K,I) = C * Z(K,I) - S * H
  180       CONTINUE
  200    CONTINUE
         E(L) = S * P
         D(L) = C * P
         IF (DABS(E(L)) .GT. B) GO TO 130
  220    D(L) = D(L) + F
  240 CONTINUE
C     ********** ORDER EIGENVALUES AND EIGENVECTORS **********
      DO 300 II = 2, N
         I = II - 1
         K = I
         P = D(I)
C
         DO 260 J = II, N
            IF (D(J) .GE. P) GO TO 260
            K = J
            P = D(J)
  260    CONTINUE
C
         IF (K .EQ. I) GO TO 300
         D(K) = D(I)
         D(I) = P
         DO 280 J = 1, N
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P
  280    CONTINUE
  300 CONTINUE
      GO TO 1001
C     ********** SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS **********
 1000 IERR = L
 1001 RETURN
C     ********** LAST CARD OF TQL2 **********
      END

max_stars_repo_name: jfecroft/DOS

max_issues_repo_name: jfecroft/DOS

max_forks_repo_name: jfecroft/DOS


text:
      subroutine coefsdk(s12,mtsq,ct,cv,c1)
      implicit none
      include 'constants.f'
      include 'epinv.f'
      include 'epinv2.f'
      include 'scheme.f'
      include 'scale.f'
      include 'alfacut.f'
      include 'includect.f'
      double precision cv,ct,c1,s12,mtsq,ddilog,rsq,omrsq,eta,Kfun,
     & wlog,rlog,mulog

      if (scheme .eq.'dred') then
C------        eta=0 4d-hel
         eta=0d0
      elseif (scheme .eq. 'tH-V') then
C------       eta=1 t'Hooft Veltman
         eta=1d0
      endif

      rsq=s12/mtsq
      omrsq=1d0-rsq
      wlog=dlog(omrsq)
      rlog=dlog(rsq)

c--- epinv stands for (4*pi)^ep/Gamma(1-ep)/ep  (as usual)
c----  ct is the integrated counter-term, including alpha-dependence,
c----  see Eq. (10) of arXiv:1102.1967
      mulog=dlog(musq/mtsq)
      if (includect) then
      ct=(epinv2*epinv+epinv*mulog+0.5d0*mulog**2)
     & +(epinv+mulog)*(2.5d0-2d0*wlog)
     & +25d0/4d0+0.5d0*(1d0/omrsq**2-8d0/omrsq+7d0)*rlog
     & +0.5d0/omrsq+2d0*ddilog(omrsq)-5d0*pisqo6
     & -5d0*wlog+2d0*wlog**2+eta/2d0
     & -2d0*dlog(aff)**2-(3.5d0-4d0*aff+aff**2/2d0)*dlog(aff)
     & +2d0*(1d0-aff)*rsq/omrsq*dlog(rsq/(1d0-aff+rsq*aff))
      else
      ct=0d0
      endif

C---- this routine has been constructed from
C---- %\cite{Gottschalk:1980rv}
C---- \bibitem{Gottschalk:1980rv}
C---- T.~Gottschalk,
C---- %``Chromodynamic Corrections To Neutrino Production Of Heavy Quarks,''
C---- Phys.\ Rev.\ D {\bf 23}, 56 (1981).
C---- %%CITATION = PHRVA,D23,56;%%
C----- Adapted from Eqs.(A8,A9)

      Kfun=1d0/rsq*wlog
      c1=2d0*Kfun
      cv=-(epinv2*epinv+epinv*mulog+0.5d0*mulog**2)
     & -(epinv+mulog)*(2.5d0-2d0*wlog)
     & -0.5d0*(11d0+eta)-pisqo6-2d0*ddilog(rsq)
     &  +3d0*wlog-2d0*wlog**2-Kfun

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      IVALUE = 0
   10 IVALUE=IVALUE+1
      WRITE 301,IVALUE
  301 FORMAT(I5)
      IF(IVALUE-IVALUE/6*6) 10,20,10
   20 STOP
      END

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
c-------------------------------------------------------------------------------
c  vk 06/08 There are 3 options of carry over states treatment
c  1. irdco = 0: there is no need to input snow states from input deck or grids,
c                zero snow states will be assumed; this is default option.
c                User should provide another IRDCO value if
c                grids or input deck states are available and should be used.
c  2. irdco = 1: all snow state are available from grids or input deck.
c  3. irdco = 2: only 7 snow states (we,neghs,liqw,tindex,accmax,sndpt,sntmp)
c                are available from grids or input deck; other snow states will
c                be generated by program. Use this option with caushion because
c                it will not guaranty time consistent states.
c-------------------------------------------------------------------------------
c
C MEMBER CKCO19
C  (from old member FCPIN19)
C
      SUBROUTINE CKCO19(IRDCO,DTHR,PLWHC,SI,ADC,WE,NEGHS,LIQW,
     +                  TINDEX,ACCMAX,SB,SBAESC,SBWS,STORGE,AEADJ,
     +                  NEXLAG,EXLAG,SNDPT,SNTMP)
C.......................................
C     THIS SUBROUTINE CHECKS CARRYOVER VALUES FOR THE 'SNOW-17 '
C        OPERATION.
C.......................................
C     INITIALLY WRITTEN BY... ERIC ANDERSON - HRL   MAY |980
CVK   MODIFIED 4/00 BY V. KOREN TO ADD TWO MORE STATES
C.......................................

      REAL NEGHS,LIQW
      DIMENSION ADC(11), EXLAG(7)
C
C.......................................
C     CHECK IF A SNOW COVER EXISTS.
      IF(IRDCO.EQ.0) GO TO 100
      IF(WE.EQ.0.0) GO TO 100
      GO TO 110
C.......................................
C     SET TO NO SNOW CONDITIONS.
  100 WE=0.0
      NEGHS=0.0
      LIQW=0.0
      TINDEX=0.0
      ACCMAX=0.0
      SB=0.0
      SBAESC=0.0
      SBWS=0.0

CVK  TWO MORE STATES
      SNDPT=0.0
      SNTMP=0.0

  105 STORGE=0.0
      AEADJ=0.0
      DO 101 N=1,NEXLAG
  101 EXLAG(N)=0.0
      RETURN
C.......................................
C     SNOW COVER EXISTS -- MAKE CHECKS
  110 IF(NEGHS.LT.0.0) NEGHS=0.0
      xxs=PLWHC*WE
      IF(LIQW.GT.xxs) LIQW=xxs
cck      IF(LIQW.GT.PLWHC*WE) LIQW=PLWHC*WE
c      IF(TINDEX.GT.0.0) TINDEX=0.0
      xxs=WE+LIQW
      IF(ACCMAX.LT.xxs) ACCMAX=xxs
cck      IF(ACCMAX.LT.(WE+LIQW)) ACCMAX=WE+LIQW

CVK    IF NO SNDPT & SNTMP STATES ENTERED, CALCULATE
CVK    THEM BASED ON OTHER AVAILABLE STATES
      IF(SNDPT .LE. 0.0) SNDPT=0.1*WE/0.2
      IF(SNTMP .EQ. 0. .AND. NEGHS .GT. 0.) SNTMP=TINDEX
      IF(SNTMP .GT. 0.0 .OR. SNTMP .EQ. -999.) SNTMP=TINDEX

c vk 06/08      IF(IRDCO.EQ.2) GO TO 120
      IF(IRDCO.EQ.1) GO TO 120
C
C     COMPUTE CARD 10 VALUES.
      TWE=WE+LIQW
      SBWS=TWE
      AI=ACCMAX
      IF(ACCMAX.GT.SI) AI=SI
      IF(TWE.LT.AI) GO TO 111
      SB=TWE
      SBAESC=1.0
      GO TO 105
  111 R=(TWE/AI)*10.0+1.0
      N=R
      FN=N
      R=R-FN
      SBAESC=ADC(N)+(ADC(N+1)-ADC(N))*R
      IF(SBAESC.GT.1.0) SBAESC=1.0
      SB=TWE+0.2*DTHR
      GO TO 105
C
C     CARD 10 VALUES READ IN - NO CHECKS
  120 CONTINUE
      RETURN
      END

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


text:
C$Procedure      ZZHLP028 ( private help text )

      SUBROUTINE ZZHLP028 ( BEGIN, FINISH, TEXT )

C$ Abstract
C
C     Fill out a portion of the help text needed by percy.
C
C     Private routine intended solely for the support of Inspekt
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     PRIVATE
C
C$ Declarations

      IMPLICIT NONE
      INTEGER               BEGIN ( * )
      INTEGER               FINISH( * )
      CHARACTER*(*)         TEXT  ( * )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     BEGIN      O   Indexes of begins of text help
C     FINISH     O   Indexes of ends of text help
C     TEXT       O   A block of text help.
C
C$ Exceptions
C
C     Error free.
C
C$ Particulars
C
C     This routine simply fills begin and end markers as well
C     as actual text for a block of help text for percy.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber      (JPL)
C
C$ Literature_References
C
C     None.
C
C$ Version
C
C-    Inspekt Version 1.0.0, 1-AUG-1997 (WLT)
C
C
C-&

      INTEGER               I
      INTEGER               J

      J = FINISH ( 1 )
      I = BEGIN  ( 1 )

      FINISH(1) = J
      BEGIN (1) = I

      CALL REPMC ( TEXT(1), '*', '*', TEXT(1) )
      TEXT( 2603 ) = 'Quit Help'
      TEXT( 2604 ) = 'Help'
      FINISH( 66 ) = 2604

      BEGIN ( 67 ) = 2605
      TEXT( 2605 ) = 'Listed here is the formal syntax for ev'
     .//             'ery'
      TEXT( 2606 ) = 'command recognized by Inspekt.'
      TEXT( 2607 ) = ' '
      TEXT( 2608 ) = '@exliteral'
      TEXT( 2609 ) = 'DEFINE   @name'
      TEXT( 2610 ) = 'DEFINE   @name (1:)@word'
      TEXT( 2611 ) = 'UNDEFINE @name'
      TEXT( 2612 ) = ' '
      TEXT( 2613 ) = 'DISCARD'
      TEXT( 2614 ) = ' '
      TEXT( 2615 ) = 'DO @word'
      TEXT( 2616 ) = ' '
      TEXT( 2617 ) = 'ECHO'
      TEXT( 2618 ) = 'NO ECHO'
      TEXT( 2619 ) = ' '
      TEXT( 2620 ) = 'EDIT @word'
      TEXT( 2621 ) = ' '
      TEXT( 2622 ) = 'EXIT'
      TEXT( 2623 ) = 'HELP'
      TEXT( 2624 ) = 'LOAD EK @word'
      TEXT( 2625 ) = 'LOAD LEAPSECONDS @word'
      TEXT( 2626 ) = 'LOAD SCLK KERNEL @word'
      TEXT( 2627 ) = ' '
      TEXT( 2628 ) = 'RECALL'
      TEXT( 2629 ) = 'RECALL @word'
      TEXT( 2630 ) = ' '
      TEXT( 2631 ) = 'SAVE TO @word'
      TEXT( 2632 ) = ' '
      TEXT( 2633 ) = 'SAMPLE @int(1:)  SELECT (1:100)@word'
      TEXT( 2634 ) = '                 FROM   (1:)@word'
      TEXT( 2635 ) = '                (0:1){ WHERE    (1:)@wo'
     .//             'rd    }'
      TEXT( 2636 ) = '                (0:1){ ORDER BY (1:100)'
     .//             '@word }'
      TEXT( 2637 ) = ' '
      TEXT( 2638 ) = ' '
      TEXT( 2639 ) = 'SAMPLE          (1:1){ FIRST @int(1:)'
      TEXT( 2640 ) = '                     | LAST  @int(1:)'
      TEXT( 2641 ) = '                     }'
      TEXT( 2642 ) = '                 SELECT (1:100)@word'
      TEXT( 2643 ) = '                 FROM   (1:)@word'
      TEXT( 2644 ) = '                (0:1){ WHERE    (1:)@wo'
     .//             'rd    }'
      TEXT( 2645 ) = '                (0:1){ ORDER BY (1:100)'
     .//             '@word }'
      TEXT( 2646 ) = ' '
      TEXT( 2647 ) = 'SAMPLE @int(1:) (1:1){ UP TO       @int'
     .//             '(0:100) EVERY @int(1:)'
      TEXT( 2648 ) = '                     | UP TO       @int'
     .//             '(0:100)'
      TEXT( 2649 ) = '                     | STARTING AT @int'
     .//             '(0:100) EVERY @int(1:)'
      TEXT( 2650 ) = '                     | STARTING AT @int'
     .//             '(0:100)'
      TEXT( 2651 ) = '                     | CENTER   AT @int'
     .//             '(0:100) EVERY @int(1:)'
      TEXT( 2652 ) = '                     | CENTER   AT @int'
     .//             '(0:100)'
      TEXT( 2653 ) = '                     | FROM @int(0:100)'
     .//             ' TO   @int(0:100)'
      TEXT( 2654 ) = '                     }'
      TEXT( 2655 ) = '                 SELECT (1:100)@word'
      TEXT( 2656 ) = '                 FROM   (1:)@word'
      TEXT( 2657 ) = '                (0:1){ WHERE    (1:)@wo'
     .//             'rd    }'
      TEXT( 2658 ) = '                (0:1){ ORDER BY (1:100)'
     .//             '@word }'
      TEXT( 2659 ) = ' '
      TEXT( 2660 ) = 'SELECT (1:100)@word FROM (1:)@word'
      TEXT( 2661 ) = '                    (0:1){ WHERE    (1:'
     .//             ')@word    }'
      TEXT( 2662 ) = '                    (0:1){ ORDER BY (1:'
     .//             '100)@word }'
      TEXT( 2663 ) = ' '
      TEXT( 2664 ) = 'SET AUTOADJUST (1:1){ OFF | ASK | ON }'
      TEXT( 2665 ) = ' '
      TEXT( 2666 ) = 'SET COLUMN @word FORMAT (1:)@word }'
      TEXT( 2667 ) = 'SET COLUMN @word HEADING (1:)@word'
      TEXT( 2668 ) = 'SET COLUMN @word JUSTIFICATION (1:1){ L'
     .//             'EFT  | JUSTIFICATION }'
      TEXT( 2669 ) = 'SET COLUMN @word WIDTH @int'
      TEXT( 2670 ) = ' '
      TEXT( 2671 ) = 'SET DEFAULT TIME     FORMAT (1:)@word'
      TEXT( 2672 ) = 'SET DEFAULT INTEGER  FORMAT @word'
      TEXT( 2673 ) = 'SET DEFAULT FLOATING FORMAT @word'
      TEXT( 2674 ) = ' '
      TEXT( 2675 ) = 'SET DELUGE WARNING @int(1:)'
      TEXT( 2676 ) = ' '
      TEXT( 2677 ) = 'SET EDITOR (1:)@word'
      TEXT( 2678 ) = ' '
      TEXT( 2679 ) = 'SET FORMAT (0:1){ SPACED | MARKED } TAB'
     .//             'ULAR (0:1){ PRESERVED }'
      TEXT( 2680 ) = 'SET FORMAT FLAGGED (0:1){ PRESERVED }'
      TEXT( 2681 ) = 'SET FORMAT MARK @word(%)'
      TEXT( 2682 ) = 'SET FORMAT VERBATIM'
      TEXT( 2683 ) = ' '
      TEXT( 2684 ) = 'SET FORMAT DELIMITED (0:1){ PRESERVED }'
      TEXT( 2685 ) = '                     (0:2){ DELIMITER @'
     .//             'word(SPACE|%)'
      TEXT( 2686 ) = '                          | QUOTE @word'
     .//             '(%) }'
      TEXT( 2687 ) = ' '
      TEXT( 2688 ) = 'SET HEADER FREQUENCY (1:1){ 0'
      TEXT( 2689 ) = '                          | FIRST | 1ST'
      TEXT( 2690 ) = '                          | ALL'
      TEXT( 2691 ) = '                          | EVERY @int('
     .//             '2:) }'
      TEXT( 2692 ) = ' '
      TEXT( 2693 ) = 'SET HELP WAIT'
      TEXT( 2694 ) = 'SET HELP NO WAIT'
      TEXT( 2695 ) = ' '
      TEXT( 2696 ) = 'SET PAGE HEIGHT @int'
      TEXT( 2697 ) = 'SET PAGE TITLE (1:)@word'
      TEXT( 2698 ) = 'SET PAGE WIDTH @int'
      TEXT( 2699 ) = ' '

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!----------------------------------------------------------------------------
! This file is part of UCLALES.
!
! UCLALES is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3 of the License, or
! (at your option) any later version.
!
! UCLALES is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1999-2008, Bjorn B. Stevens, Dep't Atmos and Ocean Sci, UCLA
!----------------------------------------------------------------------------
!
module forc

  use defs, only      : cp
  use radiation, only : d4stream
  use stat, only : sflg
  implicit none

  ! these are now all namelist parameters
  character (len=10) :: case_name = 'none'
  character (len=50) :: radsounding = 'datafiles/dsrt.lay'  ! Juha: Added so the radiation background sounding can be given
                                                            ! from the NAMELIST
  REAL :: sfc_albedo = 0.05
  REAL  :: div = 0.
  LOGICAL :: useMcICA = .TRUE.
  LOGICAL :: RadConstPress = .FALSE. ! Keep constant pressure levels
  INTEGER :: RadPrecipBins = 0 ! Add precipitation bins to cloud water (for level 3 and up)

contains
  !
  ! -------------------------------------------------------------------
  ! subroutine forcings:  calls the appropriate large-scale forcings
  !
  subroutine forcings(time_in, cntlat, sst)

    use grid, only: nxp, nyp, nzp, zm, zt, dzt, dzm, dn0, iradtyp, a_rc     &
         , a_rflx, a_sflx, albedo, a_tt, a_tp, a_rt, a_rp, a_pexnr, a_temp  &
         , a_rv, a_rpp, a_npp, CCN, pi0, pi1, level, a_ut, a_up, a_vt, a_vp, &
         a_ncloudp, a_nprecpp, a_mprecpp, a_ri, a_nicep, a_nsnowp, &
         a_fus, a_fds, a_fuir, a_fdir

    USE mo_submctl, ONLY : nspec,nprc,ira,fra

    use mpi_interface, only : myid, appl_abort

    real, optional, intent (in) :: time_in, cntlat, sst

    real :: xka, fr0, fr1, xref1, xref2
    REAL :: znc(nzp,nxp,nyp), zrc(nzp,nxp,nyp), zni(nzp,nxp,nyp), zri(nzp,nxp,nyp)

    ! DIVERGENCE GIVEN FROM NAMELIST
    if (trim(case_name) == 'atex') then
       xka = 130.
       fr0 = 74.
       fr1 = 0.
       div = 0.
    else
       xka = 85.
       fr0 = 70.
       fr1 = 22.
       !div = 3.75e-6
    end if

    if (trim(case_name)=='ascos') THEN
        ! Full radiation calculations when saving data (stat/sflg=.TRUE. when saving)
        useMcICA=.NOT.sflg
    endif

    select case(iradtyp)
    case (1)
       ! No radiation, just large-scale forcing.
       ! Note, there's a slight discrepancy between lev 1-3 and lev 4 with a_rp
       ! (total water vs vapour): it perhaps doesn't make sence to change the
       ! tendency of condesated water due to subsidence in level4 and for level
       ! 1-3 total mixing ratio is the only prognostic water variable.
       ! -------------------------------------------------
       IF ( case_name /= 'none' ) THEN
          call case_forcing(nzp,nxp,nyp,zt,dzt,dzm,div,a_tp,a_rp,a_tt,a_rt)
       END IF

    case (2)
       ! Some original special cases....
       ! ---------------------------------
       IF ( level >= 4) THEN
          IF(myid == 0) WRITE(*,*) 'FORCING: selection not implemented for level 4 or 5, iradtyp = ',iradtyp
          CALL appl_abort(0)
       END IF

       select case(level)
       case(1)
          call smoke_rad(nzp, nxp, nyp, dn0, a_rflx, zm, dzt,a_tt,a_rp)
       case(2)
          call gcss_rad(nzp, nxp, nyp, xka, fr0, fr1, div, a_rc, dn0,     &
               a_rflx, zt, zm, dzt, a_tt, a_tp, a_rt, a_rp)

       end select
       if (trim(case_name) == 'atex') call case_forcing(nzp, nxp, nyp,    &
            zt, dzt, dzm, div, a_tp, a_rp, a_tt, a_rt)
    case (3)
       ! Radiation + large-scale forcing
       ! -------------------------------------
       if (present(time_in) .and. present(cntlat) .and. present(sst)) then

          IF (level <= 3) THEN
             znc(:,:,:) = CCN
             zrc(:,:,:) = a_rc(:,:,:) ! Cloud water only
             IF (level == 3 .AND. RadPrecipBins > 0) THEN ! Add precipitation (all or nothing)
                znc(:,:,:) = znc(:,:,:) + a_npp(:,:,:)
                zrc(:,:,:) = zrc(:,:,:) + a_rpp(:,:,:)
             ENDIF
             call d4stream(nzp, nxp, nyp, cntlat, time_in, sst, sfc_albedo, &
                  dn0, pi0, pi1, dzt, a_pexnr, a_temp, a_rv, zrc, znc, a_tt,  &
                  a_rflx, a_sflx, a_fus, a_fds, a_fuir, a_fdir, albedo, radsounding=radsounding, &
                  useMcICA=useMcICA, ConstPrs=RadConstPress)

          ELSE IF (level == 4) THEN
             znc(:,:,:) = SUM(a_ncloudp(:,:,:,:),DIM=4) ! Cloud droplets
             zrc(:,:,:) = a_rc(:,:,:) ! Cloud and aerosol water
             IF (RadPrecipBins>0) THEN ! Add precipitation bins
                ! Water is the last species (nspec+1)
                zrc(:,:,:) = zrc(:,:,:) + SUM(a_mprecpp(:,:,:,nspec*nprc+ira:nspec*nprc+min(RadPrecipBins,fra)),DIM=4)
                znc(:,:,:) = znc(:,:,:) + SUM(a_nprecpp(:,:,:,ira:min(RadPrecipBins,fra)),DIM=4)
             ENDIF
             CALL d4stream(nzp, nxp, nyp, cntlat, time_in, sst, sfc_albedo, &
                  dn0, pi0, pi1, dzt, a_pexnr, a_temp, a_rp, zrc, znc, a_tt,  &
                  a_rflx, a_sflx, a_fus, a_fds, a_fuir, a_fdir, albedo, radsounding=radsounding, &
                  useMcICA=useMcICA, ConstPrs=RadConstPress)

          ELSE IF (level == 5) THEN
             znc(:,:,:) = SUM(a_ncloudp(:,:,:,:),DIM=4) ! Cloud droplets
             zrc(:,:,:) = a_rc(:,:,:) ! Cloud and aerosol water
             IF (RadPrecipBins>0) THEN ! Add precipitation bins
                ! Water is the last species (nspec+1)
                zrc(:,:,:) = zrc(:,:,:) + SUM(a_mprecpp(:,:,:,nspec*nprc+ira:nspec*nprc+min(RadPrecipBins,fra)),DIM=4)
                znc(:,:,:) = znc(:,:,:) + SUM(a_nprecpp(:,:,:,ira:min(RadPrecipBins,fra)),DIM=4)
             ENDIF
             zni(:,:,:) = SUM(a_nicep(:,:,:,:),DIM=4) ! Ice
             zri(:,:,:) = a_ri(:,:,:) ! Ice (no aerosol ice?)
             CALL d4stream(nzp, nxp, nyp, cntlat, time_in, sst, sfc_albedo, &
                  dn0, pi0, pi1, dzt, a_pexnr, a_temp, a_rp, zrc, znc, a_tt,  &
                  a_rflx, a_sflx, a_fus, a_fds, a_fuir, a_fdir, albedo, ice=zri,nice=zni,radsounding=radsounding, &
                  useMcICA=useMcICA, ConstPrs=RadConstPress)

          END IF

          IF ( case_name /= 'none') THEN
             CALL case_forcing(nzp,nxp,nyp,zt,dzt,dzm,div,a_tp,a_rp,a_tt,a_rt)
          END IF

       else
          if (myid == 0) print *, '  ABORTING: inproper call to radiation'
          call appl_abort(0)
       end if
    case (4)
       call bellon(nzp, nxp, nyp, a_rflx, a_sflx, zt, dzt, dzm, a_tt, a_tp&
            ,a_rt, a_rp, a_ut, a_up, a_vt, a_vp)
    end select

  end subroutine forcings

  !
  ! -------------------------------------------------------------------
  ! subroutine gcss_rad:  call simple radiative parameterization and
  ! simultaneously update fields due to vertical motion as given by div
  !
  subroutine gcss_rad(n1,n2,n3,xka,fr0,fr1,div,rc,dn0,flx,zt,zm,dzt,   &
       tt,tl,rtt,rt)

    integer, intent (in):: n1,n2, n3
    real, intent (in)   :: xka, fr0, fr1, div
    real, intent (in)   :: zt(n1),zm(n1),dzt(n1),dn0(n1),rc(n1,n2,n3),   &
         tl(n1,n2,n3),rt(n1,n2,n3)
    real, intent (inout):: tt(n1,n2,n3),rtt(n1,n2,n3)
    real, intent (out)  :: flx(n1,n2,n3)

    integer :: i, j, k, km1, kp1, ki
    real    :: lwp(n2,n3), fact

    lwp=0.
    do j=3,n3-2
       do i=3,n2-2
          ki = n1
          do k=1,n1
             km1=max(1,k-1)
             lwp(i,j)=lwp(i,j)+max(0.,rc(k,i,j)*dn0(k)*(zm(k)-zm(km1)))
             flx(k,i,j)=fr1*exp(-1.*xka*lwp(i,j))
             if ( (rc(k,i,j) > 0.01e-3) .and. (rt(k,i,j) >= 0.008) ) ki=k
          enddo

          fact = div*cp*dn0(ki)
          do k=2,n1
             km1=max(2,k-1)
             lwp(i,j)=lwp(i,j)-max(0.,rc(k,i,j)*dn0(k)*(zm(k)-zm(k-1)))
             flx(k,i,j)=flx(k,i,j)+fr0*exp(-1.*xka*lwp(i,j))
             if (zm(k) > zm(ki) .and. ki > 1 .and. fact > 0.) then
                flx(k,i,j)=flx(k,i,j) + fact*(0.25*(zm(k)-zm(ki))**1.333 + &
                  zm(ki)*(zm(k)-zm(ki))**0.333333)
             end if
             tt(k,i,j) =tt(k,i,j)-(flx(k,i,j)-flx(km1,i,j))*dzt(k)/(dn0(k)*cp)
          enddo
          !
          ! subsidence
          !
          if (div /= 0.) then
             do k=2,n1-2
                kp1 = k+1
                tt(k,i,j) = tt(k,i,j) + &
                        div*zt(k)*(tl(kp1,i,j)-tl(k,i,j))*dzt(k)
                rtt(k,i,j)=rtt(k,i,j) + &
                        div*zt(k)*(rt(kp1,i,j)-rt(k,i,j))*dzt(k)
             end do
          end if
       enddo
    enddo

  end subroutine gcss_rad
  !
  ! -------------------------------------------------------------------
  ! subroutine smoke_rad:  call simple radiative parameterization for
  ! the smoke cloud
  !
  subroutine smoke_rad(n1,n2,n3,dn0,flx,zm,dzt,tt,rt)

    integer, intent (in):: n1,n2, n3
    real, intent (in)   :: zm(n1),dzt(n1),dn0(n1),rt(n1,n2,n3)
    real, intent (inout):: tt(n1,n2,n3)
    real, intent (out)  :: flx(n1,n2,n3)
    real, parameter     :: xka= 50.0, fr0=60.0

    integer :: i,j,k, km1, ki
    real    :: smoke(n2,n3)

    smoke=0.
    do j=3,n3-2
       do i=3,n2-2
          ki = n1
          do k=1,n1
             km1=max(1,k-1)
             smoke(i,j)=smoke(i,j)+max(0.,rt(k,i,j)*dn0(k)*(zm(k)-zm(km1)))
          enddo

          do k=2,n1
             km1=max(2,k-1)
             smoke(i,j)=smoke(i,j)-max(0.,rt(k,i,j)*dn0(k)*(zm(k)-zm(k-1)))
             flx(k,i,j)=fr0*exp(-1.*xka*smoke(i,j))
             tt(k,i,j) =tt(k,i,j)-(flx(k,i,j)-flx(km1,i,j))*dzt(k)/(dn0(k)*cp)
          enddo
       enddo
    enddo

  end subroutine smoke_rad
  !
  ! -------------------------------------------------------------------
  ! subroutine case_forcing: adjusts tendencies according to a specified
  ! large scale forcing.  Normally case (run) specific.
  !
  subroutine case_forcing(n1,n2,n3,zt,dzt,dzm,zdiv,tl,rt,tt,rtt)

    use mpi_interface, only : pecount, double_scalar_par_sum,myid, appl_abort
    use stat, only : get_zi

    integer, intent (in):: n1,n2, n3
    real, dimension (n1), intent (in)          :: zt, dzt, dzm
    real, intent(in)                           :: zdiv
    real, dimension (n1,n2,n3), intent (in)    :: tl, rt
    real, dimension (n1,n2,n3), intent (inout) :: tt, rtt

    integer :: i,j,k,kp1
    real, dimension (n1) :: sf
    real, parameter :: zmx_sub = 2260. ! originally 2260.

    real (kind=8) :: zig, zil
    real          :: zibar

    zig = 0.0; zil = 0.0; zibar = 0.0
    kp1= 0
    select case (trim(case_name))
    case('default')
       !
       ! User specified divergence used as a simple large scle forcing for moisture and temperature fields
       ! -------------------------------------------------------------------------------------------------
       !
       DO j=3,n3-2
          DO i=3,n2-2
             DO k=2,n1-1
                kp1 = k+1
                tt(k,i,j) = tt(k,i,j) + zdiv*zt(k)*(tl(kp1,i,j)-tl(k,i,j))*dzt(k)
                rtt(k,i,j) = rtt(k,i,j) + zdiv*zt(k)*(rt(kp1,i,j)-rt(k,i,j))*dzt(k)
             END DO
          END DO
       END DO

    case('rico')
       !
       ! calculate subsidence factor (wsub / dz)
       !
       do k=2,n1-2
          if (zt(k) < zmx_sub) then
             sf(k) =  -0.005*zt(k)/zmx_sub
          else
             sf(k) =  -0.005
          end if
          sf(k) = sf(k)*dzt(k)
       end do

       do j=3,n3-2
          do i=3,n2-2
             do k=2,n1-2
                !
                ! subsidence
                !
                kp1 = k+1
                tt(k,i,j)  =  tt(k,i,j) - ( tl(kp1,i,j) - tl(k,i,j) )*sf(k)
                rtt(k,i,j) = rtt(k,i,j) - ( rt(kp1,i,j) - rt(k,i,j) )*sf(k)
                !
                ! temperature advection and radiative cooling
                !
                tt(k,i,j) = tt(k,i,j)  - 2.5/86400.
                !
                ! moisture advection
                !
                if (zt(k) <= 2980.) then
                   rtt(k,i,j) = rtt(k,i,j)  - (1. -  1.3456*zt(k)/2980.)/8.64e7
                else
                   rtt(k,i,j) = rtt(k,i,j)  + .3456/8.64e7
                end if
             enddo
          enddo
       enddo

    case ('bomex')
       !
       ! calculate subsidence factor (wsub / dz)
       !
       do k=2,n1-2
          if (zt(k) < 1500.) then
             sf(k) =  -0.0065*zt(k)/1500.
          else
             sf(k) =  min(0.,-0.0065  + 0.0065*(zt(k)-1500.)/600.)
          end if
          sf(k) = sf(k)*dzt(k)
       end do

       do j=3,n3-2
          do i=3,n2-2
             do k=2,n1-2
                !
                ! temperature advection and radiative cooling
                !
                kp1 = k+1
                if (zt(k) < 1500.) then
                   tt(k,i,j) = tt(k,i,j) - ( tl(kp1,i,j)-tl(k,i,j) )*sf(k) &
                        - 2.315e-5
                else if (zt(k) < 2000.) then
                   tt(k,i,j) = tt(k,i,j) - ( tl(kp1,i,j)-tl(k,i,j) )*sf(k) &
                        - 2.315e-5*(1.- (zt(k)-1500.)*1.e-3)
                end if
                !
                ! moisture advection
                !
                rtt(k,i,j) = rtt(k,i,j) - ( rt(kp1,i,j) - rt(k,i,j) )*sf(k)
                if (zt(k) < 300.) then
                   rtt(k,i,j) = rtt(k,i,j)  - 1.2e-8
                elseif (zt(k) < 500.) then
                   rtt(k,i,j) = rtt(k,i,j)  - 1.2e-8*(1.- (zt(k)-300.)/200.)
                end if
             enddo
          enddo
       enddo
    case ('atex')
       !
       ! calculate subsidence factor (wsub / dz)
       !
       zil = get_zi (n1, n2, n3, 2, rt, dzm, zt, 6.5e-3)
       call double_scalar_par_sum(zil,zig)
       zibar = real(zig/pecount)

       do k=2,n1-2
          if (zt(k) < zibar) then
             sf(k) =  -0.0065*zt(k)/1500.
          else
             sf(k) =  min(0.,-0.0065*(1 - (zt(k)-zibar)/300.))
          end if
          sf(k) = sf(k)*dzt(k)
       end do

       do j=3,n3-2
          do i=3,n2-2
             do k=2,n1-2
                !
                ! temperature advection and radiative cooling
                !
                kp1 = k+1
                if (zt(k) < zibar) then
                   tt(k,i,j) = tt(k,i,j) - ( tl(kp1,i,j)-tl(k,i,j) )*sf(k) &
                        - 2.315e-5*(1. + (1.- zt(k)/zibar)/2.)
                else if (zt(k) < zibar+300.) then
                   tt(k,i,j) = tt(k,i,j) - ( tl(kp1,i,j)-tl(k,i,j) )*sf(k) &
                        - 2.315e-5*(1.- (zt(k)-zibar)/300.)
                end if
                !
                ! moisture advection
                !
                rtt(k,i,j) = rtt(k,i,j) - ( rt(kp1,i,j) - rt(k,i,j) )*sf(k)
                if (zt(k) < zibar) rtt(k,i,j) = rtt(k,i,j)  - 1.5e-8
             enddo
          enddo
       enddo
        !
    case ('ascos')
        ! ASCOS
        ! ---------
        !
        do k=2,n1-2
            ! calculate subsidence factor (wsub / dz)
            sf(k) = -5.0e-6*min(2000.0,zt(k))*dzt(k)
        end do
        !
        do j=3,n3-2
            do i=3,n2-2
                do k=2,n1-2
                    !
                    ! Temperature and humidity advection due to subsidence
                    !
                    kp1 = k+1
                    tt(k,i,j)  =  tt(k,i,j) - ( tl(kp1,i,j) - tl(k,i,j) )*sf(k)
                    rtt(k,i,j) = rtt(k,i,j) - ( rt(kp1,i,j) - rt(k,i,j) )*sf(k)
                enddo
            enddo
        enddo
        !
    case ('barba')
        ! Barbados
        ! -----------
        ! Large scale subsidence: w(z)=w0*(1-exp(z/H)), where w0=7.5 mm/s and H=1000 m.
        ! Radiative cooling rate: 2.5 K/day
        ! No temperature or humidity advection
        !
        ! calculate subsidence factor (wsub / dz)
        do k=2,n1-2
            sf(k) = -7.5e-3*(1.0-exp(-zt(k)/1000.0))*dzt(k)
        end do
        !
        do j=3,n3-2
            do i=3,n2-2
                do k=2,n1-2
                    ! Subsidence
                    kp1 = k+1
                    tt(k,i,j)  =  tt(k,i,j) - ( tl(kp1,i,j) - tl(k,i,j) )*sf(k)
                    rtt(k,i,j) = rtt(k,i,j) - ( rt(kp1,i,j) - rt(k,i,j) )*sf(k)
                    !
                    ! Radiative cooling: 2.5 K/day
                    tt(k,i,j) = tt(k,i,j)  - 2.5/86400.
                enddo
            enddo
        enddo
        !
    CASE ('amazon')
        ! Amazon
        ! --------
        ! - to be added -
        !
    case default
       if (myid == 0) print *, '  ABORTING: inproper call to radiation'
       call appl_abort(0)
    end select

  end subroutine case_forcing
  !
  ! -------------------------------------------------------------------
  ! subroutine bellon_rad:  call simple radiative parameterization
  !
  subroutine bellon(n1,n2,n3,flx,sflx,zt,dzt,dzm,tt,tl,rtt,rt, ut,u,vt,v)

    integer, intent (in) :: n1,n2, n3

    real, dimension (n1), intent (in)            :: zt, dzt, dzm
    real, dimension (n1, n2, n3), intent (inout) :: tt, tl, rtt, rt, ut,u,vt,v
    real,  dimension (n1, n2, n3), intent (out)  :: flx, sflx
    real, parameter      :: w0= 7.5e-3, H=1000., Qrate = 2.5/86400.

    integer :: i,j,k,kp1
    real    :: grad,wk

    do j=3,n3-2
       do i=3,n2-2
          !
          ! subsidence
          !
          flx(1,i,j)  = 0.
          sflx(1,i,j) = 0.
          do k=2,n1-2
             kp1 = k+1
             wk = w0*(1.-exp(-zt(k)/H))
             grad = Qrate/wk
             flx(k,i,j)  = wk*((tl(kp1,i,j)-tl(k,i,j))*dzt(k)-grad)
             sflx(k,i,j) = wk*((rt(kp1,i,j)-rt(k,i,j))*dzt(k)-grad)
             tt(k,i,j) = tt(k,i,j) + flx(k,i,j)
             rtt(k,i,j)=rtt(k,i,j) + &
                  wk*(rt(kp1,i,j)-rt(k,i,j))*dzt(k)
             ut(k,i,j) =  ut(k,i,j) + &
                  wk*(u(kp1,i,j)-u(k,i,j))*dzm(k)
             vt(k,i,j) =  vt(k,i,j) + &
                  wk*(v(kp1,i,j)-v(k,i,j))*dzm(k)
          end do
          flx(n1,  i,j)  = 0.
          flx(n1-1,i,j)  = 0.
          sflx(n1,  i,j) = 0.
          sflx(n1-1,i,j) = 0.
       enddo
    enddo

  end subroutine bellon

end module forc

max_stars_repo_name: huangynj/UCLALES-SALSA

max_issues_repo_name: huangynj/UCLALES-SALSA

max_forks_repo_name: huangynj/UCLALES-SALSA


text:
C NCLFORTSTART
      subroutine mocloops(nyaux,mlon,nlat,kdep,nrx,tlat,lat_aux_grid
     +                   ,rmlak,work1,work2,work3,wmsg,tmp1,tmp2,tmp3)
      implicit none
      integer  nyaux,mlon,nlat,kdep,nrx

      double precision     work1(mlon,nlat,kdep), work2(mlon,nlat,kdep)
     +                    ,work3(mlon,nlat,kdep), wmsg
      double precision     tmp1(nyaux,kdep,2)   , tmp2(nyaux,kdep,2)
     +                    ,tmp3(nyaux,kdep,2)
      double precision     lat_aux_grid(nyaux)  , tlat(mlon,nlat)

      integer  rmlak(mlon,nlat,2)
C NCLEND
      integer  ny, kd, nl, ml, nr
c c c logical  section(mlon,nlat,kdep)

c initilize
      do nr=1,nrx
        do kd=1,kdep
          do ny=1,nyaux
             tmp1(ny,kd,nr) = 0.0d0
             tmp2(ny,kd,nr) = 0.0d0
             tmp3(ny,kd,nr) = 0.0d0
          end do
        end do
      end do

c globe [note: rmlak(:,:,1)]

      do ny=2,nyaux

        do nl=1,nlat
          do ml=1,mlon
             if (tlat(ml,nl).ge.lat_aux_grid(ny-1) .and.
     +           tlat(ml,nl).lt.lat_aux_grid(ny)   .and.
     +           rmlak(ml,nl,1).eq.1) then

                 do kd=1,kdep
                    if (work1(ml,nl,kd).ne.wmsg)  then
                        tmp1(ny,kd,1) = tmp1(ny,kd,1) + work1(ml,nl,kd)
                        tmp2(ny,kd,1) = tmp2(ny,kd,1) + work2(ml,nl,kd)
                        tmp3(ny,kd,1) = tmp3(ny,kd,1) + work3(ml,nl,kd)
                    end if
                 end do

             end if
          end do
        end do

      end do

c atlantic [note: rmlak(:,:,2)]

      do ny=2,nyaux

        do nl=1,nlat
          do ml=1,mlon
             if (tlat(ml,nl).ge.lat_aux_grid(ny-1) .and.
     +           tlat(ml,nl).lt.lat_aux_grid(ny)   .and.
     +           rmlak(ml,nl,2).eq.1) then

                 do kd=1,kdep
                    if (work1(ml,nl,kd).ne.wmsg)  then
                        tmp1(ny,kd,2) = tmp1(ny,kd,2) + work1(ml,nl,kd)
                        tmp2(ny,kd,2) = tmp2(ny,kd,2) + work2(ml,nl,kd)
                        tmp3(ny,kd,2) = tmp3(ny,kd,2) + work3(ml,nl,kd)
                    end if
                 end do

             end if
          end do
        end do

      end do

      return
      end

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text: 20100423 14:53:42 nbsp Go check your org charts. http://iet.ucdavis.edu/orgchart.cfm Here. Mediaworks is STILL used to refer to the arm of ATS that produces media. Users/WilliamLewis, who worked as a tech in the video group until February 2010.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine velocityfactors(R,z,vfacr,vfacz)

      common /blackhole/ bhmass
      real pot

      rsph = sqrt(R*R+z*z)
      call force(R,z,fr,fz,pot)
      frbh = -bhmass*R/(rsph**3.)
      fzbh = -bhmass*z/(rsph**3.)
      vfacr = sqrt(abs((fr+frbh)/fr))
      vfacz = sqrt(abs((fz+fzbh)/fz))
      if(vfacr.lt.1.or.vfacz.lt.1) then
         write(0,*) vfacr,vfacz,fr,fz
         stop
      endif
      return
      end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
      MODULE alclla_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.3E  06:22:09  12/28/06
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
      SUBROUTINE ALCLLA(IBEG, ILAB, ILAST, ILEFT, IPTCS, IRIGH, LBLIN, &
         LLDIM, IMODE)
      INTEGER, INTENT(INOUT) :: LLDIM
      INTEGER, INTENT(IN) :: IMODE
      INTEGER, DIMENSION(:), pointer :: IBEG, ILAB, ILAST, ILEFT, &
                                        IPTCS, IRIGH, LBLIN
!...This routine performs I/O.
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
      REAL FUNCTION tfreez(salinity)
C     Constants taken from Gill, 1982.
C     Author: Robert Grumbine
C     LAST MODIFIED: 1 February 2002.
C          MODIFIED: 21 September 1994.

      IMPLICIT none

      REAL salinity
      REAL a1, a2, a3
      PARAMETER (a1 = -0.0575)
      PARAMETER (a2 =  1.710523E-3)
      PARAMETER (a3 = -2.154996E-4)

      IF (salinity .LT. 0.) THEN
        tfreez = 0.
CD        PRINT *,'tfreez was passed a negative salinity',salinity
        salinity = 1.e-5
      ELSE
        tfreez = salinity*(a1+a2*SQRT(salinity)+a3*salinity)
      ENDIF

      RETURN
      END

max_stars_repo_name: rgrumbine/mmablib

max_issues_repo_name: rgrumbine/mmablib

max_forks_repo_name: rgrumbine/mmablib


text:
The Davis Daytime Toastmasters Community Organizations club is a safe environment for you to become the Public Speaking speaker and leader you want to be.  As of September 2013, we are a mix of beginning and intermediate Toastmasters, with about 8 to 10 members in attendance each week.

The Davis Daytime Toastmasters club is part of http://www.toastmasters.org/ Toastmasters International, a NonProfit Organizations 501(c)(3) nonprofit educational organization based in Rancho Santa Margarita, California. Through club activities and conferences, Toastmasters helps its members improve their communication and leadership skills while fostering selfconfidence and personal growth.

The Daytimers club was chartered by http://www.westyost.com/ West Yost Associates, a local water resources engineering firm, as a training tool for the employees to improve their communication and better themselves.  As such, the meetings are held at the West Yost office.  However, we are a public club, and more than half of our members are from the community at large.  All are welcome to visit and check us out!

The Daytimers would like to thank the Davis Toastmasters Town and Gown Toastmasters for all their help in getting started as a club.  The Town and Gown meets on Tuesday evenings.

Here are some photos of the building we meet in:

20131223 16:09:07 nbsp The Davis Daytime Toastmasters club is moving the January 1, 2014 meeting to January 8, 2014 in observation of the New Years Day holiday.  The club will return to its normal 1st and 3rd Wednesday of each month schedule in February 2014. Users/TTT


text:
      SUBROUTINE PASSF2(IDO,L1,CC,CH,WA1)
      DOUBLE PRECISION CC
      DOUBLE PRECISION CH
      DOUBLE PRECISION WA1
      DOUBLE PRECISION TR2
      DOUBLE PRECISION TI2
      DIMENSION CC(IDO,2,L1),CH(IDO,L1,2),WA1(*)

      IF (IDO.GT.2) GO TO 102
      DO 101 K = 1,L1
          CH(1,K,1) = CC(1,1,K) + CC(1,2,K)
          CH(1,K,2) = CC(1,1,K) - CC(1,2,K)
          CH(2,K,1) = CC(2,1,K) + CC(2,2,K)
          CH(2,K,2) = CC(2,1,K) - CC(2,2,K)
  101 CONTINUE
      RETURN
  102 DO 104 K = 1,L1
          DO 103 I = 2,IDO,2
              CH(I-1,K,1) = CC(I-1,1,K) + CC(I-1,2,K)
              TR2 = CC(I-1,1,K) - CC(I-1,2,K)
              CH(I,K,1) = CC(I,1,K) + CC(I,2,K)
              TI2 = CC(I,1,K) - CC(I,2,K)
              CH(I,K,2) = WA1(I-1)*TI2 - WA1(I)*TR2
              CH(I-1,K,2) = WA1(I-1)*TR2 + WA1(I)*TI2
  103     CONTINUE
  104 CONTINUE
      RETURN
      END


text:
Diagnostic imaging, nuclear medicine, and radiation oncology. CT, open MRI, ultrasound, mammography, routine, and diagnostic radiology were  available at Davis Imaging.

The building was originally occupied by The Wherehouse which provided Video Stores movie rentals as well as dvds, cds and video games. After The Wherehouse closed, it was briefly occupied by Rubios.

Trader Joes wanted to build on the land that Davis Imaging occupied. However, RAS had a 14year lease and had spent more than $1 million in building upgrades. Therefore, it was unwilling to move.  Eventually, they did, however, in part because they lost their Sutter contracts when Sutter decided to do its imaging in house.  (See the TJs page for further details on the history).

Demolition on the building began in April 2010.

See also Medical Services.

20070603 20:04:12 nbsp (Edited) This building used to be The Wherehouse. You could rent videos and buy cds and videos here. Then it became a restaurant. I remember the display in front with cold bottles of beer sitting in a giant display of ice. The restaurant closed soon after it opened. Users/Jedron
 It was a Rubios, actually, but same premise. Was pretty happy when they opened, bummed when they closed, but have since Chuys Taqueria learned to cope. Users/KevinChin

20100318 13:00:10 nbsp This business seems to have closed. Dunno exactly when, but Id say in/around February/late January.  Users/HarrisonM

20100322 00:45:09 nbsp are they closed or just remodeling, there seems to be a lot of fencing around, can anyone verify? Users/theyangster


text:
C    @(#)swpagc.f	20.3 2/13/96
        subroutine swpagc (m,n)
C
      include 'ipfinc/parametr.inc'
      include 'ipfinc/agc.inc'
C
        do 100 i = 1, 14
           itemp = kagc(i,m)
           kagc(i,m) = kagc(i,n)
           kagc(i,n) = itemp
  100   continue

        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c     THE SCALING FUNCTIONS                                            c
c                                                                      c
c     NS_GGH(MASSIN,C)                                                 c
c     NS_GGOS(MASSIN,C)                                                c
c                                                                      c
c     INPUT :                                                          c
c                                                                      c
c       MASSIN(1)  = s                                                 c
c       MASSIN(2)  = tg                                                c
c       MASSIN(3)  = s4                                                c
c       MASSIN(6)  = m1                                                c
c       MASSIN(7)  = mg                                                c
c       MASSIN(9)  = mt                                                c
c       MASSIN(11) = ms                                                c
c       MASSIN(12) = qr                                                c
c       MASSIN(13) = qf                                                c
c       MASSIN(25) = gamma_s                                           c
c       MASSIN(26) = yesno_s4                                          c
c       MASSIN(27) = yesno_s3                                          c
c       MASSIN(30) = kinematic factor s4/(s4-m1**2)                    c
c       MASSIN(4-5,8,10,13-30) not needed                              c
c                                                                      c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c --------------------------------------------------------------------
      real*8 function NS_GGOS(massin,Ar)

      implicit none

      integer    n
      real*8     massin(1:30),Pi,CF,CA,Sn,gs,theta_s3
     &          ,s,s4,m1,m2,ms,s3j,hardfac,mjs2,coegg
     &          ,subs3fac,ggHsub,NSANG(0:9,0:9,-2:2,-2:2)
      complex*16 Ar(4),Arc(4)

      pi = 4.D0*atan(1.D0)
      CF = 4.D0/3.D0
      CA = 3.D0
      Sn = 1.D0/(16.D0 * Pi**2)

c               here the neutralino mass has to be positive
      s   = massin(1)
      s4  = massin(3)
      s3j = massin(4)
      m1  = massin(6)
      m2  = massin(7)
      ms  = massin(11)

c               the s3 regularization
      theta_s3 = 0.D0
      if ((s.gt.(ms+abs(m1))**2).and.(ms.gt.abs(m2))) theta_s3 = 1.D0

c               real kinematics built in
      subs3fac = 1.D0 + m2**2/s3j
      hardfac  = 1.D0 + m1**2/s4

      mjs2 = m2**2 - m1**2

c               the angular functions
      call ANGULAR_ARRAY_NS_GGOS(massin,theta_s3,NSANG)

      do n=1,4
         Arc(n) =  conjg( Ar(n) )
      end do

c               gs**2=4*pi*alpha_s is cut, re-appears as nlo later
      gs = 1.D0

c               coupling matrices, taking into account CC = AA/2
      coegg = gs**4 * real( Ar(1)*Arc(1) )/2.D0 * Sn**2/4.D0

c               wim's form output
c               the prefactor for the scaling functions
      NS_GGOS = 0.D0

      end

c --------------------------------------------------------------------
      real*8 function NS_GGH(massin,Ar)

      implicit none

      integer    n
      real*8     massin(1:30),Pi,CF,CA,Sn,gs,coegg
     &          ,QF,s,s4,m1,m2,mjs2,ms,tj,t1,uj,u1,t1hat,u1hat
     &          ,sptj,spuj,tjpuj,hardfac
     &          ,ggH,NSANG(0:9,0:9,-2:2,-2:2)
      complex*16 Ar(4),Arc(4)

      pi = 4.D0*atan(1.D0)
      CF = 4.D0/3.D0
      CA = 3.D0
      Sn = 1.D0/(16.D0 * Pi**2)

c               here the neutralino mass has to be positive?!?!
      s   = massin(1)
      tj  = massin(2)
      s4  = massin(3)
      m1  = massin(6)
      m2  = massin(7)
      ms  = massin(11)

c               real kinematics built in
      t1 = tj + m2**2 - m1**2

      u1 = s4 - s - tj
      uj = u1 + m1**2 - m2**2

      hardfac = 1.D0 + m1**2/s4
      mjs2 = m2**2 - m1**2

c               some denominators following wim's notes
      t1hat = -( tj*uj - mjs2*s )/(s+tj)
      u1hat = -( tj*uj - mjs2*s )/(s+uj)
      sptj  = s + tj
      spuj  = s + uj
      tjpuj = tj + uj

c               the angular functions
      call ANGULAR_ARRAY_NS_GG(massin,NSANG)

c               the factorization/renormalization scale
      QF = massin(13)

      do n=1,4
         Arc(n) =  conjg( Ar(n) )
      end do

c               gs**2=4*pi*alpha_s is cut, re-appears as nlo later
      gs = 1.D0

c               coupling matrices, taking into account CC = AA/2
      coegg = gs**4 * real( Ar(1)*Arc(1) )/2.D0 * Sn**2/4.D0

c               wim's form output

c               the prefactor for the scaling functions
      NS_GGH = 0.D0

      end

max_stars_repo_name: sliem/docker-hep

max_issues_repo_name: sliem/docker-hep

max_forks_repo_name: sliem/docker-hep


text:
      SUBROUTINE SOFTOC
C
C     SOF TABLE OF CONTENTS ROUTINE
C
C
C     THE CURRENT SUBSTRUCTURE TYPE BIT POSITIONS ARE -
C
C        NO BIT - BASIC SUBSTRUCTURE (EXCEPT IMAGE BIT)
C        BIT 30 - IMAGE SUBSTRUCTURE
C            29 - COMBINED SUBSTRUCTURE
C            28 - GUYAN REDUCTUION SUBSTRUCTURE
C            27 - MODAL REDUCTION SUBSTRUCTURE
C            26 - COMPLEX MODAL REDUCTION SUBSTRUCTURE
C
C     TO ADD A NEW SUBSTRUCTURE TYPE BIT THE FOLLOWING UPDATES ARE
C     REQUIRED.
C
C        1) INCREASE THE DEMENSION OF TYPE.
C        2) INCREASE THE VALUE OF NTYPE IN THE DATA STATEMENT.
C        3) ADD A NEW BCD TYPE VALUE TO THE DATA STATEMENT.
C
C
C     THIS ROUTINE IS CURRENTLY CODED TO HANDLE UP TO 27 SOF ITEMS
C     AUTOMATICALLY.
C     TO INCREASE THIS TO 40 ITEMS PERFORM THE FOLLOWING UPDATES.
C
C        1) CHANGE THE DIMENSION OF HDR TO (40,4)
C        2) CHANGE THE DIMENSION OF ITM TO (40)
C        3) CHANGE THE VALUE OF MAXITM IN THE DATA STATEMENT TO 40
C        4) CHANGE THE INNER GROUPS ON FORMAT 80 TO 39(A1,1X),A1
C        5) CHANGE THE INNER GROUP ON FORMAT 100 TO 39(A1,1X),A1
C
      EXTERNAL        LSHIFT,RSHIFT,ANDF
      INTEGER         AVBLKS,BLANK,DITNSB,BUF,SSNAME(2),ANDF,SS,PS,CS,
     1                HL,RSHIFT,DIRSIZ,SOFSIZ,DITSIZ,NUM(10),BLKSIZ,
     2                HIBLK,FILSIZ,TYPE(5),ITM(27),HDR(27,4)
      CHARACTER       UFM*23,UWM*25,UIM*29,SFM*25,SWM*27
      COMMON /XMSSG / UFM,UWM,UIM,SFM,SWM
      COMMON /MACHIN/ MACH,IHALF
      COMMON /ZZZZZZ/ BUF(1)
      COMMON /SOF   / DIT,DITPBN,DITLBN,DITSIZ,DITNSB,DITBL
      COMMON /SYS   / BLKSIZ,DIRSIZ,SUPSIZ,AVBLKS,HIBLK,IFRST
      COMMON /SOFCOM/ NFILES,FILNAM(10),FILSIZ(10)
      COMMON /SYSTEM/ SYSBUF,NOUT,Z1(6),NLPP,Y(2),LINE,Z2(26),NBPC,NBPW
      COMMON /ITEMDT/ NITEM,ITEM(7,1)
      DATA    TYPE  / 2HB , 2HC , 2HR , 2HM , 2HCM /
      DATA    NUM   / 1H1, 1H2, 1H3, 1H4, 1H5, 1H6 ,1H7, 1H8, 1H9, 1H0 /
      DATA    BLANK / 4H     /
      DATA    IMAGE / 4HI    /
      DATA    NTYPE / 6      /
      DATA    MAXITM/ 27     /
C
      NITM = NITEM
      IF (NITM .LE. MAXITM) GO TO 10
      NITM = MAXITM
      WRITE  (NOUT,6237) SWM,MAXITM
 6237 FORMAT (A27,' 6237, THE SOFTOC ROUTINE CAN HANDLE ONLY',I4,
     1       ' ITEMS.', /34X,'ADDITIONAL ITEMS WILL NOT BE SHOWN')
C
C     SET UP HEADINGS AND MASKS
C
   10 NSHFT = 0
      DO 30 I = 1,4
      DO 20 J = 1,NITM
   20 HDR(J,I) = KLSHFT(ITEM(1,J),NSHFT/NBPC)
      K = NITM + 1
      IF (K .GT. MAXITM) GO TO 30
      DO 22 J = K,MAXITM
   22 HDR(J,I) = BLANK
   30 NSHFT = NSHFT + NBPC
C
      LINE  = NLPP + 1
      M0009 = 1023
      M1019 = LSHIFT(1023,10)
      M2029 = LSHIFT(1023,20)
      IMASK = LSHIFT(1,30)
C
C     LOOP THROUGH DIT
C
      DO 110 JMKN = 1,DITSIZ,2
      I = (JMKN-1)/2 + 1
      CALL FDIT (I,K)
      SSNAME(1) = BUF(K  )
      SSNAME(2) = BUF(K+1)
      IF (SSNAME(1).EQ.BLANK .AND. SSNAME(2).EQ.BLANK) GO TO 110
      CALL FMDI (I,K)
C
C     TEST TYPE BITS IN MDI
C
      DO 40 IT = 2,NTYPE
      IBIT = ANDF(BUF(K+1),LSHIFT(1,31-IT))
      IF (IBIT .NE. 0) GO TO 50
   40 CONTINUE
      IT = 1
   50 IS = ANDF(BUF(K+1),IMASK)
      IM = BLANK
      IF (IS .NE. 0) IM = IMAGE
      SS = RSHIFT(ANDF(BUF(K+1),M1019),10)
      PS = ANDF(BUF(K+1),M0009)
      LL = RSHIFT(ANDF(BUF(K+2),M2029),20)
      CS = RSHIFT(ANDF(BUF(K+2),M1019),10)
      HL = ANDF(BUF(K+2),M0009)
C
C     LOOP THROUGH MDI ENTRY FOR THIS SUBSTRUCTURE DETERMINING THE
C     SIZE OF EACH EXISTING ITEM.
C
      DO 70 J = 1,NITM
      JJ = J + IFRST - 1
      IF (BUF(K+JJ) .EQ. 0) GO TO 60
      INUM = RSHIFT(BUF(K+JJ),IHALF)*BLKSIZ
      INUM = ALOG10(FLOAT(INUM)) + .3
      ITM(J) = NUM(INUM)
      IF (IS.NE.0 .AND. ITEM(4,J).EQ.0) ITM(J) = NUM(10)
      IF (PS.NE.0 .AND. IS.EQ.0 .AND. ITEM(5,J).EQ.0) ITM(J) = NUM(10)
      GO TO 70
   60 ITM(J) = BLANK
   70 CONTINUE
C
      LINE = LINE + 1
      IF (LINE .LE. NLPP) GO TO 90
      CALL PAGE1
      LINE = LINE + 9 - 4
      WRITE  (NOUT,80) HDR
   80 FORMAT (//,26X,90HS U B S T R U C T U R E   O P E R A T I N G   F
     1I L E   T A B L E   O F   C O N T E N T S , //,
     1 1H ,51X,26(A1,2X),A1,/1H ,51X,26(A1,2X),A1,/1H ,51X,26(A1,2X),A1,
     2 /,1H ,4X,12HSUBSTRUCTURE,35X,26(A1,2X),A1, /1H ,4X,3HNO.,3X,4HNAM
     3E,4X,4HTYPE,3X,2HSS,3X,2HPS,3X,2HLL,3X,2HCS,3X,2HHL,4X,80(1H-)/)
C
   90 WRITE  (NOUT,100) I,SSNAME,IM,TYPE(IT),SS,PS,LL,CS,HL,
     1                  (ITM(L),L=1,NITM)
  100 FORMAT (2X,I6,2X,2A4,2X,A1,A2,5(1X,I4),4X,26(A1,2X),A1)
  110 CONTINUE
C
C     PRINT SOF SPACE UTILIZATION MESSAGE
C
      LINE = LINE + 8
      IF (LINE .GT. NLPP) CALL PAGE1
      K    = SOFSIZ(K)
      NBLK = 0
      DO 115 I = 1,NFILES
  115 NBLK = NBLK + FILSIZ(I)
      IPER = (AVBLKS*100)/NBLK
      WRITE  (NOUT,120) K,AVBLKS,IPER,HIBLK
  120 FORMAT (//,51X,80HSIZE OF ITEM IS GIVEN IN POWERS OF TEN   (0 INDI
     1CATES DATA IS STORED IN PRIMARY) ,/,
     2        26H0*** UNUSED SPACE ON SOF = ,I9,7H WORDS.  ,/,
     3        22X,                   4HOR = ,I9,8H BLOCKS. ,/,
     4        22X,                   4HOR = ,I9,9H PERCENT.,/,
     5        26H0*** HIGHEST BLOCK USED  = ,I9)
      LINE = NLPP
      RETURN
      END


text:
      SUBROUTINE TB01UX( COMPZ, N, M, P, A, LDA, B, LDB, C, LDC, Z, LDZ,
     $                   NOBSV, NLBLCK, CTAU, TOL, IWORK, DWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute an orthogonal transformation matrix Z which reduces the
C     N-th order system (A,B,C) to the form
C
C                ( Ano  * )               ( Bno )
C       Z'*A*Z = (        ) ,      Z'*B = (     ) ,
C                ( 0   Ao )               ( Bo  )
C
C          C*Z = ( 0   Co ) ,
C
C     where the NOBSV-th order system (Ao,Bo,Co) is observable.
C     The matrix Ano of order N-NOBSV contains the unobservable
C     eigenvalues of A.
C
C     The pencil ( Ao-lambda*I ) has full column rank NOBSV for all
C                (      Co     )
C     lambda, and is in a staircase form, with
C                     _      _            _      _
C                   ( Ak,k   Ak,k-1   ... Ak,2   Ak,1   )
C                   ( _      _            _      _      )
C       ( Ao ) =    ( Ak-1,k Ak-1,k-1 ... Ak-1,2 Ak-1,1 ) ,          (1)
C       ( Co )      (  :       :      ... _ :    _ :    )
C                   (  0       0      ... A1,2   A1,1   )
C                   (                            _      )
C                   (  0       0      ... 0      A0,1   )
C           _
C     where Ai-1,i is a CTAU(i-1)-by-CTAU(i) full column rank matrix
C     (with CTAU(0) = P).
C
C     The orthogonal transformation Z, performed to reduce the system
C     matrices, can be optionally accumulated.
C
C     The reduced order system (Ao,Bo,Co) has the same transfer-function
C     matrix as the original system (A,B,C).
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     COMPZ   CHARACTER*1
C             = 'N':  do not compute Z;
C             = 'I':  Z is initialized to the unit matrix, and the
C                     orthogonal matrix Z is returned.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The dimension of the system state vector; also the order
C             of the square matrix A, the number of rows of the matrix B
C             and the number of columns of the matrix C.  N >= 0.
C
C     M       (input) INTEGER
C             The dimension of system input vector; also the number of
C             columns of the matrix B.  M >= 0.
C
C     P       (input) INTEGER
C             The dimension of system output vector; also the number of
C             rows of the matrix C.  P >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array must
C             contain the N-by-N state matrix A.
C             On exit, the leading N-by-N part of this array contains
C             the transformed state matrix Z'*A*Z,
C
C                                ( Ano  *  )
C                       Z'*A*Z = (         ) ,
C                                ( 0    Ao )
C
C             where Ao is NOBSV-by-NOBSV and Ano is
C             (N-NOBSV)-by-(N-NOBSV).
C             The matrix ( Ao ) is in the observability staircase
C                        ( Co )
C             form (1).
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension
C             (LDB,MAX(M,P))
C             On entry, the leading N-by-M part of this array must
C             contain the N-by-M input matrix B.
C             On exit, the leading N-by-M part of this array contains
C             the transformed input matrix Z'*B.
C
C     LDB     INTEGER
C             The leading dimension of the array B.
C             LDB >= MAX(1,N) if M > 0 or  P > 0;
C             LDB >= 1        if M = 0 and P = 0.
C
C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
C             On entry, the leading P-by-N part of this array must
C             contain the state/output matrix C.
C             On exit, the leading P-by-N part of this array contains
C             the transformed matrix
C
C                     C*Z = (  0   Co ) ,
C
C             where Co is P-by-NOBSV.
C             The matrix ( Ao ) is in the observability staircase
C                        ( Co )
C             form (1).
C
C     LDC     INTEGER
C             The leading dimension of the array C.
C             LDC >= MAX(1,M,P) if N > 0;
C             LDC >= 1          if N = 0.
C
C     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,*)
C             If COMPZ = 'N': Z is not referenced.
C             If COMPZ = 'I': on entry, Z need not be set;
C                             on exit, the leading N-by-N part of this
C                             array contains the orthogonal matrix Z,
C                             i.e., the product of the transformations
C                             applied to A and C on the right.
C
C     LDZ     INTEGER
C             The leading dimension of the array Z.
C             LDZ >= 1,        if COMPZ = 'N';
C             LDZ >= MAX(1,N), if COMPZ = 'I'.
C
C     NOBSV   (output) INTEGER
C             The order of the reduced matrix Ao, and the number of
C             columns of the reduced matrix Co; also, the order of the
C             observable part of the pair (C, A-lambda*I).
C
C     NLBLCK  (output) INTEGER                         _
C             The number k, of full column rank blocks Ai-1,i in the
C             staircase form of the pencil (Ao-lambda*I) (see (1)).
C                                          (     Co    )
C
C     CTAU    (output) INTEGER array, dimension (N)
C             CTAU(i), for i = 1, ..., NLBLCK, is the column dimension
C                                           _
C             of the full column rank block Ai-1,i in the staircase
C             form (1).
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION
C             The tolerance to be used in rank determinations when
C             transforming the pair (A,C). If the user sets TOL > 0,
C             then the given value of TOL is used as a lower bound for
C             reciprocal condition numbers in rank determinations; a
C             (sub)matrix whose estimated condition number is less than
C             1/TOL is considered to be of full rank.  If the user sets
C             TOL <= 0, then an implicitly computed, default tolerance,
C             defined by  TOLDEF = N*N*EPS,  is used instead, where EPS
C             is the machine precision (see LAPACK Library routine
C             DLAMCH).  TOL < 1.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (P)
C
C     DWORK   DOUBLE PRECISION array, dimension (N+MAX(1, N, 3*P, M))
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The subroutine is based on the dual of the reduction
C     algorithms of [1].
C
C     REFERENCES
C
C     [1] Varga, A.
C         Computation of Irreducible Generalized State-Space
C         Realizations.
C         Kybernetika, vol. 26, pp. 89-106, 1990.
C
C     NUMERICAL ASPECTS
C
C     The algorithm is numerically backward stable and requires
C     0( N**3 )  floating point operations.
C
C     FURTHER COMMENTS
C
C     If the system matrices A and C are badly scaled, it is
C     generally recommendable to scale them with the SLICOT routine
C     TB01ID, before calling TG01UX.
C
C     CONTRIBUTOR
C
C     A. Varga, German Aerospace Center, DLR Oberpfaffenhofen.
C     March 2002.
C
C     REVISIONS
C
C     V. Sima, Dec. 2016.
C
C     KEYWORDS
C
C     Observability, minimal realization, orthogonal canonical form,
C     orthogonal transformation.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0D0 )
C     .. Scalar Arguments ..
      CHARACTER          COMPZ
      INTEGER            INFO, LDA, LDB, LDC, LDZ, M, N, NLBLCK, NOBSV,
     $                   P
      DOUBLE PRECISION   TOL
C     .. Array Arguments ..
      INTEGER            CTAU( * ), IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, *  ),
     $                   DWORK( * ), Z( LDZ, * )
C     .. Local Scalars ..
      LOGICAL            ILZ
      INTEGER            LBA, LDWORK
C     .. Local Arrays ..
      DOUBLE PRECISION   DUM(1)
C     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
C     .. External Subroutines ..
      EXTERNAL           AB07MD, MA02BD, TB01UD, TB01XD, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          MAX
C
C     .. Executable Statements ..
C
C     Decode COMPZ.
C
      ILZ = LSAME( COMPZ, 'I' )
C
C     Test the input scalar parameters.
C
      INFO = 0
      IF( .NOT.ILZ .AND. .NOT.LSAME( COMPZ, 'N' ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( M.LT.0 ) THEN
         INFO = -3
      ELSE IF( P.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDB.LT.1 .OR. ( MAX( M, P ).GT.0 .AND. LDB.LT.N ) ) THEN
         INFO = -8
      ELSE IF( LDC.LT.1 .OR. ( LDC.LT.MAX( M, P ) .AND. N.GT.0 ) ) THEN
         INFO = -10
      ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
         INFO = -12
      ELSE IF( TOL.GE.ONE ) THEN
         INFO = -16
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'TB01UX', -INFO )
         RETURN
      END IF
C
C     Build the dual system.
C
      CALL AB07MD( 'Zero D', N, M, P, A, LDA, B, LDB, C, LDC, DUM, 1,
     $             INFO )
C
      LDWORK = MAX( 1, N, 3*P, M )
      CALL TB01UD( COMPZ, N, P, M, A, LDA, B, LDB, C, LDC, NOBSV,
     $             NLBLCK, CTAU, Z, LDZ, DWORK, TOL, IWORK, DWORK(N+1),
     $             LDWORK, INFO )
      IF( NLBLCK.GT.1 ) THEN
         LBA = CTAU(1) + CTAU(2) - 1
      ELSE IF( NLBLCK.EQ.1 ) THEN
         LBA = CTAU(1) - 1
      ELSE
         LBA = 0
      END IF
C
C     Compute the pertransposed dual system exploiting matrix shapes.
C
      LBA = MAX( LBA, N-NOBSV-1 )
      CALL TB01XD( 'Zero D', N, P, M, LBA, MAX( 0, N-1 ), A, LDA, B,
     $             LDB, C, LDC, DUM, 1, INFO )
      IF( ILZ )
     $   CALL MA02BD( 'Right', N, N, Z, LDZ )
      DWORK(1) = DWORK(N+1)
      RETURN
C *** Last line of TB01UX ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


program fancy
use, intrinsic :: iso_fortran_env, only : stderr=>ERROR_UNIT,stdin=>INPUT_UNIT,stdout=>OUTPUT_UNIT
use M_attr, only : attr, attr_update, attr_mode
implicit none
integer :: i

character(len=256),parameter ::   help_list_nodash(*)=[character(len=256) :: &
!'<clear>', &
'<E><w><bo>F</bo>ortran <bo>P</bo>ackage <bo>M</bo>anager:', &
'<E><w>                                                                        ', &
'<E><w><bo>USAGE: <r><un>fpm</un></r><w> [ <r><un>SUBCOMMAND</un></r> <w>[<m>SUBCOMMAND_OPTIONS</m><w>] ]|&
&<bo>[<g>--list</g><w>|<g>--help</g><w>|<g>--version</g><w>]', &
'<E><w>                                                                        ', &
'<E><w> where <r><un><bo>SUBCOMMAND</bo></un></r> <w>is commonly one of &
&<bo><g>new</g><w>|<g>build</g><w>|<g>run</g><w>|<g>test</g></bo> <w>but may be any of', &
'<E><w>                                                                        ', &
'<E><w>  <bo><g>build</g></bo>     <w>Compile the package placing results in the "build" directory', &
'<E><w>  <bo><g>help</g></bo>      <w>Display help                                                ', &
'<E><w>  <bo><g>list</g></bo>      <w>Display this list of subcommand descriptions                ', &
'<E><w>  <bo><g>new</g></bo>       <w>Create a new Fortran package directory with sample files    ', &
'<E><w>  <bo><g>run</g></bo>       <w>Run the local package application programs                  ', &
'<E><w>  <bo><g>test</g></bo>      <w>Run the test programs                                       ', &
'<E><w>  <bo><g>update</g></bo>    <w>Update and manage project dependencies                      ', &
'<E><w>  <bo><g>install</g></bo>   <w>Install project                                             ', &
'<E><w>                                                                        ', &
'<E><w> <bo>Enter "<r><un>fpm</un></r> <g>--list</g><w></bo>" for a brief list of subcommand options. Enter', &
'<E><w> "<bo><r><un>fpm</un></r> <g>--help</g></bo><w>" or &
&"<bo><r><un>fpm</un> <un>SUBCOMMAND</un></r> <g>--help</g><w></bo>" for detailed descriptions.', &
'<E><w>                                                                        ', &
'<E><w> <bo>Note:</bo> color mode is controlled by the environment variable FPM_COLOR. The', &
'<E><w>       set of allowable values is {<g><bo>always,never,auto</bo></g><w>}. The default is "<g><bo>auto</bo></g><w>". ', &
'<E><w> ']

character(len=256),parameter ::   help_list_dash(*) = [character(len=256) :: &
!'<clear>', &
'<E><w><bo>F</bo>ortran </bo>P</bo>ackage <bo>M</bo>anager:', &
'<E><w>                                                           ', &
'<E><w> <bo><g>build</g><w> [<g>--compiler</g> <m>COMPILER_NAME</m><w>] &
&[<g>--profile</g> <m>PROF</m><w>] [<g>--flag</g> <m>FFLAGS</m><w>] [<g>--list</g><w>]', &
'<E><w>                                                           ', &
'<E><w> <bo><g>help</g><w> [<r><un>NAME(s)</un></r><w>]                       ', &
'<E><w>                                                           ', &
'<E><w> <bo><g>new</g><w> <r><un>NAME</un></r><w> [[<g>--lib</g><w>|<g>--src</g><w>] [<g>--app</g><w>] &
&[<g>--test</g><w>] [<g>--example</g><w>]]|', &
'<E><w> <bo>         [<g>--full</g><w>|<g>--bare</g><w>][<g>--backfill</g><w>]', &
'<E><w>                                                           ', &
'<E><w> <bo><g>update</g><w> <w>[<r><un>NAME(s)</un></r><w>] [<g>--fetch-only</g><w>] [<g>--clean</g><w>] [<g>--verbose</g><w>]', &
'<E><w>                                                           ', &
'<E><w> <bo><g>list</g><w> [<g>--list</g><w>]', &
'<E><w>                                                           ', &
'<E><w> <bo><g>run</g><w>  [[<g>--target</g><w>] <r><un>NAME(s)</un></r> <w>[<g>--example</g><w>]&
&[<g>--profile</g><w> <m>PROF</m><w>] [<g>--flag</g><w> <m>FFLAGS</m><w>] [<g>--all</g><w>]', &
'<E><w> <bo>     [<g>--runner</g><w> <m>"CMD"</m><w>] [<g>--compiler</g><w> <m>COMPILER_NAME</m><w>] &
&[<g>--list</g><w>] [-- <m>ARGS</m><w>]', &
'<E><w>                                                           ', &
'<E><w> <bo><g>test</g><w> [[<g>--target</g><w>] <r><un>NAME(s)</un></r><w>] [<g>--profile</g><w> <m>PROF</m><w>] &
&[<g>--flag</g><w> <m>FFLAGS</m><w>] [<g>--runner</g><w> <m>"CMD"</m><w>]', &
'<E><w> <bo>     [<g>--list</g><w>] [<g>--compiler</g><w> <m>COMPILER_NAME</m><w>] [-- <m>ARGS</m><w>]', &
'<E><w>                                                           ', &
'<E><w> <bo><g>install</g><w> [<g>--profile</g><w> <m>PROF</m>] [<g>--flag</g><w> <m>FFLAGS</m>] &
&[<g>--no-rebuild</g><w>] [<g>--prefix</g><w> <m>PATH</m><w>] ', &
'<E><w> <bo>        <w>[<m>OPTIONS</m><w>]', &
'<E><w> ']

! shows a limitation if periods are replaced with spaces the CHAR= option does not
! count the colored blanks so if no text after the blanks too much padding is appended to the line
! need to count last character with a background attribute or something different
!nvfortran bug!character(len=*),parameter :: text(*)= [character(len=132) :: &
character(len=*),parameter :: text(22)= [character(len=132) :: &
'<E><e>                                                     ',&
'<E><w><bo>       <C>                                   </C><E></bo>',&
'<E><w><bo>       <C>      LLL   </C><bo><w> F</bo>ortran   <C>            </bo><W> <e></C><E>',&
'<E><w><bo>       <C>     LL LL  </C><bo><w>  P</bo>ackage  <C>            </bo><W> <e></C><E>',&
'<E><w><bo>       <C>     LL     </C><bo><w>   M</bo>anager <C>            </bo><W> <e></C><E>',&
'<E><w><bo>       <C>   LLLLLL                          <W> <E></bo>',&
'<E><w><bo>       <C>     LL                            <W> <E></bo>',&
'<E><w><bo>       <C>     LL               LLLLL LLL    <W> <E></bo>',&
'<E><w><bo>       <C>     LL     LLLLLL     LL  L  LL   <W> <E></bo>',&
'<E><w><bo>       <C>     LL    LL     L    LL  L  LL   <W> <E></bo>',&
'<E><w><bo>       <C>     LL    LL     L    LL  L  LL   <W> <E></bo>',&
'<E><w><bo>       <C>     LL    LLLLLLL     LL  L  LL   <W> <E></bo>',&
'<E><w><bo>       <C>           LL                      <W> <E></bo>',&
'<E><w><bo>       <C>           LL                      <W> <E></bo>',&
'<E><w><bo>       <C>           LL                      <W> <E></bo>',&
'<E><w><bo>       <C>                                   <W> <E></bo>',&
'<E><w><bo>        <W>                                   <E></bo><e>   ',&
'<E><bo>Program:<w>     fpm(1)                                     ',&
'<E><bo>Description:<w> package manager and build system for Fortran',&
'<E><bo>Version:<w>     0.3.0, alpha                               ',&
'<E><bo>License:<w>     MIT                                        ',&
'<E><bo>Home Page:<w>   https://github.com/fortran-lang/fpm        ']
   call attr_mode(manner='color')

   write(*,'(a)')(attr(trim(help_list_nodash(i)),chars=80),i=1,size(help_list_nodash))
   call paws()

   write(*,'(a)')(attr(trim(help_list_dash(i)),chars=80),i=1,size(help_list_dash))
   call paws()

   ! add custom keywords
   call attr_mode(manner='plain')
   write(*,'(a)')(attr(trim(text(i)),chars=80),i=1,size(text))
   call paws()

   call attr_mode(manner='color')
   write(*,'(a)')(attr(trim(text(i)),chars=80),i=1,size(text))
   call paws()
contains

subroutine paws()
character(len=1) :: letter
integer :: ios

   write(stdout,'(a)',advance='no',iostat=ios)'Enter [RETURN] to continue ...'
   read(*,'(a)',iostat=ios)letter
   write(stdout,'(a)')attr('<clear>')
   flush(unit=stdout,iostat=ios)

end subroutine paws

end program fancy