text:
C
C     file genbun.f
C
      SUBROUTINE GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,IERROR,W)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 1999 by UCAR                   *
C     *                                                               *
C     *       UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                      FISHPACK version 4.1                     *
C     *                                                               *
C     *     A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF      *
C     *                                                               *
C     *      SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS        *
C     *                                                               *
C     *                             BY                                *
C     *                                                               *
C     *        JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET         *
C     *                                                               *
C     *                             OF                                *
C     *                                                               *
C     *         THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH          *
C     *                                                               *
C     *                BOULDER, COLORADO  (80307)  U.S.A.             *
C     *                                                               *
C     *                   WHICH IS SPONSORED BY                       *
C     *                                                               *
C     *              THE NATIONAL SCIENCE FOUNDATION                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C
C DIMENSION OF           A(M),B(M),C(M),Y(IDIMY,N),
C                        W(SEE PARAMETER LIST)
C ARGUMENTS
C
C LATEST REVISION        NOVEMBER 1988
C
C PURPOSE                THE NAME OF THIS PACKAGE IS A MNEMONIC FOR THE
C                        GENERALIZED BUNEMAN ALGORITHM.
C
C                        IT SOLVES THE REAL LINEAR SYSTEM OF EQUATIONS
C
C                        A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J)
C                        + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J)
C
C                        FOR I = 1,2,...,M  AND  J = 1,2,...,N.
C
C                        INDICES I+1 AND I-1 ARE EVALUATED MODULO M,
C                        I.E., X(0,J) = X(M,J) AND X(M+1,J) = X(1,J),
C                        AND X(I,0) MAY EQUAL 0, X(I,2), OR X(I,N),
C                        AND X(I,N+1) MAY EQUAL 0, X(I,N-1), OR X(I,1)
C                        DEPENDING ON AN INPUT PARAMETER.
C
C USAGE                  CALL GENBUN (NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,
C                                     IERROR,W)
C
C ARGUMENTS
C
C ON INPUT               NPEROD
C
C                          INDICATES THE VALUES THAT X(I,0) AND
C                          X(I,N+1) ARE ASSUMED TO HAVE.
C
C                          = 0  IF X(I,0) = X(I,N) AND X(I,N+1) =
C                               X(I,1).
C                          = 1  IF X(I,0) = X(I,N+1) = 0  .
C                          = 2  IF X(I,0) = 0 AND X(I,N+1) = X(I,N-1).
C                          = 3  IF X(I,0) = X(I,2) AND X(I,N+1) =
C                               X(I,N-1).
C                          = 4  IF X(I,0) = X(I,2) AND X(I,N+1) = 0.
C
C                        N
C                          THE NUMBER OF UNKNOWNS IN THE J-DIRECTION.
C                          N MUST BE GREATER THAN 2.
C
C                        MPEROD
C                          = 0 IF A(1) AND C(M) ARE NOT ZERO
C                          = 1 IF A(1) = C(M) = 0
C
C                        M
C                          THE NUMBER OF UNKNOWNS IN THE I-DIRECTION.
C                          N MUST BE GREATER THAN 2.
C
C                        A,B,C
C                          ONE-DIMENSIONAL ARRAYS OF LENGTH M THAT
C                          SPECIFY THE COEFFICIENTS IN THE LINEAR
C                          EQUATIONS GIVEN ABOVE.  IF MPEROD = 0
C                          THE ARRAY ELEMENTS MUST NOT DEPEND UPON
C                          THE INDEX I, BUT MUST BE CONSTANT.
C                          SPECIFICALLY, THE SUBROUTINE CHECKS THE
C                          FOLLOWING CONDITION .
C
C                            A(I) = C(1)
C                            C(I) = C(1)
C                            B(I) = B(1)
C
C                          FOR I=1,2,...,M.
C
C                        IDIMY
C                          THE ROW (OR FIRST) DIMENSION OF THE
C                          TWO-DIMENSIONAL ARRAY Y AS IT APPEARS
C                          IN THE PROGRAM CALLING GENBUN.
C                          THIS PARAMETER IS USED TO SPECIFY THE
C                          VARIABLE DIMENSION OF Y.
C                          IDIMY MUST BE AT LEAST M.
C
C                        Y
C                          A TWO-DIMENSIONAL COMPLEX ARRAY THAT
C                          SPECIFIES THE VALUES OF THE RIGHT SIDE
C                          OF THE LINEAR SYSTEM OF EQUATIONS GIVEN
C                          ABOVE.
C                          Y MUST BE DIMENSIONED AT LEAST M*N.
C
C                        W
C                          A ONE-DIMENSIONAL ARRAY THAT MUST
C                          BE PROVIDED BY THE USER FOR WORK
C                          SPACE.  W MAY REQUIRE UP TO 4*N +
C                          (10 + INT(LOG2(N)))*M LOCATIONS.
C                          THE ACTUAL NUMBER OF LOCATIONS USED IS
C                          COMPUTED BY GENBUN AND IS RETURNED IN
C                          LOCATION W(1).
C
C
C  ON OUTPUT             Y
C
C                          CONTAINS THE SOLUTION X.
C
C                        IERROR
C                          AN ERROR FLAG WHICH INDICATES INVALID
C                          INPUT PARAMETERS  EXCEPT FOR NUMBER
C                          ZERO, A SOLUTION IS NOT ATTEMPTED.
C
C                          = 0  NO ERROR.
C                          = 1  M .LE. 2  .
C                          = 2  N .LE. 2
C                          = 3  IDIMY .LT. M
C                          = 4  NPEROD .LT. 0 OR NPEROD .GT. 4
C                          = 5  MPEROD .LT. 0 OR MPEROD .GT. 1
C                          = 6  A(I) .NE. C(1) OR C(I) .NE. C(1) OR
C                               B(I) .NE. B(1) FOR
C                               SOME I=1,2,...,M.
C                          = 7  A(1) .NE. 0 OR C(M) .NE. 0 AND
C                                 MPEROD = 1
C
C                        W
C                          W(1) CONTAINS THE REQUIRED LENGTH OF W.
C
C SPECIAL CONDITONS      NONE
C
C I/O                    NONE
C
C PRECISION              SINGLE
C
C REQUIRED LIBRARY       COMF AND GNBNAUX FROM FISHPACK
C FILES
C
C LANGUAGE               FORTRAN
C
C HISTORY                WRITTEN IN 1979 BY ROLAND SWEET OF NCAR'S
C                        SCIENTIFIC COMPUTING DIVISION.  MADE AVAILABLE
C                        ON NCAR'S PUBLIC LIBRARIES IN JANUARY, 1980.
C
C ALGORITHM              THE LINEAR SYSTEM IS SOLVED BY A CYCLIC
C                        REDUCTION ALGORITHM DESCRIBED IN THE
C                        REFERENCE.
C
C PORTABILITY            FORTRAN 77 --
C                        THE MACHINE DEPENDENT CONSTANT PI IS
C                        DEFINED IN FUNCTION PIMACH.
C
C REFERENCES             SWEET, R., "A CYCLIC REDUCTION ALGORITHM FOR
C                        SOLVING BLOCK TRIDIAGONAL SYSTEMS OF ARBITRARY
C                        DIMENSIONS," SIAM J. ON NUMER. ANAL., 14 (1977)
C                        PP. 706-720.
C
C ACCURACY               THIS TEST WAS PERFORMED ON A CDC 7600:
C
C                        A UNIFORM RANDOM NUMBER GENERATOR WAS USED
C                        TO CREATE A SOLUTION ARRAY X FOR THE SYSTEM
C                        GIVEN IN THE 'PURPOSE' DESCRIPTION ABOVE
C                        WITH
C                          A(I) = C(I) = -0.5*B(I) = 1, I=1,2,...,M
C
C                        AND, WHEN MPEROD = 1
C
C                          A(1) = C(M) = 0
C                          A(M) = C(1) = 2.
C
C                        THE SOLUTION X WAS SUBSTITUTED INTO THE
C                        GIVEN SYSTEM  AND, USING DOUBLE PRECISION
C                        A RIGHT SIDE Y WAS COMPUTED.
C                        USING THIS ARRAY Y, SUBROUTINE GENBUN
C                        WAS CALLED TO PRODUCE APPROXIMATE
C                        SOLUTION Z.  THEN RELATIVE ERROR
C                          E = MAX(ABS(Z(I,J)-X(I,J)))/
C                              MAX(ABS(X(I,J)))
C                        WAS COMPUTED, WHERE THE TWO MAXIMA ARE TAKEN
C                        OVER I=1,2,...,M AND J=1,...,N.
C
C                        THE VALUE OF E IS GIVEN IN THE TABLE
C                        BELOW FOR SOME TYPICAL VALUES OF M AND N.
C
C                   M (=N)    MPEROD    NPEROD    T(MSECS)    E
C                   ------    ------    ------    --------  ------
C
C                     31        0         0          36     6.E-14
C                     31        1         1          21     4.E-13
C                     31        1         3          41     3.E-13
C                     32        0         0          29     9.E-14
C                     32        1         1          32     3.E-13
C                     32        1         3          48     1.E-13
C                     33        0         0          36     9.E-14
C                     33        1         1          30     4.E-13
C                     33        1         3          34     1.E-13
C                     63        0         0         150     1.E-13
C                     63        1         1          91     1.E-12
C                     63        1         3         173     2.E-13
C                     64        0         0         122     1.E-13
C                     64        1         1         128     1.E-12
C                     64        1         3         199     6.E-13
C                     65        0         0         143     2.E-13
C                     65        1         1         120     1.E-12
C                     65        1         3         138     4.E-13
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
      DIMENSION       Y(IDIMY,1)
      DIMENSION       W(*)       ,B(*)       ,A(*)       ,C(*)
C
      IERROR = 0
      IF (M .LE. 2) IERROR = 1
      IF (N .LE. 2) IERROR = 2
      IF (IDIMY .LT. M) IERROR = 3
      IF (NPEROD.LT.0 .OR. NPEROD.GT.4) IERROR = 4
      IF (MPEROD.LT.0 .OR. MPEROD.GT.1) IERROR = 5
      IF (MPEROD .EQ. 1) GO TO 102
      DO 101 I=2,M
         IF (A(I) .NE. C(1)) GO TO 103
         IF (C(I) .NE. C(1)) GO TO 103
         IF (B(I) .NE. B(1)) GO TO 103
  101 CONTINUE
      GO TO 104
  102 IF (A(1).NE.0. .OR. C(M).NE.0.) IERROR = 7
      GO TO 104
  103 IERROR = 6
  104 IF (IERROR .NE. 0) RETURN
      MP1 = M+1
      IWBA = MP1
      IWBB = IWBA+M
      IWBC = IWBB+M
      IWB2 = IWBC+M
      IWB3 = IWB2+M
      IWW1 = IWB3+M
      IWW2 = IWW1+M
      IWW3 = IWW2+M
      IWD = IWW3+M
      IWTCOS = IWD+M
      IWP = IWTCOS+4*N
      DO 106 I=1,M
         K = IWBA+I-1
         W(K) = -A(I)
         K = IWBC+I-1
         W(K) = -C(I)
         K = IWBB+I-1
         W(K) = 2.-B(I)
         DO 105 J=1,N
            Y(I,J) = -Y(I,J)
  105    CONTINUE
  106 CONTINUE
      MP = MPEROD+1
      NP = NPEROD+1
      GO TO (114,107),MP
  107 GO TO (108,109,110,111,123),NP
  108 CALL POISP2 (M,N,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
     1             W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
     2             W(IWP))
      GO TO 112
  109 CALL POISD2 (M,N,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWW1),
     1             W(IWD),W(IWTCOS),W(IWP))
      GO TO 112
  110 CALL POISN2 (M,N,1,2,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
     1             W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
     2             W(IWP))
      GO TO 112
  111 CALL POISN2 (M,N,1,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
     1             W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
     2             W(IWP))
  112 IPSTOR = W(IWW1)
      IREV = 2
      IF (NPEROD .EQ. 4) GO TO 124
  113 GO TO (127,133),MP
  114 CONTINUE
C
C     REORDER UNKNOWNS WHEN MP =0
C
      MH = (M+1)/2
      MHM1 = MH-1
      MODD = 1
      IF (MH*2 .EQ. M) MODD = 2
      DO 119 J=1,N
         DO 115 I=1,MHM1
            MHPI = MH+I
            MHMI = MH-I
            W(I) = Y(MHMI,J)-Y(MHPI,J)
            W(MHPI) = Y(MHMI,J)+Y(MHPI,J)
  115    CONTINUE
         W(MH) = 2.*Y(MH,J)
         GO TO (117,116),MODD
  116    W(M) = 2.*Y(M,J)
  117    CONTINUE
         DO 118 I=1,M
            Y(I,J) = W(I)
  118    CONTINUE
  119 CONTINUE
      K = IWBC+MHM1-1
      I = IWBA+MHM1
      W(K) = 0.
      W(I) = 0.
      W(K+1) = 2.*W(K+1)
      GO TO (120,121),MODD
  120 CONTINUE
      K = IWBB+MHM1-1
      W(K) = W(K)-W(I-1)
      W(IWBC-1) = W(IWBC-1)+W(IWBB-1)
      GO TO 122
  121 W(IWBB-1) = W(K+1)
  122 CONTINUE
      GO TO 107
C
C     REVERSE COLUMNS WHEN NPEROD = 4.
C
  123 IREV = 1
      NBY2 = N/2
  124 DO 126 J=1,NBY2
         MSKIP = N+1-J
         DO 125 I=1,M
            A1 = Y(I,J)
            Y(I,J) = Y(I,MSKIP)
            Y(I,MSKIP) = A1
  125    CONTINUE
  126 CONTINUE
      GO TO (110,113),IREV
  127 CONTINUE
      DO 132 J=1,N
         DO 128 I=1,MHM1
            MHMI = MH-I
            MHPI = MH+I
            W(MHMI) = .5*(Y(MHPI,J)+Y(I,J))
            W(MHPI) = .5*(Y(MHPI,J)-Y(I,J))
  128    CONTINUE
         W(MH) = .5*Y(MH,J)
         GO TO (130,129),MODD
  129    W(M) = .5*Y(M,J)
  130    CONTINUE
         DO 131 I=1,M
            Y(I,J) = W(I)
  131    CONTINUE
  132 CONTINUE
  133 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR W ARRAY.
