!
!################################################################
!################################################################
!##### #####
!##### SUBROUTINE RETRPOIS #####
!##### #####
!##### Copyright (c) 1993 #####
!###### 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;
! =0, iterates have converged, RETURN to RETRPTPR;
! =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
! add them to pc.
!
!-----------------------------------------------------------------------
!
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