| 1 | !-------------------------------------------------------------- |
|---|
| 2 | ! |
|---|
| 3 | ! BLAS/LAPACK-like subroutines used by the integration algorithms |
|---|
| 4 | ! It is recommended to replace them by calls to the optimized |
|---|
| 5 | ! BLAS/LAPACK library for your machine |
|---|
| 6 | ! |
|---|
| 7 | ! (C) Adrian Sandu, Aug. 2004 |
|---|
| 8 | ! Virginia Polytechnic Institute and State University |
|---|
| 9 | !-------------------------------------------------------------- |
|---|
| 10 | |
|---|
| 11 | |
|---|
| 12 | !-------------------------------------------------------------- |
|---|
| 13 | SUBROUTINE WCOPY(N,X,incX,Y,incY) |
|---|
| 14 | !-------------------------------------------------------------- |
|---|
| 15 | ! copies a vector, x, to a vector, y: y <- x |
|---|
| 16 | ! only for incX=incY=1 |
|---|
| 17 | ! after BLAS |
|---|
| 18 | ! replace this by the function from the optimized BLAS implementation: |
|---|
| 19 | ! CALL SCOPY(N,X,1,Y,1) or CALL DCOPY(N,X,1,Y,1) |
|---|
| 20 | !-------------------------------------------------------------- |
|---|
| 21 | ! USE KPP_ROOT_Precision |
|---|
| 22 | |
|---|
| 23 | INTEGER :: i,incX,incY,M,MP1,N |
|---|
| 24 | KPP_REAL :: X(N),Y(N) |
|---|
| 25 | |
|---|
| 26 | IF (N.LE.0) RETURN |
|---|
| 27 | |
|---|
| 28 | M = MOD(N,8) |
|---|
| 29 | IF( M .NE. 0 ) THEN |
|---|
| 30 | DO i = 1,M |
|---|
| 31 | Y(i) = X(i) |
|---|
| 32 | END DO |
|---|
| 33 | IF( N .LT. 8 ) RETURN |
|---|
| 34 | END IF |
|---|
| 35 | MP1 = M+1 |
|---|
| 36 | DO i = MP1,N,8 |
|---|
| 37 | Y(i) = X(i) |
|---|
| 38 | Y(i + 1) = X(i + 1) |
|---|
| 39 | Y(i + 2) = X(i + 2) |
|---|
| 40 | Y(i + 3) = X(i + 3) |
|---|
| 41 | Y(i + 4) = X(i + 4) |
|---|
| 42 | Y(i + 5) = X(i + 5) |
|---|
| 43 | Y(i + 6) = X(i + 6) |
|---|
| 44 | Y(i + 7) = X(i + 7) |
|---|
| 45 | END DO |
|---|
| 46 | |
|---|
| 47 | END SUBROUTINE WCOPY |
|---|
| 48 | |
|---|
| 49 | |
|---|
| 50 | !-------------------------------------------------------------- |
|---|
| 51 | SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY) |
|---|
| 52 | !-------------------------------------------------------------- |
|---|
| 53 | ! constant times a vector plus a vector: y <- y + Alpha*x |
|---|
| 54 | ! only for incX=incY=1 |
|---|
| 55 | ! after BLAS |
|---|
| 56 | ! replace this by the function from the optimized BLAS implementation: |
|---|
| 57 | ! CALL SAXPY(N,Alpha,X,1,Y,1) or CALL DAXPY(N,Alpha,X,1,Y,1) |
|---|
| 58 | !-------------------------------------------------------------- |
|---|
| 59 | |
|---|
| 60 | INTEGER :: i,incX,incY,M,MP1,N |
|---|
| 61 | KPP_REAL :: X(N),Y(N),Alpha |
|---|
| 62 | KPP_REAL, PARAMETER :: ZERO = 0.0_dp |
|---|
| 63 | |
|---|
| 64 | IF (Alpha .EQ. ZERO) RETURN |
|---|
| 65 | IF (N .LE. 0) RETURN |
|---|
| 66 | |
|---|
| 67 | M = MOD(N,4) |
|---|
| 68 | IF( M .NE. 0 ) THEN |
|---|
| 69 | DO i = 1,M |
|---|
| 70 | Y(i) = Y(i) + Alpha*X(i) |
|---|
| 71 | END DO |
|---|
| 72 | IF( N .LT. 4 ) RETURN |
|---|
| 73 | END IF |
|---|
| 74 | MP1 = M + 1 |
|---|
| 75 | DO i = MP1,N,4 |
|---|
| 76 | Y(i) = Y(i) + Alpha*X(i) |
|---|
| 77 | Y(i + 1) = Y(i + 1) + Alpha*X(i + 1) |
|---|
| 78 | Y(i + 2) = Y(i + 2) + Alpha*X(i + 2) |
|---|
| 79 | Y(i + 3) = Y(i + 3) + Alpha*X(i + 3) |
|---|
| 80 | END DO |
|---|
| 81 | |
|---|
| 82 | END SUBROUTINE WAXPY |
|---|
| 83 | |
|---|
| 84 | |
|---|
| 85 | |
|---|
| 86 | !-------------------------------------------------------------- |
|---|
| 87 | SUBROUTINE WSCAL(N,Alpha,X,incX) |
|---|
| 88 | !-------------------------------------------------------------- |
|---|
| 89 | ! constant times a vector: x(1:N) <- Alpha*x(1:N) |
|---|
| 90 | ! only for incX=incY=1 |
|---|
| 91 | ! after BLAS |
|---|
| 92 | ! replace this by the function from the optimized BLAS implementation: |
|---|
| 93 | ! CALL SSCAL(N,Alpha,X,1) or CALL DSCAL(N,Alpha,X,1) |
|---|
| 94 | !-------------------------------------------------------------- |
|---|
| 95 | |
|---|
| 96 | INTEGER :: i,incX,M,MP1,N |
|---|
| 97 | KPP_REAL :: X(N),Alpha |
|---|
| 98 | KPP_REAL, PARAMETER :: ZERO=0.0_dp, ONE=1.0_dp |
|---|
| 99 | |
|---|
| 100 | IF (Alpha .EQ. ONE) RETURN |
|---|
| 101 | IF (N .LE. 0) RETURN |
|---|
| 102 | |
|---|
| 103 | M = MOD(N,5) |
|---|
| 104 | IF( M .NE. 0 ) THEN |
|---|
| 105 | IF (Alpha .EQ. (-ONE)) THEN |
|---|
| 106 | DO i = 1,M |
|---|
| 107 | X(i) = -X(i) |
|---|
| 108 | END DO |
|---|
| 109 | ELSEIF (Alpha .EQ. ZERO) THEN |
|---|
| 110 | DO i = 1,M |
|---|
| 111 | X(i) = ZERO |
|---|
| 112 | END DO |
|---|
| 113 | ELSE |
|---|
| 114 | DO i = 1,M |
|---|
| 115 | X(i) = Alpha*X(i) |
|---|
| 116 | END DO |
|---|
| 117 | END IF |
|---|
| 118 | IF( N .LT. 5 ) RETURN |
|---|
| 119 | END IF |
|---|
| 120 | MP1 = M + 1 |
|---|
| 121 | IF (Alpha .EQ. (-ONE)) THEN |
|---|
| 122 | DO i = MP1,N,5 |
|---|
| 123 | X(i) = -X(i) |
|---|
| 124 | X(i + 1) = -X(i + 1) |
|---|
| 125 | X(i + 2) = -X(i + 2) |
|---|
| 126 | X(i + 3) = -X(i + 3) |
|---|
| 127 | X(i + 4) = -X(i + 4) |
|---|
| 128 | END DO |
|---|
| 129 | ELSEIF (Alpha .EQ. ZERO) THEN |
|---|
| 130 | DO i = MP1,N,5 |
|---|
| 131 | X(i) = ZERO |
|---|
| 132 | X(i + 1) = ZERO |
|---|
| 133 | X(i + 2) = ZERO |
|---|
| 134 | X(i + 3) = ZERO |
|---|
| 135 | X(i + 4) = ZERO |
|---|
| 136 | END DO |
|---|
| 137 | ELSE |
|---|
| 138 | DO i = MP1,N,5 |
|---|
| 139 | X(i) = Alpha*X(i) |
|---|
| 140 | X(i + 1) = Alpha*X(i + 1) |
|---|
| 141 | X(i + 2) = Alpha*X(i + 2) |
|---|
| 142 | X(i + 3) = Alpha*X(i + 3) |
|---|
| 143 | X(i + 4) = Alpha*X(i + 4) |
|---|
| 144 | END DO |
|---|
| 145 | END IF |
|---|
| 146 | |
|---|
| 147 | END SUBROUTINE WSCAL |
|---|
| 148 | |
|---|
| 149 | !-------------------------------------------------------------- |
|---|
| 150 | KPP_REAL FUNCTION WLAMCH( C ) |
|---|
| 151 | !-------------------------------------------------------------- |
|---|
| 152 | ! returns epsilon machine |
|---|
| 153 | ! after LAPACK |
|---|
| 154 | ! replace this by the function from the optimized LAPACK implementation: |
|---|
| 155 | ! CALL SLAMCH('E') or CALL DLAMCH('E') |
|---|
| 156 | !-------------------------------------------------------------- |
|---|
| 157 | ! USE KPP_ROOT_Precision |
|---|
| 158 | |
|---|
| 159 | CHARACTER :: C |
|---|
| 160 | INTEGER :: i |
|---|
| 161 | KPP_REAL, SAVE :: Eps |
|---|
| 162 | KPP_REAL :: Suma |
|---|
| 163 | KPP_REAL, PARAMETER :: ONE=1.0_dp, HALF=0.5_dp |
|---|
| 164 | LOGICAL, SAVE :: First=.TRUE. |
|---|
| 165 | |
|---|
| 166 | IF (First) THEN |
|---|
| 167 | First = .FALSE. |
|---|
| 168 | Eps = HALF**(16) |
|---|
| 169 | DO i = 17, 80 |
|---|
| 170 | Eps = Eps*HALF |
|---|
| 171 | CALL WLAMCH_ADD(ONE,Eps,Suma) |
|---|
| 172 | IF (Suma.LE.ONE) GOTO 10 |
|---|
| 173 | END DO |
|---|
| 174 | PRINT*,'ERROR IN WLAMCH. EPS < ',Eps |
|---|
| 175 | RETURN |
|---|
| 176 | 10 Eps = Eps*2 |
|---|
| 177 | i = i-1 |
|---|
| 178 | END IF |
|---|
| 179 | |
|---|
| 180 | WLAMCH = Eps |
|---|
| 181 | |
|---|
| 182 | END FUNCTION WLAMCH |
|---|
| 183 | |
|---|
| 184 | SUBROUTINE WLAMCH_ADD( A, B, Suma ) |
|---|
| 185 | ! USE KPP_ROOT_Precision |
|---|
| 186 | |
|---|
| 187 | KPP_REAL A, B, Suma |
|---|
| 188 | Suma = A + B |
|---|
| 189 | |
|---|
| 190 | END SUBROUTINE WLAMCH_ADD |
|---|
| 191 | !-------------------------------------------------------------- |
|---|
| 192 | |
|---|
| 193 | |
|---|
| 194 | !-------------------------------------------------------------- |
|---|
| 195 | SUBROUTINE SET2ZERO(N,Y) |
|---|
| 196 | !-------------------------------------------------------------- |
|---|
| 197 | ! copies zeros into the vector y: y <- 0 |
|---|
| 198 | ! after BLAS |
|---|
| 199 | !-------------------------------------------------------------- |
|---|
| 200 | |
|---|
| 201 | INTEGER :: i,M,MP1,N |
|---|
| 202 | KPP_REAL :: Y(N) |
|---|
| 203 | KPP_REAL, PARAMETER :: ZERO = 0.0d0 |
|---|
| 204 | |
|---|
| 205 | IF (N.LE.0) RETURN |
|---|
| 206 | |
|---|
| 207 | M = MOD(N,8) |
|---|
| 208 | IF( M .NE. 0 ) THEN |
|---|
| 209 | DO i = 1,M |
|---|
| 210 | Y(i) = ZERO |
|---|
| 211 | END DO |
|---|
| 212 | IF( N .LT. 8 ) RETURN |
|---|
| 213 | END IF |
|---|
| 214 | MP1 = M+1 |
|---|
| 215 | DO i = MP1,N,8 |
|---|
| 216 | Y(i) = ZERO |
|---|
| 217 | Y(i + 1) = ZERO |
|---|
| 218 | Y(i + 2) = ZERO |
|---|
| 219 | Y(i + 3) = ZERO |
|---|
| 220 | Y(i + 4) = ZERO |
|---|
| 221 | Y(i + 5) = ZERO |
|---|
| 222 | Y(i + 6) = ZERO |
|---|
| 223 | Y(i + 7) = ZERO |
|---|
| 224 | END DO |
|---|
| 225 | |
|---|
| 226 | END SUBROUTINE SET2ZERO |
|---|
| 227 | |
|---|
| 228 | |
|---|
| 229 | !-------------------------------------------------------------- |
|---|
| 230 | KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY) |
|---|
| 231 | !-------------------------------------------------------------- |
|---|
| 232 | ! dot produce: wdot = x(1:N)*y(1:N) |
|---|
| 233 | ! only for incX=incY=1 |
|---|
| 234 | ! after BLAS |
|---|
| 235 | ! replace this by the function from the optimized BLAS implementation: |
|---|
| 236 | ! CALL SDOT(N,X,1,Y,1) or CALL DDOT(N,X,1,Y,1) |
|---|
| 237 | !-------------------------------------------------------------- |
|---|
| 238 | ! USE messy_mecca_kpp_Precision |
|---|
| 239 | !-------------------------------------------------------------- |
|---|
| 240 | IMPLICIT NONE |
|---|
| 241 | INTEGER :: N, incX, incY |
|---|
| 242 | KPP_REAL :: DX(N), DY(N) |
|---|
| 243 | |
|---|
| 244 | INTEGER :: i, IX, IY, M, MP1, NS |
|---|
| 245 | |
|---|
| 246 | WDOT = 0.0D0 |
|---|
| 247 | IF (N .LE. 0) RETURN |
|---|
| 248 | IF (incX .EQ. incY) IF (incX-1) 5,20,60 |
|---|
| 249 | ! |
|---|
| 250 | ! Code for unequal or nonpositive increments. |
|---|
| 251 | ! |
|---|
| 252 | 5 IX = 1 |
|---|
| 253 | IY = 1 |
|---|
| 254 | IF (incX .LT. 0) IX = (-N+1)*incX + 1 |
|---|
| 255 | IF (incY .LT. 0) IY = (-N+1)*incY + 1 |
|---|
| 256 | DO i = 1,N |
|---|
| 257 | WDOT = WDOT + DX(IX)*DY(IY) |
|---|
| 258 | IX = IX + incX |
|---|
| 259 | IY = IY + incY |
|---|
| 260 | END DO |
|---|
| 261 | RETURN |
|---|
| 262 | ! |
|---|
| 263 | ! Code for both increments equal to 1. |
|---|
| 264 | ! |
|---|
| 265 | ! Clean-up loop so remaining vector length is a multiple of 5. |
|---|
| 266 | ! |
|---|
| 267 | 20 M = MOD(N,5) |
|---|
| 268 | IF (M .EQ. 0) GO TO 40 |
|---|
| 269 | DO i = 1,M |
|---|
| 270 | WDOT = WDOT + DX(i)*DY(i) |
|---|
| 271 | END DO |
|---|
| 272 | IF (N .LT. 5) RETURN |
|---|
| 273 | 40 MP1 = M + 1 |
|---|
| 274 | DO i = MP1,N,5 |
|---|
| 275 | WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) + & |
|---|
| 276 | DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4) |
|---|
| 277 | END DO |
|---|
| 278 | RETURN |
|---|
| 279 | ! |
|---|
| 280 | ! Code for equal, positive, non-unit increments. |
|---|
| 281 | ! |
|---|
| 282 | 60 NS = N*incX |
|---|
| 283 | DO i = 1,NS,incX |
|---|
| 284 | WDOT = WDOT + DX(i)*DY(i) |
|---|
| 285 | END DO |
|---|
| 286 | |
|---|
| 287 | END FUNCTION WDOT |
|---|
| 288 | |
|---|
| 289 | |
|---|
| 290 | !-------------------------------------------------------------- |
|---|
| 291 | SUBROUTINE WADD(N,X,Y,Z) |
|---|
| 292 | !-------------------------------------------------------------- |
|---|
| 293 | ! adds two vectors: z <- x + y |
|---|
| 294 | ! BLAS - like |
|---|
| 295 | !-------------------------------------------------------------- |
|---|
| 296 | ! USE KPP_ROOT_Precision |
|---|
| 297 | |
|---|
| 298 | INTEGER :: i, M, MP1, N |
|---|
| 299 | KPP_REAL :: X(N),Y(N),Z(N) |
|---|
| 300 | |
|---|
| 301 | IF (N.LE.0) RETURN |
|---|
| 302 | |
|---|
| 303 | M = MOD(N,5) |
|---|
| 304 | IF( M /= 0 ) THEN |
|---|
| 305 | DO i = 1,M |
|---|
| 306 | Z(i) = X(i) + Y(i) |
|---|
| 307 | END DO |
|---|
| 308 | IF( N < 5 ) RETURN |
|---|
| 309 | END IF |
|---|
| 310 | MP1 = M+1 |
|---|
| 311 | DO i = MP1,N,5 |
|---|
| 312 | Z(i) = X(i) + Y(i) |
|---|
| 313 | Z(i + 1) = X(i + 1) + Y(i + 1) |
|---|
| 314 | Z(i + 2) = X(i + 2) + Y(i + 2) |
|---|
| 315 | Z(i + 3) = X(i + 3) + Y(i + 3) |
|---|
| 316 | Z(i + 4) = X(i + 4) + Y(i + 4) |
|---|
| 317 | END DO |
|---|
| 318 | |
|---|
| 319 | END SUBROUTINE WADD |
|---|
| 320 | |
|---|
| 321 | |
|---|
| 322 | |
|---|
| 323 | !-------------------------------------------------------------- |
|---|
| 324 | SUBROUTINE WGEFA(N,A,Ipvt,info) |
|---|
| 325 | !-------------------------------------------------------------- |
|---|
| 326 | ! WGEFA FACTORS THE MATRIX A (N,N) BY |
|---|
| 327 | ! GAUSS ELIMINATION WITH PARTIAL PIVOTING |
|---|
| 328 | ! LINPACK - LIKE |
|---|
| 329 | !-------------------------------------------------------------- |
|---|
| 330 | ! |
|---|
| 331 | INTEGER :: N,Ipvt(N),info |
|---|
| 332 | KPP_REAL :: A(N,N) |
|---|
| 333 | KPP_REAL :: t, dmax, da |
|---|
| 334 | INTEGER :: j,k,l |
|---|
| 335 | KPP_REAL, PARAMETER :: ZERO = 0.0, ONE = 1.0 |
|---|
| 336 | |
|---|
| 337 | info = 0 |
|---|
| 338 | |
|---|
| 339 | size: IF (n > 1) THEN |
|---|
| 340 | |
|---|
| 341 | col: DO k = 1, n-1 |
|---|
| 342 | |
|---|
| 343 | ! find l = pivot index |
|---|
| 344 | ! l = idamax(n-k+1,A(k,k),1) + k - 1 |
|---|
| 345 | l = k; dmax = abs(A(k,k)) |
|---|
| 346 | DO j = k+1,n |
|---|
| 347 | da = ABS(A(j,k)) |
|---|
| 348 | IF (da > dmax) THEN |
|---|
| 349 | l = j; dmax = da |
|---|
| 350 | END IF |
|---|
| 351 | END DO |
|---|
| 352 | Ipvt(k) = l |
|---|
| 353 | |
|---|
| 354 | ! zero pivot implies this column already triangularized |
|---|
| 355 | IF (ABS(A(l,k)) < TINY(ZERO)) THEN |
|---|
| 356 | info = k |
|---|
| 357 | return |
|---|
| 358 | ELSE |
|---|
| 359 | IF (l /= k) THEN |
|---|
| 360 | t = A(l,k); A(l,k) = A(k,k); A(k,k) = t |
|---|
| 361 | END IF |
|---|
| 362 | t = -ONE/A(k,k) |
|---|
| 363 | CALL WSCAL(n-k,t,A(k+1,k),1) |
|---|
| 364 | DO j = k+1, n |
|---|
| 365 | t = A(l,j) |
|---|
| 366 | IF (l /= k) THEN |
|---|
| 367 | A(l,j) = A(k,j); A(k,j) = t |
|---|
| 368 | END IF |
|---|
| 369 | CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1) |
|---|
| 370 | END DO |
|---|
| 371 | END IF |
|---|
| 372 | |
|---|
| 373 | END DO col |
|---|
| 374 | |
|---|
| 375 | END IF size |
|---|
| 376 | |
|---|
| 377 | Ipvt(N) = N |
|---|
| 378 | IF (ABS(A(N,N)) == ZERO) info = N |
|---|
| 379 | |
|---|
| 380 | END SUBROUTINE WGEFA |
|---|
| 381 | |
|---|
| 382 | |
|---|
| 383 | !-------------------------------------------------------------- |
|---|
| 384 | SUBROUTINE WGESL(Trans,N,A,Ipvt,b) |
|---|
| 385 | !-------------------------------------------------------------- |
|---|
| 386 | ! WGESL solves the system |
|---|
| 387 | ! a * x = b or trans(a) * x = b |
|---|
| 388 | ! using the factors computed by WGEFA. |
|---|
| 389 | ! |
|---|
| 390 | ! Trans = 'N' to solve A*x = b , |
|---|
| 391 | ! = 'T' to solve transpose(A)*x = b |
|---|
| 392 | ! LINPACK - LIKE |
|---|
| 393 | !-------------------------------------------------------------- |
|---|
| 394 | |
|---|
| 395 | INTEGER :: N,Ipvt(N) |
|---|
| 396 | CHARACTER :: trans |
|---|
| 397 | KPP_REAL :: A(N,N),b(N) |
|---|
| 398 | KPP_REAL :: t |
|---|
| 399 | INTEGER :: k,kb,l |
|---|
| 400 | |
|---|
| 401 | |
|---|
| 402 | SELECT CASE (Trans) |
|---|
| 403 | |
|---|
| 404 | CASE ('n','N') ! Solve A * x = b |
|---|
| 405 | |
|---|
| 406 | ! first solve L*y = b |
|---|
| 407 | IF (n >= 2) THEN |
|---|
| 408 | DO k = 1, n-1 |
|---|
| 409 | l = Ipvt(k) |
|---|
| 410 | t = b(l) |
|---|
| 411 | IF (l /= k) THEN |
|---|
| 412 | b(l) = b(k) |
|---|
| 413 | b(k) = t |
|---|
| 414 | END IF |
|---|
| 415 | CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1) |
|---|
| 416 | END DO |
|---|
| 417 | END IF |
|---|
| 418 | ! now solve U*x = y |
|---|
| 419 | DO kb = 1, n |
|---|
| 420 | k = n + 1 - kb |
|---|
| 421 | b(k) = b(k)/a(k,k) |
|---|
| 422 | t = -b(k) |
|---|
| 423 | CALL WAXPY(k-1,t,a(1,k),1,b(1),1) |
|---|
| 424 | END DO |
|---|
| 425 | |
|---|
| 426 | CASE ('t','T') ! Solve transpose(A) * x = b |
|---|
| 427 | |
|---|
| 428 | ! first solve trans(U)*y = b |
|---|
| 429 | DO k = 1, n |
|---|
| 430 | t = WDOT(k-1,a(1,k),1,b(1),1) |
|---|
| 431 | b(k) = (b(k) - t)/a(k,k) |
|---|
| 432 | END DO |
|---|
| 433 | ! now solve trans(L)*x = y |
|---|
| 434 | IF (n >= 2) THEN |
|---|
| 435 | DO kb = 1, n-1 |
|---|
| 436 | k = n - kb |
|---|
| 437 | b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1) |
|---|
| 438 | l = Ipvt(k) |
|---|
| 439 | IF (l /= k) THEN |
|---|
| 440 | t = b(l); b(l) = b(k); b(k) = t |
|---|
| 441 | END IF |
|---|
| 442 | END DO |
|---|
| 443 | END IF |
|---|
| 444 | |
|---|
| 445 | END SELECT |
|---|
| 446 | |
|---|
| 447 | END SUBROUTINE WGESL |
|---|