C
      W(1) = IPSTOR+IWP-1
      RETURN
      END
      SUBROUTINE POISD2 (MR,NR,ISTAG,BA,BB,BC,Q,IDIMQ,B,W,D,TCOS,P)
C
C     SUBROUTINE TO SOLVE POISSON'S EQUATION FOR DIRICHLET BOUNDARY
C     CONDITIONS.
C
C     ISTAG = 1 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A.
C     ISTAG = 2 IF THE LAST DIAGONAL BLOCK IS THE MATRIX A+I.
C
      DIMENSION       Q(IDIMQ,1) ,BA(*)      ,BB(*)      ,BC(*)      ,
     1                TCOS(*)    ,B(*)       ,D(*)       ,W(*)       ,
     2                P(*)
      M = MR
      N = NR
      JSH = 0
      FI = 1./FLOAT(ISTAG)
      IP = -M
      IPSTOR = 0
      GO TO (101,102),ISTAG
  101 KR = 0
      IRREG = 1
      IF (N .GT. 1) GO TO 106
      TCOS(1) = 0.
      GO TO 103
  102 KR = 1
      JSTSAV = 1
      IRREG = 2
      IF (N .GT. 1) GO TO 106
      TCOS(1) = -1.
  103 DO 104 I=1,M
         B(I) = Q(I,1)
  104 CONTINUE
      CALL TRIX (1,0,M,BA,BB,BC,B,TCOS,D,W)
      DO 105 I=1,M
         Q(I,1) = B(I)
  105 CONTINUE
      GO TO 183
  106 LR = 0
      DO 107 I=1,M
         P(I) = 0.
  107 CONTINUE
      NUN = N
      JST = 1
      JSP = N
C
C     IRREG = 1 WHEN NO IRREGULARITIES HAVE OCCURRED, OTHERWISE IT IS 2.
C
  108 L = 2*JST
      NODD = 2-2*((NUN+1)/2)+NUN
C
C     NODD = 1 WHEN NUN IS ODD, OTHERWISE IT IS 2.
C
      GO TO (110,109),NODD
  109 JSP = JSP-L
      GO TO 111
  110 JSP = JSP-JST
      IF (IRREG .NE. 1) JSP = JSP-L
  111 CONTINUE
C
C     REGULAR REDUCTION
C
      CALL COSGEN (JST,1,0.5,0.0,TCOS)
      IF (L .GT. JSP) GO TO 118
      DO 117 J=L,JSP,L
         JM1 = J-JSH
         JP1 = J+JSH
         JM2 = J-JST
         JP2 = J+JST
         JM3 = JM2-JSH
         JP3 = JP2+JSH
         IF (JST .NE. 1) GO TO 113
         DO 112 I=1,M
            B(I) = 2.*Q(I,J)
            Q(I,J) = Q(I,JM2)+Q(I,JP2)
  112    CONTINUE
         GO TO 115
  113    DO 114 I=1,M
            T = Q(I,J)-Q(I,JM1)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2)
            B(I) = T+Q(I,J)-Q(I,JM3)-Q(I,JP3)
            Q(I,J) = T
  114    CONTINUE
  115    CONTINUE
         CALL TRIX (JST,0,M,BA,BB,BC,B,TCOS,D,W)
         DO 116 I=1,M
            Q(I,J) = Q(I,J)+B(I)
  116    CONTINUE
  117 CONTINUE
C
C     REDUCTION FOR LAST UNKNOWN
C
  118 GO TO (119,136),NODD
  119 GO TO (152,120),IRREG
C
C     ODD NUMBER OF UNKNOWNS
C
  120 JSP = JSP+L
      J = JSP
      JM1 = J-JSH
      JP1 = J+JSH
      JM2 = J-JST
      JP2 = J+JST
      JM3 = JM2-JSH
      GO TO (123,121),ISTAG
  121 CONTINUE
      IF (JST .NE. 1) GO TO 123
      DO 122 I=1,M
         B(I) = Q(I,J)
         Q(I,J) = 0.
  122 CONTINUE
      GO TO 130
  123 GO TO (124,126),NODDPR
  124 DO 125 I=1,M
         IP1 = IP+I
         B(I) = .5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))+P(IP1)+Q(I,J)
  125 CONTINUE
      GO TO 128
  126 DO 127 I=1,M
         B(I) = .5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))+Q(I,JP2)-Q(I,JP1)+Q(I,J)
  127 CONTINUE
  128 DO 129 I=1,M
         Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
  129 CONTINUE
  130 CALL TRIX (JST,0,M,BA,BB,BC,B,TCOS,D,W)
      IP = IP+M
      IPSTOR = MAX0(IPSTOR,IP+M)
      DO 131 I=1,M
         IP1 = IP+I
         P(IP1) = Q(I,J)+B(I)
         B(I) = Q(I,JP2)+P(IP1)
  131 CONTINUE
      IF (LR .NE. 0) GO TO 133
      DO 132 I=1,JST
         KRPI = KR+I
         TCOS(KRPI) = TCOS(I)
  132 CONTINUE
      GO TO 134
  133 CONTINUE
      CALL COSGEN (LR,JSTSAV,0.,FI,TCOS(JST+1))
      CALL MERGE (TCOS,0,JST,JST,LR,KR)
  134 CONTINUE
      CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS)
      CALL TRIX (KR,KR,M,BA,BB,BC,B,TCOS,D,W)
      DO 135 I=1,M
         IP1 = IP+I
         Q(I,J) = Q(I,JM2)+B(I)+P(IP1)
  135 CONTINUE
      LR = KR
      KR = KR+L
      GO TO 152
C
C     EVEN NUMBER OF UNKNOWNS
C
  136 JSP = JSP+L
      J = JSP
      JM1 = J-JSH
      JP1 = J+JSH
      JM2 = J-JST
      JP2 = J+JST
      JM3 = JM2-JSH
      GO TO (137,138),IRREG
  137 CONTINUE
      JSTSAV = JST
      IDEG = JST
      KR = L
      GO TO 139
  138 CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS)
      CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1))
      IDEG = KR
      KR = KR+JST
  139 IF (JST .NE. 1) GO TO 141
      IRREG = 2
      DO 140 I=1,M
         B(I) = Q(I,J)
         Q(I,J) = Q(I,JM2)
  140 CONTINUE
      GO TO 150
  141 DO 142 I=1,M
         B(I) = Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))
  142 CONTINUE
      GO TO (143,145),IRREG
  143 DO 144 I=1,M
         Q(I,J) = Q(I,JM2)+.5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
  144 CONTINUE
      IRREG = 2
      GO TO 150
  145 CONTINUE
      GO TO (146,148),NODDPR
  146 DO 147 I=1,M
         IP1 = IP+I
         Q(I,J) = Q(I,JM2)+P(IP1)
  147 CONTINUE
      IP = IP-M
      GO TO 150
  148 DO 149 I=1,M
         Q(I,J) = Q(I,JM2)+Q(I,J)-Q(I,JM1)
  149 CONTINUE
  150 CALL TRIX (IDEG,LR,M,BA,BB,BC,B,TCOS,D,W)
      DO 151 I=1,M
         Q(I,J) = Q(I,J)+B(I)
  151 CONTINUE
  152 NUN = NUN/2
      NODDPR = NODD
      JSH = JST
      JST = 2*JST
      IF (NUN .GE. 2) GO TO 108
C
C     START SOLUTION.
C
      J = JSP
      DO 153 I=1,M
         B(I) = Q(I,J)
  153 CONTINUE
      GO TO (154,155),IRREG
  154 CONTINUE
      CALL COSGEN (JST,1,0.5,0.0,TCOS)
      IDEG = JST
      GO TO 156
  155 KR = LR+JST
      CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS)
      CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1))
      IDEG = KR
  156 CONTINUE
      CALL TRIX (IDEG,LR,M,BA,BB,BC,B,TCOS,D,W)
      JM1 = J-JSH
      JP1 = J+JSH
      GO TO (157,159),IRREG
  157 DO 158 I=1,M
         Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I)
  158 CONTINUE
      GO TO 164
  159 GO TO (160,162),NODDPR
  160 DO 161 I=1,M
         IP1 = IP+I
         Q(I,J) = P(IP1)+B(I)
  161 CONTINUE
      IP = IP-M
      GO TO 164
  162 DO 163 I=1,M
         Q(I,J) = Q(I,J)-Q(I,JM1)+B(I)
  163 CONTINUE
  164 CONTINUE
C
C     START BACK SUBSTITUTION.
C
      JST = JST/2
      JSH = JST/2
      NUN = 2*NUN
      IF (NUN .GT. N) GO TO 183
      DO 182 J=JST,N,L
         JM1 = J-JSH
         JP1 = J+JSH
         JM2 = J-JST
         JP2 = J+JST
         IF (J .GT. JST) GO TO 166
         DO 165 I=1,M
            B(I) = Q(I,J)+Q(I,JP2)
  165    CONTINUE
         GO TO 170
  166    IF (JP2 .LE. N) GO TO 168
         DO 167 I=1,M
            B(I) = Q(I,J)+Q(I,JM2)
  167    CONTINUE
         IF (JST .LT. JSTSAV) IRREG = 1
         GO TO (170,171),IRREG
  168    DO 169 I=1,M
            B(I) = Q(I,J)+Q(I,JM2)+Q(I,JP2)
  169    CONTINUE
  170    CONTINUE
         CALL COSGEN (JST,1,0.5,0.0,TCOS)
         IDEG = JST
         JDEG = 0
         GO TO 172
  171    IF (J+L .GT. N) LR = LR-JST
         KR = JST+LR
         CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS)
         CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1))
         IDEG = KR
         JDEG = LR
  172    CONTINUE
         CALL TRIX (IDEG,JDEG,M,BA,BB,BC,B,TCOS,D,W)
         IF (JST .GT. 1) GO TO 174
         DO 173 I=1,M
            Q(I,J) = B(I)
  173    CONTINUE
         GO TO 182
  174    IF (JP2 .GT. N) GO TO 177
  175    DO 176 I=1,M
            Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I)
  176    CONTINUE
         GO TO 182
  177    GO TO (175,178),IRREG
  178    IF (J+JSH .GT. N) GO TO 180
         DO 179 I=1,M
            IP1 = IP+I
            Q(I,J) = B(I)+P(IP1)
  179    CONTINUE
         IP = IP-M
         GO TO 182
  180    DO 181 I=1,M
            Q(I,J) = B(I)+Q(I,J)-Q(I,JM1)
  181    CONTINUE
  182 CONTINUE
      L = L/2
      GO TO 164
  183 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR P VECTORS.
C
      W(1) = IPSTOR
      RETURN
      END
      SUBROUTINE POISN2 (M,N,ISTAG,MIXBND,A,BB,C,Q,IDIMQ,B,B2,B3,W,W2,
     1                   W3,D,TCOS,P)
C
C     SUBROUTINE TO SOLVE POISSON'S EQUATION WITH NEUMANN BOUNDARY
C     CONDITIONS.
C
C     ISTAG = 1 IF THE LAST DIAGONAL BLOCK IS A.
C     ISTAG = 2 IF THE LAST DIAGONAL BLOCK IS A-I.
C     MIXBND = 1 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTH BOUNDARIES.
C     MIXBND = 2 IF HAVE NEUMANN BOUNDARY CONDITIONS AT BOTTOM AND
C     DIRICHLET CONDITION AT TOP.  (FOR THIS CASE, MUST HAVE ISTAG = 1.)
C
      DIMENSION       A(*)       ,BB(*)      ,C(*)       ,Q(IDIMQ,*) ,
     1                B(*)       ,B2(*)      ,B3(*)      ,W(*)       ,
     2                W2(*)      ,W3(*)      ,D(*)       ,TCOS(*)    ,
     3                K(4)       ,P(*)
      EQUIVALENCE     (K(1),K1)  ,(K(2),K2)  ,(K(3),K3)  ,(K(4),K4)
      FISTAG = 3-ISTAG
      FNUM = 1./FLOAT(ISTAG)
      FDEN = 0.5*FLOAT(ISTAG-1)
      MR = M
      IP = -MR
      IPSTOR = 0
      I2R = 1
      JR = 2
      NR = N
      NLAST = N
      KR = 1
      LR = 0
      GO TO (101,103),ISTAG
  101 CONTINUE
      DO 102 I=1,MR
         Q(I,N) = .5*Q(I,N)
  102 CONTINUE
      GO TO (103,104),MIXBND
  103 IF (N .LE. 3) GO TO 155
  104 CONTINUE
      JR = 2*I2R
      NROD = 1
      IF ((NR/2)*2 .EQ. NR) NROD = 0
      GO TO (105,106),MIXBND
  105 JSTART = 1
      GO TO 107
  106 JSTART = JR
      NROD = 1-NROD
  107 CONTINUE
      JSTOP = NLAST-JR
      IF (NROD .EQ. 0) JSTOP = JSTOP-I2R
      CALL COSGEN (I2R,1,0.5,0.0,TCOS)
      I2RBY2 = I2R/2
      IF (JSTOP .GE. JSTART) GO TO 108
      J = JR
      GO TO 116
  108 CONTINUE
C
C     REGULAR REDUCTION.
C
      DO 115 J=JSTART,JSTOP,JR
         JP1 = J+I2RBY2
         JP2 = J+I2R
         JP3 = JP2+I2RBY2
         JM1 = J-I2RBY2
         JM2 = J-I2R
         JM3 = JM2-I2RBY2
         IF (J .NE. 1) GO TO 109
         JM1 = JP1
         JM2 = JP2
         JM3 = JP3
  109    CONTINUE
         IF (I2R .NE. 1) GO TO 111
         IF (J .EQ. 1) JM2 = JP2
         DO 110 I=1,MR
            B(I) = 2.*Q(I,J)
            Q(I,J) = Q(I,JM2)+Q(I,JP2)
  110    CONTINUE
         GO TO 113
  111    CONTINUE
         DO 112 I=1,MR
            FI = Q(I,J)
            Q(I,J) = Q(I,J)-Q(I,JM1)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2)
            B(I) = FI+Q(I,J)-Q(I,JM3)-Q(I,JP3)
  112    CONTINUE
  113    CONTINUE
         CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W)
         DO 114 I=1,MR
            Q(I,J) = Q(I,J)+B(I)
  114    CONTINUE
