| 1 | C==============================================================================
|
|---|
| 2 | C
|
|---|
| 3 | C LIBRARY ROUTINES FOR SPECTRAL METHODS
|
|---|
| 4 | C
|
|---|
| 5 | C March 1989
|
|---|
| 6 | C
|
|---|
| 7 | C For questions, comments or suggestions, please contact:
|
|---|
| 8 | C
|
|---|
| 9 | C Einar Malvin Ronquist
|
|---|
| 10 | C Room 3-243
|
|---|
| 11 | C Department of Mechanical Engineering
|
|---|
| 12 | C Massachusetts Institute of Technology
|
|---|
| 13 | C 77 Massachusetts Avenue
|
|---|
| 14 | C Cambridge, MA 0299
|
|---|
| 15 | C U.S.A.
|
|---|
| 16 | C
|
|---|
| 17 | C------------------------------------------------------------------------------
|
|---|
| 18 | C
|
|---|
| 19 | C ABBRIVIATIONS:
|
|---|
| 20 | C
|
|---|
| 21 | C M - Set of mesh points
|
|---|
| 22 | C Z - Set of collocation/quadrature points
|
|---|
| 23 | C W - Set of quadrature weights
|
|---|
| 24 | C H - Lagrangian interpolant
|
|---|
| 25 | C D - Derivative operator
|
|---|
| 26 | C I - Interpolation operator
|
|---|
| 27 | C GL - Gauss Legendre
|
|---|
| 28 | C GLL - Gauss-Lobatto Legendre
|
|---|
| 29 | C GJ - Gauss Jacobi
|
|---|
| 30 | C GLJ - Gauss-Lobatto Jacobi
|
|---|
| 31 | C
|
|---|
| 32 | C
|
|---|
| 33 | C MAIN ROUTINES:
|
|---|
| 34 | C
|
|---|
| 35 | C Points and weights:
|
|---|
| 36 | C
|
|---|
| 37 | C ZWGL Compute Gauss Legendre points and weights
|
|---|
| 38 | C ZWGLL Compute Gauss-Lobatto Legendre points and weights
|
|---|
| 39 | C ZWGJ Compute Gauss Jacobi points and weights (general)
|
|---|
| 40 | C ZWGLJ Compute Gauss-Lobatto Jacobi points and weights (general)
|
|---|
| 41 | C
|
|---|
| 42 | C Lagrangian interpolants:
|
|---|
| 43 | C
|
|---|
| 44 | C HGL Compute Gauss Legendre Lagrangian interpolant
|
|---|
| 45 | C HGLL Compute Gauss-Lobatto Legendre Lagrangian interpolant
|
|---|
| 46 | C HGJ Compute Gauss Jacobi Lagrangian interpolant (general)
|
|---|
| 47 | C HGLJ Compute Gauss-Lobatto Jacobi Lagrangian interpolant (general)
|
|---|
| 48 | C
|
|---|
| 49 | C Derivative operators:
|
|---|
| 50 | C
|
|---|
| 51 | C DGLL Compute Gauss-Lobatto Legendre derivative matrix
|
|---|
| 52 | C DGLLGL Compute derivative matrix for a staggered mesh (GLL->GL)
|
|---|
| 53 | C DGJ Compute Gauss Jacobi derivative matrix (general)
|
|---|
| 54 | C DGLJ Compute Gauss-Lobatto Jacobi derivative matrix (general)
|
|---|
| 55 | C DGLJGJ Compute derivative matrix for a staggered mesh (GLJ->GJ) (general)
|
|---|
| 56 | C
|
|---|
| 57 | C Interpolation operators:
|
|---|
| 58 | C
|
|---|
| 59 | C IGLM Compute interpolation operator GL -> M
|
|---|
| 60 | C IGLLM Compute interpolation operator GLL -> M
|
|---|
| 61 | C IGJM Compute interpolation operator GJ -> M (general)
|
|---|
| 62 | C IGLJM Compute interpolation operator GLJ -> M (general)
|
|---|
| 63 | C
|
|---|
| 64 | C Other:
|
|---|
| 65 | C
|
|---|
| 66 | C PNLEG Compute Legendre polynomial of degree N
|
|---|
| 67 | C PNDLEG Compute derivative of Legendre polynomial of degree N
|
|---|
| 68 | C
|
|---|
| 69 | C Comments:
|
|---|
| 70 | C
|
|---|
| 71 | C Note that many of the above routines exist in both single and
|
|---|
| 72 | C double precision. If the name of the single precision routine is
|
|---|
| 73 | C SUB, the double precision version is called SUBD. In most cases
|
|---|
| 74 | C all the "low-level" arithmetic is done in double precision, even
|
|---|
| 75 | C for the single precsion versions.
|
|---|
| 76 | C
|
|---|
| 77 | C Useful references:
|
|---|
| 78 | C
|
|---|
| 79 | C [1] Gabor Szego: Orthogonal Polynomials, American Mathematical Society,
|
|---|
| 80 | C Providence, Rhode Island, 1939.
|
|---|
| 81 | C [2] Abramowitz & Stegun: Handbook of Mathematical Functions,
|
|---|
| 82 | C Dover, New York, 1972.
|
|---|
| 83 | C [3] Canuto, Hussaini, Quarteroni & Zang: Spectral Methods in Fluid
|
|---|
| 84 | C Dynamics, Springer-Verlag, 1988.
|
|---|
| 85 | C
|
|---|
| 86 | C
|
|---|
| 87 | C==============================================================================
|
|---|
| 88 | C
|
|---|
| 89 | C--------------------------------------------------------------------
|
|---|
| 90 | SUBROUTINE ZWGL (Z,W,NP)
|
|---|
| 91 | C--------------------------------------------------------------------
|
|---|
| 92 | C
|
|---|
| 93 | C Generate NP Gauss Legendre points (Z) and weights (W)
|
|---|
| 94 | C associated with Jacobi polynomial P(N)(alpha=0,beta=0).
|
|---|
| 95 | C The polynomial degree N=NP-1.
|
|---|
| 96 | C Z and W are in single precision, but all the arithmetic
|
|---|
| 97 | C operations are done in double precision.
|
|---|
| 98 | C
|
|---|
| 99 | C--------------------------------------------------------------------
|
|---|
| 100 | REAL Z(1),W(1)
|
|---|
| 101 | ALPHA = 0.
|
|---|
| 102 | BETA = 0.
|
|---|
| 103 | CALL ZWGJ (Z,W,NP,ALPHA,BETA)
|
|---|
| 104 | RETURN
|
|---|
| 105 | END
|
|---|
| 106 | C
|
|---|
| 107 | SUBROUTINE ZWGLL (Z,W,NP)
|
|---|
| 108 | C--------------------------------------------------------------------
|
|---|
| 109 | C
|
|---|
| 110 | C Generate NP Gauss-Lobatto Legendre points (Z) and weights (W)
|
|---|
| 111 | C associated with Jacobi polynomial P(N)(alpha=0,beta=0).
|
|---|
| 112 | C The polynomial degree N=NP-1.
|
|---|
| 113 | C Z and W are in single precision, but all the arithmetic
|
|---|
| 114 | C operations are done in double precision.
|
|---|
| 115 | C
|
|---|
| 116 | C--------------------------------------------------------------------
|
|---|
| 117 | REAL Z(1),W(1)
|
|---|
| 118 | ALPHA = 0.
|
|---|
| 119 | BETA = 0.
|
|---|
| 120 | CALL ZWGLJ (Z,W,NP,ALPHA,BETA)
|
|---|
| 121 | RETURN
|
|---|
| 122 | END
|
|---|
| 123 | C
|
|---|
| 124 | SUBROUTINE ZWGJ (Z,W,NP,ALPHA,BETA)
|
|---|
| 125 | C--------------------------------------------------------------------
|
|---|
| 126 | C
|
|---|
| 127 | C Generate NP GAUSS JACOBI points (Z) and weights (W)
|
|---|
| 128 | C associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
|
|---|
| 129 | C The polynomial degree N=NP-1.
|
|---|
| 130 | C Single precision version.
|
|---|
| 131 | C
|
|---|
| 132 | C--------------------------------------------------------------------
|
|---|
| 133 | PARAMETER (NMAX=84)
|
|---|
| 134 | PARAMETER (lzd = NMAX)
|
|---|
| 135 | REAL*8 ZD(lzd),WD(lzd),APHAD,BETAD
|
|---|
| 136 | REAL Z(1),W(1),ALPHA,BETA
|
|---|
| 137 | C
|
|---|
| 138 | NPMAX = lzd
|
|---|
| 139 | IF (NP.GT.NPMAX) THEN
|
|---|
| 140 | WRITE (6,*) 'Too large polynomial degree in ZWGJ'
|
|---|
| 141 | WRITE (6,*) 'Maximum polynomial degree is',NMAX
|
|---|
| 142 | WRITE (6,*) 'Here NP=',NP
|
|---|
| 143 | call exitt
|
|---|
| 144 | ENDIF
|
|---|
| 145 | ALPHAD = ALPHA
|
|---|
| 146 | BETAD = BETA
|
|---|
| 147 | CALL ZWGJD (ZD,WD,NP,ALPHAD,BETAD)
|
|---|
| 148 | DO 100 I=1,NP
|
|---|
| 149 | Z(I) = ZD(I)
|
|---|
| 150 | W(I) = WD(I)
|
|---|
| 151 | 100 CONTINUE
|
|---|
| 152 | RETURN
|
|---|
| 153 | END
|
|---|
| 154 | C
|
|---|
| 155 | SUBROUTINE ZWGJD (Z,W,NP,ALPHA,BETA)
|
|---|
| 156 | C--------------------------------------------------------------------
|
|---|
| 157 | C
|
|---|
| 158 | C Generate NP GAUSS JACOBI points (Z) and weights (W)
|
|---|
| 159 | C associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
|
|---|
| 160 | C The polynomial degree N=NP-1.
|
|---|
| 161 | C Double precision version.
|
|---|
| 162 | C
|
|---|
| 163 | C--------------------------------------------------------------------
|
|---|
| 164 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 165 | REAL*8 Z(1),W(1),ALPHA,BETA
|
|---|
| 166 | C
|
|---|
| 167 | N = NP-1
|
|---|
| 168 | DN = ((N))
|
|---|
| 169 | ONE = 1.
|
|---|
| 170 | TWO = 2.
|
|---|
| 171 | APB = ALPHA+BETA
|
|---|
| 172 | C
|
|---|
| 173 | IF (NP.LE.0) THEN
|
|---|
| 174 | WRITE (6,*) 'ZWGJD: Minimum number of Gauss points is 1',np
|
|---|
| 175 | call exitt
|
|---|
| 176 | ENDIF
|
|---|
| 177 | IF ((ALPHA.LE.-ONE).OR.(BETA.LE.-ONE)) THEN
|
|---|
| 178 | WRITE (6,*) 'ZWGJD: Alpha and Beta must be greater than -1'
|
|---|
| 179 | call exitt
|
|---|
| 180 | ENDIF
|
|---|
| 181 | C
|
|---|
| 182 | IF (NP.EQ.1) THEN
|
|---|
| 183 | Z(1) = (BETA-ALPHA)/(APB+TWO)
|
|---|
| 184 | W(1) = GAMMAF(ALPHA+ONE)*GAMMAF(BETA+ONE)/GAMMAF(APB+TWO)
|
|---|
| 185 | $ * TWO**(APB+ONE)
|
|---|
| 186 | RETURN
|
|---|
| 187 | ENDIF
|
|---|
| 188 | C
|
|---|
| 189 | CALL JACG (Z,NP,ALPHA,BETA)
|
|---|
| 190 | C
|
|---|
| 191 | NP1 = N+1
|
|---|
| 192 | NP2 = N+2
|
|---|
| 193 | DNP1 = ((NP1))
|
|---|
| 194 | DNP2 = ((NP2))
|
|---|
| 195 | FAC1 = DNP1+ALPHA+BETA+ONE
|
|---|
| 196 | FAC2 = FAC1+DNP1
|
|---|
| 197 | FAC3 = FAC2+ONE
|
|---|
| 198 | FNORM = PNORMJ(NP1,ALPHA,BETA)
|
|---|
| 199 | RCOEF = (FNORM*FAC2*FAC3)/(TWO*FAC1*DNP2)
|
|---|
| 200 | DO 100 I=1,NP
|
|---|
| 201 | CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,NP2,ALPHA,BETA,Z(I))
|
|---|
| 202 | W(I) = -RCOEF/(P*PDM1)
|
|---|
| 203 | 100 CONTINUE
|
|---|
| 204 | RETURN
|
|---|
| 205 | END
|
|---|
| 206 | C
|
|---|
| 207 | SUBROUTINE ZWGLJ (Z,W,NP,ALPHA,BETA)
|
|---|
| 208 | C--------------------------------------------------------------------
|
|---|
| 209 | C
|
|---|
| 210 | C Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W)
|
|---|
| 211 | C associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
|
|---|
| 212 | C The polynomial degree N=NP-1.
|
|---|
| 213 | C Single precision version.
|
|---|
| 214 | C
|
|---|
| 215 | C--------------------------------------------------------------------
|
|---|
| 216 | PARAMETER (NMAX=84)
|
|---|
| 217 | PARAMETER (lzd = NMAX)
|
|---|
| 218 | REAL*8 ZD(lzd),WD(lzd),ALPHAD,BETAD
|
|---|
| 219 | REAL Z(1),W(1),ALPHA,BETA
|
|---|
| 220 | C
|
|---|
| 221 | NPMAX = lzd
|
|---|
| 222 | IF (NP.GT.NPMAX) THEN
|
|---|
| 223 | WRITE (6,*) 'Too large polynomial degree in ZWGLJ'
|
|---|
| 224 | WRITE (6,*) 'Maximum polynomial degree is',NMAX
|
|---|
| 225 | WRITE (6,*) 'Here NP=',NP
|
|---|
| 226 | call exitt
|
|---|
| 227 | ENDIF
|
|---|
| 228 | ALPHAD = ALPHA
|
|---|
| 229 | BETAD = BETA
|
|---|
| 230 | CALL ZWGLJD (ZD,WD,NP,ALPHAD,BETAD)
|
|---|
| 231 | DO 100 I=1,NP
|
|---|
| 232 | Z(I) = ZD(I)
|
|---|
| 233 | W(I) = WD(I)
|
|---|
| 234 | 100 CONTINUE
|
|---|
| 235 | RETURN
|
|---|
| 236 | END
|
|---|
| 237 | C
|
|---|
| 238 | SUBROUTINE ZWGLJD (Z,W,NP,ALPHA,BETA)
|
|---|
| 239 | C--------------------------------------------------------------------
|
|---|
| 240 | C
|
|---|
| 241 | C Generate NP GAUSS LOBATTO JACOBI points (Z) and weights (W)
|
|---|
| 242 | C associated with Jacobi polynomial P(N)(alpha>-1,beta>-1).
|
|---|
| 243 | C The polynomial degree N=NP-1.
|
|---|
| 244 | C Double precision version.
|
|---|
| 245 | C
|
|---|
| 246 | C--------------------------------------------------------------------
|
|---|
| 247 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 248 | REAL*8 Z(NP),W(NP),ALPHA,BETA
|
|---|
| 249 | C
|
|---|
| 250 | N = NP-1
|
|---|
| 251 | NM1 = N-1
|
|---|
| 252 | ONE = 1.
|
|---|
| 253 | TWO = 2.
|
|---|
| 254 | C
|
|---|
| 255 | IF (NP.LE.1) THEN
|
|---|
| 256 | WRITE (6,*) 'ZWGLJD: Minimum number of Gauss-Lobatto points is 2'
|
|---|
| 257 | WRITE (6,*) 'ZWGLJD: alpha,beta:',alpha,beta,np
|
|---|
| 258 | call exitt
|
|---|
| 259 | ENDIF
|
|---|
| 260 | IF ((ALPHA.LE.-ONE).OR.(BETA.LE.-ONE)) THEN
|
|---|
| 261 | WRITE (6,*) 'ZWGLJD: Alpha and Beta must be greater than -1'
|
|---|
| 262 | call exitt
|
|---|
| 263 | ENDIF
|
|---|
| 264 | C
|
|---|
| 265 | IF (NM1.GT.0) THEN
|
|---|
| 266 | ALPG = ALPHA+ONE
|
|---|
| 267 | BETG = BETA+ONE
|
|---|
| 268 | CALL ZWGJD (Z(2),W(2),NM1,ALPG,BETG)
|
|---|
| 269 | ENDIF
|
|---|
| 270 | Z(1) = -ONE
|
|---|
| 271 | Z(NP) = ONE
|
|---|
| 272 | DO 100 I=2,NP-1
|
|---|
| 273 | W(I) = W(I)/(ONE-Z(I)**2)
|
|---|
| 274 | 100 CONTINUE
|
|---|
| 275 | CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z(1))
|
|---|
| 276 | W(1) = ENDW1 (N,ALPHA,BETA)/(TWO*PD)
|
|---|
| 277 | CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z(NP))
|
|---|
| 278 | W(NP) = ENDW2 (N,ALPHA,BETA)/(TWO*PD)
|
|---|
| 279 | C
|
|---|
| 280 | RETURN
|
|---|
| 281 | END
|
|---|
| 282 | C
|
|---|
| 283 | REAL*8 FUNCTION ENDW1 (N,ALPHA,BETA)
|
|---|
| 284 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 285 | REAL*8 ALPHA,BETA
|
|---|
| 286 | ZERO = 0.
|
|---|
| 287 | ONE = 1.
|
|---|
| 288 | TWO = 2.
|
|---|
| 289 | THREE = 3.
|
|---|
| 290 | FOUR = 4.
|
|---|
| 291 | APB = ALPHA+BETA
|
|---|
| 292 | IF (N.EQ.0) THEN
|
|---|
| 293 | ENDW1 = ZERO
|
|---|
| 294 | RETURN
|
|---|
| 295 | ENDIF
|
|---|
| 296 | F1 = GAMMAF(ALPHA+TWO)*GAMMAF(BETA+ONE)/GAMMAF(APB+THREE)
|
|---|
| 297 | F1 = F1*(APB+TWO)*TWO**(APB+TWO)/TWO
|
|---|
| 298 | IF (N.EQ.1) THEN
|
|---|
| 299 | ENDW1 = F1
|
|---|
| 300 | RETURN
|
|---|
| 301 | ENDIF
|
|---|
| 302 | FINT1 = GAMMAF(ALPHA+TWO)*GAMMAF(BETA+ONE)/GAMMAF(APB+THREE)
|
|---|
| 303 | FINT1 = FINT1*TWO**(APB+TWO)
|
|---|
| 304 | FINT2 = GAMMAF(ALPHA+TWO)*GAMMAF(BETA+TWO)/GAMMAF(APB+FOUR)
|
|---|
| 305 | FINT2 = FINT2*TWO**(APB+THREE)
|
|---|
| 306 | F2 = (-TWO*(BETA+TWO)*FINT1 + (APB+FOUR)*FINT2)
|
|---|
| 307 | $ * (APB+THREE)/FOUR
|
|---|
| 308 | IF (N.EQ.2) THEN
|
|---|
| 309 | ENDW1 = F2
|
|---|
| 310 | RETURN
|
|---|
| 311 | ENDIF
|
|---|
| 312 | DO 100 I=3,N
|
|---|
| 313 | DI = ((I-1))
|
|---|
| 314 | ABN = ALPHA+BETA+DI
|
|---|
| 315 | ABNN = ABN+DI
|
|---|
| 316 | A1 = -(TWO*(DI+ALPHA)*(DI+BETA))/(ABN*ABNN*(ABNN+ONE))
|
|---|
| 317 | A2 = (TWO*(ALPHA-BETA))/(ABNN*(ABNN+TWO))
|
|---|
| 318 | A3 = (TWO*(ABN+ONE))/((ABNN+TWO)*(ABNN+ONE))
|
|---|
| 319 | F3 = -(A2*F2+A1*F1)/A3
|
|---|
| 320 | F1 = F2
|
|---|
| 321 | F2 = F3
|
|---|
| 322 | 100 CONTINUE
|
|---|
| 323 | ENDW1 = F3
|
|---|
| 324 | RETURN
|
|---|
| 325 | END
|
|---|
| 326 | C
|
|---|
| 327 | REAL*8 FUNCTION ENDW2 (N,ALPHA,BETA)
|
|---|
| 328 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 329 | REAL*8 ALPHA,BETA
|
|---|
| 330 | ZERO = 0.
|
|---|
| 331 | ONE = 1.
|
|---|
| 332 | TWO = 2.
|
|---|
| 333 | THREE = 3.
|
|---|
| 334 | FOUR = 4.
|
|---|
| 335 | APB = ALPHA+BETA
|
|---|
| 336 | IF (N.EQ.0) THEN
|
|---|
| 337 | ENDW2 = ZERO
|
|---|
| 338 | RETURN
|
|---|
| 339 | ENDIF
|
|---|
| 340 | F1 = GAMMAF(ALPHA+ONE)*GAMMAF(BETA+TWO)/GAMMAF(APB+THREE)
|
|---|
| 341 | F1 = F1*(APB+TWO)*TWO**(APB+TWO)/TWO
|
|---|
| 342 | IF (N.EQ.1) THEN
|
|---|
| 343 | ENDW2 = F1
|
|---|
| 344 | RETURN
|
|---|
| 345 | ENDIF
|
|---|
| 346 | FINT1 = GAMMAF(ALPHA+ONE)*GAMMAF(BETA+TWO)/GAMMAF(APB+THREE)
|
|---|
| 347 | FINT1 = FINT1*TWO**(APB+TWO)
|
|---|
| 348 | FINT2 = GAMMAF(ALPHA+TWO)*GAMMAF(BETA+TWO)/GAMMAF(APB+FOUR)
|
|---|
| 349 | FINT2 = FINT2*TWO**(APB+THREE)
|
|---|
| 350 | F2 = (TWO*(ALPHA+TWO)*FINT1 - (APB+FOUR)*FINT2)
|
|---|
| 351 | $ * (APB+THREE)/FOUR
|
|---|
| 352 | IF (N.EQ.2) THEN
|
|---|
| 353 | ENDW2 = F2
|
|---|
| 354 | RETURN
|
|---|
| 355 | ENDIF
|
|---|
| 356 | DO 100 I=3,N
|
|---|
| 357 | DI = ((I-1))
|
|---|
| 358 | ABN = ALPHA+BETA+DI
|
|---|
| 359 | ABNN = ABN+DI
|
|---|
| 360 | A1 = -(TWO*(DI+ALPHA)*(DI+BETA))/(ABN*ABNN*(ABNN+ONE))
|
|---|
| 361 | A2 = (TWO*(ALPHA-BETA))/(ABNN*(ABNN+TWO))
|
|---|
| 362 | A3 = (TWO*(ABN+ONE))/((ABNN+TWO)*(ABNN+ONE))
|
|---|
| 363 | F3 = -(A2*F2+A1*F1)/A3
|
|---|
| 364 | F1 = F2
|
|---|
| 365 | F2 = F3
|
|---|
| 366 | 100 CONTINUE
|
|---|
| 367 | ENDW2 = F3
|
|---|
| 368 | RETURN
|
|---|
| 369 | END
|
|---|
| 370 | C
|
|---|
| 371 | REAL*8 FUNCTION GAMMAF (X)
|
|---|
| 372 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 373 | REAL*8 X
|
|---|
| 374 | ZERO = 0.0
|
|---|
| 375 | HALF = 0.5
|
|---|
| 376 | ONE = 1.0
|
|---|
| 377 | TWO = 2.0
|
|---|
| 378 | FOUR = 4.0
|
|---|
| 379 | PI = FOUR*ATAN(ONE)
|
|---|
| 380 | GAMMAF = ONE
|
|---|
| 381 | IF (X.EQ.-HALF) GAMMAF = -TWO*SQRT(PI)
|
|---|
| 382 | IF (X.EQ. HALF) GAMMAF = SQRT(PI)
|
|---|
| 383 | IF (X.EQ. ONE ) GAMMAF = ONE
|
|---|
| 384 | IF (X.EQ. TWO ) GAMMAF = ONE
|
|---|
| 385 | IF (X.EQ. 1.5 ) GAMMAF = SQRT(PI)/2.
|
|---|
| 386 | IF (X.EQ. 2.5) GAMMAF = 1.5*SQRT(PI)/2.
|
|---|
| 387 | IF (X.EQ. 3.5) GAMMAF = 0.5*(2.5*(1.5*SQRT(PI)))
|
|---|
| 388 | IF (X.EQ. 3. ) GAMMAF = 2.
|
|---|
| 389 | IF (X.EQ. 4. ) GAMMAF = 6.
|
|---|
| 390 | IF (X.EQ. 5. ) GAMMAF = 24.
|
|---|
| 391 | IF (X.EQ. 6. ) GAMMAF = 120.
|
|---|
| 392 | RETURN
|
|---|
| 393 | END
|
|---|
| 394 | C
|
|---|
| 395 | REAL*8 FUNCTION PNORMJ (N,ALPHA,BETA)
|
|---|
| 396 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 397 | REAL*8 ALPHA,BETA
|
|---|
| 398 | ONE = 1.
|
|---|
| 399 | TWO = 2.
|
|---|
| 400 | DN = ((N))
|
|---|
| 401 | CONST = ALPHA+BETA+ONE
|
|---|
| 402 | IF (N.LE.1) THEN
|
|---|
| 403 | PROD = GAMMAF(DN+ALPHA)*GAMMAF(DN+BETA)
|
|---|
| 404 | PROD = PROD/(GAMMAF(DN)*GAMMAF(DN+ALPHA+BETA))
|
|---|
| 405 | PNORMJ = PROD * TWO**CONST/(TWO*DN+CONST)
|
|---|
| 406 | RETURN
|
|---|
| 407 | ENDIF
|
|---|
| 408 | PROD = GAMMAF(ALPHA+ONE)*GAMMAF(BETA+ONE)
|
|---|
| 409 | PROD = PROD/(TWO*(ONE+CONST)*GAMMAF(CONST+ONE))
|
|---|
| 410 | PROD = PROD*(ONE+ALPHA)*(TWO+ALPHA)
|
|---|
| 411 | PROD = PROD*(ONE+BETA)*(TWO+BETA)
|
|---|
| 412 | DO 100 I=3,N
|
|---|
| 413 | DINDX = ((I))
|
|---|
| 414 | FRAC = (DINDX+ALPHA)*(DINDX+BETA)/(DINDX*(DINDX+ALPHA+BETA))
|
|---|
| 415 | PROD = PROD*FRAC
|
|---|
| 416 | 100 CONTINUE
|
|---|
| 417 | PNORMJ = PROD * TWO**CONST/(TWO*DN+CONST)
|
|---|
| 418 | RETURN
|
|---|
| 419 | END
|
|---|
| 420 | C
|
|---|
| 421 | SUBROUTINE JACG (XJAC,NP,ALPHA,BETA)
|
|---|
| 422 | C--------------------------------------------------------------------
|
|---|
| 423 | C
|
|---|
| 424 | C Compute NP Gauss points XJAC, which are the zeros of the
|
|---|
| 425 | C Jacobi polynomial J(NP) with parameters ALPHA and BETA.
|
|---|
| 426 | C ALPHA and BETA determines the specific type of Gauss points.
|
|---|
| 427 | C Examples:
|
|---|
| 428 | C ALPHA = BETA = 0.0 -> Legendre points
|
|---|
| 429 | C ALPHA = BETA = -0.5 -> Chebyshev points
|
|---|
| 430 | C
|
|---|
| 431 | C--------------------------------------------------------------------
|
|---|
| 432 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 433 | REAL*8 XJAC(1)
|
|---|
| 434 | DATA KSTOP /10/
|
|---|
| 435 | DATA EPS/1.0e-12/
|
|---|
| 436 | N = NP-1
|
|---|
| 437 | one = 1.
|
|---|
| 438 | DTH = 4.*ATAN(one)/(2.*((N))+2.)
|
|---|
| 439 | DO 40 J=1,NP
|
|---|
| 440 | IF (J.EQ.1) THEN
|
|---|
| 441 | X = COS((2.*(((J))-1.)+1.)*DTH)
|
|---|
| 442 | ELSE
|
|---|
| 443 | X1 = COS((2.*(((J))-1.)+1.)*DTH)
|
|---|
| 444 | X2 = XLAST
|
|---|
| 445 | X = (X1+X2)/2.
|
|---|
| 446 | ENDIF
|
|---|
| 447 | DO 30 K=1,KSTOP
|
|---|
| 448 | CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,NP,ALPHA,BETA,X)
|
|---|
| 449 | RECSUM = 0.
|
|---|
| 450 | JM = J-1
|
|---|
| 451 | DO 29 I=1,JM
|
|---|
| 452 | RECSUM = RECSUM+1./(X-XJAC(NP-I+1))
|
|---|
| 453 | 29 CONTINUE
|
|---|
| 454 | DELX = -P/(PD-RECSUM*P)
|
|---|
| 455 | X = X+DELX
|
|---|
| 456 | IF (ABS(DELX) .LT. EPS) GOTO 31
|
|---|
| 457 | 30 CONTINUE
|
|---|
| 458 | 31 CONTINUE
|
|---|
| 459 | XJAC(NP-J+1) = X
|
|---|
| 460 | XLAST = X
|
|---|
| 461 | 40 CONTINUE
|
|---|
| 462 | DO 200 I=1,NP
|
|---|
| 463 | XMIN = 2.
|
|---|
| 464 | DO 100 J=I,NP
|
|---|
| 465 | IF (XJAC(J).LT.XMIN) THEN
|
|---|
| 466 | XMIN = XJAC(J)
|
|---|
| 467 | JMIN = J
|
|---|
| 468 | ENDIF
|
|---|
| 469 | 100 CONTINUE
|
|---|
| 470 | IF (JMIN.NE.I) THEN
|
|---|
| 471 | SWAP = XJAC(I)
|
|---|
| 472 | XJAC(I) = XJAC(JMIN)
|
|---|
| 473 | XJAC(JMIN) = SWAP
|
|---|
| 474 | ENDIF
|
|---|
| 475 | 200 CONTINUE
|
|---|
| 476 | RETURN
|
|---|
| 477 | END
|
|---|
| 478 | C
|
|---|
| 479 | SUBROUTINE JACOBF (POLY,PDER,POLYM1,PDERM1,POLYM2,PDERM2,
|
|---|
| 480 | $ N,ALP,BET,X)
|
|---|
| 481 | C--------------------------------------------------------------------
|
|---|
| 482 | C
|
|---|
| 483 | C Computes the Jacobi polynomial (POLY) and its derivative (PDER)
|
|---|
| 484 | C of degree N at X.
|
|---|
| 485 | C
|
|---|
| 486 | C--------------------------------------------------------------------
|
|---|
| 487 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 488 | APB = ALP+BET
|
|---|
| 489 | POLY = 1.
|
|---|
| 490 | PDER = 0.
|
|---|
| 491 | IF (N .EQ. 0) RETURN
|
|---|
| 492 | POLYL = POLY
|
|---|
| 493 | PDERL = PDER
|
|---|
| 494 | POLY = (ALP-BET+(APB+2.)*X)/2.
|
|---|
| 495 | PDER = (APB+2.)/2.
|
|---|
| 496 | IF (N .EQ. 1) RETURN
|
|---|
| 497 | DO 20 K=2,N
|
|---|
| 498 | DK = ((K))
|
|---|
| 499 | A1 = 2.*DK*(DK+APB)*(2.*DK+APB-2.)
|
|---|
| 500 | A2 = (2.*DK+APB-1.)*(ALP**2-BET**2)
|
|---|
| 501 | B3 = (2.*DK+APB-2.)
|
|---|
| 502 | A3 = B3*(B3+1.)*(B3+2.)
|
|---|
| 503 | A4 = 2.*(DK+ALP-1.)*(DK+BET-1.)*(2.*DK+APB)
|
|---|
| 504 | POLYN = ((A2+A3*X)*POLY-A4*POLYL)/A1
|
|---|
| 505 | PDERN = ((A2+A3*X)*PDER-A4*PDERL+A3*POLY)/A1
|
|---|
| 506 | PSAVE = POLYL
|
|---|
| 507 | PDSAVE = PDERL
|
|---|
| 508 | POLYL = POLY
|
|---|
| 509 | POLY = POLYN
|
|---|
| 510 | PDERL = PDER
|
|---|
| 511 | PDER = PDERN
|
|---|
| 512 | 20 CONTINUE
|
|---|
| 513 | POLYM1 = POLYL
|
|---|
| 514 | PDERM1 = PDERL
|
|---|
| 515 | POLYM2 = PSAVE
|
|---|
| 516 | PDERM2 = PDSAVE
|
|---|
| 517 | RETURN
|
|---|
| 518 | END
|
|---|
| 519 | C
|
|---|
| 520 | REAL FUNCTION HGJ (II,Z,ZGJ,NP,ALPHA,BETA)
|
|---|
| 521 | C---------------------------------------------------------------------
|
|---|
| 522 | C
|
|---|
| 523 | C Compute the value of the Lagrangian interpolant HGJ through
|
|---|
| 524 | C the NP Gauss Jacobi points ZGJ at the point Z.
|
|---|
| 525 | C Single precision version.
|
|---|
| 526 | C
|
|---|
| 527 | C---------------------------------------------------------------------
|
|---|
| 528 | PARAMETER (NMAX=84)
|
|---|
| 529 | PARAMETER (lzd = NMAX)
|
|---|
| 530 | REAL*8 ZD,ZGJD(lzd),HGJD,ALPHAD,BETAD
|
|---|
| 531 | REAL Z,ZGJ(1),ALPHA,BETA
|
|---|
| 532 | NPMAX = lzd
|
|---|
| 533 | IF (NP.GT.NPMAX) THEN
|
|---|
| 534 | WRITE (6,*) 'Too large polynomial degree in HGJ'
|
|---|
| 535 | WRITE (6,*) 'Maximum polynomial degree is',NMAX
|
|---|
| 536 | WRITE (6,*) 'Here NP=',NP
|
|---|
| 537 | call exitt
|
|---|
| 538 | ENDIF
|
|---|
| 539 | ZD = Z
|
|---|
| 540 | DO 100 I=1,NP
|
|---|
| 541 | ZGJD(I) = ZGJ(I)
|
|---|
| 542 | 100 CONTINUE
|
|---|
| 543 | ALPHAD = ALPHA
|
|---|
| 544 | BETAD = BETA
|
|---|
| 545 | HGJ = HGJD (II,ZD,ZGJD,NP,ALPHAD,BETAD)
|
|---|
| 546 | RETURN
|
|---|
| 547 | END
|
|---|
| 548 | C
|
|---|
| 549 | REAL*8 FUNCTION HGJD (II,Z,ZGJ,NP,ALPHA,BETA)
|
|---|
| 550 | C---------------------------------------------------------------------
|
|---|
| 551 | C
|
|---|
| 552 | C Compute the value of the Lagrangian interpolant HGJD through
|
|---|
| 553 | C the NZ Gauss-Lobatto Jacobi points ZGJ at the point Z.
|
|---|
| 554 | C Double precision version.
|
|---|
| 555 | C
|
|---|
| 556 | C---------------------------------------------------------------------
|
|---|
| 557 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 558 | REAL*8 Z,ZGJ(1),ALPHA,BETA
|
|---|
| 559 | EPS = 1.e-5
|
|---|
| 560 | ONE = 1.
|
|---|
| 561 | ZI = ZGJ(II)
|
|---|
| 562 | DZ = Z-ZI
|
|---|
| 563 | IF (ABS(DZ).LT.EPS) THEN
|
|---|
| 564 | HGJD = ONE
|
|---|
| 565 | RETURN
|
|---|
| 566 | ENDIF
|
|---|
| 567 | CALL JACOBF (PZI,PDZI,PM1,PDM1,PM2,PDM2,NP,ALPHA,BETA,ZI)
|
|---|
| 568 | CALL JACOBF (PZ,PDZ,PM1,PDM1,PM2,PDM2,NP,ALPHA,BETA,Z)
|
|---|
| 569 | HGJD = PZ/(PDZI*(Z-ZI))
|
|---|
| 570 | RETURN
|
|---|
| 571 | END
|
|---|
| 572 | C
|
|---|
| 573 | REAL FUNCTION HGLJ (II,Z,ZGLJ,NP,ALPHA,BETA)
|
|---|
| 574 | C---------------------------------------------------------------------
|
|---|
| 575 | C
|
|---|
| 576 | C Compute the value of the Lagrangian interpolant HGLJ through
|
|---|
| 577 | C the NZ Gauss-Lobatto Jacobi points ZGLJ at the point Z.
|
|---|
| 578 | C Single precision version.
|
|---|
| 579 | C
|
|---|
| 580 | C---------------------------------------------------------------------
|
|---|
| 581 | PARAMETER (NMAX=84)
|
|---|
| 582 | PARAMETER (lzd = NMAX)
|
|---|
| 583 | REAL*8 ZD,ZGLJD(lzd),HGLJD,ALPHAD,BETAD
|
|---|
| 584 | REAL Z,ZGLJ(1),ALPHA,BETA
|
|---|
| 585 | NPMAX = lzd
|
|---|
| 586 | IF (NP.GT.NPMAX) THEN
|
|---|
| 587 | WRITE (6,*) 'Too large polynomial degree in HGLJ'
|
|---|
| 588 | WRITE (6,*) 'Maximum polynomial degree is',NMAX
|
|---|
| 589 | WRITE (6,*) 'Here NP=',NP
|
|---|
| 590 | call exitt
|
|---|
| 591 | ENDIF
|
|---|
| 592 | ZD = Z
|
|---|
| 593 | DO 100 I=1,NP
|
|---|
| 594 | ZGLJD(I) = ZGLJ(I)
|
|---|
| 595 | 100 CONTINUE
|
|---|
| 596 | ALPHAD = ALPHA
|
|---|
| 597 | BETAD = BETA
|
|---|
| 598 | HGLJ = HGLJD (II,ZD,ZGLJD,NP,ALPHAD,BETAD)
|
|---|
| 599 | RETURN
|
|---|
| 600 | END
|
|---|
| 601 | C
|
|---|
| 602 | REAL*8 FUNCTION HGLJD (I,Z,ZGLJ,NP,ALPHA,BETA)
|
|---|
| 603 | C---------------------------------------------------------------------
|
|---|
| 604 | C
|
|---|
| 605 | C Compute the value of the Lagrangian interpolant HGLJD through
|
|---|
| 606 | C the NZ Gauss-Lobatto Jacobi points ZJACL at the point Z.
|
|---|
| 607 | C Double precision version.
|
|---|
| 608 | C
|
|---|
| 609 | C---------------------------------------------------------------------
|
|---|
| 610 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 611 | REAL*8 Z,ZGLJ(1),ALPHA,BETA
|
|---|
| 612 | EPS = 1.e-5
|
|---|
| 613 | ONE = 1.
|
|---|
| 614 | ZI = ZGLJ(I)
|
|---|
| 615 | DZ = Z-ZI
|
|---|
| 616 | IF (ABS(DZ).LT.EPS) THEN
|
|---|
| 617 | HGLJD = ONE
|
|---|
| 618 | RETURN
|
|---|
| 619 | ENDIF
|
|---|
| 620 | N = NP-1
|
|---|
| 621 | DN = ((N))
|
|---|
| 622 | EIGVAL = -DN*(DN+ALPHA+BETA+ONE)
|
|---|
| 623 | CALL JACOBF (PI,PDI,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,ZI)
|
|---|
| 624 | CONST = EIGVAL*PI+ALPHA*(ONE+ZI)*PDI-BETA*(ONE-ZI)*PDI
|
|---|
| 625 | CALL JACOBF (P,PD,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z)
|
|---|
| 626 | HGLJD = (ONE-Z**2)*PD/(CONST*(Z-ZI))
|
|---|
| 627 | RETURN
|
|---|
| 628 | END
|
|---|
| 629 | C
|
|---|
| 630 | SUBROUTINE DGJ (D,DT,Z,NZ,lzd,ALPHA,BETA)
|
|---|
| 631 | C-----------------------------------------------------------------
|
|---|
| 632 | C
|
|---|
| 633 | C Compute the derivative matrix D and its transpose DT
|
|---|
| 634 | C associated with the Nth order Lagrangian interpolants
|
|---|
| 635 | C through the NZ Gauss Jacobi points Z.
|
|---|
| 636 | C Note: D and DT are square matrices.
|
|---|
| 637 | C Single precision version.
|
|---|
| 638 | C
|
|---|
| 639 | C-----------------------------------------------------------------
|
|---|
| 640 | PARAMETER (NMAX=84)
|
|---|
| 641 | PARAMETER (lzdD = NMAX)
|
|---|
| 642 | REAL*8 DD(lzdD,lzdD),DTD(lzdD,lzdD),ZD(lzdD),ALPHAD,BETAD
|
|---|
| 643 | REAL D(lzd,lzd),DT(lzd,lzd),Z(1),ALPHA,BETA
|
|---|
| 644 | C
|
|---|
| 645 | IF (NZ.LE.0) THEN
|
|---|
| 646 | WRITE (6,*) 'DGJ: Minimum number of Gauss points is 1'
|
|---|
| 647 | call exitt
|
|---|
| 648 | ENDIF
|
|---|
| 649 | IF (NZ .GT. NMAX) THEN
|
|---|
| 650 | WRITE (6,*) 'Too large polynomial degree in DGJ'
|
|---|
| 651 | WRITE (6,*) 'Maximum polynomial degree is',NMAX
|
|---|
| 652 | WRITE (6,*) 'Here Nz=',Nz
|
|---|
| 653 | call exitt
|
|---|
| 654 | ENDIF
|
|---|
| 655 | IF ((ALPHA.LE.-1.).OR.(BETA.LE.-1.)) THEN
|
|---|
| 656 | WRITE (6,*) 'DGJ: Alpha and Beta must be greater than -1'
|
|---|
| 657 | call exitt
|
|---|
| 658 | ENDIF
|
|---|
| 659 | ALPHAD = ALPHA
|
|---|
| 660 | BETAD = BETA
|
|---|
| 661 | DO 100 I=1,NZ
|
|---|
| 662 | ZD(I) = Z(I)
|
|---|
| 663 | 100 CONTINUE
|
|---|
| 664 | CALL DGJD (DD,DTD,ZD,NZ,lzdD,ALPHAD,BETAD)
|
|---|
| 665 | DO 200 I=1,NZ
|
|---|
| 666 | DO 200 J=1,NZ
|
|---|
| 667 | D(I,J) = DD(I,J)
|
|---|
| 668 | DT(I,J) = DTD(I,J)
|
|---|
| 669 | 200 CONTINUE
|
|---|
| 670 | RETURN
|
|---|
| 671 | END
|
|---|
| 672 | C
|
|---|
| 673 | SUBROUTINE DGJD (D,DT,Z,NZ,lzd,ALPHA,BETA)
|
|---|
| 674 | C-----------------------------------------------------------------
|
|---|
| 675 | C
|
|---|
| 676 | C Compute the derivative matrix D and its transpose DT
|
|---|
| 677 | C associated with the Nth order Lagrangian interpolants
|
|---|
| 678 | C through the NZ Gauss Jacobi points Z.
|
|---|
| 679 | C Note: D and DT are square matrices.
|
|---|
| 680 | C Double precision version.
|
|---|
| 681 | C
|
|---|
| 682 | C-----------------------------------------------------------------
|
|---|
| 683 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 684 | REAL*8 D(lzd,lzd),DT(lzd,lzd),Z(1),ALPHA,BETA
|
|---|
| 685 | N = NZ-1
|
|---|
| 686 | DN = ((N))
|
|---|
| 687 | ONE = 1.
|
|---|
| 688 | TWO = 2.
|
|---|
| 689 | C
|
|---|
| 690 | IF (NZ.LE.1) THEN
|
|---|
| 691 | WRITE (6,*) 'DGJD: Minimum number of Gauss-Lobatto points is 2'
|
|---|
| 692 | call exitt
|
|---|
| 693 | ENDIF
|
|---|
| 694 | IF ((ALPHA.LE.-ONE).OR.(BETA.LE.-ONE)) THEN
|
|---|
| 695 | WRITE (6,*) 'DGJD: Alpha and Beta must be greater than -1'
|
|---|
| 696 | call exitt
|
|---|
| 697 | ENDIF
|
|---|
| 698 | C
|
|---|
| 699 | DO 200 I=1,NZ
|
|---|
| 700 | DO 200 J=1,NZ
|
|---|
| 701 | CALL JACOBF (PI,PDI,PM1,PDM1,PM2,PDM2,NZ,ALPHA,BETA,Z(I))
|
|---|
| 702 | CALL JACOBF (PJ,PDJ,PM1,PDM1,PM2,PDM2,NZ,ALPHA,BETA,Z(J))
|
|---|
| 703 | IF (I.NE.J) D(I,J) = PDI/(PDJ*(Z(I)-Z(J)))
|
|---|
| 704 | IF (I.EQ.J) D(I,J) = ((ALPHA+BETA+TWO)*Z(I)+ALPHA-BETA)/
|
|---|
| 705 | $ (TWO*(ONE-Z(I)**2))
|
|---|
| 706 | DT(J,I) = D(I,J)
|
|---|
| 707 | 200 CONTINUE
|
|---|
| 708 | RETURN
|
|---|
| 709 | END
|
|---|
| 710 | C
|
|---|
| 711 | SUBROUTINE DGLJ (D,DT,Z,NZ,lzd,ALPHA,BETA)
|
|---|
| 712 | C-----------------------------------------------------------------
|
|---|
| 713 | C
|
|---|
| 714 | C Compute the derivative matrix D and its transpose DT
|
|---|
| 715 | C associated with the Nth order Lagrangian interpolants
|
|---|
| 716 | C through the NZ Gauss-Lobatto Jacobi points Z.
|
|---|
| 717 | C Note: D and DT are square matrices.
|
|---|
| 718 | C Single precision version.
|
|---|
| 719 | C
|
|---|
| 720 | C-----------------------------------------------------------------
|
|---|
| 721 | PARAMETER (NMAX=84)
|
|---|
| 722 | PARAMETER (lzdD = NMAX)
|
|---|
| 723 | REAL*8 DD(lzdD,lzdD),DTD(lzdD,lzdD),ZD(lzdD),ALPHAD,BETAD
|
|---|
| 724 | REAL D(lzd,lzd),DT(lzd,lzd),Z(1),ALPHA,BETA
|
|---|
| 725 | C
|
|---|
| 726 | IF (NZ.LE.1) THEN
|
|---|
| 727 | WRITE (6,*) 'DGLJ: Minimum number of Gauss-Lobatto points is 2'
|
|---|
| 728 | call exitt
|
|---|
| 729 | ENDIF
|
|---|
| 730 | IF (NZ .GT. NMAX) THEN
|
|---|
| 731 | WRITE (6,*) 'Too large polynomial degree in DGLJ'
|
|---|
| 732 | WRITE (6,*) 'Maximum polynomial degree is',NMAX
|
|---|
| 733 | WRITE (6,*) 'Here NZ=',NZ
|
|---|
| 734 | call exitt
|
|---|
| 735 | ENDIF
|
|---|
| 736 | IF ((ALPHA.LE.-1.).OR.(BETA.LE.-1.)) THEN
|
|---|
| 737 | WRITE (6,*) 'DGLJ: Alpha and Beta must be greater than -1'
|
|---|
| 738 | call exitt
|
|---|
| 739 | ENDIF
|
|---|
| 740 | ALPHAD = ALPHA
|
|---|
| 741 | BETAD = BETA
|
|---|
| 742 | DO 100 I=1,NZ
|
|---|
| 743 | ZD(I) = Z(I)
|
|---|
| 744 | 100 CONTINUE
|
|---|
| 745 | CALL DGLJD (DD,DTD,ZD,NZ,lzdD,ALPHAD,BETAD)
|
|---|
| 746 | DO 200 I=1,NZ
|
|---|
| 747 | DO 200 J=1,NZ
|
|---|
| 748 | D(I,J) = DD(I,J)
|
|---|
| 749 | DT(I,J) = DTD(I,J)
|
|---|
| 750 | 200 CONTINUE
|
|---|
| 751 | RETURN
|
|---|
| 752 | END
|
|---|
| 753 | C
|
|---|
| 754 | SUBROUTINE DGLJD (D,DT,Z,NZ,lzd,ALPHA,BETA)
|
|---|
| 755 | C-----------------------------------------------------------------
|
|---|
| 756 | C
|
|---|
| 757 | C Compute the derivative matrix D and its transpose DT
|
|---|
| 758 | C associated with the Nth order Lagrangian interpolants
|
|---|
| 759 | C through the NZ Gauss-Lobatto Jacobi points Z.
|
|---|
| 760 | C Note: D and DT are square matrices.
|
|---|
| 761 | C Double precision version.
|
|---|
| 762 | C
|
|---|
| 763 | C-----------------------------------------------------------------
|
|---|
| 764 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 765 | REAL*8 D(lzd,lzd),DT(lzd,lzd),Z(1),ALPHA,BETA
|
|---|
| 766 | N = NZ-1
|
|---|
| 767 | DN = ((N))
|
|---|
| 768 | ONE = 1.
|
|---|
| 769 | TWO = 2.
|
|---|
| 770 | EIGVAL = -DN*(DN+ALPHA+BETA+ONE)
|
|---|
| 771 | C
|
|---|
| 772 | IF (NZ.LE.1) THEN
|
|---|
| 773 | WRITE (6,*) 'DGLJD: Minimum number of Gauss-Lobatto points is 2'
|
|---|
| 774 | call exitt
|
|---|
| 775 | ENDIF
|
|---|
| 776 | IF ((ALPHA.LE.-ONE).OR.(BETA.LE.-ONE)) THEN
|
|---|
| 777 | WRITE (6,*) 'DGLJD: Alpha and Beta must be greater than -1'
|
|---|
| 778 | call exitt
|
|---|
| 779 | ENDIF
|
|---|
| 780 | C
|
|---|
| 781 | DO 200 I=1,NZ
|
|---|
| 782 | DO 200 J=1,NZ
|
|---|
| 783 | CALL JACOBF (PI,PDI,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z(I))
|
|---|
| 784 | CALL JACOBF (PJ,PDJ,PM1,PDM1,PM2,PDM2,N,ALPHA,BETA,Z(J))
|
|---|
| 785 | CI = EIGVAL*PI-(BETA*(ONE-Z(I))-ALPHA*(ONE+Z(I)))*PDI
|
|---|
| 786 | CJ = EIGVAL*PJ-(BETA*(ONE-Z(J))-ALPHA*(ONE+Z(J)))*PDJ
|
|---|
| 787 | IF (I.NE.J) D(I,J) = CI/(CJ*(Z(I)-Z(J)))
|
|---|
| 788 | IF ((I.EQ.J).AND.(I.NE.1).AND.(I.NE.NZ))
|
|---|
| 789 | $ D(I,J) = (ALPHA*(ONE+Z(I))-BETA*(ONE-Z(I)))/
|
|---|
| 790 | $ (TWO*(ONE-Z(I)**2))
|
|---|
| 791 | IF ((I.EQ.J).AND.(I.EQ.1))
|
|---|
| 792 | $ D(I,J) = (EIGVAL+ALPHA)/(TWO*(BETA+TWO))
|
|---|
| 793 | IF ((I.EQ.J).AND.(I.EQ.NZ))
|
|---|
| 794 | $ D(I,J) = -(EIGVAL+BETA)/(TWO*(ALPHA+TWO))
|
|---|
| 795 | DT(J,I) = D(I,J)
|
|---|
| 796 | 200 CONTINUE
|
|---|
| 797 | RETURN
|
|---|
| 798 | END
|
|---|
| 799 | C
|
|---|
| 800 | SUBROUTINE DGLL (D,DT,Z,NZ,lzd)
|
|---|
| 801 | C-----------------------------------------------------------------
|
|---|
| 802 | C
|
|---|
| 803 | C Compute the derivative matrix D and its transpose DT
|
|---|
| 804 | C associated with the Nth order Lagrangian interpolants
|
|---|
| 805 | C through the NZ Gauss-Lobatto Legendre points Z.
|
|---|
| 806 | C Note: D and DT are square matrices.
|
|---|
| 807 | C
|
|---|
| 808 | C-----------------------------------------------------------------
|
|---|
| 809 | PARAMETER (NMAX=84)
|
|---|
| 810 | REAL D(lzd,lzd),DT(lzd,lzd),Z(1)
|
|---|
| 811 | N = NZ-1
|
|---|
| 812 | IF (NZ .GT. NMAX) THEN
|
|---|
| 813 | WRITE (6,*) 'Subroutine DGLL'
|
|---|
| 814 | WRITE (6,*) 'Maximum polynomial degree =',NMAX
|
|---|
| 815 | WRITE (6,*) 'Polynomial degree =',NZ
|
|---|
| 816 | ENDIF
|
|---|
| 817 | IF (NZ .EQ. 1) THEN
|
|---|
| 818 | D(1,1) = 0.
|
|---|
| 819 | RETURN
|
|---|
| 820 | ENDIF
|
|---|
| 821 | FN = (N)
|
|---|
| 822 | d0 = FN*(FN+1.)/4.
|
|---|
| 823 | DO 200 I=1,NZ
|
|---|
| 824 | DO 200 J=1,NZ
|
|---|
| 825 | D(I,J) = 0.
|
|---|
| 826 | IF (I.NE.J) D(I,J) = PNLEG(Z(I),N)/
|
|---|
| 827 | $ (PNLEG(Z(J),N)*(Z(I)-Z(J)))
|
|---|
| 828 | IF ((I.EQ.J).AND.(I.EQ.1)) D(I,J) = -d0
|
|---|
| 829 | IF ((I.EQ.J).AND.(I.EQ.NZ)) D(I,J) = d0
|
|---|
| 830 | DT(J,I) = D(I,J)
|
|---|
| 831 | 200 CONTINUE
|
|---|
| 832 | RETURN
|
|---|
| 833 | END
|
|---|
| 834 | C
|
|---|
| 835 | REAL FUNCTION HGLL (I,Z,ZGLL,NZ)
|
|---|
| 836 | C---------------------------------------------------------------------
|
|---|
| 837 | C
|
|---|
| 838 | C Compute the value of the Lagrangian interpolant L through
|
|---|
| 839 | C the NZ Gauss-Lobatto Legendre points ZGLL at the point Z.
|
|---|
| 840 | C
|
|---|
| 841 | C---------------------------------------------------------------------
|
|---|
| 842 | REAL ZGLL(1)
|
|---|
| 843 | EPS = 1.E-5
|
|---|
| 844 | DZ = Z - ZGLL(I)
|
|---|
| 845 | IF (ABS(DZ) .LT. EPS) THEN
|
|---|
| 846 | HGLL = 1.
|
|---|
| 847 | RETURN
|
|---|
| 848 | ENDIF
|
|---|
| 849 | N = NZ - 1
|
|---|
| 850 | ALFAN = (N)*((N)+1.)
|
|---|
| 851 | HGLL = - (1.-Z*Z)*PNDLEG(Z,N)/
|
|---|
| 852 | $ (ALFAN*PNLEG(ZGLL(I),N)*(Z-ZGLL(I)))
|
|---|
| 853 | RETURN
|
|---|
| 854 | END
|
|---|
| 855 | C
|
|---|
| 856 | REAL FUNCTION HGL (I,Z,ZGL,NZ)
|
|---|
| 857 | C---------------------------------------------------------------------
|
|---|
| 858 | C
|
|---|
| 859 | C Compute the value of the Lagrangian interpolant HGL through
|
|---|
| 860 | C the NZ Gauss Legendre points ZGL at the point Z.
|
|---|
| 861 | C
|
|---|
| 862 | C---------------------------------------------------------------------
|
|---|
| 863 | REAL ZGL(1)
|
|---|
| 864 | EPS = 1.E-5
|
|---|
| 865 | DZ = Z - ZGL(I)
|
|---|
| 866 | IF (ABS(DZ) .LT. EPS) THEN
|
|---|
| 867 | HGL = 1.
|
|---|
| 868 | RETURN
|
|---|
| 869 | ENDIF
|
|---|
| 870 | N = NZ-1
|
|---|
| 871 | HGL = PNLEG(Z,NZ)/(PNDLEG(ZGL(I),NZ)*(Z-ZGL(I)))
|
|---|
| 872 | RETURN
|
|---|
| 873 | END
|
|---|
| 874 | C
|
|---|
| 875 | REAL FUNCTION PNLEG (Z,N)
|
|---|
| 876 | C---------------------------------------------------------------------
|
|---|
| 877 | C
|
|---|
| 878 | C Compute the value of the Nth order Legendre polynomial at Z.
|
|---|
| 879 | C (Simpler than JACOBF)
|
|---|
| 880 | C Based on the recursion formula for the Legendre polynomials.
|
|---|
| 881 | C
|
|---|
| 882 | C---------------------------------------------------------------------
|
|---|
| 883 | C
|
|---|
| 884 | C This next statement is to overcome the underflow bug in the i860.
|
|---|
| 885 | C It can be removed at a later date. 11 Aug 1990 pff.
|
|---|
| 886 | C
|
|---|
| 887 | IF(ABS(Z) .LT. 1.0E-25) Z = 0.0
|
|---|
| 888 | C
|
|---|
| 889 | P1 = 1.
|
|---|
| 890 | IF (N.EQ.0) THEN
|
|---|
| 891 | PNLEG = P1
|
|---|
| 892 | RETURN
|
|---|
| 893 | ENDIF
|
|---|
| 894 | P2 = Z
|
|---|
| 895 | P3 = P2
|
|---|
| 896 | DO 10 K = 1, N-1
|
|---|
| 897 | FK = (K)
|
|---|
| 898 | P3 = ((2.*FK+1.)*Z*P2 - FK*P1)/(FK+1.)
|
|---|
| 899 | P1 = P2
|
|---|
| 900 | P2 = P3
|
|---|
| 901 | 10 CONTINUE
|
|---|
| 902 | PNLEG = P3
|
|---|
| 903 | if (n.eq.0) pnleg = 1.
|
|---|
| 904 | RETURN
|
|---|
| 905 | END
|
|---|
| 906 | C
|
|---|
| 907 | REAL FUNCTION PNDLEG (Z,N)
|
|---|
| 908 | C----------------------------------------------------------------------
|
|---|
| 909 | C
|
|---|
| 910 | C Compute the derivative of the Nth order Legendre polynomial at Z.
|
|---|
| 911 | C (Simpler than JACOBF)
|
|---|
| 912 | C Based on the recursion formula for the Legendre polynomials.
|
|---|
| 913 | C
|
|---|
| 914 | C----------------------------------------------------------------------
|
|---|
| 915 | P1 = 1.
|
|---|
| 916 | P2 = Z
|
|---|
| 917 | P1D = 0.
|
|---|
| 918 | P2D = 1.
|
|---|
| 919 | P3D = 1.
|
|---|
| 920 | DO 10 K = 1, N-1
|
|---|
| 921 | FK = (K)
|
|---|
| 922 | P3 = ((2.*FK+1.)*Z*P2 - FK*P1)/(FK+1.)
|
|---|
| 923 | P3D = ((2.*FK+1.)*P2 + (2.*FK+1.)*Z*P2D - FK*P1D)/(FK+1.)
|
|---|
| 924 | P1 = P2
|
|---|
| 925 | P2 = P3
|
|---|
| 926 | P1D = P2D
|
|---|
| 927 | P2D = P3D
|
|---|
| 928 | 10 CONTINUE
|
|---|
| 929 | PNDLEG = P3D
|
|---|
| 930 | IF (N.eq.0) pndleg = 0.
|
|---|
| 931 | RETURN
|
|---|
| 932 | END
|
|---|
| 933 | C
|
|---|
| 934 | SUBROUTINE DGLLGL (D,DT,ZM1,ZM2,IM12,NZM1,NZM2,ND1,ND2)
|
|---|
| 935 | C-----------------------------------------------------------------------
|
|---|
| 936 | C
|
|---|
| 937 | C Compute the (one-dimensional) derivative matrix D and its
|
|---|
| 938 | C transpose DT associated with taking the derivative of a variable
|
|---|
| 939 | C expanded on a Gauss-Lobatto Legendre mesh (M1), and evaluate its
|
|---|
| 940 | C derivative on a Guass Legendre mesh (M2).
|
|---|
| 941 | C Need the one-dimensional interpolation operator IM12
|
|---|
| 942 | C (see subroutine IGLLGL).
|
|---|
| 943 | C Note: D and DT are rectangular matrices.
|
|---|
| 944 | C
|
|---|
| 945 | C-----------------------------------------------------------------------
|
|---|
| 946 | REAL D(ND2,ND1), DT(ND1,ND2), ZM1(ND1), ZM2(ND2), IM12(ND2,ND1)
|
|---|
| 947 | IF (NZM1.EQ.1) THEN
|
|---|
| 948 | D (1,1) = 0.
|
|---|
| 949 | DT(1,1) = 0.
|
|---|
| 950 | RETURN
|
|---|
| 951 | ENDIF
|
|---|
| 952 | EPS = 1.E-6
|
|---|
| 953 | NM1 = NZM1-1
|
|---|
| 954 | DO 10 IP = 1, NZM2
|
|---|
| 955 | DO 10 JQ = 1, NZM1
|
|---|
| 956 | ZP = ZM2(IP)
|
|---|
| 957 | ZQ = ZM1(JQ)
|
|---|
| 958 | IF ((ABS(ZP) .LT. EPS).AND.(ABS(ZQ) .LT. EPS)) THEN
|
|---|
| 959 | D(IP,JQ) = 0.
|
|---|
| 960 | ELSE
|
|---|
| 961 | D(IP,JQ) = (PNLEG(ZP,NM1)/PNLEG(ZQ,NM1)
|
|---|
| 962 | $ -IM12(IP,JQ))/(ZP-ZQ)
|
|---|
| 963 | ENDIF
|
|---|
| 964 | DT(JQ,IP) = D(IP,JQ)
|
|---|
| 965 | 10 CONTINUE
|
|---|
| 966 | RETURN
|
|---|
| 967 | END
|
|---|
| 968 | C
|
|---|
| 969 | SUBROUTINE DGLJGJ (D,DT,ZGL,ZG,IGLG,NPGL,NPG,ND1,ND2,ALPHA,BETA)
|
|---|
| 970 | C-----------------------------------------------------------------------
|
|---|
| 971 | C
|
|---|
| 972 | C Compute the (one-dimensional) derivative matrix D and its
|
|---|
| 973 | C transpose DT associated with taking the derivative of a variable
|
|---|
| 974 | C expanded on a Gauss-Lobatto Jacobi mesh (M1), and evaluate its
|
|---|
| 975 | C derivative on a Guass Jacobi mesh (M2).
|
|---|
| 976 | C Need the one-dimensional interpolation operator IM12
|
|---|
| 977 | C (see subroutine IGLJGJ).
|
|---|
| 978 | C Note: D and DT are rectangular matrices.
|
|---|
| 979 | C Single precision version.
|
|---|
| 980 | C
|
|---|
| 981 | C-----------------------------------------------------------------------
|
|---|
| 982 | REAL D(ND2,ND1), DT(ND1,ND2), ZGL(ND1), ZG(ND2), IGLG(ND2,ND1)
|
|---|
| 983 | PARAMETER (NMAX=84)
|
|---|
| 984 | PARAMETER (NDD = NMAX)
|
|---|
| 985 | REAL*8 DD(NDD,NDD), DTD(NDD,NDD)
|
|---|
| 986 | REAL*8 ZGD(NDD), ZGLD(NDD), IGLGD(NDD,NDD)
|
|---|
| 987 | REAL*8 ALPHAD, BETAD
|
|---|
| 988 | C
|
|---|
| 989 | IF (NPGL.LE.1) THEN
|
|---|
| 990 | WRITE(6,*) 'DGLJGJ: Minimum number of Gauss-Lobatto points is 2'
|
|---|
| 991 | call exitt
|
|---|
| 992 | ENDIF
|
|---|
| 993 | IF (NPGL.GT.NMAX) THEN
|
|---|
| 994 | WRITE(6,*) 'Polynomial degree too high in DGLJGJ'
|
|---|
| 995 | WRITE(6,*) 'Maximum polynomial degree is',NMAX
|
|---|
| 996 | WRITE(6,*) 'Here NPGL=',NPGL
|
|---|
| 997 | call exitt
|
|---|
| 998 | ENDIF
|
|---|
| 999 | IF ((ALPHA.LE.-1.).OR.(BETA.LE.-1.)) THEN
|
|---|
| 1000 | WRITE(6,*) 'DGLJGJ: Alpha and Beta must be greater than -1'
|
|---|
| 1001 | call exitt
|
|---|
| 1002 | ENDIF
|
|---|
| 1003 | C
|
|---|
| 1004 | ALPHAD = ALPHA
|
|---|
| 1005 | BETAD = BETA
|
|---|
| 1006 | DO 100 I=1,NPG
|
|---|
| 1007 | ZGD(I) = ZG(I)
|
|---|
| 1008 | DO 100 J=1,NPGL
|
|---|
| 1009 | IGLGD(I,J) = IGLG(I,J)
|
|---|
| 1010 | 100 CONTINUE
|
|---|
| 1011 | DO 200 I=1,NPGL
|
|---|
| 1012 | ZGLD(I) = ZGL(I)
|
|---|
| 1013 | 200 CONTINUE
|
|---|
| 1014 | CALL DGLJGJD (DD,DTD,ZGLD,ZGD,IGLGD,NPGL,NPG,NDD,NDD,ALPHAD,BETAD)
|
|---|
| 1015 | DO 300 I=1,NPG
|
|---|
| 1016 | DO 300 J=1,NPGL
|
|---|
| 1017 | D(I,J) = DD(I,J)
|
|---|
| 1018 | DT(J,I) = DTD(J,I)
|
|---|
| 1019 | 300 CONTINUE
|
|---|
| 1020 | RETURN
|
|---|
| 1021 | END
|
|---|
| 1022 | C
|
|---|
| 1023 | SUBROUTINE DGLJGJD (D,DT,ZGL,ZG,IGLG,NPGL,NPG,ND1,ND2,ALPHA,BETA)
|
|---|
| 1024 | C-----------------------------------------------------------------------
|
|---|
| 1025 | C
|
|---|
| 1026 | C Compute the (one-dimensional) derivative matrix D and its
|
|---|
| 1027 | C transpose DT associated with taking the derivative of a variable
|
|---|
| 1028 | C expanded on a Gauss-Lobatto Jacobi mesh (M1), and evaluate its
|
|---|
| 1029 | C derivative on a Guass Jacobi mesh (M2).
|
|---|
| 1030 | C Need the one-dimensional interpolation operator IM12
|
|---|
| 1031 | C (see subroutine IGLJGJ).
|
|---|
| 1032 | C Note: D and DT are rectangular matrices.
|
|---|
| 1033 | C Double precision version.
|
|---|
| 1034 | C
|
|---|
| 1035 | C-----------------------------------------------------------------------
|
|---|
| 1036 | IMPLICIT REAL*8 (A-H,O-Z)
|
|---|
| 1037 | REAL*8 D(ND2,ND1), DT(ND1,ND2), ZGL(ND1), ZG(ND2)
|
|---|
| 1038 | REAL*8 IGLG(ND2,ND1), ALPHA, BETA
|
|---|
| 1039 | C
|
|---|
| 1040 | IF (NPGL.LE.1) THEN
|
|---|
| 1041 | WRITE(6,*) 'DGLJGJD: Minimum number of Gauss-Lobatto points is 2'
|
|---|
| 1042 | call exitt
|
|---|
| 1043 | ENDIF
|
|---|
| 1044 | IF ((ALPHA.LE.-1.).OR.(BETA.LE.-1.)) THEN
|
|---|
| 1045 | WRITE(6,*) 'DGLJGJD: Alpha and Beta must be greater than -1'
|
|---|
| 1046 | call exitt
|
|---|
| 1047 | ENDIF
|
|---|
| 1048 | C
|
|---|
| 1049 | EPS = 1.e-6
|
|---|
| 1050 | ONE = 1.
|
|---|
| 1051 | TWO = 2.
|
|---|
| 1052 | NGL = NPGL-1
|
|---|
| 1053 | DN = ((NGL))
|
|---|
| 1054 | EIGVAL = -DN*(DN+ALPHA+BETA+ONE)
|
|---|
| 1055 | C
|
|---|
| 1056 | DO 100 I=1,NPG
|
|---|
| 1057 | DO 100 J=1,NPGL
|
|---|
| 1058 | DZ = ABS(ZG(I)-ZGL(J))
|
|---|
| 1059 | IF (DZ.LT.EPS) THEN
|
|---|
| 1060 | D(I,J) = (ALPHA*(ONE+ZG(I))-BETA*(ONE-ZG(I)))/
|
|---|
| 1061 | $ (TWO*(ONE-ZG(I)**2))
|
|---|
| 1062 | ELSE
|
|---|
| 1063 | CALL JACOBF (PI,PDI,PM1,PDM1,PM2,PDM2,NGL,ALPHA,BETA,ZG(I))
|
|---|
| 1064 | CALL JACOBF (PJ,PDJ,PM1,PDM1,PM2,PDM2,NGL,ALPHA,BETA,ZGL(J))
|
|---|
| 1065 | FACI = ALPHA*(ONE+ZG(I))-BETA*(ONE-ZG(I))
|
|---|
| 1066 | FACJ = ALPHA*(ONE+ZGL(J))-BETA*(ONE-ZGL(J))
|
|---|
| 1067 | CONST = EIGVAL*PJ+FACJ*PDJ
|
|---|
| 1068 | D(I,J) = ((EIGVAL*PI+FACI*PDI)*(ZG(I)-ZGL(J))
|
|---|
| 1069 | $ -(ONE-ZG(I)**2)*PDI)/(CONST*(ZG(I)-ZGL(J))**2)
|
|---|
| 1070 | ENDIF
|
|---|
| 1071 | DT(J,I) = D(I,J)
|
|---|
| 1072 | 100 CONTINUE
|
|---|
| 1073 | RETURN
|
|---|
| 1074 | END
|
|---|
| 1075 | C
|
|---|
| 1076 | SUBROUTINE IGLM (I12,IT12,Z1,Z2,lz1,lz2,ND1,ND2)
|
|---|
| 1077 | C----------------------------------------------------------------------
|
|---|
| 1078 | C
|
|---|
| 1079 | C Compute the one-dimensional interpolation operator (matrix) I12
|
|---|
| 1080 | C ands its transpose IT12 for interpolating a variable from a
|
|---|
| 1081 | C Gauss Legendre mesh (1) to a another mesh M (2).
|
|---|
| 1082 | C Z1 : lz1 Gauss Legendre points.
|
|---|
| 1083 | C Z2 : lz2 points on mesh M.
|
|---|
| 1084 | C
|
|---|
| 1085 | C--------------------------------------------------------------------
|
|---|
| 1086 | REAL I12(ND2,ND1),IT12(ND1,ND2),Z1(ND1),Z2(ND2)
|
|---|
| 1087 | IF (lz1 .EQ. 1) THEN
|
|---|
| 1088 | I12 (1,1) = 1.
|
|---|
| 1089 | IT12(1,1) = 1.
|
|---|
| 1090 | RETURN
|
|---|
| 1091 | ENDIF
|
|---|
| 1092 | DO 10 I=1,lz2
|
|---|
| 1093 | ZI = Z2(I)
|
|---|
| 1094 | DO 10 J=1,lz1
|
|---|
| 1095 | I12 (I,J) = HGL(J,ZI,Z1,lz1)
|
|---|
| 1096 | IT12(J,I) = I12(I,J)
|
|---|
| 1097 | 10 CONTINUE
|
|---|
| 1098 | RETURN
|
|---|
| 1099 | END
|
|---|
| 1100 | c
|
|---|
| 1101 | SUBROUTINE IGLLM (I12,IT12,Z1,Z2,lz1,lz2,ND1,ND2)
|
|---|
| 1102 | C----------------------------------------------------------------------
|
|---|
| 1103 | C
|
|---|
| 1104 | C Compute the one-dimensional interpolation operator (matrix) I12
|
|---|
| 1105 | C ands its transpose IT12 for interpolating a variable from a
|
|---|
| 1106 | C Gauss-Lobatto Legendre mesh (1) to a another mesh M (2).
|
|---|
| 1107 | C Z1 : lz1 Gauss-Lobatto Legendre points.
|
|---|
| 1108 | C Z2 : lz2 points on mesh M.
|
|---|
| 1109 | C
|
|---|
| 1110 | C--------------------------------------------------------------------
|
|---|
| 1111 | REAL I12(ND2,ND1),IT12(ND1,ND2),Z1(ND1),Z2(ND2)
|
|---|
| 1112 | IF (lz1 .EQ. 1) THEN
|
|---|
| 1113 | I12 (1,1) = 1.
|
|---|
| 1114 | IT12(1,1) = 1.
|
|---|
| 1115 | RETURN
|
|---|
| 1116 | ENDIF
|
|---|
| 1117 | DO 10 I=1,lz2
|
|---|
| 1118 | ZI = Z2(I)
|
|---|
| 1119 | DO 10 J=1,lz1
|
|---|
| 1120 | I12 (I,J) = HGLL(J,ZI,Z1,lz1)
|
|---|
| 1121 | IT12(J,I) = I12(I,J)
|
|---|
| 1122 | 10 CONTINUE
|
|---|
| 1123 | RETURN
|
|---|
| 1124 | END
|
|---|
| 1125 | C
|
|---|
| 1126 | SUBROUTINE IGJM (I12,IT12,Z1,Z2,lz1,lz2,ND1,ND2,ALPHA,BETA)
|
|---|
| 1127 | C----------------------------------------------------------------------
|
|---|
| 1128 | C
|
|---|
| 1129 | C Compute the one-dimensional interpolation operator (matrix) I12
|
|---|
| 1130 | C ands its transpose IT12 for interpolating a variable from a
|
|---|
| 1131 | C Gauss Jacobi mesh (1) to a another mesh M (2).
|
|---|
| 1132 | C Z1 : lz1 Gauss Jacobi points.
|
|---|
| 1133 | C Z2 : lz2 points on mesh M.
|
|---|
| 1134 | C Single precision version.
|
|---|
| 1135 | C
|
|---|
| 1136 | C--------------------------------------------------------------------
|
|---|
| 1137 | REAL I12(ND2,ND1),IT12(ND1,ND2),Z1(ND1),Z2(ND2)
|
|---|
| 1138 | IF (lz1 .EQ. 1) THEN
|
|---|
| 1139 | I12 (1,1) = 1.
|
|---|
| 1140 | IT12(1,1) = 1.
|
|---|
| 1141 | RETURN
|
|---|
| 1142 | ENDIF
|
|---|
| 1143 | DO 10 I=1,lz2
|
|---|
| 1144 | ZI = Z2(I)
|
|---|
| 1145 | DO 10 J=1,lz1
|
|---|
| 1146 | I12 (I,J) = HGJ(J,ZI,Z1,lz1,ALPHA,BETA)
|
|---|
| 1147 | IT12(J,I) = I12(I,J)
|
|---|
| 1148 | 10 CONTINUE
|
|---|
| 1149 | RETURN
|
|---|
| 1150 | END
|
|---|
| 1151 | c
|
|---|
| 1152 | SUBROUTINE IGLJM (I12,IT12,Z1,Z2,lz1,lz2,ND1,ND2,ALPHA,BETA)
|
|---|
| 1153 | C----------------------------------------------------------------------
|
|---|
| 1154 | C
|
|---|
| 1155 | C Compute the one-dimensional interpolation operator (matrix) I12
|
|---|
| 1156 | C ands its transpose IT12 for interpolating a variable from a
|
|---|
| 1157 | C Gauss-Lobatto Jacobi mesh (1) to a another mesh M (2).
|
|---|
| 1158 | C Z1 : lz1 Gauss-Lobatto Jacobi points.
|
|---|
| 1159 | C Z2 : lz2 points on mesh M.
|
|---|
| 1160 | C Single precision version.
|
|---|
| 1161 | C
|
|---|
| 1162 | C--------------------------------------------------------------------
|
|---|
| 1163 | REAL I12(ND2,ND1),IT12(ND1,ND2),Z1(ND1),Z2(ND2)
|
|---|
| 1164 | IF (lz1 .EQ. 1) THEN
|
|---|
| 1165 | I12 (1,1) = 1.
|
|---|
| 1166 | IT12(1,1) = 1.
|
|---|
| 1167 | RETURN
|
|---|
| 1168 | ENDIF
|
|---|
| 1169 | DO 10 I=1,lz2
|
|---|
| 1170 | ZI = Z2(I)
|
|---|
| 1171 | DO 10 J=1,lz1
|
|---|
| 1172 | I12 (I,J) = HGLJ(J,ZI,Z1,lz1,ALPHA,BETA)
|
|---|
| 1173 | IT12(J,I) = I12(I,J)
|
|---|
| 1174 | 10 CONTINUE
|
|---|
| 1175 | RETURN
|
|---|
| 1176 | END
|
|---|
