!

!     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