C
C     END OF REDUCTION FOR REGULAR UNKNOWNS.
C
  115 CONTINUE
C
C     BEGIN SPECIAL REDUCTION FOR LAST UNKNOWN.
C
      J = JSTOP+JR
  116 NLAST = J
      JM1 = J-I2RBY2
      JM2 = J-I2R
      JM3 = JM2-I2RBY2
      IF (NROD .EQ. 0) GO TO 128
C
C     ODD NUMBER OF UNKNOWNS
C
      IF (I2R .NE. 1) GO TO 118
      DO 117 I=1,MR
         B(I) = FISTAG*Q(I,J)
         Q(I,J) = Q(I,JM2)
  117 CONTINUE
      GO TO 126
  118 DO 119 I=1,MR
         B(I) = Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))
  119 CONTINUE
      IF (NRODPR .NE. 0) GO TO 121
      DO 120 I=1,MR
         II = IP+I
         Q(I,J) = Q(I,JM2)+P(II)
  120 CONTINUE
      IP = IP-MR
      GO TO 123
  121 CONTINUE
      DO 122 I=1,MR
         Q(I,J) = Q(I,J)-Q(I,JM1)+Q(I,JM2)
  122 CONTINUE
  123 IF (LR .EQ. 0) GO TO 124
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(KR+1))
      GO TO 126
  124 CONTINUE
      DO 125 I=1,MR
         B(I) = FISTAG*B(I)
  125 CONTINUE
  126 CONTINUE
      CALL COSGEN (KR,1,0.5,FDEN,TCOS)
      CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W)
      DO 127 I=1,MR
         Q(I,J) = Q(I,J)+B(I)
  127 CONTINUE
      KR = KR+I2R
      GO TO 151
  128 CONTINUE
C
C     EVEN NUMBER OF UNKNOWNS
C
      JP1 = J+I2RBY2
      JP2 = J+I2R
      IF (I2R .NE. 1) GO TO 135
      DO 129 I=1,MR
         B(I) = Q(I,J)
  129 CONTINUE
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      IP = 0
      IPSTOR = MR
      GO TO (133,130),ISTAG
  130 DO 131 I=1,MR
         P(I) = B(I)
         B(I) = B(I)+Q(I,N)
  131 CONTINUE
      TCOS(1) = 1.
      TCOS(2) = 0.
      CALL TRIX (1,1,MR,A,BB,C,B,TCOS,D,W)
      DO 132 I=1,MR
         Q(I,J) = Q(I,JM2)+P(I)+B(I)
  132 CONTINUE
      GO TO 150
  133 CONTINUE
      DO 134 I=1,MR
         P(I) = B(I)
         Q(I,J) = Q(I,JM2)+2.*Q(I,JP2)+3.*B(I)
  134 CONTINUE
      GO TO 150
  135 CONTINUE
      DO 136 I=1,MR
         B(I) = Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))
  136 CONTINUE
      IF (NRODPR .NE. 0) GO TO 138
      DO 137 I=1,MR
         II = IP+I
         B(I) = B(I)+P(II)
  137 CONTINUE
      GO TO 140
  138 CONTINUE
      DO 139 I=1,MR
         B(I) = B(I)+Q(I,JP2)-Q(I,JP1)
  139 CONTINUE
  140 CONTINUE
      CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W)
      IP = IP+MR
      IPSTOR = MAX0(IPSTOR,IP+MR)
      DO 141 I=1,MR
         II = IP+I
         P(II) = B(I)+.5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
         B(I) = P(II)+Q(I,JP2)
  141 CONTINUE
      IF (LR .EQ. 0) GO TO 142
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(I2R+1))
      CALL MERGE (TCOS,0,I2R,I2R,LR,KR)
      GO TO 144
  142 DO 143 I=1,I2R
         II = KR+I
         TCOS(II) = TCOS(I)
  143 CONTINUE
  144 CALL COSGEN (KR,1,0.5,FDEN,TCOS)
      IF (LR .NE. 0) GO TO 145
      GO TO (146,145),ISTAG
  145 CONTINUE
      CALL TRIX (KR,KR,MR,A,BB,C,B,TCOS,D,W)
      GO TO 148
  146 CONTINUE
      DO 147 I=1,MR
         B(I) = FISTAG*B(I)
  147 CONTINUE
  148 CONTINUE
      DO 149 I=1,MR
         II = IP+I
         Q(I,J) = Q(I,JM2)+P(II)+B(I)
  149 CONTINUE
  150 CONTINUE
      LR = KR
      KR = KR+JR
  151 CONTINUE
      GO TO (152,153),MIXBND
  152 NR = (NLAST-1)/JR+1
      IF (NR .LE. 3) GO TO 155
      GO TO 154
  153 NR = NLAST/JR
      IF (NR .LE. 1) GO TO 192
  154 I2R = JR
      NRODPR = NROD
      GO TO 104
  155 CONTINUE
C
C      BEGIN SOLUTION
C
      J = 1+JR
      JM1 = J-I2R
      JP1 = J+I2R
      JM2 = NLAST-I2R
      IF (NR .EQ. 2) GO TO 184
      IF (LR .NE. 0) GO TO 170
      IF (N .NE. 3) GO TO 161
C
C     CASE N = 3.
C
      GO TO (156,168),ISTAG
  156 CONTINUE
      DO 157 I=1,MR
         B(I) = Q(I,2)
  157 CONTINUE
      TCOS(1) = 0.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 158 I=1,MR
         Q(I,2) = B(I)
         B(I) = 4.*B(I)+Q(I,1)+2.*Q(I,3)
  158 CONTINUE
      TCOS(1) = -2.
      TCOS(2) = 2.
      I1 = 2
      I2 = 0
      CALL TRIX (I1,I2,MR,A,BB,C,B,TCOS,D,W)
      DO 159 I=1,MR
         Q(I,2) = Q(I,2)+B(I)
         B(I) = Q(I,1)+2.*Q(I,2)
  159 CONTINUE
      TCOS(1) = 0.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 160 I=1,MR
         Q(I,1) = B(I)
  160 CONTINUE
      JR = 1
      I2R = 0
      GO TO 194
C
C     CASE N = 2**P+1
C
  161 CONTINUE
      GO TO (162,170),ISTAG
  162 CONTINUE
      DO 163 I=1,MR
         B(I) = Q(I,J)+.5*Q(I,1)-Q(I,JM1)+Q(I,NLAST)-Q(I,JM2)
  163 CONTINUE
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
      DO 164 I=1,MR
         Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I)
         B(I) = Q(I,1)+2.*Q(I,NLAST)+4.*Q(I,J)
  164 CONTINUE
      JR2 = 2*JR
      CALL COSGEN (JR,1,0.0,0.0,TCOS)
      DO 165 I=1,JR
         I1 = JR+I
         I2 = JR+1-I
         TCOS(I1) = -TCOS(I2)
  165 CONTINUE
      CALL TRIX (JR2,0,MR,A,BB,C,B,TCOS,D,W)
      DO 166 I=1,MR
         Q(I,J) = Q(I,J)+B(I)
         B(I) = Q(I,1)+2.*Q(I,J)
  166 CONTINUE
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
      DO 167 I=1,MR
         Q(I,1) = .5*Q(I,1)-Q(I,JM1)+B(I)
  167 CONTINUE
      GO TO 194
C
C     CASE OF GENERAL N WITH NR = 3 .
C
  168 DO 169 I=1,MR
         B(I) = Q(I,2)
         Q(I,2) = 0.
         B2(I) = Q(I,3)
         B3(I) = Q(I,1)
  169 CONTINUE
      JR = 1
      I2R = 0
      J = 2
      GO TO 177
  170 CONTINUE
      DO 171 I=1,MR
         B(I) = .5*Q(I,1)-Q(I,JM1)+Q(I,J)
  171 CONTINUE
      IF (NROD .NE. 0) GO TO 173
      DO 172 I=1,MR
         II = IP+I
         B(I) = B(I)+P(II)
  172 CONTINUE
      GO TO 175
  173 DO 174 I=1,MR
         B(I) = B(I)+Q(I,NLAST)-Q(I,JM2)
  174 CONTINUE
  175 CONTINUE
      DO 176 I=1,MR
         T = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
         Q(I,J) = T
         B2(I) = Q(I,NLAST)+T
         B3(I) = Q(I,1)+2.*T
  176 CONTINUE
  177 CONTINUE
      K1 = KR+2*JR-1
      K2 = KR+JR
      TCOS(K1+1) = -2.
      K4 = K1+3-ISTAG
      CALL COSGEN (K2+ISTAG-2,1,0.0,FNUM,TCOS(K4))
      K4 = K1+K2+1
      CALL COSGEN (JR-1,1,0.0,1.0,TCOS(K4))
      CALL MERGE (TCOS,K1,K2,K1+K2,JR-1,0)
      K3 = K1+K2+LR
      CALL COSGEN (JR,1,0.5,0.0,TCOS(K3+1))
      K4 = K3+JR+1
      CALL COSGEN (KR,1,0.5,FDEN,TCOS(K4))
      CALL MERGE (TCOS,K3,JR,K3+JR,KR,K1)
      IF (LR .EQ. 0) GO TO 178
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(K4))
      CALL MERGE (TCOS,K3,JR,K3+JR,LR,K3-LR)
      CALL COSGEN (KR,1,0.5,FDEN,TCOS(K4))
  178 K3 = KR
      K4 = KR
      CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3)
      DO 179 I=1,MR
         B(I) = B(I)+B2(I)+B3(I)
  179 CONTINUE
      TCOS(1) = 2.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 180 I=1,MR
         Q(I,J) = Q(I,J)+B(I)
         B(I) = Q(I,1)+2.*Q(I,J)
  180 CONTINUE
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
      IF (JR .NE. 1) GO TO 182
      DO 181 I=1,MR
         Q(I,1) = B(I)
  181 CONTINUE
      GO TO 194
  182 CONTINUE
      DO 183 I=1,MR
         Q(I,1) = .5*Q(I,1)-Q(I,JM1)+B(I)
  183 CONTINUE
      GO TO 194
  184 CONTINUE
      IF (N .NE. 2) GO TO 188
C
C     CASE  N = 2
C
      DO 185 I=1,MR
         B(I) = Q(I,1)
  185 CONTINUE
      TCOS(1) = 0.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 186 I=1,MR
         Q(I,1) = B(I)
         B(I) = 2.*(Q(I,2)+B(I))*FISTAG
  186 CONTINUE
      TCOS(1) = -FISTAG
      TCOS(2) = 2.
      CALL TRIX (2,0,MR,A,BB,C,B,TCOS,D,W)
      DO 187 I=1,MR
         Q(I,1) = Q(I,1)+B(I)
  187 CONTINUE
      JR = 1
      I2R = 0
      GO TO 194
  188 CONTINUE
C
C     CASE OF GENERAL N AND NR = 2 .
C
      DO 189 I=1,MR
         II = IP+I
         B3(I) = 0.
         B(I) = Q(I,1)+2.*P(II)
         Q(I,1) = .5*Q(I,1)-Q(I,JM1)
         B2(I) = 2.*(Q(I,1)+Q(I,NLAST))
  189 CONTINUE
      K1 = KR+JR-1
      TCOS(K1+1) = -2.
      K4 = K1+3-ISTAG
      CALL COSGEN (KR+ISTAG-2,1,0.0,FNUM,TCOS(K4))
      K4 = K1+KR+1
      CALL COSGEN (JR-1,1,0.0,1.0,TCOS(K4))
      CALL MERGE (TCOS,K1,KR,K1+KR,JR-1,0)
      CALL COSGEN (KR,1,0.5,FDEN,TCOS(K1+1))
      K2 = KR
      K4 = K1+K2+1
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(K4))
      K3 = LR
      K4 = 0
      CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3)
      DO 190 I=1,MR
         B(I) = B(I)+B2(I)
  190 CONTINUE
      TCOS(1) = 2.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 191 I=1,MR
         Q(I,1) = Q(I,1)+B(I)
  191 CONTINUE
      GO TO 194
  192 DO 193 I=1,MR
         B(I) = Q(I,NLAST)
  193 CONTINUE
      GO TO 196
  194 CONTINUE
C
C     START BACK SUBSTITUTION.
C
      J = NLAST-JR
      DO 195 I=1,MR
         B(I) = Q(I,NLAST)+Q(I,J)
  195 CONTINUE
  196 JM2 = NLAST-I2R
      IF (JR .NE. 1) GO TO 198
      DO 197 I=1,MR
         Q(I,NLAST) = 0.
  197 CONTINUE
      GO TO 202
  198 CONTINUE
      IF (NROD .NE. 0) GO TO 200
      DO 199 I=1,MR
         II = IP+I
         Q(I,NLAST) = P(II)
  199 CONTINUE
      IP = IP-MR
      GO TO 202
  200 DO 201 I=1,MR
         Q(I,NLAST) = Q(I,NLAST)-Q(I,JM2)
  201 CONTINUE
  202 CONTINUE
      CALL COSGEN (KR,1,0.5,FDEN,TCOS)
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(KR+1))
      IF (LR .NE. 0) GO TO 204
      DO 203 I=1,MR
         B(I) = FISTAG*B(I)
  203 CONTINUE
  204 CONTINUE
      CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W)
      DO 205 I=1,MR
         Q(I,NLAST) = Q(I,NLAST)+B(I)
  205 CONTINUE
      NLASTP = NLAST
  206 CONTINUE
      JSTEP = JR
      JR = I2R
      I2R = I2R/2
      IF (JR .EQ. 0) GO TO 222
      GO TO (207,208),MIXBND
  207 JSTART = 1+JR
      GO TO 209
  208 JSTART = JR
  209 CONTINUE
      KR = KR-JR
      IF (NLAST+JR .GT. N) GO TO 210
      KR = KR-JR
      NLAST = NLAST+JR
      JSTOP = NLAST-JSTEP
      GO TO 211
  210 CONTINUE
      JSTOP = NLAST-JR
  211 CONTINUE
      LR = KR-JR
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      DO 221 J=JSTART,JSTOP,JSTEP
         JM2 = J-JR
         JP2 = J+JR
         IF (J .NE. JR) GO TO 213
         DO 212 I=1,MR
            B(I) = Q(I,J)+Q(I,JP2)
  212    CONTINUE
         GO TO 215
  213    CONTINUE
         DO 214 I=1,MR
            B(I) = Q(I,J)+Q(I,JM2)+Q(I,JP2)
  214    CONTINUE
  215    CONTINUE
         IF (JR .NE. 1) GO TO 217
         DO 216 I=1,MR
            Q(I,J) = 0.
  216    CONTINUE
         GO TO 219
  217    CONTINUE
         JM1 = J-I2R
         JP1 = J+I2R
         DO 218 I=1,MR
            Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
  218    CONTINUE
  219    CONTINUE
         CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
         DO 220 I=1,MR
            Q(I,J) = Q(I,J)+B(I)
  220    CONTINUE
  221 CONTINUE
      NROD = 1
      IF (NLAST+I2R .LE. N) NROD = 0
      IF (NLASTP .NE. NLAST) GO TO 194
      GO TO 206
  222 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR P VECTORS.
