diff --git a/fish90/src/hst_procs.f90 b/fish90/src/hst_procs.f90 new file mode 100644 index 0000000000000000000000000000000000000000..f438ec0416c8da45c55da0e5911711327e42946b --- /dev/null +++ b/fish90/src/hst_procs.f90 @@ -0,0 +1,1650 @@ +! +! file hstcrt.f +! +! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . +! . . +! . copyright (c) 2004 by UCAR . +! . . +! . UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH . +! . . +! . all rights reserved . +! . . +! . . +! . FISHPACK version 5.0 . +! . . +! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . +! +! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +! * * +! * F I S H P A C K * +! * * +! * * +! * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * +! * * +! * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * +! * * +! * (Version 5.0 , JUNE 2004) * +! * * +! * BY * +! * * +! * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * +! * * +! * OF * +! * * +! * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * +! * * +! * BOULDER, COLORADO (80307) U.S.A. * +! * * +! * WHICH IS SPONSORED BY * +! * * +! * THE NATIONAL SCIENCE FOUNDATION * +! * * +! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +! +! SUBROUTINE HSTCRT (A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD, +! + ELMBDA,F,IDIMF,PERTRB,IERROR) +! +! DIMENSION OF BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N) +! ARGUMENTS +! +! LATEST REVISION June 2004 +! +! PURPOSE SOLVES THE STANDARD FIVE-POINT FINITE +! DIFFERENCE APPROXIMATION TO THE HELMHOLTZ +! EQUATION +! (D/DX)(DU/DX) + (D/DY)(DU/DY) + LAMBDA*U +! = F(X,Y) +! ON A STAGGERED GRID IN CARTESIAN COORDINATES. +! +! USAGE CALL HSTCRT (A,B,M,MBDCND,BDA,BDB,C,D +! N,NBDCND,BDC,BDD,ELMBDA, +! F,IDIMF,PERTRB,IERROR) +! +! ARGUMENTS +! ON INPUT +! +! A,B +! THE RANGE OF X, I.E. A .LE. X .LE. B. +! A MUST BE LESS THAN B. +! +! M +! THE NUMBER OF GRID POINTS IN THE +! INTERVAL (A,B). THE GRID POINTS +! IN THE X-DIRECTION ARE GIVEN BY +! X(I) = A + (I-0.5)DX FOR I=1,2,...,M +! WHERE DX =(B-A)/M. M MUST BE GREATER +! THAN 2. +! +! MBDCND +! INDICATES THE TYPE OF BOUNDARY CONDITIONS +! AT X = A AND X = B. +! +! = 0 IF THE SOLUTION IS PERIODIC IN X, +! U(M+I,J) = U(I,J). +! +! = 1 IF THE SOLUTION IS SPECIFIED AT +! X = A AND X = B. +! +! = 2 IF THE SOLUTION IS SPECIFIED AT +! X = A AND THE DERIVATIVE +! OF THE SOLUTION WITH RESPECT TO X +! IS SPECIFIED AT X = B. +! +! = 3 IF THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO X IS SPECIFIED +! AT X = A AND X = B. +! +! = 4 IF THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO X IS SPECIFIED +! AT X = A AND THE SOLUTION IS +! SPECIFIED AT X = B. +! +! BDA +! A ONE-DIMENSIONAL ARRAY OF LENGTH N +! THAT SPECIFIES THE BOUNDARY VALUES +! (IF ANY) OF THE SOLUTION AT X = A. +! +! WHEN MBDCND = 1 OR 2, +! BDA(J) = U(A,Y(J)) , J=1,2,...,N. +! +! WHEN MBDCND = 3 OR 4, +! BDA(J) = (D/DX)U(A,Y(J)) , J=1,2,...,N. +! +! BDB +! A ONE-DIMENSIONAL ARRAY OF LENGTH N +! THAT SPECIFIES THE BOUNDARY VALUES +! OF THE SOLUTION AT X = B. +! +! WHEN MBDCND = 1 OR 4 +! BDB(J) = U(B,Y(J)) , J=1,2,...,N. +! +! WHEN MBDCND = 2 OR 3 +! BDB(J) = (D/DX)U(B,Y(J)) , J=1,2,...,N. +! +! C,D +! THE RANGE OF Y, I.E. C .LE. Y .LE. D. +! C MUST BE LESS THAN D. +! +! +! N +! THE NUMBER OF UNKNOWNS IN THE INTERVAL +! (C,D). THE UNKNOWNS IN THE Y-DIRECTION +! ARE GIVEN BY Y(J) = C + (J-0.5)DY, +! J=1,2,...,N, WHERE DY = (D-C)/N. +! N MUST BE GREATER THAN 2. +! +! NBDCND +! INDICATES THE TYPE OF BOUNDARY CONDITIONS +! AT Y = C AND Y = D. +! +! +! = 0 IF THE SOLUTION IS PERIODIC IN Y, I.E. +! U(I,J) = U(I,N+J). +! +! = 1 IF THE SOLUTION IS SPECIFIED AT Y = C +! AND Y = D. +! +! = 2 IF THE SOLUTION IS SPECIFIED AT Y = C +! AND THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO Y IS SPECIFIED AT +! Y = D. +! +! = 3 IF THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO Y IS SPECIFIED AT +! Y = C AND Y = D. +! +! = 4 IF THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO Y IS SPECIFIED AT +! Y = C AND THE SOLUTION IS SPECIFIED +! AT Y = D. +! +! BDC +! A ONE DIMENSIONAL ARRAY OF LENGTH M THAT +! SPECIFIES THE BOUNDARY VALUES OF THE +! SOLUTION AT Y = C. +! +! WHEN NBDCND = 1 OR 2, +! BDC(I) = U(X(I),C) , I=1,2,...,M. +! +! WHEN NBDCND = 3 OR 4, +! BDC(I) = (D/DY)U(X(I),C), I=1,2,...,M. +! +! WHEN NBDCND = 0, BDC IS A DUMMY VARIABLE. +! +! BDD +! A ONE-DIMENSIONAL ARRAY OF LENGTH M THAT +! SPECIFIES THE BOUNDARY VALUES OF THE +! SOLUTION AT Y = D. +! +! WHEN NBDCND = 1 OR 4, +! BDD(I) = U(X(I),D) , I=1,2,...,M. +! +! WHEN NBDCND = 2 OR 3, +! BDD(I) = (D/DY)U(X(I),D) , I=1,2,...,M. +! +! WHEN NBDCND = 0, BDD IS A DUMMY VARIABLE. +! +! ELMBDA +! THE CONSTANT LAMBDA IN THE HELMHOLTZ +! EQUATION. IF LAMBDA IS GREATER THAN 0, +! A SOLUTION MAY NOT EXIST. HOWEVER, +! HSTCRT WILL ATTEMPT TO FIND A SOLUTION. +! +! F +! A TWO-DIMENSIONAL ARRAY THAT SPECIFIES +! THE VALUES OF THE RIGHT SIDE OF THE +! HELMHOLTZ EQUATION. FOR I=1,2,...,M +! AND J=1,2,...,N +! +! F(I,J) = F(X(I),Y(J)) . +! +! F MUST BE DIMENSIONED AT LEAST M X N. +! +! IDIMF +! THE ROW (OR FIRST) DIMENSION OF THE ARRAY +! F AS IT APPEARS IN THE PROGRAM CALLING +! HSTCRT. THIS PARAMETER IS USED TO SPECIFY +! THE VARIABLE DIMENSION OF F. +! IDIMF MUST BE AT LEAST M. +! +! +! ON OUTPUT F +! CONTAINS THE SOLUTION U(I,J) OF THE FINITE +! DIFFERENCE APPROXIMATION FOR THE GRID POINT +! (X(I),Y(J)) FOR I=1,2,...,M, J=1,2,...,N. +! +! PERTRB +! IF A COMBINATION OF PERIODIC OR DERIVATIVE +! BOUNDARY CONDITIONS IS SPECIFIED FOR A +! POISSON EQUATION (LAMBDA = 0), A SOLUTION +! MAY NOT EXIST. PERTRB IS A CONSTANT, +! CALCULATED AND SUBTRACTED FROM F, WHICH +! ENSURES THAT A SOLUTION EXISTS. HSTCRT +! THEN COMPUTES THIS SOLUTION, WHICH IS A +! LEAST SQUARES SOLUTION TO THE ORIGINAL +! APPROXIMATION. THIS SOLUTION PLUS ANY +! CONSTANT IS ALSO A SOLUTION; HENCE, THE +! SOLUTION IS NOT UNIQUE. THE VALUE OF +! PERTRB SHOULD BE SMALL COMPARED TO THE +! RIGHT SIDE F. OTHERWISE, A SOLUTION IS +! OBTAINED TO AN ESSENTIALLY DIFFERENT PROBLEM. +! THIS COMPARISON SHOULD ALWAYS BE MADE TO +! INSURE THAT A MEANINGFUL SOLUTION HAS BEEN +! OBTAINED. +! +! IERROR +! AN ERROR FLAG THAT INDICATES INVALID INPUT +! PARAMETERS. EXCEPT TO NUMBERS 0 AND 6, +! A SOLUTION IS NOT ATTEMPTED. +! +! = 0 NO ERROR +! +! = 1 A .GE. B +! +! = 2 MBDCND .LT. 0 OR MBDCND .GT. 4 +! +! = 3 C .GE. D +! +! = 4 N .LE. 2 +! +! = 5 NBDCND .LT. 0 OR NBDCND .GT. 4 +! +! = 6 LAMBDA .GT. 0 +! +! = 7 IDIMF .LT. M +! +! = 8 M .LE. 2 +! +! SINCE THIS IS THE ONLY MEANS OF INDICATING +! A POSSIBLY INCORRECT CALL TO HSTCRT, THE +! USER SHOULD TEST IERROR AFTER THE CALL. +! +! = 20 If the dynamic allocation of real and +! complex work space required for solution +! fails (for example if N,M are too large +! for your computer) +! +! +! I/O NONE +! +! PRECISION SINGLE +! +! REQUIRED LIBRARY fish.f,comf.f,genbun.f,gnbnaux.f,poistg.f +! FILES +! +! LANGUAGE FORTRAN 90 +! +! HISTORY WRITTEN BY ROLAND SWEET AT NCAR IN 1977. +! RELEASED ON NCAR'S PUBLIC SOFTWARE LIBRARIES +! IN JANUARY 1980. +! Revised in June 2004 by John Adams using +! Fortran 90 dynamically allocated work space. +! +! PORTABILITY FORTRAN 90 +! +! ALGORITHM THIS SUBROUTINE DEFINES THE FINITE-DIFFERENCE +! EQUATIONS, INCORPORATES BOUNDARY DATA, ADJUSTS +! THE RIGHT SIDE WHEN THE SYSTEM IS SINGULAR +! AND CALLS EITHER POISTG OR GENBUN WHICH SOLVES +! THE LINEAR SYSTEM OF EQUATIONS. +! +! TIMING FOR LARGE M AND N, THE OPERATION COUNT +! IS ROUGHLY PROPORTIONAL TO M*N*LOG2(N). +! +! ACCURACY THE SOLUTION PROCESS EMPLOYED RESULTS IN A +! LOSS OF NO MORE THAN FOUR SIGNIFICANT DIGITS +! FOR N AND M AS LARGE AS 64. MORE DETAILED +! INFORMATION ABOUT ACCURACY CAN BE FOUND IN +! THE DOCUMENTATION FOR PACKAGE POISTG WHICH +! SOLVES THE FINITE DIFFERENCE EQUATIONS. +! +! REFERENCES U. SCHUMANN AND R. SWEET,"A DIRECT METHOD +! FOR THE SOLUTION OF POISSON'S EQUATION WITH +! BOUNDARY CONDITIONS ON A STAGGERED GRID OF +! ARBITRARY SIZE," J. COMP. PHYS. 20(1976), +! PP. 171-182. +!*********************************************************************** + SUBROUTINE HSTCRT(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, NBDCND, & + vecBDC, vecBDD, ELMBDA, matF, IDIMF, PERTRB, IERROR) + USE fish + implicit none + + DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & + ONE = 1.0D0, TWO = 2.0D0, & + FOUR = 4.0D0 + +!----------------------------------------------- +! D u m m y A r g u m e n t s +!----------------------------------------------- + INTEGER,INTENT(IN) :: M,MBDCND,N,NBDCND,IDIMF + INTEGER,INTENT(OUT) :: IERROR + DOUBLE PRECISION,INTENT(IN) :: A,B,C,D,ELMBDA + DOUBLE PRECISION,INTENT(OUT) :: PERTRB + DOUBLE PRECISION,DIMENSION(N),INTENT(IN) :: vecBDA,vecBDB + DOUBLE PRECISION,DIMENSION(M),INTENT(IN) :: vecBDC,vecBDD + DOUBLE PRECISION,DIMENSION(IDIMF,N),INTENT(INOUT) :: matF +!----------------------------------------------- +! Allocatable arrays +!----------------------------------------------- + DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: work +!----------------------------------------------- +! L o c a l V a r i a b l e s +!----------------------------------------------- + INTEGER :: irwk, istatus + + ! JANNIS: add interface + interface + SUBROUTINE HSTCRTT(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, NBDCND, & + vecBDC, vecBDD, ELMBDA, matF, IDIMF, PERTRB, IERROR, vecW, IW) + INTEGER,INTENT(IN) :: M,MBDCND,N,NBDCND,IDIMF,IW + INTEGER,INTENT(OUT) :: IERROR + DOUBLE PRECISION,INTENT(IN) :: A,B,C,D,ELMBDA + DOUBLE PRECISION,INTENT(OUT) :: PERTRB + DOUBLE PRECISION,DIMENSION(N),INTENT(IN) :: vecBDA,vecBDB + DOUBLE PRECISION,DIMENSION(M),INTENT(IN) :: vecBDC,vecBDD + DOUBLE PRECISION,DIMENSION(IDIMF,N),INTENT(INOUT) :: matF + DOUBLE PRECISION,DIMENSION(IW),INTENT(INOUT) :: vecW + end subroutine HSTCRTT + end interface +!----------------------------------------------- +! +! CHECK FOR INVALID PARAMETERS. +! + IERROR = 0 + + IF (A >= B) IERROR = 1 + IF (MBDCND<0 .OR. MBDCND>4) IERROR = 2 + IF (C >= D) IERROR = 3 + IF (N <= 2) IERROR = 4 + IF (NBDCND<0 .OR. NBDCND>4) IERROR = 5 + IF (IDIMF < M) IERROR = 7 + IF (M <= 2) IERROR = 8 + IF (IERROR /= 0) RETURN +!C! write(*,*) 'hstcrt: vecBDC:',vecBDC(1:size(vecBDC)) +!C! write(*,*) 'hstcrt: BDD:',BDD(:) +!C! write(*,*) 'hstcrt: ELMBDA:',ELMBDA +!C! write(*,*) 'hstcrt: matF:',matF(1:IDIMF,:) +!C! write(*,*) 'hstcrt: IDIMF:',IDIMF +! compute and allocate required real work space + CALL GEN_SPACE (N, M, IRWK) + irwk = irwk + 3*M + allocate(work(irwK),STAT=istatus) +!@! write(*,*) 'HSTCRT: allocate work(irwk); irwk=',irwk +! return if allocation failed (e.g., if n,m are too large) + IF (istatus > 0) THEN + write(*,*) 'HSTCRT: error allocate work(irwk); irwk=',irwk + RETURN + END IF + +! check that allocation was successful + call hstcrtt(a,b,m,mbdcnd,vecbda,vecbdb,c,d,n,nbdcnd, & + vecbdc,vecbdd,elmbda,matf,idimf,pertrb,ierror, & + work,size(work)) +! +!@! write(*,*) 'HSTCRT: deallocate work(irwk); irwk=',irwk + deallocate(work,STAT=istatus) + IF (istatus > 0) THEN + write(*,*) 'HSTCRT: error deallocate work' + RETURN + END IF +! + END SUBROUTINE HSTCRT + + SUBROUTINE HSTCRTT(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, NBDCND, & + vecBDC, vecBDD, ELMBDA, matF, IDIMF, PERTRB, & + IERROR, vecW, IW) + +! USE genbunal +! USE poisson + + implicit none + + DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & + ONE = 1.0D0, TWO = 2.0D0, & + FOUR = 4.0D0 + +!----------------------------------------------- +! D u m m y A r g u m e n t s +!----------------------------------------------- + INTEGER,INTENT(IN) :: M,MBDCND,N,NBDCND,IDIMF,IW + INTEGER,INTENT(OUT) :: IERROR + DOUBLE PRECISION,INTENT(IN) :: A,B,C,D,ELMBDA + DOUBLE PRECISION,INTENT(OUT) :: PERTRB + DOUBLE PRECISION,DIMENSION(N),INTENT(IN) :: vecBDA,vecBDB + DOUBLE PRECISION,DIMENSION(M),INTENT(IN) :: vecBDC,vecBDD + DOUBLE PRECISION,DIMENSION(IDIMF,N),INTENT(INOUT) :: matF + DOUBLE PRECISION,DIMENSION(IW),INTENT(INOUT) :: vecW +!----------------------------------------------- +! L o c a l V a r i a b l e s +!----------------------------------------------- + INTEGER :: NPEROD, MPEROD, NP, MP, ID2, ID3, ID4, & + J, IERR1, IW2 + DOUBLE PRECISION :: DELTAX,TWDELX,DELXSQ,DELTAY,TWDELY,DELYSQ,TWDYSQ,S,ST2 +!----------------------------------------------- + + NPEROD = NBDCND + MPEROD = 0 + IF (MBDCND > 0) MPEROD = 1 + DELTAX = (B - A)/DBLE(M) + TWDELX = ONE/DELTAX + DELXSQ = TWO/DELTAX**2 + DELTAY = (D - C)/DBLE(N) + TWDELY = ONE/DELTAY + DELYSQ = DELTAY**2 + TWDYSQ = TWO/DELYSQ + NP = NBDCND + 1 + MP = MBDCND + 1 +! +! DEFINE THE A,B,C COEFFICIENTS IN W-ARRAY. +! + ID2 = M + ID3 = ID2 + M + ID4 = ID3 + M + S = (DELTAY/DELTAX)**2 + ST2 = TWO*S + vecW(:M) = S + vecW(ID2+1:M+ID2) = (-ST2) + ELMBDA*DELYSQ + vecW(ID3+1:M+ID3) = S +! +! ENTER BOUNDARY DATA FOR X-BOUNDARIES. +! + GO TO (111,102,102,104,104) MP + 102 CONTINUE + matF(1,:N) = matF(1,:N) - vecBDA(:N)*DELXSQ + vecW(ID2+1) = vecW(ID2+1) - vecW(1) + GO TO 106 + 104 CONTINUE + matF(1,:N) = matF(1,:N) + vecBDA(:N)*TWDELX + vecW(ID2+1) = vecW(ID2+1) + vecW(1) + 106 CONTINUE + GO TO (111,107,109,109,107) MP + 107 CONTINUE + matF(M,:N) = matF(M,:N) - vecBDB(:N)*DELXSQ + vecW(ID3) = vecW(ID3) - vecW(1) + GO TO 111 + 109 CONTINUE + matF(M,:N) = matF(M,:N) - vecBDB(:N)*TWDELX + vecW(ID3) = vecW(ID3) + vecW(1) + 111 CONTINUE + GO TO (121,112,112,114,114) NP + 112 CONTINUE + matF(:M,1) = matF(:M,1) - vecBDC(:M)*TWDYSQ + GO TO 116 + 114 CONTINUE + matF(:M,1) = matF(:M,1) + vecBDC(:M)*TWDELY + 116 CONTINUE + GO TO (121,117,119,119,117) NP + 117 CONTINUE + matF(:M,N) = matF(:M,N) - vecBDD(:M)*TWDYSQ + GO TO 121 + 119 CONTINUE + matF(:M,N) = matF(:M,N) - vecBDD(:M)*TWDELY + 121 CONTINUE + matF(:M,:N) = matF(:M,:N)*DELYSQ + IF (MPEROD /= 0) THEN + vecW(1) = ZERO + vecW(ID4) = ZERO + ENDIF + PERTRB = ZERO + IF (ELMBDA >= ZERO) THEN + IF (ELMBDA /= ZERO) THEN + IERROR = 6 + ELSE + GO TO (127,133,133,127,133) MP + 127 CONTINUE + GO TO (128,133,133,128,133) NP +! +! FOR SINGULAR PROBLEMS MUST ADJUST DATA TO INSURE THAT A SOLUTION +! WILL EXIST. +! + 128 CONTINUE + S = ZERO + DO J = 1, N + S = S + SUM(matF(:M,J)) + END DO + PERTRB = S/DBLE(M*N) + matF(:M,:N) = matF(:M,:N) - PERTRB + PERTRB = PERTRB/DELYSQ +! +! SOLVE THE EQUATION. +! + ENDIF + ENDIF + 133 CONTINUE + IERR1 = 0 +! Workarray splitted into parts of length M and IW + IW2 = size(vecW)-ID4 + IF (NPEROD /= 0) THEN + CALL POISTGG (NPEROD, N, MPEROD, M, vecW(1:M), vecW(ID2+1:ID2+M), & + vecW(ID3+1:ID3+M),IDIMF, matF, IERR1, vecW(ID4+1:),IW2) + ELSE + CALL GENBUNN (NPEROD, N, MPEROD, M, vecW(1:M), vecW(ID2+1:ID2+M), & + vecW(ID3+1:ID3+M), IDIMF, matF, IERR1, vecW(ID4+1:),IW2) + ENDIF +! +! REVISION HISTORY--- +! +! SEPTEMBER 1973 VERSION 1 +! APRIL 1976 VERSION 2 +! JANUARY 1978 VERSION 3 +! DECEMBER 1979 VERSION 3.1 +! FEBRUARY 1985 DOCUMENTATION UPGRADE +! NOVEMBER 1988 VERSION 3.2, FORTRAN 77 CHANGES +! June 2004 Version 5.0, Fortran 90 changes +!----------------------------------------------------------------------- + END SUBROUTINE HSTCRTT + + +! file hstcyl.f +! +! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . +! . . +! . copyright (c) 2004 by UCAR . +! . . +! . UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH . +! . . +! . all rights reserved . +! . . +! . . +! . FISHPACK version 5.0 . +! . . +! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . +! +! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +! * * +! * F I S H P A C K * +! * * +! * * +! * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * +! * * +! * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * +! * * +! * (Version 5.0 , JUNE 2004) * +! * * +! * BY * +! * * +! * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * +! * * +! * OF * +! * * +! * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * +! * * +! * BOULDER, COLORADO (80307) U.S.A. * +! * * +! * WHICH IS SPONSORED BY * +! * * +! * THE NATIONAL SCIENCE FOUNDATION * +! * * +! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +! +! SUBROUTINE HSTCYL (A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD, +! + ELMBDA,F,IDIMF,PERTRB,IERROR) +! +! DIMENSION OF BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N) +! ARGUMENTS +! +! LATEST REVISION June 2004 +! +! PURPOSE SOLVES THE STANDARD FIVE-POINT FINITE +! DIFFERENCE APPROXIMATION ON A STAGGERED +! GRID TO THE MODIFIED HELMHOLTZ EQUATION +! IN CYLINDRICAL COORDINATES. THIS EQUATION +! +! (1/R)(D/DR)(R(DU/DR)) + (D/DZ)(DU/DZ) +! +! + LAMBDA*(1/R**2)*U = F(R,Z) +! +! IS A TWO-DIMENSIONAL MODIFIED HELMHOLTZ +! EQUATION RESULTING FROM THE FOURIER TRANSFORM +! OF A THREE-DIMENSIONAL POISSON EQUATION. +! +! USAGE CALL HSTCYL (A,B,M,MBDCND,BDA,BDB,C,D,N, +! NBDCND,BDC,BDD,ELMBDA,F,IDIMF, +! PERTRB,IERROR) +! +! ARGUMENTS +! ON INPUT A,B +! +! THE RANGE OF R, I.E. A .LE. R .LE. B. +! A MUST BE LESS THAN B AND A MUST BE +! BE NON-NEGATIVE. +! +! M +! THE NUMBER OF GRID POINTS IN THE INTERVAL +! (A,B). THE GRID POINTS IN THE R-DIRECTION +! R-DIRECTION ARE GIVEN BY +! R(I) = A + (I-0.5)DR FOR I=1,2,...,M +! WHERE DR =(B-A)/M. +! M MUST BE GREATER THAN 2. +! +! MBDCND +! INDICATES THE TYPE OF BOUNDARY CONDITIONS +! AT R = A AND R = B. +! +! = 1 IF THE SOLUTION IS SPECIFIED AT R = A +! (SEE NOTE BELOW) AND R = B. +! +! = 2 IF THE SOLUTION IS SPECIFIED AT R = A +! (SEE NOTE BELOW) AND THE DERIVATIVE +! OF THE SOLUTION WITH RESPECT TO R IS +! SPECIFIED AT R = B. +! +! = 3 IF THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO R IS SPECIFIED AT +! R = A (SEE NOTE BELOW) AND R = B. +! +! = 4 IF THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO R IS SPECIFIED AT +! R = A (SEE NOTE BELOW) AND THE +! SOLUTION IS SPECIFIED AT R = B. +! +! = 5 IF THE SOLUTION IS UNSPECIFIED AT +! R = A = 0 AND THE SOLUTION IS +! SPECIFIED AT R = B. +! +! = 6 IF THE SOLUTION IS UNSPECIFIED AT +! R = A = 0 AND THE DERIVATIVE OF THE +! SOLUTION WITH RESPECT TO R IS SPECIFIED +! AT R = B. +! +! NOTE: +! IF A = 0, DO NOT USE MBDCND = 1,2,3, OR 4, +! BUT INSTEAD USE MBDCND = 5 OR 6. +! THE RESULTING APPROXIMATION GIVES THE ONLY +! MEANINGFUL BOUNDARY CONDITION, +! I.E. DU/DR = 0. +! (SEE D. GREENSPAN, 'INTRODUCTORY NUMERICAL +! ANALYSIS OF ELLIPTIC BOUNDARY VALUE +! PROBLEMS,' HARPER AND ROW, 1965, CHAPTER 5.) +! +! BDA +! A ONE-DIMENSIONAL ARRAY OF LENGTH N THAT +! SPECIFIES THE BOUNDARY VALUES (IF ANY) +! OF THE SOLUTION AT R = A. +! +! WHEN MBDCND = 1 OR 2, +! BDA(J) = U(A,Z(J)) , J=1,2,...,N. +! +! WHEN MBDCND = 3 OR 4, +! BDA(J) = (D/DR)U(A,Z(J)) , J=1,2,...,N. +! +! WHEN MBDCND = 5 OR 6, BDA IS A DUMMY +! VARIABLE. +! +! BDB +! A ONE-DIMENSIONAL ARRAY OF LENGTH N THAT +! SPECIFIES THE BOUNDARY VALUES OF THE +! SOLUTION AT R = B. +! +! WHEN MBDCND = 1,4,OR 5, +! BDB(J) = U(B,Z(J)) , J=1,2,...,N. +! +! WHEN MBDCND = 2,3, OR 6, +! BDB(J) = (D/DR)U(B,Z(J)) , J=1,2,...,N. +! +! C,D +! THE RANGE OF Z, I.E. C .LE. Z .LE. D. +! C MUST BE LESS THAN D. +! +! N +! THE NUMBER OF UNKNOWNS IN THE INTERVAL +! (C,D). THE UNKNOWNS IN THE Z-DIRECTION +! ARE GIVEN BY Z(J) = C + (J-0.5)DZ, +! J=1,2,...,N, WHERE DZ = (D-C)/N. +! N MUST BE GREATER THAN 2. +! +! NBDCND +! INDICATES THE TYPE OF BOUNDARY CONDITIONS +! AT Z = C AND Z = D. +! +! = 0 IF THE SOLUTION IS PERIODIC IN Z, I.E. +! U(I,J) = U(I,N+J). +! +! = 1 IF THE SOLUTION IS SPECIFIED AT Z = C +! AND Z = D. +! +! = 2 IF THE SOLUTION IS SPECIFIED AT Z = C +! AND THE DERIVATIVE OF THE SOLUTION WITH +! RESPECT TO Z IS SPECIFIED AT Z = D. +! +! = 3 IF THE DERIVATIVE OF THE SOLUTION WITH +! RESPECT TO Z IS SPECIFIED AT Z = C +! AND Z = D. +! +! = 4 IF THE DERIVATIVE OF THE SOLUTION WITH +! RESPECT TO Z IS SPECIFIED AT Z = C AND +! THE SOLUTION IS SPECIFIED AT Z = D. +! +! BDC +! A ONE DIMENSIONAL ARRAY OF LENGTH M THAT +! SPECIFIES THE BOUNDARY VALUES OF THE +! SOLUTION AT Z = C. +! +! WHEN NBDCND = 1 OR 2, +! BDC(I) = U(R(I),C) , I=1,2,...,M. +! +! WHEN NBDCND = 3 OR 4, +! BDC(I) = (D/DZ)U(R(I),C), I=1,2,...,M. +! +! WHEN NBDCND = 0, BDC IS A DUMMY VARIABLE. +! +! BDD +! A ONE-DIMENSIONAL ARRAY OF LENGTH M THAT +! SPECIFIES THE BOUNDARY VALUES OF THE +! SOLUTION AT Z = D. +! +! WHEN NBDCND = 1 OR 4, +! BDD(I) = U(R(I),D) , I=1,2,...,M. +! +! WHEN NBDCND = 2 OR 3, +! BDD(I) = (D/DZ)U(R(I),D) , I=1,2,...,M. +! +! WHEN NBDCND = 0, BDD IS A DUMMY VARIABLE. +! +! ELMBDA +! THE CONSTANT LAMBDA IN THE MODIFIED +! HELMHOLTZ EQUATION. IF LAMBDA IS GREATER +! THAN 0, A SOLUTION MAY NOT EXIST. +! HOWEVER, HSTCYL WILL ATTEMPT TO FIND A +! SOLUTION. LAMBDA MUST BE ZERO WHEN +! MBDCND = 5 OR 6. +! +! F +! A TWO-DIMENSIONAL ARRAY THAT SPECIFIES +! THE VALUES OF THE RIGHT SIDE OF THE +! MODIFIED HELMHOLTZ EQUATION. +! FOR I=1,2,...,M AND J=1,2,...,N +! F(I,J) = F(R(I),Z(J)) . +! F MUST BE DIMENSIONED AT LEAST M X N. +! +! IDIMF +! THE ROW (OR FIRST) DIMENSION OF THE ARRAY +! F AS IT APPEARS IN THE PROGRAM CALLING +! HSTCYL. THIS PARAMETER IS USED TO SPECIFY +! THE VARIABLE DIMENSION OF F. IDIMF MUST +! BE AT LEAST M. +! +! ON OUTPUT +! +! F +! CONTAINS THE SOLUTION U(I,J) OF THE FINITE +! DIFFERENCE APPROXIMATION FOR THE GRID POINT +! (R(I),Z(J)) FOR I=1,2,...,M, J=1,2,...,N. +! +! PERTRB +! IF A COMBINATION OF PERIODIC, DERIVATIVE, +! OR UNSPECIFIED BOUNDARY CONDITIONS IS +! SPECIFIED FOR A POISSON EQUATION +! (LAMBDA = 0), A SOLUTION MAY NOT EXIST. +! PERTRB IS A CONSTANT, CALCULATED AND +! SUBTRACTED FROM F, WHICH ENSURES THAT A +! SOLUTION EXISTS. HSTCYL THEN COMPUTES +! THIS SOLUTION, WHICH IS A LEAST SQUARES +! SOLUTION TO THE ORIGINAL APPROXIMATION. +! THIS SOLUTION PLUS ANY CONSTANT IS ALSO +! A SOLUTION; HENCE, THE SOLUTION IS NOT +! UNIQUE. THE VALUE OF PERTRB SHOULD BE +! SMALL COMPARED TO THE RIGHT SIDE F. +! OTHERWISE, A SOLUTION IS OBTAINED TO AN +! ESSENTIALLY DIFFERENT PROBLEM. +! THIS COMPARISON SHOULD ALWAYS BE MADE TO +! INSURE THAT A MEANINGFUL SOLUTION HAS BEEN +! OBTAINED. +! +! IERROR +! AN ERROR FLAG THAT INDICATES INVALID INPUT +! PARAMETERS. EXCEPT TO NUMBERS 0 AND 11, +! A SOLUTION IS NOT ATTEMPTED. +! +! = 0 NO ERROR +! +! = 1 A .LT. 0 +! +! = 2 A .GE. B +! +! = 3 MBDCND .LT. 1 OR MBDCND .GT. 6 +! +! = 4 C .GE. D +! +! = 5 N .LE. 2 +! +! = 6 NBDCND .LT. 0 OR NBDCND .GT. 4 +! +! = 7 A = 0 AND MBDCND = 1,2,3, OR 4 +! +! = 8 A .GT. 0 AND MBDCND .GE. 5 +! +! = 9 M .LE. 2 +! +! = 10 IDIMF .LT. M +! +! = 11 LAMBDA .GT. 0 +! +! = 12 A=0, MBDCND .GE. 5, ELMBDA .NE. 0 +! +! SINCE THIS IS THE ONLY MEANS OF INDICATING +! A POSSIBLY INCORRECT CALL TO HSTCYL, THE +! USER SHOULD TEST IERROR AFTER THE CALL. +! +! = 20 If the dynamic allocation of real and +! complex work space required for solution +! fails (for example if N,M are too large +! for your computer) +! I/O NONE +! +! PRECISION SINGLE +! +! REQUIRED LIBRARY fish.f,comf.f,genbun.f,gnbnaux.f,poistg.f +! FILES +! +! LANGUAGE FORTRAN 90 +! +! HISTORY WRITTEN BY ROLAND SWEET AT NCAR IN 1977. +! RELEASED ON NCAR'S PUBLIC SOFTWARE LIBRARIES +! IN JANUARY 1980. +! Revised in June 2004 by John Adams using +! Fortran 90 dynamically allocated work space. +! +! PORTABILITY FORTRAN 90 +! +! ALGORITHM THIS SUBROUTINE DEFINES THE FINITE-DIFFERENCE +! EQUATIONS, INCORPORATES BOUNDARY DATA, ADJUSTS +! THE RIGHT SIDE WHEN THE SYSTEM IS SINGULAR AND +! CALLS EITHER POISTG OR GENBUN WHICH SOLVES THE +! LINEAR SYSTEM OF EQUATIONS. +! +! TIMING FOR LARGE M AND N, THE OPERATION COUNT +! IS ROUGHLY PROPORTIONAL TO M*N*LOG2(N). +! +! ACCURACY THE SOLUTION PROCESS RESULTS IN A LOSS +! OF NO MORE THAN FOUR SIGNIFICANT DIGITS +! FOR N AND M AS LARGE AS 64. +! MORE DETAILED INFORMATION ABOUT ACCURACY +! CAN BE FOUND IN THE DOCUMENTATION FOR +! SUBROUTINE POISTG WHICH IS THE ROUTINE THAT +! ACTUALLY SOLVES THE FINITE DIFFERENCE +! EQUATIONS. +! +! REFERENCES U. SCHUMANN AND R. SWEET, "A DIRECT METHOD FOR +! THE SOLUTION OF POISSON'S EQUATION WITH NEUMANN +! BOUNDARY CONDITIONS ON A STAGGERED GRID OF +! ARBITRARY SIZE," J. COMP. PHYS. 20(1976), +! PP. 171-182. +!*********************************************************************** + SUBROUTINE HSTCYL(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, & + NBDCND, vecBDC, vecBDD, ELMBDA, matF, IDIMF, & + PERTRB, IERROR) + USE fish + implicit none + + DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & + ONE = 1.0D0, TWO = 2.0D0, & + FOUR = 4.0D0 + +!----------------------------------------------- +! D u m m y A r g u m e n t s +!----------------------------------------------- + INTEGER, INTENT(IN) :: M, MBDCND, N, NBDCND, IDIMF + INTEGER, INTENT(OUT) :: IERROR + DOUBLE PRECISION, INTENT(IN) :: A, B, C, D, ELMBDA + DOUBLE PRECISION, INTENT(OUT) :: PERTRB + DOUBLE PRECISION, DIMENSION(:), INTENT(IN) :: vecBDA, vecBDB, vecBDC, & + vecBDD + DOUBLE PRECISION, DIMENSION(IDIMF,N), INTENT(INOUT) :: matF +!----------------------------------------------- +! Allocatable arrays +!----------------------------------------------- + DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: work +!----------------------------------------------- +! L o c a l V a r i a b l e s +!----------------------------------------------- + INTEGER :: irwk, istatus +!----------------------------------------------- + ! JANNIS: add interface + interface + SUBROUTINE HSTCYLL(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, NBDCND, & + vecBDC, vecBDD, ELMBDA, matF, IDIMF, PERTRB, IERROR, W) + INTEGER, INTENT(IN) :: M, MBDCND, N, NBDCND, IDIMF + INTEGER, INTENT(OUT) :: IERROR + DOUBLE PRECISION, INTENT(IN) :: A, B, C, D, ELMBDA + DOUBLE PRECISION, INTENT(OUT) :: PERTRB + DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecBDA, vecBDB,vecBDC, & + vecBDD + DOUBLE PRECISION, DIMENSION(IDIMF,N), INTENT(INOUT) :: matF + DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: W + end subroutine hstcyll + end interface + + IERROR = 0 + IF (A < ZERO) IERROR = 1 + IF (A >= B) IERROR = 2 + IF (MBDCND<=0 .OR. MBDCND>=7) IERROR = 3 + IF (C >= D) IERROR = 4 + IF (N <= 2) IERROR = 5 + IF (NBDCND<0 .OR. NBDCND>=5) IERROR = 6 + IF (A==ZERO .AND. MBDCND/=5 .AND. MBDCND/=6) IERROR = 7 + IF (A>ZERO .AND. MBDCND>=5) IERROR = 8 + IF (IDIMF < M) IERROR = 10 + IF (M <= 2) IERROR = 9 + IF (A==ZERO .AND. MBDCND>=5 .AND. ELMBDA/=ZERO) IERROR = 12 + IF (IERROR /= 0) RETURN +! allocate real work space +! compute and allocate required real work space + CALL GEN_SPACE (N, M, IRWK) + IRWK = IRWK + 3*M + allocate(work(irwK),STAT=istatus) +!@! write(*,*) 'HSTCYL: allocate work(irwk); irwk=',irwk +! return if allocation failed (e.g., if n,m are too large) + IF (istatus > 0) THEN + write(*,*) 'HSTCYL: error allocate work(irwk); irwk=',irwk + RETURN + END IF +! check that allocation was successful + call HSTCYLL(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, NBDCND, & + vecBDC, vecBDD, ELMBDA, matF, IDIMF, PERTRB, IERROR, work) +! release allocated work space +!@! write(*,*) 'HSTCYL: deallocate work(irwk); irwk=',irwk + deallocate(work,STAT=istatus) + IF (istatus > 0) THEN + write(*,*) 'HSTCYL: error deallocate work' + RETURN + END IF +! +! RETURN + END SUBROUTINE HSTCYL + + SUBROUTINE HSTCYLL(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, NBDCND, & + vecBDC, vecBDD, ELMBDA, matF, IDIMF, PERTRB, IERROR, W) + +! USE poisson +! USE genbunal + + implicit none + + DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & + ONE = 1.0D0, TWO = 2.0D0, & + FOUR = 4.0D0 + +!----------------------------------------------- +! D u m m y A r g u m e n t s +!----------------------------------------------- + INTEGER, INTENT(IN) :: M, MBDCND, N, NBDCND, IDIMF + INTEGER, INTENT(OUT) :: IERROR + DOUBLE PRECISION, INTENT(IN) :: A, B, C, D, ELMBDA + DOUBLE PRECISION, INTENT(OUT) :: PERTRB + DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecBDA, vecBDB,vecBDC, & + vecBDD + DOUBLE PRECISION, DIMENSION(IDIMF,N), INTENT(INOUT) :: matF + DOUBLE PRECISION, DIMENSION(:), INTENT(OUT) :: W +!----------------------------------------------- +! L o c a l V a r i a b l e s +!----------------------------------------------- + INTEGER :: NP, IWB, IWC, IWR, I, J, K, LP, IERR1 + DOUBLE PRECISION :: DELTAR, DLRSQ, DELTHT, DLTHSQ, A1 +!----------------------------------------------- + DELTAR = (B - A)/DBLE(M) + DLRSQ = DELTAR**2 + DELTHT = (D - C)/DBLE(N) + DLTHSQ = DELTHT**2 + NP = NBDCND + 1 +! +! DEFINE A,B,C COEFFICIENTS IN W-ARRAY. +! + IWB = M + IWC = IWB + M + IWR = IWC + M + DO I = 1, M + J = IWR + I + W(J) = A + (DBLE(I) - HALF)*DELTAR + W(I) = (A + DBLE(I - 1)*DELTAR)/(DLRSQ*W(J)) + K = IWC + I + W(K) = (A + DBLE(I)*DELTAR)/(DLRSQ*W(J)) + K = IWB + I + W(K) = ELMBDA/W(J)**2 - TWO/DLRSQ + END DO +! +! ENTER BOUNDARY DATA FOR R-BOUNDARIES. +! + GO TO (102,102,104,104,106,106) MBDCND + 102 CONTINUE + A1 = TWO*W(1) + W(IWB+1) = W(IWB+1) - W(1) + matF(1,:N) = matF(1,:N) - A1*vecBDA(:N) + GO TO 106 + 104 CONTINUE + A1 = DELTAR*W(1) + W(IWB+1) = W(IWB+1) + W(1) + matF(1,:N) = matF(1,:N) + A1*vecBDA(:N) + 106 CONTINUE + GO TO (107,109,109,107,107,109) MBDCND + 107 CONTINUE + W(IWC) = W(IWC) - W(IWR) + A1 = TWO*W(IWR) + matF(M,:N) = matF(M,:N) - A1*vecBDB(:N) + GO TO 111 + 109 CONTINUE + W(IWC) = W(IWC) + W(IWR) + A1 = DELTAR*W(IWR) + matF(M,:N) = matF(M,:N) - A1*vecBDB(:N) +! +! ENTER BOUNDARY DATA FOR THETA-BOUNDARIES. +! + 111 CONTINUE + A1 = TWO/DLTHSQ + GO TO (121,112,112,114,114) NP + 112 CONTINUE + matF(:M,1) = matF(:M,1) - A1*vecBDC(:M) + GO TO 116 + 114 CONTINUE + A1 = ONE/DELTHT + matF(:M,1) = matF(:M,1) + A1*vecBDC(:M) + 116 CONTINUE + A1 = TWO/DLTHSQ + GO TO (121,117,119,119,117) NP + 117 CONTINUE + matF(:M,N) = matF(:M,N) - A1*vecBDD(:M) + GO TO 121 + 119 CONTINUE + A1 = ONE/DELTHT + matF(:M,N) = matF(:M,N) - A1*vecBDD(:M) + 121 CONTINUE + PERTRB = ZERO + IF (ELMBDA >= ZERO) THEN + IF (ELMBDA /= ZERO) THEN + IERROR = 11 + ELSE + GO TO (130,130,124,130,130,124) MBDCND + 124 CONTINUE + GO TO (125,130,130,125,130) NP + 125 CONTINUE + DO I = 1, M + A1 = ZERO + A1 = SUM(matF(I,:N)) + J = IWR + I + PERTRB = PERTRB + A1*W(J) + END DO + PERTRB = PERTRB/(DBLE(M*N)*HALF*(A + B)) + matF(:M,:N) = matF(:M,:N) - PERTRB + ENDIF + ENDIF + 130 CONTINUE + W(:M) = W(:M)*DLTHSQ + W(IWC+1:M+IWC) = W(IWC+1:M+IWC)*DLTHSQ + W(IWB+1:M+IWB) = W(IWB+1:M+IWB)*DLTHSQ + matF(:M,:N) = matF(:M,:N)*DLTHSQ + LP = NBDCND + W(1) = ZERO + W(IWR) = ZERO +! +! SOLVE THE SYSTEM OF EQUATIONS. +! + IERR1 = 0 + IF (NBDCND /= 0) THEN + CALL POISTGG (LP, N, 1, M, W(1:M), W(IWB+1:IWB+M), W(IWC+1:IWC+M), & + IDIMF, matF, IERR1, W(IWR+1:), SIZE(W)-IWR) + ELSE + CALL GENBUNN (LP, N, 1, M, W(1:M), W(IWB+1:IWB+M), W(IWC+1:IWC+M), & + IDIMF, matF, IERR1, W(IWR+1:), SIZE(W)-IWR) + ENDIF +! RETURN +! +! REVISION HISTORY--- +! +! SEPTEMBER 1973 VERSION 1 +! APRIL 1976 VERSION 2 +! JANUARY 1978 VERSION 3 +! DECEMBER 1979 VERSION 3.1 +! FEBRUARY 1985 DOCUMENTATION UPGRADE +! NOVEMBER 1988 VERSION 3.2, FORTRAN 77 CHANGES +! June 2004 Version 5.0, Fortran 90 changes +!----------------------------------------------------------------------- + END SUBROUTINE HSTCYLL +! +! file mhstcyl.f +! +! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . +! . . +! . copyright (c) 2004 by UCAR . +! . . +! . UNIVERSITY CORPORATION for ATMOSPHERIC RESEARCH . +! . . +! . all rights reserved . +! . . +! . . +! . FISHPACK version 5.0 . +! . . +! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . +! +! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +! * * +! * F I S H P A C K * +! * * +! * * +! * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * +! * * +! * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * +! * * +! * (Version 5.0 , JUNE 2004) * +! * * +! * BY * +! * * +! * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * +! * * +! * OF * +! * * +! * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * +! * * +! * BOULDER, COLORADO (80307) U.S.A. * +! * * +! * WHICH IS SPONSORED BY * +! * * +! * THE NATIONAL SCIENCE FOUNDATION * +! * * +! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * +! +! SUBROUTINE MHSTCYL (A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD, +! + ELMBDA,ES,F,IDIMF,PERTRB,IERROR) +! +! DIMENSION OF BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N) +! ARGUMENTS +! +! LATEST REVISION June 2004 +! +! PURPOSE SOLVES THE STANDARD FIVE-POINT FINITE +! DIFFERENCE APPROXIMATION ON A STAGGERED +! GRID TO THE MODIFIED HELMHOLTZ EQUATION +! IN CYLINDRICAL COORDINATES. THIS EQUATION +! +! (D2U/DR2) + (D/DZ)(DU/DZ) +! +! + LAMBDA*(1/R**2)*U + ES * U = F(R,Z) +! +! IS A TWO-DIMENSIONAL MODIFIED HELMHOLTZ +! EQUATION RESULTING FROM THE FOURIER TRANSFORM +! OF A THREE-DIMENSIONAL POISSON EQUATION. +! +! USAGE CALL MHSTCYL (A,B,M,MBDCND,BDA,BDB,C,D,N, +! NBDCND,BDC,BDD,ELMBDA,ES,F,IDIMF, +! PERTRB,IERROR) +! +! ARGUMENTS +! ON INPUT A,B +! +! THE RANGE OF R, I.E. A .LE. R .LE. B. +! A MUST BE LESS THAN B AND A MUST BE +! BE NON-NEGATIVE. +! +! M +! THE NUMBER OF GRID POINTS IN THE INTERVAL +! (A,B). THE GRID POINTS IN THE R-DIRECTION +! R-DIRECTION ARE GIVEN BY +! R(I) = A + (I-0.5)DR FOR I=1,2,...,M +! WHERE DR =(B-A)/M. +! M MUST BE GREATER THAN 2. +! +! MBDCND +! INDICATES THE TYPE OF BOUNDARY CONDITIONS +! AT R = A AND R = B. +! +! = 1 IF THE SOLUTION IS SPECIFIED AT R = A +! (SEE NOTE BELOW) AND R = B. +! +! = 2 IF THE SOLUTION IS SPECIFIED AT R = A +! (SEE NOTE BELOW) AND THE DERIVATIVE +! OF THE SOLUTION WITH RESPECT TO R IS +! SPECIFIED AT R = B. +! +! = 3 IF THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO R IS SPECIFIED AT +! R = A (SEE NOTE BELOW) AND R = B. +! +! = 4 IF THE DERIVATIVE OF THE SOLUTION +! WITH RESPECT TO R IS SPECIFIED AT +! R = A (SEE NOTE BELOW) AND THE +! SOLUTION IS SPECIFIED AT R = B. +! +! = 5 IF THE SOLUTION IS UNSPECIFIED AT +! R = A = 0 AND THE SOLUTION IS +! SPECIFIED AT R = B. +! +! = 6 IF THE SOLUTION IS UNSPECIFIED AT +! R = A = 0 AND THE DERIVATIVE OF THE +! SOLUTION WITH RESPECT TO R IS SPECIFIED +! AT R = B. +! +! NOTE: +! IF A = 0, DO NOT USE MBDCND = 1,2,3, OR 4, +! BUT INSTEAD USE MBDCND = 5 OR 6. +! THE RESULTING APPROXIMATION GIVES THE ONLY +! MEANINGFUL BOUNDARY CONDITION, +! I.E. DU/DR = 0. +! (SEE D. GREENSPAN, 'INTRODUCTORY NUMERICAL +! ANALYSIS OF ELLIPTIC BOUNDARY VALUE +! PROBLEMS,' HARPER AND ROW, 1965, CHAPTER 5.) +! +! BDA +! A ONE-DIMENSIONAL ARRAY OF LENGTH N THAT +! SPECIFIES THE BOUNDARY VALUES (IF ANY) +! OF THE SOLUTION AT R = A. +! +! WHEN MBDCND = 1 OR 2, +! BDA(J) = U(A,Z(J)) , J=1,2,...,N. +! +! WHEN MBDCND = 3 OR 4, +! BDA(J) = (D/DR)U(A,Z(J)) , J=1,2,...,N. +! +! WHEN MBDCND = 5 OR 6, BDA IS A DUMMY +! VARIABLE. +! +! BDB +! A ONE-DIMENSIONAL ARRAY OF LENGTH N THAT +! SPECIFIES THE BOUNDARY VALUES OF THE +! SOLUTION AT R = B. +! +! WHEN MBDCND = 1,4,OR 5, +! BDB(J) = U(B,Z(J)) , J=1,2,...,N. +! +! WHEN MBDCND = 2,3, OR 6, +! BDB(J) = (D/DR)U(B,Z(J)) , J=1,2,...,N. +! +! C,D +! THE RANGE OF Z, I.E. C .LE. Z .LE. D. +! C MUST BE LESS THAN D. +! +! N +! THE NUMBER OF UNKNOWNS IN THE INTERVAL +! (C,D). THE UNKNOWNS IN THE Z-DIRECTION +! ARE GIVEN BY Z(J) = C + (J-0.5)DZ, +! J=1,2,...,N, WHERE DZ = (D-C)/N. +! N MUST BE GREATER THAN 2. +! +! NBDCND +! INDICATES THE TYPE OF BOUNDARY CONDITIONS +! AT Z = C AND Z = D. +! +! = 0 IF THE SOLUTION IS PERIODIC IN Z, I.E. +! U(I,J) = U(I,N+J). +! +! = 1 IF THE SOLUTION IS SPECIFIED AT Z = C +! AND Z = D. +! +! = 2 IF THE SOLUTION IS SPECIFIED AT Z = C +! AND THE DERIVATIVE OF THE SOLUTION WITH +! RESPECT TO Z IS SPECIFIED AT Z = D. +! +! = 3 IF THE DERIVATIVE OF THE SOLUTION WITH +! RESPECT TO Z IS SPECIFIED AT Z = C +! AND Z = D. +! +! = 4 IF THE DERIVATIVE OF THE SOLUTION WITH +! RESPECT TO Z IS SPECIFIED AT Z = C AND +! THE SOLUTION IS SPECIFIED AT Z = D. +! +! BDC +! A ONE DIMENSIONAL ARRAY OF LENGTH M THAT +! SPECIFIES THE BOUNDARY VALUES OF THE +! SOLUTION AT Z = C. +! +! WHEN NBDCND = 1 OR 2, +! BDC(I) = U(R(I),C) , I=1,2,...,M. +! +! WHEN NBDCND = 3 OR 4, +! BDC(I) = (D/DZ)U(R(I),C), I=1,2,...,M. +! +! WHEN NBDCND = 0, BDC IS A DUMMY VARIABLE. +! +! BDD +! A ONE-DIMENSIONAL ARRAY OF LENGTH M THAT +! SPECIFIES THE BOUNDARY VALUES OF THE +! SOLUTION AT Z = D. +! +! WHEN NBDCND = 1 OR 4, +! BDD(I) = U(R(I),D) , I=1,2,...,M. +! +! WHEN NBDCND = 2 OR 3, +! BDD(I) = (D/DZ)U(R(I),D) , I=1,2,...,M. +! +! WHEN NBDCND = 0, BDD IS A DUMMY VARIABLE. +! +! ELMBDA +! THE CONSTANT LAMBDA IN THE MODIFIED +! HELMHOLTZ EQUATION. IF LAMBDA IS GREATER +! THAN 0, A SOLUTION MAY NOT EXIST. +! HOWEVER, MHSTCYL WILL ATTEMPT TO FIND A +! SOLUTION. LAMBDA MUST BE ZERO WHEN +! MBDCND = 5 OR 6. +! +! F +! A TWO-DIMENSIONAL ARRAY THAT SPECIFIES +! THE VALUES OF THE RIGHT SIDE OF THE +! MODIFIED HELMHOLTZ EQUATION. +! FOR I=1,2,...,M AND J=1,2,...,N +! F(I,J) = F(R(I),Z(J)) . +! F MUST BE DIMENSIONED AT LEAST M X N. +! +! IDIMF +! THE ROW (OR FIRST) DIMENSION OF THE ARRAY +! F AS IT APPEARS IN THE PROGRAM CALLING +! MHSTCYL. THIS PARAMETER IS USED TO SPECIFY +! THE VARIABLE DIMENSION OF F. IDIMF MUST +! BE AT LEAST M. +! +! ON OUTPUT +! +! F +! CONTAINS THE SOLUTION U(I,J) OF THE FINITE +! DIFFERENCE APPROXIMATION FOR THE GRID POINT +! (R(I),Z(J)) FOR I=1,2,...,M, J=1,2,...,N. +! +! PERTRB +! IF A COMBINATION OF PERIODIC, DERIVATIVE, +! OR UNSPECIFIED BOUNDARY CONDITIONS IS +! SPECIFIED FOR A POISSON EQUATION +! (LAMBDA = 0), A SOLUTION MAY NOT EXIST. +! PERTRB IS A CONSTANT, CALCULATED AND +! SUBTRACTED FROM F, WHICH ENSURES THAT A +! SOLUTION EXISTS. MHSTCYL THEN COMPUTES +! THIS SOLUTION, WHICH IS A LEAST SQUARES +! SOLUTION TO THE ORIGINAL APPROXIMATION. +! THIS SOLUTION PLUS ANY CONSTANT IS ALSO +! A SOLUTION; HENCE, THE SOLUTION IS NOT +! UNIQUE. THE VALUE OF PERTRB SHOULD BE +! SMALL COMPARED TO THE RIGHT SIDE F. +! OTHERWISE, A SOLUTION IS OBTAINED TO AN +! ESSENTIALLY DIFFERENT PROBLEM. +! THIS COMPARISON SHOULD ALWAYS BE MADE TO +! INSURE THAT A MEANINGFUL SOLUTION HAS BEEN +! OBTAINED. +! +! IERROR +! AN ERROR FLAG THAT INDICATES INVALID INPUT +! PARAMETERS. EXCEPT TO NUMBERS 0 AND 11, +! A SOLUTION IS NOT ATTEMPTED. +! +! = 0 NO ERROR +! +! = 1 A .LT. 0 +! +! = 2 A .GE. B +! +! = 3 MBDCND .LT. 1 OR MBDCND .GT. 6 +! +! = 4 C .GE. D +! +! = 5 N .LE. 2 +! +! = 6 NBDCND .LT. 0 OR NBDCND .GT. 4 +! +! = 7 A = 0 AND MBDCND = 1,2,3, OR 4 +! +! = 8 A .GT. 0 AND MBDCND .GE. 5 +! +! = 9 M .LE. 2 +! +! = 10 IDIMF .LT. M +! +! = 11 LAMBDA .GT. 0 +! +! = 12 A=0, MBDCND .GE. 5, ELMBDA .NE. 0 +! +! SINCE THIS IS THE ONLY MEANS OF INDICATING +! A POSSIBLY INCORRECT CALL TO MHSTCYL, THE +! USER SHOULD TEST IERROR AFTER THE CALL. +! +! = 20 If the dynamic allocation of real and +! complex work space required for solution +! fails (for example if N,M are too large +! for your computer) +! I/O NONE +! +! PRECISION SINGLE +! +! REQUIRED LIBRARY fish.f,comf.f,genbun.f,gnbnaux.f,poistg.f +! FILES +! +! LANGUAGE FORTRAN 90 +! +! HISTORY WRITTEN BY ROLAND SWEET AT NCAR IN 1977. +! RELEASED ON NCAR'S PUBLIC SOFTWARE LIBRARIES +! IN JANUARY 1980. +! Revised in June 2004 by John Adams using +! Fortran 90 dynamically allocated work space. +! +! PORTABILITY FORTRAN 90 +! +! ALGORITHM THIS SUBROUTINE DEFINES THE FINITE-DIFFERENCE +! EQUATIONS, INCORPORATES BOUNDARY DATA, ADJUSTS +! THE RIGHT SIDE WHEN THE SYSTEM IS SINGULAR AND +! CALLS EITHER POISTG OR GENBUN WHICH SOLVES THE +! LINEAR SYSTEM OF EQUATIONS. +! +! TIMING FOR LARGE M AND N, THE OPERATION COUNT +! IS ROUGHLY PROPORTIONAL TO M*N*LOG2(N). +! +! ACCURACY THE SOLUTION PROCESS RESULTS IN A LOSS +! OF NO MORE THAN FOUR SIGNIFICANT DIGITS +! FOR N AND M AS LARGE AS 64. +! MORE DETAILED INFORMATION ABOUT ACCURACY +! CAN BE FOUND IN THE DOCUMENTATION FOR +! SUBROUTINE POISTG WHICH IS THE ROUTINE THAT +! ACTUALLY SOLVES THE FINITE DIFFERENCE +! EQUATIONS. +! +! REFERENCES U. SCHUMANN AND R. SWEET, "A DIRECT METHOD FOR +! THE SOLUTION OF POISSON'S EQUATION WITH NEUMANN +! BOUNDARY CONDITIONS ON A STAGGERED GRID OF +! ARBITRARY SIZE," J. COMP. PHYS. 20(1976), +! PP. 171-182. +!*********************************************************************** + SUBROUTINE MHSTCYL(A, B, M, MBDCND, vecBDA, vecBDB, C, D, & + N, NBDCND, vecBDC, vecBDD, ELMBDA, ES, & + matF, IDIMF, PERTRB, IERROR) + USE fish + implicit none + + DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & + ONE = 1.0D0, TWO = 2.0D0, & + FOUR = 4.0D0 + +!----------------------------------------------- +! D u m m y A r g u m e n t s +!----------------------------------------------- + INTEGER, INTENT(IN) :: M, MBDCND, N, NBDCND, IDIMF + INTEGER, INTENT(OUT) :: IERROR + DOUBLE PRECISION, INTENT(IN) :: A, B, C, D, ELMBDA, ES + DOUBLE PRECISION, INTENT(OUT) :: PERTRB + DOUBLE PRECISION, DIMENSION(N), INTENT(IN) :: vecBDA, vecBDB + DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecBDC, vecBDD + DOUBLE PRECISION, DIMENSION(IDIMF,N), INTENT(OUT) :: matF +!----------------------------------------------- +! Allocatable arrays +!----------------------------------------------- + DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE :: work +!----------------------------------------------- +! L o c a l V a r i a b l e s +!----------------------------------------------- + INTEGER :: irwk, istatus +!----------------------------------------------- + IERROR = 0 + IF (A < ZERO) IERROR = 1 + IF (A >= B) IERROR = 2 + IF (MBDCND<=0 .OR. MBDCND>=7) IERROR = 3 + IF (C >= D) IERROR = 4 + IF (N <= 2) IERROR = 5 + IF (NBDCND<0 .OR. NBDCND>=5) IERROR = 6 +!Mod IF (A==ZERO .AND. MBDCND/=5 .AND. MBDCND/=6) IERROR = 7 + IF (A>ZERO .AND. MBDCND>=5) IERROR = 8 + IF (IDIMF < M) IERROR = 10 + IF (M <= 2) IERROR = 9 +!Mod IF (A==ZERO .AND. MBDCND>=5 .AND. ELMBDA/=ZERO) IERROR = 12 + IF (IERROR /= 0) RETURN +! allocate real work space +! compute and allocate required real work space + CALL GEN_SPACE (N, M, IRWK) + IRWK = IRWK + 3*M +!@! write(*,*) 'MHSTCYL: allocate work(irwk); irwk=',irwk + allocate(work(irwk),STAT=istatus) +! return if allocation failed (e.g., if n,m are too large) + IF (istatus > 0) THEN + write(*,*) 'HSTCRT: error allocate work(irwk); irwk=',irwk + RETURN + END IF + call MHSTCYLL(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, NBDCND, & + vecBDC, vecBDD, ELMBDA, ES, matF, IDIMF, PERTRB, & + IERROR, work, IRWK ) +! release allocated work space +!@! write(*,*) 'MHSTCYL: deallocate work(irwk); irwk=',irwk + deallocate(work,STAT=istatus) + IF (istatus > 0) THEN + write(*,*) 'MHSTCYL: error deallocate work' + RETURN + END IF +! + END SUBROUTINE MHSTCYL + + SUBROUTINE MHSTCYLL(A, B, M, MBDCND, vecBDA, vecBDB, C, D, N, NBDCND, & + vecBDC, vecBDD, ELMBDA, ES, matF, IDIMF, & + PERTRB, IERROR, W, IW) +! USE genbunal +! USE poisson + + implicit none + + DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, & + ONE = 1.0D0, TWO = 2.0D0, & + FOUR = 4.0D0 + +!----------------------------------------------- +! D u m m y A r g u m e n t s +!----------------------------------------------- + INTEGER, INTENT(IN) :: M, MBDCND, N, NBDCND, IDIMF, IW + INTEGER, INTENT(OUT) :: IERROR + DOUBLE PRECISION, INTENT(IN) :: A, B, C, D, ELMBDA, ES + DOUBLE PRECISION, INTENT(OUT) :: PERTRB + DOUBLE PRECISION, DIMENSION(N), INTENT(IN) :: vecBDA, vecBDB + DOUBLE PRECISION, DIMENSION(M), INTENT(IN) :: vecBDC, vecBDD + DOUBLE PRECISION, DIMENSION(IDIMF,N), INTENT(OUT) :: matF + DOUBLE PRECISION, DIMENSION(IW), INTENT(INOUT) :: W +!----------------------------------------------- +! L o c a l V a r i a b l e s +!----------------------------------------------- + INTEGER :: NP, IWB, IWC, IWR, I, J, K, LP, IERR1 + DOUBLE PRECISION :: DELTAR, DLRSQ, DELTHT, DLTHSQ, A1 +!----------------------------------------------- + DELTAR = (B - A)/DBLE(M) + DLRSQ = DELTAR**2 + DELTHT = (D - C)/DBLE(N) + DLTHSQ = DELTHT**2 + NP = NBDCND + 1 +! +! DEFINE A,B,C COEFFICIENTS IN W-ARRAY. +! + IWB = M + IWC = IWB + M + IWR = IWC + M + DO I = 1, M + J = IWR + I + W(J) = A + (DBLE(I) - HALF)*DELTAR + W(I) = (W(J) - HALF * DELTAR)/(DLRSQ*W(J)) + K = IWC + I + W(K) = (W(J) + HALF * DELTAR)/(DLRSQ*W(J)) + K = IWB + I + W(K) = ELMBDA/W(J)**2 + ES - TWO/DLRSQ + END DO +! +! ENTER BOUNDARY DATA FOR R-BOUNDARIES. +! + GO TO (102,102,104,104,106,106) MBDCND + 102 CONTINUE + A1 = TWO*W(1) + W(IWB+1) = W(IWB+1) - W(1) + matF(1,:N) = matF(1,:N) - A1*vecBDA(:N) + GO TO 106 + 104 CONTINUE + A1 = DELTAR*W(1) + W(IWB+1) = W(IWB+1) + W(1) + matF(1,:N) = matF(1,:N) + A1*vecBDA(:N) + 106 CONTINUE + GO TO (107,109,109,107,107,109) MBDCND + 107 CONTINUE + W(IWC) = W(IWC) - W(IWR) + A1 = TWO*W(IWR) + matF(M,:N) = matF(M,:N) - A1*vecBDB(:N) + GO TO 111 + 109 CONTINUE + W(IWC) = W(IWC) + W(IWR) + A1 = DELTAR*W(IWR) + matF(M,:N) = matF(M,:N) - A1*vecBDB(:N) +! +! ENTER BOUNDARY DATA FOR THETA-BOUNDARIES. +! + 111 CONTINUE + A1 = TWO/DLTHSQ + GO TO (121,112,112,114,114) NP + 112 CONTINUE + matF(:M,1) = matF(:M,1) - A1*vecBDC(:M) + GO TO 116 + 114 CONTINUE + A1 = ONE/DELTHT + matF(:M,1) = matF(:M,1) + A1*vecBDC(:M) + 116 CONTINUE + A1 = TWO/DLTHSQ + GO TO (121,117,119,119,117) NP + 117 CONTINUE + matF(:M,N) = matF(:M,N) - A1*vecBDD(:M) + GO TO 121 + 119 CONTINUE + A1 = ONE/DELTHT + matF(:M,N) = matF(:M,N) - A1*vecBDD(:M) + 121 CONTINUE + PERTRB = ZERO + IF (ELMBDA >= ZERO) THEN + IF (ELMBDA /= ZERO) THEN + IERROR = 11 + ELSE + GO TO (130,130,124,130,130,124) MBDCND + 124 CONTINUE + GO TO (125,130,130,125,130) NP + 125 CONTINUE + DO I = 1, M + A1 = ZERO + A1 = SUM(matF(I,:N)) + J = IWR + I + PERTRB = PERTRB + A1*W(J) + END DO + PERTRB = PERTRB/(DBLE(M*N)*HALF*(A + B)) + matF(:M,:N) = matF(:M,:N) - PERTRB + ENDIF + ENDIF + 130 CONTINUE + W(:M) = W(:M)*DLTHSQ + W(IWC+1:M+IWC) = W(IWC+1:M+IWC)*DLTHSQ + W(IWB+1:M+IWB) = W(IWB+1:M+IWB)*DLTHSQ + matF(:M,:N) = matF(:M,:N)*DLTHSQ + LP = NBDCND + W(1) = ZERO + W(IWR) = ZERO +! +! SOLVE THE SYSTEM OF EQUATIONS. +! + IERR1 = 0 + IF (NBDCND /= 0) THEN + CALL POISTGG (LP, N, 1, M, W, W(IWB+1:IWB+M), W(IWC+1:IWC+M), & + IDIMF, matF,IERR1, W(IWR+1:), size(W)-IWR) + ELSE + CALL GENBUNN (LP, N, 1, M, W, W(IWB+1:IWB+M), W(IWC+1:IWC+M), & + IDIMF, matF,IERR1, W(IWR+1:), size(W)-IWR) + ENDIF +! +! REVISION HISTORY--- +! +! SEPTEMBER 1973 VERSION 1 +! APRIL 1976 VERSION 2 +! JANUARY 1978 VERSION 3 +! DECEMBER 1979 VERSION 3.1 +! FEBRUARY 1985 DOCUMENTATION UPGRADE +! NOVEMBER 1988 VERSION 3.2, FORTRAN 77 CHANGES +! June 2004 Version 5.0, Fortran 90 changes +! February 2009 Modified version of HSTCYL +!----------------------------------------------------------------------- + END SUBROUTINE MHSTCYLL