```
!
!################################################################
!################################################################
!#####                                                      #####
!#####               SUBROUTINE RETRPOIS                    #####
!#####                                                      #####
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!#####                                                      #####
!################################################################
!################################################################
! SUBROUTINE retrpois(nx,ny,nz,k,                                         & 1
biga,bigb,pc,                                                &
tem4,pnew,d,alpha,beta,pprt)
!
!---------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine solves a 2-D Poisson equation:
!
!  d2p/dx2  +  d2p/dy2  =  dbiga/dx + dbigb/dy
!
!  subject to the Dirichlet conditions:
!
!  p = pbkg on the i = 2 and i = nx-1 physical boundaries
!
!  p = pbkg on the j = 2 and j = ny-1 physical boundaries
!
!
!  The solution is obtained by successively applying an alternating
!  direction implicit (ADI) scheme to solve tri-diagonal matrix
!  equations in the i and j directions.  This scheme is known as the
!  "sweeping method" or Thomas algorithm.  (see P. J. Roache, 1982,
!  "Computational Fluid Dynamics").
!
!---------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Alan Shapiro and Steve Lazarus
!  2/2/93.
!
!  MODIFICATION HISTORY:
!
!  6/18/96 (Limin Zhao and Alan Shapiro)
!  Added the option for Dirichlet boundary condition, which turns
!  out very good for properply retrieving pressure and potential
!  temperature outside of rainwater area.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the z-direction (vertical)
!
!    biga     Sum of terms in x-momentum equation (except for dp/dx)
!    bigb     Sum of terms in y-momentum equation (except for dp/dy)
!    pc       First guess for the solution of the Poisson equation.
!
!
!---------------------------------------------------------------------
!
!  OUTPUT:
!
!  pc        Solution of the Poisson equation.  pc differs from the
!            perturbation pressure by an arbitrary function of height.
!
!---------------------------------------------------------------------
!
!  WORK ARRAYS:
!
!    alpha    Coefficient in solution to a tri-diagonal matrix equation.
!    beta     Coefficient in solution to a tri-diagonal matrix equation.
!    tem4     Temporary work array.
!
!   (These arrays are defined and used locally (i.e. inside this
!    subroutine), they may also be passed into routines called by
!    this one. Exiting the call to this subroutine, these temporary
!    work arrays may be used for other purposes, and therefore their
!    contents may be overwritten. Please examine the usage of work
!    arrays before you alter the code.)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE     ! Force explicit declarations
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'     ! Model control constants
INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
!
INTEGER :: nx, ny, nz        ! Number of grid points in x, y and z directions

REAL :: biga(nx,ny,nz)       ! Sum of terms in x-momentum equation
! (except for dp/dx).

REAL :: bigb(nx,ny,nz)       ! Sum of terms in y-momentum equation
! (except for dp/dy)

REAL :: pc(nx,ny,nz)         ! Solution of the Poisson equation.
! pc differs from the perturbation pressure
! by an arbitrary function of height.
REAL :: tem4(nx,ny,nz)
REAL :: pprt(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
INTEGER :: i, j, k
INTEGER :: iter      ! Iteration number
INTEGER :: itmax     ! Maximum number of iterations
INTEGER :: icon      ! Used in error tolerance check;
! =1, iterates have not converged, continue iterating

INTEGER :: inum      ! Number of grid points on a horizontal plane

REAL :: rdenom       ! Reciprocal denominator, used to simplify a calculation
REAL :: prms         ! Root mean square value of pc along an i or j line
REAL :: tay          ! Over-relaxation coefficient
REAL :: tol          ! An error tolerance parameter
REAL :: resid        ! Discrepancy between updated  pc and old pc at a point
REAL :: compar       ! Compare resid with the quantity compar (=tol*prms)
REAL :: pave         ! Horizontal average of pc

PARAMETER(itmax=500, tol=0.00001)

REAL :: pnew(nx,ny)  ! Updated value of pc

REAL :: d(nx,ny)     ! Inhomogeneous term in a tri-diagonal matrix equation

REAL :: a, b, c      ! Coefficients in a tri-diagonal matrix equation

REAL :: alpha(nx*ny) ! Coefficient in recursive solution to a tri-diagonal
! matrix equation

REAL :: beta(nx*ny)  ! Coefficient in recursive solution to a tri-diagonal
! matrix equation

REAL :: rdx2         ! Reciprocal of (dx)**2
REAL :: rdy2         ! Reciprocal of (dy)**2

REAL :: bc_opt
INTEGER :: ips
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code ...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!----------------------------------------------------------------------
!
!  Estimate the optimum over-relaxation coefficient, tay, from
!  the formula (cf G.J. Haltiner and R.T. Williams, 1980: "Numerical
!  Prediction and Dynamic Meteorology, second edition":
!
!  tay = 2. - pi*sqrt(2/nx**2 + 2/ny**2)
!
!
!----------------------------------------------------------------------
!

tay = 2. - 3.14159*SQRT((2./(nx*nx) + 2./(ny*ny)))

PRINT *, 'estimated optimum over-relaxation coeff. =', tay

rdx2 = dxinv*dxinv        ! reciprocal of (dx)**2
rdy2 = dyinv*dyinv        ! reciprocal of (dy)**2

!
!-----------------------------------------------------------------------
!
!  Begin the big do-loop for iterating back and forth between the
!  solutions of the two tri-diagonal matrix equations.
!
!-----------------------------------------------------------------------
!
DO iter = 1, itmax

IF (iter == itmax) THEN
PRINT 2000, k
STOP 1
END IF

2000  FORMAT(/, 1X, 'Warning& ! Poisson solver did not converge.',  &
&        /, 1X, 'klevel =', i5, /)

IF (MOD(iter,20) == 0) PRINT 3000, k, iter

3000  FORMAT(2X, 'klevel =', i5, 5X, 'iter = ', i5)
!
!-----------------------------------------------------------------------
!
!  First sweep in the i direction (put in an outer j do loop since we will
!  sweep in the i-direction for each value of j).  We are thus solving
!  a tri-diagonal matrix equation of the form:
!
!          a*pnew(i+1,j) + b*pnew(i,j) + c*pnew(i-1,j) = d(i,j)
!
!    where a, b, c and d are of the form:
!
!    a = c = 1/(dx*dx),    b = -2*(1/(dx*dx) + 1/(dy*dy)),
!
!    d = tem4 - (pc(i,j+1) + pc(i,j-1))/(dy*dy)
!
!
!-----------------------------------------------------------------------
!
a = rdx2
c = a
b = -2.*(rdx2 + rdy2)

DO j = 2, ny-2

DO i = 2, nx-2
d(i,j) = tem4(i,j,k) - rdy2*(pc(i,j+1,k) + pc(i,j-1,k))
END DO
!
!-----------------------------------------------------------------------
!
!  The solution of the tri-diagonal matrix equation satisfies the
!  recursion relation:
!
!            pnew(i-1,j) = alpha(i)*pnew(i,j) + beta(i)
!
!  where alpha and beta satisfy the recursion formulae:
!
!            alpha(i) = -a(i)/(b + c*alpha(i-1));
!
!            beta(i) = (d(i-1,j) - c*beta(i-1))/(b + c*alpha(i-1))
!
!  Write the boundary values of the alpha and beta coefficients at
!  i=2 from the Neumann boundary condition on the west boundary.  Then
!  compute the alpha and beta coefficients along i-lines from the
!  recursion formulae.
!
!-----------------------------------------------------------------------
!
bc_opt = 0

IF (bc_opt == 1) THEN        !Neumann B.C.
alpha(2) = 1.
beta(2) = - biga(2,j,k)*dx
ELSE
alpha(2) = 0.0
beta(2)  = pprt(1,j,k)
END IF

DO i = 3, nx-1
rdenom = 1./(b + c*alpha(i-1))
alpha(i) = - a*rdenom
beta(i) = (d(i-1,j) - c*beta(i-1))*rdenom
END DO

!
!-----------------------------------------------------------------------
!
!  Compute the updated pc (pnew), at nx-1 from the Neumann boundary
!  condition at the east boundary.  Sweep downwards in the decreasing
!  i-direction to get the updated pc (pnew) from the recursion formula.
!
!-----------------------------------------------------------------------
!
IF (bc_opt == 1) THEN
pnew(nx-1,j) = (biga(nx-1,j,k)*dx                               &
+ beta(nx-1))/(1. - alpha(nx-1))
ELSE
pnew(nx-1,j) = pprt(nx-1,j,k)
END IF

DO i = nx-1, 2, -1
pnew(i-1,j) = alpha(i)*pnew(i,j) + beta(i)
END DO
!
!-----------------------------------------------------------------------
!
!  Compute the over-relaxed residuals and add them to pc.
!
!-----------------------------------------------------------------------
!

DO i = 1, nx-1
resid = tay*(pnew(i,j) - pc(i,j,k))
pc(i,j,k) = pc(i,j,k) + resid
END DO

END DO

!
!-----------------------------------------------------------------------
!
!  Now sweep in the j direction (put in an outer i-do loop since we will
!  sweep in the j-direction for each value of i).  We are thus solving
!  a tri-diagonal matrix equation of the form:
!
!          a*pnew(i,j+1) + b*pnew(i,j) + c*pnew(i,j-1) = d(i,j)
!
!    where a, b, c and d are of the form:
!
!    a = c = 1/(dy*dy),    b = -2*(1/(dx*dx) + 1/(dy*dy)),
!
!    d = tem4 - (pc(i+1,j) + pc(i-1,j))/(dx*dx)
!
!-----------------------------------------------------------------------
!
a = rdy2
c = a
b = -2.*(rdx2 + rdy2)

DO i = 2, nx-2

DO j = 2, ny-2
d(i,j) = tem4(i,j,k) - rdx2*(pc(i+1,j,k) + pc(i-1,j,k))
END DO
!
!-----------------------------------------------------------------------
!
!  The solution of the tri-diagonal matrix equation satisfies the
!  recursion relation:
!
!            pnew(i,j-1) = alpha(j)*pnew(i,j) + beta(j)
!
!  where alpha and beta satisfy the recursion formulae:
!
!            alpha(j) = -a(j)/(b + c*alpha(j-1));
!
!            beta(j) = (d(i,j-1) - c*beta(j-1))/(b + c*alpha(j-1))
!
!  Write the boundary values of the alpha and beta coefficients at
!  j=2 from the Neumann boundary condition on the south boundary.  Then
!  compute the alpha and beta coefficients along j-lines from the
!  recursion formulae.
!
!-----------------------------------------------------------------------
!
IF (bc_opt == 1) THEN
alpha(2) = 1.
beta(2) = - bigb(i,2,k)*dy
ELSE
alpha(2) = 0.0
beta(2)  = pprt(i,1,k)
END IF

DO j = 3, ny-1
rdenom = 1./(b + c*alpha(j-1))
alpha(j) = - a*rdenom
beta(j) = (d(i,j-1) - c*beta(j-1))*rdenom
END DO
!
!-----------------------------------------------------------------------
!
!  Compute the updated pc (pnew), at ny-1 from the Neumann boundary
!  condition at the north boundary.  Sweep downwards in the decreasing
!  j-direction to get the updated pc (pnew) from the recursion formula.
!
!-----------------------------------------------------------------------
!
IF (bc_opt == 1) THEN
pnew(i,ny-1) = (bigb(i,ny-1,k)*dy                               &
+ beta(ny-1))/(1. - alpha(ny-1))
ELSE
pnew(i,ny-1) = pprt(i,ny-1,k)
END IF

DO j = ny-1, 2, -1
pnew(i,j-1) = alpha(j)*pnew(i,j) + beta(j)
END DO
!
!-----------------------------------------------------------------------
!
!  Compute the rms value of the new pc along a j-line for the error
!  tolerance comparison.  Also compute the over-relaxed residuals and
!
!-----------------------------------------------------------------------
!
prms = 0.
DO j = 1, ny-1
prms = prms + pnew(i,j)*pnew(i,j)
END DO
prms = SQRT(prms/(ny-1))
compar = tol*prms

icon = 0
DO j = 1, ny-1
resid = tay*(pnew(i,j) - pc(i,j,k))
icon = icon + IFIX(ABS(resid/compar))
pc(i,j,k) = pc(i,j,k) + resid
END DO

END DO

IF (icon == 0) RETURN
!
!-----------------------------------------------------------------------
!
!  Since the solution of Poisson's equation with Neumann boundary
!  conditions isn't unique (there's an arbitrary constant in the
!  solution), subtract off the average value of pc to avoid a potential
!  drift that could affect the convergence.
!
!-----------------------------------------------------------------------
!
IF (bc_opt /= 1) CYCLE
WRITE(6,*)'subtract mean pressue'
pave = 0.
inum = 0
DO i = 1, nx-1
DO j = 1, ny-1
pave = pave + pc(i,j,k)
inum = inum + 1
END DO
END DO
pave = pave/FLOAT(inum)

DO i = 1, nx-1
DO j = 1, ny-1
pc(i,j,k) = pc(i,j,k) - pave
END DO
END DO

END DO

RETURN
END SUBROUTINE retrpois
```