C
      W(1) = IPSTOR
      RETURN
      END
      SUBROUTINE POISP2 (M,N,A,BB,C,Q,IDIMQ,B,B2,B3,W,W2,W3,D,TCOS,P)
C
C     SUBROUTINE TO SOLVE POISSON EQUATION WITH PERIODIC BOUNDARY
C     CONDITIONS.
C
      DIMENSION       A(*)       ,BB(*)      ,C(*)       ,Q(IDIMQ,1) ,
     1                B(*)       ,B2(*)      ,B3(*)      ,W(*)       ,
     2                W2(*)      ,W3(*)      ,D(*)       ,TCOS(*)    ,
     3                P(*)
      MR = M
      NR = (N+1)/2
      NRM1 = NR-1
      IF (2*NR .NE. N) GO TO 107
C
C     EVEN NUMBER OF UNKNOWNS
C
      DO 102 J=1,NRM1
         NRMJ = NR-J
         NRPJ = NR+J
         DO 101 I=1,MR
            S = Q(I,NRMJ)-Q(I,NRPJ)
            T = Q(I,NRMJ)+Q(I,NRPJ)
            Q(I,NRMJ) = S
            Q(I,NRPJ) = T
  101    CONTINUE
  102 CONTINUE
      DO 103 I=1,MR
         Q(I,NR) = 2.*Q(I,NR)
         Q(I,N) = 2.*Q(I,N)
  103 CONTINUE
      CALL POISD2 (MR,NRM1,1,A,BB,C,Q,IDIMQ,B,W,D,TCOS,P)
      IPSTOR = W(1)
      CALL POISN2 (MR,NR+1,1,1,A,BB,C,Q(1,NR),IDIMQ,B,B2,B3,W,W2,W3,D,
     1             TCOS,P)
      IPSTOR = MAX0(IPSTOR,INT(W(1)))
      DO 105 J=1,NRM1
         NRMJ = NR-J
         NRPJ = NR+J
         DO 104 I=1,MR
            S = .5*(Q(I,NRPJ)+Q(I,NRMJ))
            T = .5*(Q(I,NRPJ)-Q(I,NRMJ))
            Q(I,NRMJ) = S
            Q(I,NRPJ) = T
  104    CONTINUE
  105 CONTINUE
      DO 106 I=1,MR
         Q(I,NR) = .5*Q(I,NR)
         Q(I,N) = .5*Q(I,N)
  106 CONTINUE
      GO TO 118
  107 CONTINUE
C
C     ODD  NUMBER OF UNKNOWNS
C
      DO 109 J=1,NRM1
         NRPJ = N+1-J
         DO 108 I=1,MR
            S = Q(I,J)-Q(I,NRPJ)
            T = Q(I,J)+Q(I,NRPJ)
            Q(I,J) = S
            Q(I,NRPJ) = T
  108    CONTINUE
  109 CONTINUE
      DO 110 I=1,MR
         Q(I,NR) = 2.*Q(I,NR)
  110 CONTINUE
      LH = NRM1/2
      DO 112 J=1,LH
         NRMJ = NR-J
         DO 111 I=1,MR
            S = Q(I,J)
            Q(I,J) = Q(I,NRMJ)
            Q(I,NRMJ) = S
  111    CONTINUE
  112 CONTINUE
      CALL POISD2 (MR,NRM1,2,A,BB,C,Q,IDIMQ,B,W,D,TCOS,P)
      IPSTOR = W(1)
      CALL POISN2 (MR,NR,2,1,A,BB,C,Q(1,NR),IDIMQ,B,B2,B3,W,W2,W3,D,
     1             TCOS,P)
      IPSTOR = MAX0(IPSTOR,INT(W(1)))
      DO 114 J=1,NRM1
         NRPJ = NR+J
         DO 113 I=1,MR
            S = .5*(Q(I,NRPJ)+Q(I,J))
            T = .5*(Q(I,NRPJ)-Q(I,J))
            Q(I,NRPJ) = T
            Q(I,J) = S
  113    CONTINUE
  114 CONTINUE
      DO 115 I=1,MR
         Q(I,NR) = .5*Q(I,NR)
  115 CONTINUE
      DO 117 J=1,LH
         NRMJ = NR-J
         DO 116 I=1,MR
            S = Q(I,J)
            Q(I,J) = Q(I,NRMJ)
            Q(I,NRMJ) = S
  116    CONTINUE
  117 CONTINUE
  118 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR P VECTORS.
C
      W(1) = IPSTOR
      RETURN
C
C REVISION HISTORY---
C
C SEPTEMBER 1973    VERSION 1
C APRIL     1976    VERSION 2
C JANUARY   1978    VERSION 3
C DECEMBER  1979    VERSION 3.1
C FEBRUARY  1985    DOCUMENTATION UPGRADE
C NOVEMBER  1988    VERSION 3.2, FORTRAN 77 CHANGES
C-----------------------------------------------------------------------
      END

max_stars_repo_name: ipelupessy/omuse

max_issues_repo_name: ipelupessy/omuse

max_forks_repo_name: ipelupessy/omuse


text:
C
C       This routine gives an example of how to use the IBIS-2
C       subroutine package using the FORTRAN bridges.
C
C	Uses: Ported VICAR runtime library,
C       IBIS-2 subroutine library.
C

	include 'VICMAIN_FOR'

	subroutine main44
	! Call the C routine from here so that
	! The file-I/O will work ok
	call main44_c
	return
	end

	subroutine fortran_test
	logical xvptst
	if (xvptst('A4TEST')) then
		call a4test
	else
		call copytest
	endif
	return
	end

c
c  a4test tests the fortran a4 conversion
c   for old IBIS files.
c

	subroutine a4test

	integer i,unit,ibis,status
	character*4  buf(10)

	call xvunit(unit,'inp',1,status,0)
	if (status.ne.1) call xvsignal(unit,status,1)

	call  ibis_file_open(unit,ibis,'read',0,0,
     +                                ' ',' ',status)
	if (status .ne. 1) call ibis_signal_u(unit,status,1)

	call  ibis_column_set(ibis,'FORMAT','A4',1,status)
	if (status .ne. 1) call ibis_signal(ibis,status,1)
	call  ibis_column_set(ibis,'U_FORMAT','A4',1,status)
	if (status .ne. 1) call ibis_signal(ibis,status,1)

	call  ibis_column_read(ibis,buf,1,1,3,status)
	if (status .ne. 1) call ibis_signal(ibis,status,1)

	do i=1,3
		call xvmessage(buf(i),' ')
	enddo

	call ibis_file_close(ibis,' ',status)
	return
	end

C	copytest: Copies one general new IBIS file to another,
C	possibly with a change in organization (row/column).
C
	subroutine copytest
	integer inunit,outunit,status
	integer count
	integer  ibis_out,ibis_in,ibufsize
	parameter (ibufsize=250)
	integer nc, nr, col,colsize
	integer row, row_inc, nrow
	character*16  org
	character*6 format(1000)
	character*6 cfmt
	logical ascii
	character*64 cbuffer(ibufsize*8)
	real*8 buffer(ibufsize*64)
	! The equivalence is just to conserve memory:
	equivalence (buffer,cbuffer)

C   open the input IBIS file, abort on error. The four '0' arguments
C   in the ibis_file_open are only used for creating output files

	call xvunit( inunit, 'inp', 1, status, ' ')
	call ibis_file_open( inunit, ibis_in, 'read', 0,0,' ',' ', status)
	if (status.ne.1) call ibis_signal_u( inunit, status, 1)

C  get the # rows & columns and the column-format
C  of the input IBIS file:

	count = ibis_file_get( ibis_in,'nc', nc, 1,1 )
	if (count.lt.0) call ibis_signal( ibis_in, count, 1)
	count = ibis_file_get( ibis_in,'nr', nr, 1,1 )
	if (count.lt.0) call ibis_signal( ibis_in, count, 1)
	count = ibis_file_get( ibis_in,'formats', format, 1, nc )
	if (count.lt.0) call ibis_signal( ibis_in, count, 1)

C open the output file with identical structure. The routine returns
C in variable 'ibis_out' the handle to the IBIS file, which should be used
C in all IBIS file-io operations

	call xvpone( 'org', org, 1, 0 )
	call xvunit( outunit, 'out', 1,status, ' ')
	call ibis_file_open( outunit, ibis_out, 'write', nc,
     +  	             nr, format, org, status )
	if (status.ne.1) call ibis_signal_u( outunit, status, 1)

C   copy the column data from input to output

	do col=1,nc

		call ibis_column_get( ibis_in, 'format', cfmt, col, status )
		if (status .ne.1) call ibis_signal( ibis_in, status, 1)
		ascii = (cfmt(1:1).eq.'a'.or.cfmt(1:1).eq.'A')
		if (ascii) then
			call ibis_column_set( ibis_in, 'u_format','a64',col,status)
			if (status .ne.1) call ibis_signal( ibis_in, status, 1)
			call ibis_column_set( ibis_out, 'u_format','a64',col,status)
			if (status .ne.1) call ibis_signal( ibis_out, status, 1)
		endif

		! Get the size of a column's element, in bytes
		call ibis_column_get( ibis_in, 'size', colsize, col, status )
		if (status .ne.1) call ibis_signal( ibis_in, status, 1)

		! Compute #elements that will fit into buffer
		row_inc = (ibufsize + colsize  - 1)/colsize

		! Actually do some writing. The ibis_column_read/write calls
		! read/write <nrow> rows of column #<col> from buffer,
		! starting at row #<row>. The data is automatically converted
		! into the native host format of the given column.

		do row = 1,nr,row_inc
			nrow = MIN( nr+1-row , row_inc)

			! If this is an ASCII column, need to use the A*64 buffer:
			! Note that we could have used the same buffer if we
			! just turned the U_FORMAT to NONE.

			if (ascii) then
			   call ibis_column_read( ibis_in, cbuffer,
     +				col, row, nrow,status )
			else
			   call ibis_column_read( ibis_in, buffer,
     +				col, row, nrow,status )
			endif

     			! We know that ibis_in is a valid IBIS handle here
     			! so in case of error, we can call ibis_signal

			if (status.ne.1) call ibis_signal( ibis_in, status, 1)

			! If this is an ASCII column, need to use CHAR:
			if (ascii) then
			   call ibis_column_write( ibis_out, cbuffer,
     +				col, row, nrow,status )
			else
			   call ibis_column_write( ibis_out, buffer,
     +				col, row, nrow,status )
			endif
			if (status.ne.1) call ibis_signal( ibis_out, status, 1)
		enddo
	enddo

C close up shop

	call ibis_file_close( ibis_in, ' ', status )
	if (status.ne.1) call ibis_signal( ibis_in, status, 1)
	call ibis_file_close( ibis_out,' ', status )
	if (status.ne.1) call ibis_signal( ibis_out, status, 1)

	return
	end

ext: f

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
! SPDX-Identifer: MIT
module test_string_operator
    use testdrive, only : new_unittest, unittest_type, error_type, check
    use stdlib_string_type, only : string_type, assignment(=), len, &
        operator(>), operator(<), operator(>=), operator(<=), &
        operator(/=), operator(==), operator(//)
    implicit none

contains

    !> Collect all exported unit tests
    subroutine collect_string_operator(testsuite)
        !> Collection of tests
        type(unittest_type), allocatable, intent(out) :: testsuite(:)

        testsuite = [ &
            new_unittest("gt", test_gt), &
            new_unittest("lt", test_lt), &
            new_unittest("ge", test_ge), &
            new_unittest("le", test_le), &
            new_unittest("eq", test_eq), &
            new_unittest("ne", test_ne), &
            new_unittest("concat", test_concat) &
            ]
    end subroutine collect_string_operator

    subroutine test_gt(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(string_type) :: string
        logical :: res

        string = "bcd"
        res = string > "abc"
        call check(error, res .eqv. .true.)
        if (allocated(error)) return

        res = string > "bcd"
        call check(error, res .eqv. .false.)
        if (allocated(error)) return

        res = string > "cde"
        call check(error, res .eqv. .false.)
    end subroutine test_gt

    subroutine test_lt(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(string_type) :: string
        logical :: res

        string = "bcd"
        res = string < "abc"
        call check(error, res .eqv. .false.)
        if (allocated(error)) return

        res = string < "bcd"
        call check(error, res .eqv. .false.)
        if (allocated(error)) return

        res = string < "cde"
        call check(error, res .eqv. .true.)
    end subroutine test_lt

    subroutine test_ge(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(string_type) :: string
        logical :: res

        string = "bcd"
        res = string >= "abc"
        call check(error, res .eqv. .true.)
        if (allocated(error)) return

        res = string >= "bcd"
        call check(error, res .eqv. .true.)
        if (allocated(error)) return

        res = string >= "cde"
        call check(error, res .eqv. .false.)
    end subroutine test_ge

    subroutine test_le(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(string_type) :: string
        logical :: res

        string = "bcd"
        res = string <= "abc"
        call check(error, res .eqv. .false.)
        if (allocated(error)) return

        res = string <= "bcd"
        call check(error, res .eqv. .true.)
        if (allocated(error)) return

        res = string <= "cde"
        call check(error, res .eqv. .true.)
    end subroutine test_le

    subroutine test_eq(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(string_type) :: string
        logical :: res

        string = "bcd"
        res = string == "abc"
        call check(error, res .eqv. .false.)
        if (allocated(error)) return

        res = string == "bcd"
        call check(error, res .eqv. .true.)
        if (allocated(error)) return

        res = string == "cde"
        call check(error, res .eqv. .false.)
    end subroutine test_eq

    subroutine test_ne(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(string_type) :: string
        logical :: res

        string = "bcd"
        res = string /= "abc"
        call check(error, res .eqv. .true.)
        if (allocated(error)) return

        res = string /= "bcd"
        call check(error, res .eqv. .false.)
        if (allocated(error)) return

        res = string /= "cde"
        call check(error, res .eqv. .true.)
    end subroutine test_ne

    subroutine test_concat(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        type(string_type) :: string

        string = "Hello, "
        string = string // "World!"
        call check(error, len(string) == 13)
    end subroutine test_concat

end module test_string_operator

program tester
    use, intrinsic :: iso_fortran_env, only : error_unit
    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
    use test_string_operator, only : collect_string_operator
    implicit none
    integer :: stat, is
    type(testsuite_type), allocatable :: testsuites(:)
    character(len=*), parameter :: fmt = '("#", *(1x, a))'

    stat = 0

    testsuites = [ &
        new_testsuite("string-operator", collect_string_operator) &
        ]

    do is = 1, size(testsuites)
        write(error_unit, fmt) "Testing:", testsuites(is)%name
        call run_testsuite(testsuites(is)%collect, error_unit, stat)
    end do

    if (stat > 0) then
        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
        error stop
    end if
end program

max_stars_repo_name: freevryheid/stdlib

max_issues_repo_name: seanwallawalla-forks/stdlib

max_forks_repo_name: seanwallawalla-forks/stdlib


text:
! { dg-do "run" }
! Transformational intrinsic MATMUL as initialization expression.

  REAL, PARAMETER :: PI = 3.141592654, theta = PI/6.0

  REAL, PARAMETER :: unity(2,2) = RESHAPE([1.0, 0.0, 0.0, 1.0], [2, 2])
  REAL, PARAMETER :: m1(2,2)    = RESHAPE([COS(theta), SIN(theta), -SIN(theta), COS(theta)], [2, 2])
  REAL, PARAMETER :: m2(2,2)    = RESHAPE([COS(theta), -SIN(theta), SIN(theta), COS(theta)], [2, 2])
  REAL, PARAMETER :: m(2,2)     = MATMUL(m1, m2)

  IF (ANY(ABS(m - unity) > EPSILON(0.0))) CALL abort()
END

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE LMG1ALG(ISUM,ISUMI,LCA,LCIA,LCJA,LCU1,LCFRHS,LCIG,
     1                   ISIZ1,ISIZ2,ISIZ3,ISIZ4,ICG,NCOL,NROW,NLAY,
     2                   IN,IOUT,IFREFM)
C
C     ******************************************************************
C     LMG / AMG can no longer be distributed.  Write message to this
C     effect and stop execution.
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
C
C1------PRINT A MESSAGE IDENTIFYING AMG PACKAGE
      WRITE(IOUT,500)
  500 FORMAT(1X,/,1X,'LMG -- ALGEBRAIC MULTI-GRID SOLUTION PACKAGE:',/,
     &    1X,'BECAUSE OF LICENSING RESTRICTIONS, THE ALGEBRAIC',
     &    ' MULTI-GRID SOLVER CAN NO',/,
     &    1X,'LONGER BE DISTRIBUTED BY THE USGS.  PLEASE CHOOSE',
     &    ' ANOTHER SOLVER PACKAGE.')
      CALL USTOP(' ')
C
      RETURN
      END
C
C***********************************************************************
      SUBROUTINE LMG1RPG(IN,MXITER,MXCYC,BCLOSE,DAMP,IOUTAMG,IOUT,
     &                  IFREFM,ICG,IADAMP,DUP,DLOW,HCLOSE)
C
      END
C
C***********************************************************************
      SUBROUTINE LMG1AP(HNEW,IBOUND,CR,CC,CV,HCOF,RHS,A,IA,JA,U,FRHS,IG,
     &                  ISIZ1,ISIZ2,ISIZ3,ISIZ4,KITER,BCLOSE,DAMP,ICNVG,
     &                  KSTP,KPER,MXITER,MXCYC,NCOL,NROW,NLAY,NODES,
     &                  HNOFLO,IOUT,IOUTAMG,ICG,IADAMP,DUP,DLOW)
C
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      DOUBLEPRECISION HNEW
      DOUBLEPRECISION A(ISIZ1),U(ISIZ4),FRHS(ISIZ4)
      DIMENSION  IBOUND(NODES), HNEW(NODES),CR(NODES), CC(NODES),
     1   CV(NODES), HCOF(NODES), RHS(NODES), IA(ISIZ2),
     2   JA(ISIZ1),IG(ISIZ3)
C
      END

max_stars_repo_name: quanpands/wflow

max_issues_repo_name: quanpands/wflow

max_forks_repo_name: quanpands/wflow


text:
      SUBROUTINE READ_LWF
     &          (lu_lwf,print_lwf,
     &           bearng,rhomx,rlat,rlon,rrho,
     &           mxps,nrps,mxsgmnt,nrsgmnt,sgmnt,
     &           mxprm,nrprm,param,nrcmp,nrlwf,
     &           mxpts,nrpts,xy,amp_lwf,phs_lwf,
     &           begin_file,end_file)

c***********************************************************************
c                         subroutine read_lwf
c***********************************************************************

c  Program Source:  Naval Ocean Systems Center - Code 542

c  Date:
c     10 Apr 1990

c  Function:
c     Reads mode sum data.

c  Parameters passed:
c     lu_lwf           [i] number of the logical unit; must already be
c                          openned by the calling routine
c     print_lwf        [i] >0 prints parameters from passed data

c     mxps             [i] dimension of array SGMNT in calling routine
c     mxsgmnt          [i] dimension of array SGMNT in calling routine
c     mxprm            [i] dimension of array PARAM in calling routine

c     mxpts            [i] dimension of arrays XY, AMP_LWF and PHS_LWF

c     begin_file       [l] indicates beginning of file; must be set to
c                          .TRUE. at  by the calling routine

c  Parameters returned:
c     bearng           [r] bearing of path; degrees E of N
c     rhomx            [r] maximum range of path; km
c     rlat             [r] receiver  latitude; degrees N
c     rlon             [r] receiver longitude; degrees W
c     rrho             [r] range to receiver; km

c     nrps             [i] number of parameters passed in SGMNT
c     nrsgmnt          [i] number of segments passed in SGMNT
c     sgmnt            [r] path segmentation parameters (see below)

c     nrprm            [i] number of elements in the array PARAM
c     param            [r] parameters (see below)
c     nrcmp            [i] number of components to be read
c     nrlwf            [i] number of parametric cases to be read

c     nrpts            [i] number of points in XY, AMP_LWF and PHS_LWF
c     xy               [r] XY values for AMP_LWF and PHS_LWF
c     amp_lwf          [r] array of amplitude of signal; dB
c     phs_lwf          [r] array of relative phase; degrees

c     end_file         [l] indicates end of the input file; set to
c                          .TRUE. by this routine as required

c  Common blocks referenced:

c  Functions and subroutines referenced:
c     min0
c     print
c     read

c  References:

c  Change History:
c     14 Mar 91     The parameter list may contain the date and time if
c                   LWPM was run for a specific date and time.

c     04 Dec 91     Added additional levels of print.

c     20 Feb 92     Added profile index to PARAM.

c     21 Oct 95     Changed to get the LOG unit from LWPC_LUN.CMN.

c  Notes:
c     Wherever indicated, ranges are distances from the transmitter.

c     The file begins with a header record which identifies the program
c     which generated the data, the date the program was run,
c     transmitter parameters and the parameters of the paths. This
c     record is processed by the routine READ_HDR.

c     Subsequent records contain data for the individual mode sums.
c     These mode sums are generated parametrically. To allow for future
c     expansion the parametric parameters are passed in a arrays along
c     with the number of elements in the array.

c     The first of these records contains information on the specific
c     path and a summary of the parameters of the segments used to
c     compute the mode sum. To allow for future expansion, the path
c     parameters are passed in an array along with the number of
c     elements in the array. This program reads these parameters into
c     a two dimensional array, SGMNT, which is assumed to contain the
c     following:

c      1)rho        range to the segment; km
c      2)lat        latitude;  degrees North
c      3)lon        longitude; degrees West
c      4)azim       geomagnetic azimuth; degrees E of path direction
c      5)dip        geomagnetic dip; degrees below horizontal
c      6)bfield     geomagnetic intensity; Gauss
c      7)sigma      ground conductivity; Seimens
c      8)epsr       ground dielectric constant
c      9)beta       exponential ionospheric conductivity slope; km-1
c     10)hprime     exponential ionospheric conductivity height; km

Change: 20 Feb 92   This parameter added
c     11)pindex     profile index for non-exponential profiles

Change: 06 Mar 92   This parameter modified by LW_FULL_MC
c     11)pindex     profile index+65536*(mc step index)

c     This record also contains the values of the parameters for which
c     the mode sums were generated. See LWFLDS for more details
c     regarding the following:

c     nrcmp         number of components
c     nrpts         number of points in primary output
c     nrlwf         number of parameteric variations

c     The next record contains the parameters of the mode sum and the
c     mode sum data. The parameters are read into the array, PARAM,
c     which is assumed to contain the following:

c      1)power      power; Kw
c      2)dist       distance; km
c      3)incl       inclination of the antenna; degrees from vertical
c      4)headng     heading     of the antenna; degrees E of N
c      5)talt       altitude of the transmitter; km
c      6)ralt       altitude of the    reciever; km

Change: 14 Mar 91   These parameters may be added by LW_VS_D
c      7)month      month number
c      8)day        day of the month
c      9)year       year
c     10)UT         Universal time; hours and minutes

c     The mode sum data are:

c     xy            value of the parameter which varies
c     amp           amplitude     ; dB above 1uv/m
c     phs           relative phase; degrees

c*******************!***************************************************

c     LWPC parameters
      include      'lwpc_lun.cmn'

      character, parameter :: label='zyx'
      logical       begin_file,end_file
      integer       lu_lwf,print_lwf,
     &              mxprm,nrprm,nrcmp,nrpts,nrlwf,
     &              mxps,nrps,mxsgmnt,nrsgmnt,mxpts
      real     *  4 bearng,rhomx,rlat,rlon,rrho,
     &              sgmnt(mxps,mxsgmnt),param(mxprm),
     &              xy(mxpts),amp_lwf(mxpts),phs_lwf(mxpts)

      if (begin_file) then

c        At the start of the file
         begin_file=.false.
         end_file=.false.
         ncmp=0
         if (print_lwf .eq. 1) then

c           Heading for summary output
            write(lwpcLOG_lun,
     &          '(''nc nrpt bearng  rhomx  rlat   rlon   rrho    '',
     &            ''pwr    dist incl headng talt ralt'')')
         end if
      end if

      if (ncmp .eq. 0) then

c        Read path parameters
         read (lu_lwf,end=99) bearng,rhomx,rlat,rlon,rrho,
     &                        nrps,nrsgmnt,nrprm,nrpts,nrcmp,nrlwf,
     &                      ((sgmnt(i,j),i=1,nrps),j=1,nrsgmnt)
         nlwf=0
      end if

      do i=1,mxprm
         param(i)=0.
      end do
      read (lu_lwf) (param(i),i=1,nrprm),
     &              (xy(i),amp_lwf(i),phs_lwf(i),i=1,nrpts)
      ncmp=ncmp+1

      if (print_lwf .eq. 1) then

c        Summary output
         if (ncmp .eq. 1) then
            write(lwpcLOG_lun,
     &          '(i2,i5,f7.1,f7.0,f6.1,f7.1,f7.0,
     &            f7.0,f8.0,f5.0,f7.1,2f5.1)')
     &            nrcmp,nrpts,bearng,rhomx,rlat,rlon,rrho,
     &           (param(i),i=1,MIN0(nrprm,6))
         end if
      else
     &if (print_lwf .gt. 1) then

c        Heading for detailed outputs
         write(lwpcLOG_lun,
     &       '(''nc nrpt bearng  rhomx  rlat   rlon   rrho    '',
     &         ''pwr    dist incl headng talt ralt'')')
         write(lwpcLOG_lun,
     &       '(i2,i5,f7.1,f7.0,f6.1,f7.1,f7.0,
     &         f7.0,f8.0,f5.0,f7.1,2f5.1)')
     &         nrcmp,nrpts,bearng,rhomx,rlat,rlon,rrho,
     &        (param(i),i=1,MIN0(nrprm,6))

         if (print_lwf .eq. 2) then

c           Output vs range
            write(lwpcLOG_lun,
     &          '(a,'' component''/
     &          1(''  dist   amplitude  phase  ''))')
     &            label(ncmp:ncmp)
            nl=(nrpts-1)/1+1
            do i1=1,nl
               write(lwpcLOG_lun,
     &             '(1(f7.0,2f10.4))')
     &              (xy(i),amp_lwf(i),phs_lwf(i),i=i1,nrpts,nl)
            end do
         else
     &   if (print_lwf .eq. 3) then

c           Output at receiver only
            if (ncmp .eq. 1) then
              write(lwpcLOG_lun,'('' E  dist   amplitude  phase  ''))')
            end if
            write(lwpcLOG_lun,
     &          '(a,f7.0,2f10.4)')
     &            label(ncmp:ncmp),
     &            xy(nrpts),amp_lwf(nrpts),phs_lwf(nrpts)
         end if
      end if

      if (ncmp .eq. nrcmp) then
         ncmp=0
         nlwf=nlwf+1
         if (nlwf .ge. nrlwf) nlwf=0
      end if

c     End of path
      RETURN

c     End of file
99    end_file=.true.
      nlwf=0

      END      ! READ_LWF

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


text:
      SUBROUTINE POISX0(A, C, LA, LC, MODEL, N, P, QTR, X, YN)
      INTEGER LA, LC, MODEL, N, P
      DOUBLE PRECISION A(LA,N), C(LC), QTR(P), X(P), YN(2,N)
      EXTERNAL DL7ITV, DL7SVX, DL7SVN,DQ7ADR, DR7MDC, DV7SCL, DV7SCP
      DOUBLE PRECISION DL7SVX, DL7SVN, DR7MDC
      INTEGER I
      DOUBLE PRECISION SX, W, WRT, WY, YN1
      DOUBLE PRECISION HALF, ONE, ZERO
      DATA HALF/0.5D+0/, ONE/1.D+0/, ZERO/0.D+0/
C
C *** BODY ***
C
      CALL DV7SCP(LC, C, ZERO)
      CALL DV7SCP(P, QTR, ZERO)
      DO 30 I = 1, N
         W = YN(2,I)
         IF (W .LE. ZERO) GO TO 40
         WRT =  SQRT(W)
         YN1 = YN(1,I) / YN(2,I)
         IF (MODEL .EQ. 2) GO TO 10
            WY = WRT * YN1
            GO TO 20
 10      WY = WRT * DLOG(  MAX(YN1, HALF/W))
 20      CALL DV7SCL(P, X, WRT, A(1,I))
         CALL DQ7ADR(P, QTR, C, X, WY)
 30      CONTINUE
      SX = DL7SVX(P, C, X, X)
      IF (SX .LE. ZERO) GO TO 40
      IF (DL7SVN(P, C, X, X)/SX .LE. DR7MDC(3)) GO TO 40
      CALL DL7ITV(P, X, C, QTR)
      GO TO 999
 40   W = ONE
      IF (MODEL .EQ. 2) W = ZERO
      CALL DV7SCP(P, X, W)
C
 999  RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
!Akash Rajak -"Hello, World!" in Fortran

program hello
  print *, 'Hello, World!'
end program hello


text:
      SUBROUTINE GP3
C
C     GP3 IS THE MAIN CONTROL PROGRAM FOR MODULE GP3.
C     IF PLOAD2 CARDS ARE PRESENT, GP3C IS EXECUTED TO BUILD PLOAD DATA
C     ON SCRATCH FILE 2 (SCR2). GP3A IS EXECUTED TO BUILD THE STATIC
C     LOADS TABLE (SLT). GP3B IS EXECUTED TO BUILD THE GRID POINT
C     TEMPERATURE TABLE (GPTT).
C     GP3D IS EXECUTED TO BUILD THE ELEMENT TEMPERATURE TABLE (ETT) FROM
C     THE GPTT AND ANY TEMPP1,TEMPP2,TEMPP3, AND TEMPRB DATA PRESENT.
C
      INTEGER         BUF1   ,BUF2  ,BUF  ,SYSBUF,PLOAD2,TWO   ,SLT   ,
     1                GPTT  ,GEOM3  ,BUF3 ,STATUS,SPERLK
      COMMON /BLANK / NOGRAV ,NOLOAD,NOTEMP
      COMMON /GP3COM/ GEOM3 ,EQEXIN,GEOM2 ,SLT   ,GPTT  ,SCR1  ,SCR2  ,
     1                BUF1  ,BUF2  ,BUF(50)      ,CARDID(60)   ,IDNO(30)
     2,               CARDDT(60)   ,MASK(60)     ,STATUS(60)   ,NTYPES,
     3                IPLOAD,IGRAV ,PLOAD2(2)    ,LOAD(2)      ,NOPLD2,
     4                TEMP(2)      ,TEMPD(2)     ,TEMPP1(2)    ,
     5                TEMPP2(2)    ,TEMPP3(2)    ,TEMPRB(2)    ,BUF3  ,
     6                PLOAD3(2)    ,IPLD3
      COMMON /SYSTEM/ SYSBUF,SY(93),SPERLK
      COMMON /ZZZZZZ/ Z(1)
      COMMON /TWO   / TWO(32)
C
C     TURN PARAMETERS ON. INITIALIZE BUFFER POINTERS.
C     READ TRAILER ON GEOM3. IF PURGED, EXIT.
C
      CALL DELSET
C
      IF (SPERLK .EQ. 0) GO TO 20
      DO 10 I = 1,60,2
      STATUS(I  ) =-1
   10 STATUS(I+1) = 0
   20 NOLOAD = -1
      NOGRAV = -1
      NOTEMP = -1
      BUF1   = KORSZ(Z) - SYSBUF - 2
      BUF2   = BUF1 - SYSBUF
      BUF3   = BUF2 - SYSBUF - 2
      BUF(1) = GEOM3
      CALL RDTRL (BUF)
      IF (BUF(1) .NE. GEOM3) RETURN
C
C     IF THE SLT IS PURGED, BYPASS THE SLT PHASE OF GP3.
C     OTHERWISE, IF PLOAD2 CARDS PRESENT, EXECUTE GP3C.
C     EXECUTE GP3A TO COMPLETE SLT PHASE.
C
      BUF(7) = SLT
      CALL RDTRL (BUF(7))
      IF (BUF(7) .NE. SLT) GO TO 30
      CALL GP3C
      CALL GP3A
C
C     IF THE GPTT IS NOT PURGED, EXECUTE GP3B TO BUILD IT.
C
   30 BUF(7) = GPTT
      CALL RDTRL (BUF(7))
      IF (BUF(7) .NE. GPTT) RETURN
C
C     GP3B WILL FORM A GPTT ON SCR1 AND THEN GP3D WILL READ SCR1 AND
C     THE TEMPP1,TEMPP2,TEMPP3, AND TEMPRB DATA FROM GEOM3 TO FORM THE
C     ETT (ELEMENT TEMPERATURE TABLE) ON THE OUTPUT FILE GPTT.
C
      CALL GP3B
      CALL GP3D
      RETURN
      END


text:
      PRI
15   &NT 15
     END

max_stars_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler

max_issues_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler

max_forks_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler


text:
      subroutine get_laps_lat_lon(dir,ext,nx,ny,lat,lon,istatus)
c
c *** Modified from code generated by gennet (B. Schwartz, FSL)
c
      implicit none
c
      include 'netcdf.inc'
c
      integer nx,ny
      real   lat(nx,ny),
     .       lon(nx,ny)
c
      integer len,elen,i,istatus
      integer start(10),count(10)
     .       ,ncid,ntp,nvdim,nvs,lenstr,ndsize,rcode
      integer vdims(10)
c
      character*(*) dir
      character*(*) ext
      character*31  dummy
c
      logical exists
c_______________________________________________________________________________
c
c *** Create netcdf file name and check for existence.
c
      len=index(dir,' ')-1
      if (dir(len:len) .ne. '/') then
          dir(len+1:len+1)='/'
          len=len+1
      endif
      elen=index(ext,' ')-1
      inquire(file=dir(1:len)//'static.'//ext(1:elen),exist=exists)
      if (.not. exists) then
         print *,' *** LAPS static file does not exist:  ',
     .           dir(1:len)//'static.'//ext(1:elen)
         istatus=0
         return
      endif
c
c *** Open netcdf file.
c
      ncid=ncopn(dir(1:len)//'static.'//ext(1:elen),ncnowrit,rcode)
c
c *** Statements to fill lat, lon.
c
      call ncvinq(ncid,1,dummy,ntp,nvdim,vdims,nvs,rcode)
      lenstr=1
      do i=1,nvdim
         call ncdinq(ncid,vdims(i),dummy,ndsize,rcode)
         lenstr=lenstr*ndsize
         start(i)=1
         count(i)=ndsize
      enddo
      call ncvgt(ncid,1,start,count,lat,rcode)
c
      call ncvinq(ncid,2,dummy,ntp,nvdim,vdims,nvs,rcode)
      lenstr=1
      do i=1,nvdim
         call ncdinq(ncid,vdims(i),dummy,ndsize,rcode)
         lenstr=lenstr*ndsize
         start(i)=1
         count(i)=ndsize
      enddo
      call ncvgt(ncid,2,start,count,lon,rcode)
c
c *** Close netcdf file.
c
      call ncclos(ncid,rcode)
c
      istatus=1
      return
      end
c
c===============================================================================
c
      subroutine get_laps_wind(dir,i4time,ext,nx,ny,nz,uw,vw,istatus)
c
c *** Modified from code generated by gennet (B. Schwartz, FSL)
c
      implicit none
c
      include 'netcdf.inc'
c
      integer nx,ny,nz
      real   uw(nx,ny,nz)   !LAPS 3d u-wind (m/s)
     .      ,vw(nx,ny,nz)   !LAPS 3d v-wind (m/s)
     .      ,gd(nx,ny,nz)
c
      integer len,elen,i4time,i,j,k,istatus
      integer start(10),count(10),
     .        ncid,ntp,nvdim,nvs,lenstr,ndsize,rcode
      integer vdims(10)
c
      character*(*) dir,ext
      character*500 ldir,lext,fname
      character*31  dummy
      character*9   atime
c
      logical exists
c_______________________________________________________________________________
c
c *** Create netcdf file name and check for existence.
c
      call make_fnam_lp(i4time,atime,istatus)
      ldir=dir
      lext=ext
      len=index(ldir,' ')-1
      if (ldir(len:len) .ne. '/') then
         ldir(len+1:len+1)='/'
         len=len+1
      endif
      fname=ldir(1:len)//atime
      len=len+9
c
      elen=index(lext//' ',' ')-1
      if (elen .gt. 0) then
         fname=fname(1:len)//'.'//lext(1:elen)
         len=len+elen+1
      endif
      inquire(file=fname(1:len),exist=exists)
      if (.not. exists) then
         print *,' *** LAPS 3d wind file does not exist:  ',
     .           fname(1:len)
         istatus=0
         return
      endif
c
c *** Open netcdf file.
c
      ncid=ncopn(fname(1:len),ncnowrit,rcode)
c
c *** Statements to fill 3d u-wind (uw).
c
      call ncvinq(ncid,1,dummy,ntp,nvdim,vdims,nvs,rcode)
      lenstr=1
      do i=1,nvdim
         call ncdinq(ncid,vdims(i),dummy,ndsize,rcode)
         lenstr=lenstr*ndsize
         start(i)=1
         count(i)=ndsize
      enddo
      call ncvgt(ncid,1,start,count,gd,rcode)
      do k=1,nz
      do j=1,ny
      do i=1,nx
         uw(i,j,k)=gd(i,j,nz+1-k)
      enddo
      enddo
      enddo
c
c *** Statements to fill 3d v-wind (vw).
c
      call ncvinq(ncid,2,dummy,ntp,nvdim,vdims,nvs,rcode)
      lenstr=1
      do i=1,nvdim
         call ncdinq(ncid,vdims(i),dummy,ndsize,rcode)
         lenstr=lenstr*ndsize
         start(i)=1
         count(i)=ndsize
      enddo
      call ncvgt(ncid,2,start,count,gd,rcode)
      do k=1,nz
      do j=1,ny
      do i=1,nx
         vw(i,j,k)=gd(i,j,nz+1-k)
      enddo
      enddo
      enddo
c
c *** Close netcdf file.
c
      call ncclos(ncid,rcode)
c
      istatus=1
      return
      end
c
c===============================================================================
c
      subroutine get_laps_2d_sfc(dir,i4time,ext,nx,ny,nfld,fldid
     .                      ,grid,istatus)
c
c *** Modified from code generated by gennet (B. Schwartz, FSL)
c
      implicit none
c
      include 'netcdf.inc'
c
      integer nx,ny,nfld
      real   grid(nx,ny,nfld)
      integer fldid(nfld)
c
      integer len,elen,i4time,i,n,istatus
      integer start(10),count(10)
     .       ,ncid,ntp,nvdim,nvs,lenstr,ndsize,rcode
      integer vdims(10)
c
      character*(*) dir,ext
      character*500 ldir,lext,fname
      character*31  dummy
      character*9   atime
c
      logical exists
c_______________________________________________________________________________
c
c *** Create netcdf file name and check for existence.
c
      call make_fnam_lp(i4time,atime,istatus)
      ldir=dir
      len=index(ldir,' ')-1
      if (ldir(len:len) .ne. '/') then
         ldir(len+1:len+1)='/'
         len=len+1
      endif
      fname=ldir(1:len)//atime
      len=len+9
c
      lext=ext
      elen=index(lext,' ')-1
      if (elen .gt. 0) then
         fname=fname(1:len)//'.'//lext(1:elen)
         len=len+elen+1
      endif
      inquire(file=fname(1:len),exist=exists)
      if (.not. exists) then
         print *,' *** LAPS surface file does not exist:  '
     .          ,fname(1:len)
         istatus=0
         return
      endif
c
c *** Open netcdf file.
c
      ncid=ncopn(fname(1:len),ncnowrit,rcode)
c
      do n=1,nfld
c
c *** Statements to fill requested 2d field.
c
         call ncvinq(ncid,fldid(n),dummy,ntp,nvdim,vdims,nvs,rcode)
         lenstr=1
         do i=1,nvdim
            call ncdinq(ncid,vdims(i),dummy,ndsize,rcode)
            lenstr=lenstr*ndsize
            start(i)=1
            count(i)=ndsize
         enddo
         call ncvgt(ncid,fldid(n),start,count,grid(1,1,n),rcode)
c
      enddo
c
c *** Close netcdf file.
c
      call ncclos(ncid,rcode)
c
      istatus=1
      return
      end
c
c===============================================================================
c
      subroutine get_laps_sfc_elev(dir,ext,nx,ny,sht,istatus)
c
c *** Modified from code generated by gennet (B. Schwartz, FSL)
c
      implicit none
c
      include 'netcdf.inc'
c
      integer nx,ny
      real   sht(nx,ny)
c
      integer len,elen,i,istatus
      integer start(10),count(10)
     .       ,ncid,ntp,nvdim,nvs,lenstr,ndsize,rcode
      integer vdims(10)
c
      character*(*) dir
      character*500 ldir
      character*(*) ext
      character*500 lext
      character*31  dummy
c
      logical exists
c_______________________________________________________________________________
c
c *** Create netcdf file name and check for existence.
c
      ldir=dir
      len=index(ldir,' ')-1
      if (ldir(len:len) .ne. '/') then
         ldir(len+1:len+1)='/'
         len=len+1
      endif
      lext=ext
      elen=index(lext,' ')-1
      inquire(file=ldir(1:len)//'static.'//lext(1:elen),exist=exists)
      if (.not. exists) then
         print *,' *** LAPS static file does not exist:  '
     .          ,ldir(1:len)//'static.'//lext(1:elen)
         istatus=0
         return
      endif
c
c *** Open netcdf file.
c
      ncid=ncopn(ldir(1:len)//'static.'//lext(1:elen),ncnowrit,rcode)
c
c *** Statements to fill sht.
c
      call ncvinq(ncid,3,dummy,ntp,nvdim,vdims,nvs,rcode)
      lenstr=1
      do i=1,nvdim
         call ncdinq(ncid,vdims(i),dummy,ndsize,rcode)
         lenstr=lenstr*ndsize
         start(i)=1
         count(i)=ndsize
      enddo
      call ncvgt(ncid,3,start,count,sht,rcode)
c
c *** Close netcdf file.
c
      call ncclos(ncid,rcode)
c
      istatus=1
      return
      end
c
c ============================================================
c
      subroutine get_laps_3d_analysis_data(i4time,nx,ny,nz
     +,phi,t,u,v,sh,omo,istatus)
c
      implicit none

      integer   nx,ny,nz
      integer   i4time
      integer   istatus
      integer   lendlco
      integer   lends
      integer   lendt
      integer   lendw
      integer   lendsh
      integer   i,j,k
      real    r_missing_data
      real    phi(nx,ny,nz),t(nx,ny,nz)
     .       ,u(nx,ny,nz),v(nx,ny,nz),sh(nx,ny,nz)
     .       ,omo(nx,ny,nz)

      real  , allocatable :: om(:,:,:)

      character*255 tempdir,winddir,sfcdir,shdir,lcodir
      character*125 comment
      character*31  tempext,windext,sfcext,shext,lcoext
      character*10  units
      logical found_lowest

      call get_r_missing_data(r_missing_data,istatus)
      if (istatus .ne. 1) then
         print *,'Error getting r_missing_data...Abort.'
         return
      endif

      shext='lq3'
      tempext='lt1'
      windext='lw3'
      sfcext='lsx'
      lcoext='lco'

      call get_directory(tempext,tempdir,lendt)
      call get_directory(windext,winddir,lendw)
      call get_directory(sfcext,sfcdir,lends)
      call get_directory(shext,shdir,lendsh)
      call get_directory(lcoext,lcodir,lendlco)

      call get_laps_3d(i4time,nx,ny,nz
     1  ,tempext,'ht ',units,comment,phi,istatus)

      if (istatus .ne. 1) then
         print *,'Error getting LAPS height data...Abort.'
         return
      endif
c
c *** Get laps temps
c
      call get_laps_3d(i4time,nx,ny,nz
     1  ,tempext,'t3 ',units,comment,t,istatus)

      if (istatus .ne. 1) then
         print *,'Error getting LAPS temp data...Abort.'
         return
      endif
c
c *** Get laps spec humidity
c
      call get_laps_3d(i4time,nx,ny,nz
     1  ,shext,'sh ',units,comment,sh,istatus)

      if(istatus .ne. 1)then
         print*,'Error getting LAPS sh  data ... Abort.'
         return
      endif
C    The specific humidity field uses the missing value
c    for below ground points, so we need to fill those in by
c    replicating the lowest valid value downward (upward in array
c    space).

      DO j = 1, ny
        DO i = 1, nx
          k = 1
          found_lowest = .false.

          DO WHILE (.NOT. found_lowest)
            IF (sh(i,j,k) .lt. 1.e37) THEN
              found_lowest = .true.
              sh(i,j,1:k) = sh(i,j,k)
            ELSE
              k = k + 1
              IF (k .ge. nz) THEN
                PRINT *, 'No valid SH found in column!'
                PRINT *, 'I/J = ', i,j
                STOP
              ENDIF
            ENDIF
          ENDDO
        ENDDO
      ENDDO

c    Make sure we got it right!
c      print *, 'Min/Max/Center values of SH:'
c      DO k = 1, nz
c        print *, minval(sh(:,:,k)),maxval(sh(:,:,k)),
c     +           sh(nx/2,ny/2,k)
c      ENDDO
c
c *** Get laps cloud omega
c
      call get_laps_3d(i4time,nx,ny,nz
     1  ,lcoext,'com',units,comment,omo,istatus)

      if(istatus .ne. 1)then
         print*,'No LAPS Cld Omega data ....'
         print*,'Initializing omo array with zero'
         call zero3d(omo,nx,ny,nz)
      endif
c
c *** Get laps wind data.
c
      call get_laps_3d(i4time,nx,ny,nz
     1  ,windext,'u3 ',units,comment,u,istatus)

      if (istatus .ne. 1) then
         print *,'Error getting LAPS time 0 u3 data...Abort.'
         return
      endif

      call get_laps_3d(i4time,nx,ny,nz
     1  ,windext,'v3 ',units,comment,v,istatus)

      if (istatus .ne. 1) then
         print *,'Error getting LAPS time 0 v3 data...Abort.'
         return
      endif

      allocate (om(nx,ny,nz))

      call get_laps_3d(i4time,nx,ny,nz
     1	,windext,'om ',units,comment,om,istatus)

      if (istatus .ne. 1) then
	 print *,'Error getting LAPS time 0 v3 data...Abort.'
	 return
      endif
      ! BLS commented this out as a test on 11/19/02
      !where(omo .eq. r_missing_data)omo=om

      deallocate(om)
      return
      end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
!-----subroutine---------------------------------------------g.guidati--
!
     subroutine SPL_EX (xa, n, x, khi, klo)
!
!-----------------------------------------------------------------------
!
!     franco guidati    IAG
!                       Uni Stuttgart
!
!
!     scope             prepare execution of spline interpolation
!                       (natural splines)
!
!.......................................................................
!               declarations
!.......................................................................
USE TIPrecision

IMPLICIT                        NONE

   ! Local variables.

INTEGER(4)                   :: k,klo,khi,n

REAL(DbKi)                   :: a,b,h,x
REAL(DbKi)                   :: xa(n)

!.......................................................................
!               execute spline interpolation
!.......................................................................

      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif

!.......................................................................
!               end of subroutine
!.......................................................................
      return
!***********************************************************************
      end

max_stars_repo_name: CarloSucameli/Openfast-noise

max_issues_repo_name: CarloSucameli/Openfast-noise

max_forks_repo_name: CarloSucameli/Openfast-noise


!  > using minimization algorithm to search spin-orbital strength
   subroutine fitsoc

      use para
      implicit none

      integer :: i
      integer :: j
      integer :: l
      integer :: n
      integer :: ik
      real(dp) :: k(3)
      real(dp), allocatable :: x(:)

      integer :: npara

      integer  :: maxfun

      ! output frequency
      integer  :: iprint

      integer  :: iexit

      integer  :: Nitermax
      integer  :: iter

      real(dp) :: ftol

      real(dp), allocatable :: W(:)
      complex(dp), allocatable :: Hamk_soc(:, :)

      real(Dp),allocatable :: y(:)
      real(Dp),allocatable :: p(:, :)

      interface
         function func(n, x)
         implicit none
         integer, parameter :: dp= kind(1d0)
         integer, intent(in) :: n
         real(dp), dimension(n), intent(in) :: x
         real(dp) :: func
         end function func
      end interface

      allocate(W(num_wann_soc))
      allocate(Hamk_soc(num_wann_soc, num_wann_soc))
      W=0d0
      Hamk_soc= 0d0

      n= 3

      knv3= Nk*Nk*Nk
      allocate(kpoints(3, knv3))
      kpoints= 0d0

      ik= 0
      do i=1, Nk
         do j=1, Nk
            do l=1, Nk
               ik= ik+ 1
               kpoints(1, ik)= dble(i)/Nk
               kpoints(2, ik)= dble(j)/Nk
               kpoints(3, ik)= dble(l)/Nk
            enddo !l
         enddo !j
      enddo !i

      !> using kline mode
      deallocate(kpoints)
      allocate(kpoints(3, nk3_band))
      kpoints= k3points

      knv3= nk3_band
      num_bands_DFT= 82

      allocate(eigval_soc(num_wann_soc, knv3))
      allocate(eigval_soc_DFT(num_bands_DFT, knv3))
      allocate(eigval_nsoc(num_wann_soc, knv3))
      eigval_soc= 0d0
      eigval_soc_DFT= 0d0

      !> get eigenvalues for soc hamiltonian
     !do ik=1, knv3
     !   k= kpoints(:, ik)
     !   W=0d0
     !   Hamk_soc= 0d0
     !   call Hamk_bulk_soc(k, Hamk_soc)
     !   call eigensystem_c('N', 'U', num_wann_soc, Hamk_soc, W)
     !   eigval_soc(:, ik)= W
     !enddo

      !> get eigenvalues for soc hamiltonian from DFT
      call EigvalFromDFT
     !eigval_soc=eigval_soc_DFT(19:18+num_wann_soc, :)
      eigval_soc(1:num_bands_DFT-17+1, :)=eigval_soc_DFT(17:num_bands_DFT, :)

      npara= Num_atom_type*2+ 1

      allocate(x(npara))
      allocate(y(npara+1))
      allocate(p(npara+1, npara))
      x=0d0
      y=0d0
      p=0d0

      !> initial value for spin-orbital coupling
      x(1:Num_atom_type)= lambda_p(1:Num_atom_type)
      x(1+Num_atom_type:2*Num_atom_type)= lambda_d(1:Num_atom_type)
      x(Num_atom_type*2+ 1)= E_fermi

      Nitermax=10000
      do i=1, npara+1
         p(i, :)= x
      enddo
      do i=1, npara
         p(i+1, i)= p(i+1, i)+0.02d0*i
      enddo

      do i=1, npara+1
         y(i)= func(npara, p(i, :))
      enddo

      ftol=1d-6
      call amoeba(npara, p, y, ftol, func, iter)

      x= p(1, :)
      write(*,*)'  Iprint',iprint, 'iter', iter
      write(*, '(3f10.4)')x
      write(*, *)' '

      open(unit=14, file='bulkek-soc.dat')
     !do i=1, num_wann_soc
      do i=1, num_bands_DFT-16
         do ik=1, knv3
            write(14, '(1000f19.9)')k3len(ik),eigval_soc(i, ik)
         enddo
         write(14, *)' '
      enddo
      close(14)

      open(unit=15, file='bulkek-nsoc.dat')
      do i=1, num_wann_soc
         do ik=1, knv3
            write(15, '(1000f19.9)')k3len(ik),eigval_nsoc(i, ik)- X(3)
         enddo
         write(15, *)' '
      enddo
      close(15)

      return
   end subroutine fitsoc

   subroutine fitsoc_p

      use para
      implicit none

      integer :: i
      integer :: j
      integer :: l
      integer :: n
      integer :: ik
      real(dp) :: k(3)
      real(dp), allocatable :: x(:)

      integer :: npara

      integer  :: maxfun

      ! output frequency
      integer  :: iprint

      integer  :: iexit

      integer  :: Nitermax
      integer  :: iter

      real(dp) :: ftol

      real(dp), allocatable :: W(:)
      complex(dp), allocatable :: Hamk_soc(:, :)

      real(Dp),allocatable :: y(:)
      real(Dp),allocatable :: p(:, :)

      interface
         function func(n, x)
         implicit none
         integer, parameter :: dp= kind(1d0)
         integer, intent(in) :: n
         real(dp), dimension(n), intent(in) :: x
         real(dp) :: func
         end function func
      end interface

      allocate(W(num_wann_soc))
      allocate(Hamk_soc(num_wann_soc, num_wann_soc))
      W=0d0
      Hamk_soc= 0d0

      knv3= Nk*Nk*Nk
      allocate(kpoints(3, knv3))
      kpoints= 0d0

      ik= 0
      do i=1, Nk
         do j=1, Nk
            do l=1, Nk
               ik= ik+ 1
               kpoints(1, ik)= dble(i)/Nk
               kpoints(2, ik)= dble(j)/Nk
               kpoints(3, ik)= dble(l)/Nk
            enddo !l
         enddo !j
      enddo !i

      !> using kline mode
      deallocate(kpoints)
      allocate(kpoints(3, nk3_band))
      kpoints= k3points

      knv3= nk3_band
      num_bands_DFT= 82
      allocate(weight(num_wann_soc, knv3))
      weight=1d0

      allocate(eigval_soc(num_wann_soc, knv3))
      allocate(eigval_soc_DFT(num_bands_DFT, knv3))
      allocate(eigval_nsoc(num_wann_soc, knv3))
      eigval_soc= 0d0
      eigval_soc_DFT= 0d0

      !> get eigenvalues for soc hamiltonian
      do ik=1, knv3
         k= kpoints(:, ik)
         if (abs(k(1))<0.2d0 .and. abs(k(2))<0.2d0 .and. abs(k(3))<0.2d0)then
            weight(:, ik)=30d0
         endif
         if (abs(k(1))<0.1d0 .and. abs(k(2))<0.1d0 .and. abs(k(3))<0.1d0)then
            weight(:, ik)=300d0
         endif
         if (abs(k(1))<0.02d0 .and. abs(k(2))<0.02d0 .and. abs(k(3))<0.02d0)then
            weight(:, ik)=3000d0
         endif
         W=0d0
         Hamk_soc= 0d0
         call Hamk_bulk_soc(k, Hamk_soc)
         call eigensystem_c('N', 'U', num_wann_soc, Hamk_soc, W)
         eigval_soc(:, ik)= W
      enddo

      !> get eigenvalues for soc hamiltonian from DFT
     !call EigvalFromDFT
     !eigval_soc=eigval_soc_DFT(19:18+num_wann_soc, :)
     !eigval_soc(1:num_bands_DFT-17+1, :)=eigval_soc_DFT(17:num_bands_DFT, :)

      npara= 2*Num_atom_type+ 1
      allocate(x(npara))
      allocate(y(npara+1))
      allocate(p(npara+1, npara))
      x= 0d0
      y=0d0
      p=0d0

      !> initial value for spin-orbital coupling
      x(1:Num_atom_type)= lambda_p(1:Num_atom_type)
      x(1+Num_atom_type:2*Num_atom_type)= lambda_d(1:Num_atom_type)
      x(Num_atom_type*2+ 1)= E_fermi

      Nitermax=10000
      do i=1, npara+1
         p(i, :)= x
      enddo
      do i=1, npara
         p(i+1, i)= p(i+1, i)+0.02d0*i
      enddo

      do i=1, npara+1
         y(i)= func(npara, p(i, :))
      enddo

      ftol=1d-6
      call amoeba(npara, p, y, ftol, func, iter)

      x= p(1, :)
      write(*,*)'iter', iter
      write(*,*)'p(atoms)'
      write(*, '(30f10.4)')lambda_p
      write(*,*)'d(atoms)'
      write(*, '(30f10.4)')lambda_d
      write(*,*)'efermi'
      write(*, *)x(Num_atom_type*2+ 1)

      open(unit=14, file='bulkek-soc.dat')
      do i=1, num_wann_soc
     !do i=1, num_bands_DFT-16
         do ik=1, knv3
            write(14, '(1000f19.9)')k3len(ik),eigval_soc(i, ik)
         enddo
         write(14, *)' '
      enddo
      close(14)

      open(unit=15, file='bulkek-nsoc.dat')
      do i=1, num_wann_soc
         do ik=1, knv3
            write(15, '(1000f19.9)')k3len(ik),eigval_nsoc(i, ik)- X(npara)
         enddo
         write(15, *)' '
      enddo
      close(15)

      return
   end subroutine fitsoc_p

   !> take eigenvalues from DFT
   !> read ek_nxy.dat
   subroutine EigvalFromDFT

      use para
      implicit none

      integer :: ik
      real(dp) :: r1

      real(dp), allocatable :: r2(:, :)

      !> only for WTe2
      allocate(r2(num_bands_DFT,180))
      r2=0d0

      open(unit=140, file='ek_nxy.dat')

      do ik=1, 100
         read(140, *)r1, r2(:, ik)
      enddo

      do ik=1, nk3_band
         eigval_soc_DFT(:, ik)= r2(:, ik)
      enddo

      close(140)

      return
   end subroutine EigvalFromDFT


text:
fl
\
\ fswrite DoubleMath.f
\ 100 fwrite DoubleMath.f
\
\
1 wconstant build_DoubleMath
\
[ifndef dum*
lockdict create dum* forthentry
$C_a_lxasm w, h137  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2Wi]ZB l, z1SyLI[ l, z2Wi]RB l, z1SyLI[ l, z2Wy]b0 l, z2Wy]j0 l, z1SyLIZ l, z2WyPb0 l,
z2WyPj0 l, z2WyPr0 l, z2WyQ00 l, zeyPW1 l, zkyPO1 l, z26QPO0 l, z1SJ04] l, z21iPfm l,
z39iPnn l, z39iPvo l, z38iQ4p l, zgy]O1 l, zoy]W1 l, zoy]b1 l, zny]j1 l, z1SL04U l,
z2WiPRD l, z1SyJQL l, z2WiPRE l, z1SyJQL l, z2WiPRF l, z1SyJQL l, z2WiPRG l, z1SV01X l,
0 l, 0 l, 0 l, 0 l,
freedict
]

[ifndef dum/mod
lockdict create dum/mod forthentry
$C_a_lxasm w, h13E  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2WiPuB l, z1SyLI[ l, z2WiaBB l, z1SyLI[ l, z2WiPmB l, z1SyLI[ l, z2WiPeB l, z1SyLI[ l,
z1SyLIZ l, z2WyQA0 l, z2WyaG0 l, z2WyaO0 l, zgyPO1 l, zoyPW1 l, zoyPb1 l, zoyPj1 l,
zoyaG1 l, zoyaO1 l, z1SS04d l, z27FaKs l, z3DFaRF l, z2WvQ00 l, z1SS04g l, z27iaKs l,
z3CiaRF l, z2WyQ01 l, zkyQ01 l, zoyaW1 l, znyab1 l, z3[yQCV l, z2WiPSt l, z1SyJQL l,
z2WiPSu l, z1SyJQL l, z2WiPSv l, z1SyJQL l, z2WiPSw l, z1SV01X l, 0 l, 0 l,
0 l, 0 l, 0 l,
freedict
]

[ifndef d+
lockdict create d+ forthentry
$C_a_lxasm w, h11D  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2WiPmB l, z1SyLI[ l, z2WiPeB l, z1SyLI[ l, z1SyLIZ l, z21iPRD l, z3OiPZE l, z1SyJQL l,
z2WiPRC l, z1SV01X l,
freedict
]

[ifndef d-
lockdict create d- forthentry
$C_a_lxasm w, h11D  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2WiPmB l, z1SyLI[ l, z2WiPeB l, z1SyLI[ l, z1SyLIZ l, z27iPRD l, z3SiPZE l, z1SyJQL l,
z2WiPRC l, z1SV01X l,
freedict
]

[ifndef du>
lockdict create du> forthentry
$C_a_lxasm w, h11D  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2WiPmB l, z1SyLI[ l, z2WiPeB l, z1SyLI[ l, z1SyLIZ l, z27FPRD l, z3FFPZE l, z2WyPO0 l,
z2WXPR6 l, z1SV01X l,
freedict
]

[ifndef du<
lockdict create du< forthentry
$C_a_lxasm w, h11D  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2WiPmB l, z1SyLI[ l, z2WiPeB l, z1SyLI[ l, z1SyLIZ l, z27FPRD l, z3FFPZE l, z2WyPO0 l,
z2WfPR6 l, z1SV01X l,
freedict
]

[ifndef d=
lockdict create d= forthentry
$C_a_lxasm w, h11D  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2WiPmB l, z1SyLI[ l, z2WiPeB l, z1SyLI[ l, z1SyLIZ l, z27FPRD l, z3FFPZE l, z2WyPO0 l,
z2WdPR6 l, z1SV01X l,
freedict
]

[ifndef du>=
lockdict create du>= forthentry
$C_a_lxasm w, h11D  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2WiPmB l, z1SyLI[ l, z2WiPeB l, z1SyLI[ l, z1SyLIZ l, z27FPRD l, z3FFPZE l, z2WyPO0 l,
z2WZPR6 l, z1SV01X l,
freedict
]

[ifndef du<=
lockdict create du<= forthentry
$C_a_lxasm w, h11D  h113  1- tuck - h9 lshift or here W@ alignl h10 lshift or l,
z2WiPmB l, z1SyLI[ l, z2WiPeB l, z1SyLI[ l, z1SyLIZ l, z27FPRD l, z3FFPZE l, z2WyPO0 l,
z2WhPR6 l, z1SV01X l,
freedict
]

\ du* ( u1lo u1hi u2lo u2hi -- u1*u2lo u1*u2hi ) u1 multiplied by u2
[ifndef du*
: du*
	dum* 2drop
;
]

\ du/mod ( u1lo u1hi u2lo u2hi -- remainderlo remainderhi quotientlo quotienthi) \ unsigned divide & mod  u1 divided by u2
[ifndef du/mod
: du/mod
	0 rot2 0 rot2 dum/mod
;
]

\
\
\ du/ ( u1lo u1hi u2lo u2hi -- u1/u2lo u1/u2hi) u1 divided by u2
[ifndef du/
: du/
	du/mod rot drop rot drop
;
]

\ du*/mod ( u1lo u1hi u2lo u2hi u3lo u3hi -- u4lo u4hi u5lo u5hi ) u5 = (u1*u2)/u3, u4 is the remainder.
\         Uses a 128bit intermediate result.
[ifndef du*/mod
: du*/mod
	>r >r dum* r> r> dum/mod
;
]

\
\ du*/ ( u1lo u1hi u2lo u2hi u3lo u3hi -- u4lo u4hi ) u4 = (u1*u2)/u3. Uses a 128bit intermediate result.
[ifndef du*/
: du*/
	>r >r dum* r> r> dum/mod rot drop rot drop
;
]

\
\
\ d# ( n1lo n1hi -- n2lo n2hi ) divide n1 by base and convert the remainder to a char and append to the output
[ifndef d#
: d#
	base W@ 0 du/mod rot drop rot tochar -1 >out W+! pad>out C!
;
]

\
\
\ d#s ( n1lo n1hi -- 0 ) execute # until the remainder is 0
[ifndef d#s
: d#s
	begin
		d# 2dup 0= swap 0= and
	until
	drop
;
]

\
\ du. ( n1lo n1hi -- ) prints the unsigned number on the top of the stack
[ifndef du.
: du.
	<# d#s #> .cstr space
;
]

\
\ d. ( n1lo n1hi -- ) prints the signed number on the top of the stack
[ifndef d.
: d.
	dup 0<
	if
		h2D emit negate
	then
	du.
;
]

\
\ dL! ( nlo nhi addr -- )
[ifndef dL!
: dL!
	tuck 4+ L! L!
;
]
\
\ dL@ ( addr -- nlo nhi)
[ifndef dL@
: dL@
	dup L@ swap 4+ L@
;
]

\ dswap ( n1lo n1hi n2lo n2hi -- n2lo n2hi n1lo n1hi)
[ifndef dswap
: dswap
	 0 ST@ 3 ST@ 1 ST! 2 ST!
	 1 ST@ 4 ST@ 2 ST! 3 ST!
;
]
\ drot ( n1lo n1hi n2lo n2hi n3lo n3hi -- n2lo n2hi n3lo n3hi  n1lo n1hi)
[ifndef drot
: drot
	 0 ST@ 3 ST@ 6 ST@ 2 ST! 5 ST! 2 ST!
	 1 ST@ 4 ST@ 7 ST@ 3 ST! 6 ST! 3 ST!
;
]
\ d>u ( u1lo u1hi -- u1lo )
[ifndef d>u
: d>u
	drop
;
]
\ u>d ( u1lo -- u1lo u1hi )
[ifndef u>d
: u>d
	0
;
]
\ ddup ( n1lo n1hi --	n1lo nihi n1lo n1hi)
[ifndef ddup
: ddup
	2dup
;
]
\ ddrop ( n1lo n1hi --	)
[ifndef ddrop
: ddrop
	2drop
;
]

\ dover ( n1lo n1hi n2lo n2hi -- n1lo n1hi n2lo n2hi n1lo n1hi)
[ifndef dover
: dover
	3 ST@ 3 ST@
;
]

\ d2dup ( n1lo nihi n2lo n2hi -- n1lo n1hi n2lo n2hi n1lo nihi n2lo n2hi)
[ifndef d2dup
: d2dup
	dover dover
;
]

...


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2015 Guido Dhondt
!
!     This program 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(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine ts_calc(xflow,Tt,Pt,kappa,r,A,Ts,icase)
!
!     Calculation of the static temperature Ts from the total
!     temperature Tt, the total pressure Pt and the mass flow xflow
!
!     this subroutine solves the implicit equation
!     xflow*dsqrt(Tt)/(A*Pt)-C*(Tt/Ts)**expon*(Tt/Ts-1)**0.5d0=0.d0
!
!          expon=-0.5d0*(kappa+1.d0)/(kappa-1.d0)
!          C=dsqrt(2.d0/r*kappa/(kappa-1.d0))
!
!     xflow is the mass flow for 360 degrees, sign is irrelevant
!     for this routine
!
!     author: Yannick Muller
!
      implicit none
!
      integer icase,i
!
      real*8 xflow,Tt,Pt,Ts,kappa,r,f,df,A,expon,Ts_old,C,TtzTs,
     &     deltaTs,TtzTs_crit, Qred_crit,Qred,h1,h2,h3
!
      intent(in) xflow,Tt,Pt,kappa,r,A,icase
!
      intent(inout) Ts
!
      expon=-0.5d0*(kappa+1.d0)/(kappa-1.d0)
!
      C=dsqrt(2.d0/r*kappa/(kappa-1.d0))
!
!     f=xflow*dsqrt(Tt)/(A*Pt)-C*(Tt/Ts)**expon*(Tt/Ts-1)**0.5d0
!
!     df=-C*Tt/Ts**expon*(expon/Ts*(Tt/Ts-1)**0.5d0
!     &     -0.5d0*Tt/Ts/Ts*(Tt/Ts-1.d0)**(-0.5d0))
!
      Ts_old=Tt
!
      if(dabs(xflow).le.1e-9) then
         Ts=Tt
         return
      endif
!
      Qred=abs(xflow)*dsqrt(Tt)/(A*Pt)
!
!     optimised estimate of T static
!
      Ts=Tt/(1+(Qred**2/C**2))
!
!     adiabatic
!
      if(icase.eq.0) then
         TtzTs_crit=(kappa+1.d0)/2.d0
!
!     isothermal
!
      else
         TtzTs_crit=(1d0+(kappa-1.d0)/(2.d0*kappa))
      endif
!
      Qred_crit=C*(TtzTs_crit)**expon*(Ttzts_crit-1.d0)**0.5d0
!
      if(Qred.ge.Qred_crit) then
         Ts=Tt/TtzTs_crit
         return
      endif
!
      i=0
!
!     start of the Newton-Raphson-Procedure to solve the nonlinear
!     equation
!
      do
         i=i+1
         Ttzts=Tt/Ts
         h1=Ttzts-1.d0
         h2=dsqrt(h1)
         h3=Ttzts**expon
!
         f=C*h2*h3
!
         df=f*(expon+0.5d0*Ttzts/h1)/Ts
!
         f=Qred-f
         deltaTs=-f/df
!
         Ts=Ts+deltaTs
!
         if((((dabs(Ts-Ts_old)/ts_old).le.1.E-8))
     &        .or.((dabs(Ts-Ts_old)).le.1.E-10)) then
            exit
         else if(i.gt.20) then
            Ts=0.9*Tt
            exit
         endif
         Ts_old=Ts
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
!>
!! @file
!! @brief Fortran interface to `thermal_funcs`
!!
!! Build with `make fortran` to create `./bin/fortran_example`

!>
!! @brief Program that prints bosonic thermal function
program thermal_funcs

  use iso_c_binding

  implicit none

  !>
  !! @brief Fortran interface to `thermal_funcs`
  interface
    function j_b(y_squared) bind (c)
      use iso_c_binding
      real(c_double) :: y_squared
      real(c_double) :: j_b
    end function
  end interface

  real(c_double) :: y_squared = 10.
  real(c_double) :: r

  r = j_b(y_squared)  ! This calls the c interface
  write(*, *) "j_b(", y_squared, ") = ", r

end program