!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module POP_SolversMod 4,16
!BOP
! !MODULE: POP_SolversMod
!
! !DESCRIPTION:
! This module contains routines and operators for solving the elliptic
! system for surface pressure in the barotropic mode.
!
! !REVISION HISTORY:
! SVN:$Id: $
! !USES:
use POP_KindsMod
use POP_ErrorMod
use POP_CommMod
use POP_ConfigMod
use POP_IOUnitsMod
use POP_BlocksMod
use POP_DistributionMod
use POP_GridHorzMod
use POP_ReductionsMod
use POP_RedistributeMod
use POP_HaloMod
use POP_GridHorzMod
use POP_FieldMod
use POP_DomainSizeMod
use domain
use grid
implicit none
private
save
! !PUBLIC MEMBER FUNCTIONS:
public :: POP_SolversInit, &
POP_SolversRun, &
POP_SolversDiagonal, &
POP_SolversGetDiagnostics
! !PUBLIC DATA MEMBERS:
!-----------------------------------------------------------------------
!
! other operator and preconditioner weights for barotropic operator
!
!-----------------------------------------------------------------------
real (POP_r8), dimension (:,:,:), allocatable, public :: &
mMaskTropic, &! land mask in barotropic distribution
btropWgtCenter, &! barotropic operater center coefficient
btropWgtNorth, &! barotropic operater north coefficient
btropWgtEast, &! barotropic operater east coefficient
btropWgtNE ! barotropic operater northeast coefficient
real (POP_r8), dimension (:,:,:), allocatable, public :: &
centerWgtClinicIndep, &! time indep center wgt on clinic distrb
centerWgtClinic ! time depend center wgt on clinic distrb
!EOP
!BOC
!*** preconditioner operator coefficients (ninept operator)
real (POP_r8), dimension (:,:,:), allocatable :: &
precondCenter, &
precondNorth, precondSouth, &
precondEast, precondWest, &
precondNE, precondSE, &
precondNW, precondSW
!-----------------------------------------------------------------------
!
! supported solvers and preconditioners
!
!-----------------------------------------------------------------------
character (POP_charLength) :: &
solverChoice
character (3), parameter :: &
solverChoicePCG = 'pcg'
character (9), parameter :: &
solverChoiceChronGear = 'ChronGear'
character (POP_charLength) :: &
preconditionerChoice
character (8), parameter :: &
precondChoiceDiag = 'diagonal'
character (4), parameter :: &
precondChoiceFile = 'file'
logical (POP_logical) :: &
Preconditioner
!-----------------------------------------------------------------------
!
! scalar convergence-related variables
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
maxIterations, &! max number of solver iterations
convergenceCheckFreq ! check convergence every freq steps
real (POP_r8) :: &
convergenceCriterion, &! convergence error criterion
residualNorm ! residual normalization
!*** convergence diagnostics
integer (POP_i4), public :: &
numIterations ! accumulated no of iterations (diagnostic)
real (POP_r8) :: &
rmsResidual ! residual (also a diagnostic)
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: POP_SolversRun
! !INTERFACE:
subroutine POP_SolversRun(sfcPressure, rhsClinic, errorCode) 1,14
! !DESCRIPTION:
! Solves the elliptic equation for surface pressure by calling
! the requested solver routine. Also redistributes necessary
! array to the barotropic distribution of blocks for better performance
! of the solver.
! The elliptic equation is
! \begin{equation}
! AF = B
! \end{equation}
! where $F$ is a field (eg surface pressure), $B$ is the right hand side
! and $A$ is the operator defined as
! \begin{equation}
! AF = a \nabla\cdot(H \nabla F)
! \end{equation}
! where $a$ is the cell area.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(in) :: &
rhsClinic ! right-hand-side of linear system
! for blocks in baroclinic distribution
! !INPUT/OUTPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(inout) :: &
sfcPressure ! on input, initial guess in baroclinic distrb
! on output, final solution for sfc pressure
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
numProcs ! number of processors in barotropic distrib
real (POP_r8), dimension(size(sfcPressure,dim=1), &
size(sfcPressure,dim=2), &
size(sfcPressure,dim=3)) :: &
pressTropic, &! surface pressure in barotropic distribution
rhsTropic ! right hand side in barotropic distribution
!-----------------------------------------------------------------------
!
! switch to the barotropic distribution for iterative solvers
!
!-----------------------------------------------------------------------
errorCode = POP_Success
call POP_RedistributeBlocks(btropWgtCenter, POP_distrbTropic, &
centerWgtClinic, POP_distrbClinic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversRun: error redistributing center operator weight')
return
endif
call POP_RedistributeBlocks(pressTropic, POP_distrbTropic, &
sfcPressure, POP_distrbClinic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversRun: error redistributing pressure')
return
endif
call POP_RedistributeBlocks(rhsTropic, POP_distrbTropic, &
rhsClinic, POP_distrbClinic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversRun: error redistributing right hand side')
return
endif
!-----------------------------------------------------------------------
!
! call proper routine based on user choice of solver
!
!-----------------------------------------------------------------------
call POP_DistributionGet(POP_distrbTropic, errorCode, &
numProcs = numProcs)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversRun: error getting num procs')
return
endif
if (POP_myTask < numProcs) then
select case(trim(solverChoice))
case (solverChoicePCG)
call pcg(pressTropic, rhsTropic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolverRun: error in PCG')
return
endif
case (solverChoiceChronGear)
call ChronGear(pressTropic, rhsTropic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolverRun: error in ChronGear')
return
endif
end select
endif
!-----------------------------------------------------------------------
!
! switch solution back to the baroclinic distribution
!
!-----------------------------------------------------------------------
call POP_RedistributeBlocks(sfcPressure, POP_distrbClinic, &
pressTropic, POP_distrbTropic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversRun: error redistributing pressure back')
return
endif
!-----------------------------------------------------------------------
!EOC
end subroutine POP_SolversRun
!***********************************************************************
!BOP
! !IROUTINE: POP_SolversInit
! !INTERFACE:
subroutine POP_SolversInit(errorCode) 1,37
! !DESCRIPTION:
! This routine initializes choice of solver, calculates the
! coefficients of the 9-point stencils for the barotropic operator and
! reads in a preconditioner if requested.
!
! !REVISION HISTORY:
! same as module
!
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables:
!
! {X,Y}{NE,SE,NW,SW} = contribution to {ne,se,nw,sw} coefficients
! from {x,y} components of divergence
! HU = depth at U points
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
i,j, &! dummy loop counters
configUnit, &! unit for configuration file
numBlocksClinic, &! num local blocks in baroclinic distribution
numBlocksTropic, &! num local blocks in barotropic distribution
iblock, &! block counter
istat ! status flag for allocates
character (POP_charLength) :: &
preconditionerFile ! file containing preconditioner
real (POP_r8) :: &
xne,xse,xnw,xsw, &! contribution to coefficients from x,y
yne,yse,ynw,ysw, &! components of divergence
ase,anw,asw
real (POP_r8), dimension(:,:,:), allocatable :: &
work0, &! temp space for computing barotropic
workNorth, &!
workEast, &!
workNE, &!
mMaskTmp ! land mask in barotropic distribution
!-----------------------------------------------------------------------
!
! read solver choice and solver constants from configuration file
!
!-----------------------------------------------------------------------
errorCode = POP_Success
if (POP_myTask == POP_masterTask) then
write(POP_stdout,POP_blankFormat)
write(POP_stdout,POP_delimFormat)
write(POP_stdout,'(a35)') ' Solver options (barotropic solver)'
write(POP_stdout,POP_delimFormat)
write(POP_stdout,POP_blankFormat)
endif
call POP_ConfigOpen(configUnit, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error opening config file')
return
endif
call POP_ConfigRead(configUnit, 'solvers', 'solverChoice', &
solverChoice, solverChoicePCG, errorCode, &
outStringBefore = 'Solver choice: ')
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error reading solver choice')
return
endif
call POP_ConfigRead(configUnit, 'solvers', 'convergenceCriterion', &
convergenceCriterion, 1.e-12_POP_r8, errorCode, &
outStringBefore = 'Solver converged for err < ')
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error reading solver convergence criterion')
return
endif
call POP_ConfigRead(configUnit, 'solvers', 'maxIterations', &
maxIterations, 1000, errorCode, &
outStringBefore = 'Solver maximum iterations: ')
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error reading solver max iterations')
return
endif
call POP_ConfigRead(configUnit, 'solvers', 'convergenceCheckFreq', &
convergenceCheckFreq, 10, errorCode, &
outStringBefore = 'Check convergence every ', &
outStringAfter = ' iterations')
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error reading convergence check frequency')
return
endif
call POP_ConfigRead(configUnit, 'solvers', 'preconditionerChoice', &
preconditionerChoice, precondChoiceDiag, &
errorCode, &
outStringBefore = 'Preconditioner choice: ')
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error reading preconditioner choice')
return
endif
if (trim(preconditionerChoice) == precondChoiceFile) then
call POP_ConfigRead(configUnit, 'solvers', 'preconditionerFile', &
preconditionerFile, 'UnknownPrecondFile', errorCode, &
outStringBefore = 'Reading preconditioner from file: ')
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error reading preconditioner file name')
return
endif
endif
call POP_ConfigClose(configUnit, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error closing config file')
return
endif
!-----------------------------------------------------------------------
!
! check inputs
!
!-----------------------------------------------------------------------
select case(trim(solverChoice))
case(solverChoicePCG)
case(solverChoiceChronGear)
case default
call POP_ErrorSet(errorCode, &
'POP_SolversInit: unknown solver - must be pcg, ChronGear')
return
end select
select case (trim(preconditionerChoice))
case(precondChoiceDiag)
usePreconditioner = .false. ! default is diagonal
case(precondChoiceFile)
usePreconditioner = .true.
case default
call POP_ErrorSet(errorCode, &
'POP_SolversInit: unknown preconditioner choice')
return
end select
!-----------------------------------------------------------------------
!
! compute nine point operator coefficients: compute on baroclinic
! decomposition first where grid info defined and redistribute
! to barotropic distribution
! leave center coefficients in baroclinic distribution to facilitate
! easy time-dependent changes in barotropic routine
!
!-----------------------------------------------------------------------
call POP_DistributionGet(POP_distrbClinic, errorCode, &
numLocalBlocks = numBlocksClinic)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error retrieving clinic local block count')
return
endif
allocate(work0 (POP_nxBlock,POP_nyBlock,numBlocksClinic), &
workNorth (POP_nxBlock,POP_nyBlock,numBlocksClinic), &
workEast (POP_nxBlock,POP_nyBlock,numBlocksClinic), &
workNE (POP_nxBlock,POP_nyBlock,numBlocksClinic), &
mMaskTmp (POP_nxBlock,POP_nyBlock,numBlocksClinic), &
centerWgtClinicIndep(POP_nxBlock,POP_nyBlock,numBlocksClinic), &
centerWgtClinic (POP_nxBlock,POP_nyBlock,numBlocksClinic), &
stat=istat)
if (istat > 0) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error allocating temporary arrays')
return
endif
!$OMP PARALLEL DO PRIVATE(iblock,i,j,xne,xse,xnw,xsw,yne,yse,ynw,ysw, &
!$OMP ase,anw,asw)
do iblock = 1,numBlocksClinic
work0 (:,:,iblock) = 0.0_POP_r8
workNorth (:,:,iblock) = 0.0_POP_r8
workEast (:,:,iblock) = 0.0_POP_r8
workNE (:,:,iblock) = 0.0_POP_r8
mMaskTmp (:,:,iblock) = 0.0_POP_r8
centerWgtClinicIndep (:,:,iblock) = 0.0_POP_r8
do j=2,POP_nyBlock
do i=2,POP_nxBlock
xne = 0.25_POP_r8*HU(i ,j ,iblock)*DXUR(i ,j ,iblock)* &
DYU (i ,j ,iblock)
xse = 0.25_POP_r8*HU(i ,j-1,iblock)*DXUR(i ,j-1,iblock)* &
DYU (i ,j-1,iblock)
xnw = 0.25_POP_r8*HU(i-1,j ,iblock)*DXUR(i-1,j ,iblock)* &
DYU (i-1,j ,iblock)
xsw = 0.25_POP_r8*HU(i-1,j-1,iblock)*DXUR(i-1,j-1,iblock)* &
DYU (i-1,j-1,iblock)
yne = 0.25_POP_r8*HU(i ,j ,iblock)*DYUR(i ,j ,iblock)* &
DXU (i ,j ,iblock)
yse = 0.25_POP_r8*HU(i ,j-1,iblock)*DYUR(i ,j-1,iblock)* &
DXU (i ,j-1,iblock)
ynw = 0.25_POP_r8*HU(i-1,j ,iblock)*DYUR(i-1,j ,iblock)* &
DXU (i-1,j ,iblock)
ysw = 0.25_POP_r8*HU(i-1,j-1,iblock)*DYUR(i-1,j-1,iblock)* &
DXU (i-1,j-1,iblock)
workNE(i,j,iblock) = xne + yne
ase = xse + yse
anw = xnw + ynw
asw = xsw + ysw
workEast (i,j,iblock) = xne + xse - yne - yse
workNorth(i,j,iblock) = yne + ynw - xne - xnw
centerWgtClinicIndep(i,j,iblock) = &
-(workNE(i,j,iblock) + ase + anw + asw)
work0 (i,j,iblock) = TAREA (i,j,iblock)**2
mMaskTmp(i,j,iblock) = RCALCT(i,j,iblock)
end do
end do
end do
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!
! redistribute operator weights and mask to barotropic distribution
!
!-----------------------------------------------------------------------
call POP_DistributionGet(POP_distrbTropic, errorCode, &
numLocalBlocks = numBlocksTropic)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error retrieving tropic local block count')
return
endif
allocate(btropWgtCenter (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
btropWgtNorth (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
btropWgtEast (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
btropWgtNE (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
mMaskTropic (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
stat=istat)
if (istat > 0) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error allocating operator weights')
return
endif
call POP_RedistributeBlocks(btropWgtNorth, POP_distrbTropic, &
workNorth, POP_distrbClinic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error redistributing north operator weight')
return
endif
call POP_RedistributeBlocks(btropWgtEast, POP_distrbTropic, &
workEast, POP_distrbClinic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error redistributing east operator weight')
return
endif
call POP_RedistributeBlocks(btropWgtNE, POP_distrbTropic, &
workNE, POP_distrbClinic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error redistributing NE operator weight')
return
endif
call POP_RedistributeBlocks(mMaskTropic, POP_distrbTropic, &
mMaskTmp, POP_distrbClinic, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error redistributing mask')
return
endif
!-----------------------------------------------------------------------
!
! calculate normalization constant (darea,darea) for rmsResidual
! in cgr routine.
!
!-----------------------------------------------------------------------
residualNorm = 1.0_POP_r8/POP_GlobalSum(work0, POP_distrbClinic, &
POP_gridHorzLocCenter, &
errorCode, &
mMask = mMaskTmp)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error computing normalization factor')
return
endif
convergenceCriterion = convergenceCriterion**2/residualNorm
deallocate(mMaskTmp, stat=istat)
if (istat > 0) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error deallocating temp mask')
return
endif
!-----------------------------------------------------------------------
!
! setup preconditioner if required
!
!-----------------------------------------------------------------------
if (usePreconditioner) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: preconditioner not supported')
return
allocate( &
precondCenter (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
precondNorth (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
precondSouth (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
precondEast (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
precondWest (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
precondNE (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
precondSE (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
precondNW (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
precondSW (POP_nxBlock,POP_nyBlock,numBlocksTropic), &
stat = istat)
if (istat > 0) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error allocating preconditioner')
return
endif
!****
!**** OLD CODE FOR READING PRECONDITIONER - NEEDS TO BE REVISED
!****
! if (POP_myTask == POP_masterTask) then
! write(stdout,*) ' Preconditioner read from file: ', &
! trim(precond_file)
! endif
!
! allocate(WORKS (POP_nxBlock,POP_nyBlock,nblocks_clinic), &
! WORKW (POP_nxBlock,POP_nyBlock,nblocks_clinic), &
! WORKSE(POP_nxBlock,POP_nyBlock,nblocks_clinic), &
! WORKSW(POP_nxBlock,POP_nyBlock,nblocks_clinic))
!
! allocate(icheck(nblocks_clinic))
!
!-----------------------------------------------------------------------
!
! read preconditioner and check that it is consistent with
! KMU field
!
!-----------------------------------------------------------------------
!
! call open_parallel_file(nu,precond_file,recl_dbl)
! call read_array(nu,workC)
! call read_array(nu,workN)
! call read_array(nu,WORKS)
! call read_array(nu,workE)
! call read_array(nu,WORKW)
! call read_array(nu,workNE)
! call read_array(nu,workNW)
! call read_array(nu,WORKSE)
! call read_array(nu,WORKSW)
! call close_parallel_file(nu)
!
! if (POP_myTask == POP_masterTask) then
! write(stdout,blank_fmt)
! write(stdout,*) ' file read: ', trim(precond_file)
! endif
!
!-----------------------------------------------------------------------
!
! check that PC is consistent with KMU field
!
!-----------------------------------------------------------------------
!
! do iblock = 1,nblocks_clinic
!
! icheck(iblock) = 0
!
! do j=1,POP_nyBlock
! do i=1,POP_nxBlock
!
! mlandne = .false.
! mlandnw = .false.
! mlandse = .false.
! mlandsw = .false.
! if (KMU(i ,j ,iblock) == 0) mlandne = .true.
! if (KMU(i-1,j ,iblock) == 0) mlandnw = .true.
! if (KMU(i ,j-1,iblock) == 0) mlandse = .true.
! if (KMU(i-1,j-1,iblock) == 0) mlandsw = .true.
!
! if (mlandne .and. workNE(i,j,iblock) /= 0.0_POP_r8) &
! icheck(iblock) = icheck(iblock) + 1
!
! if (mlandnw .and. workNW(i,j,iblock) /= 0.0_POP_r8) &
! icheck(iblock) = icheck(iblock) + 1
!
! if (mlandse .and. WORKSE(i,j,iblock) /= 0.0_POP_r8) &
! icheck(iblock) = icheck(iblock) + 1
!
! if (mlandsw .and. WORKSW(i,j,iblock) /= 0.0_POP_r8) &
! icheck(iblock) = icheck(iblock) + 1
!
! if (mlandne .and. mlandnw .and. (workN(i,j,iblock) /= 0.0_POP_r8)) &
! icheck(iblock) = icheck(iblock) + 1
! if (mlandne .and. mlandse .and. (workE(i,j,iblock) /= 0.0_POP_r8)) &
! icheck(iblock) = icheck(iblock) + 1
! if (mlandnw .and. mlandsw .and. (WORKW(i,j,iblock) /= 0.0_POP_r8)) &
! icheck(iblock) = icheck(iblock) + 1
! if (mlandse .and. mlandsw .and. (WORKS(i,j,iblock) /= 0.0_POP_r8)) &
! icheck(iblock) = icheck(iblock) + 1
! if (mlandne .and. mlandse .and. &
! mlandnw .and. mlandsw .and. (workC(i,j,iblock) /= 0.0_POP_r8)) &
! icheck(iblock) = icheck(iblock) + 1
! end do
! end do
! end do
!
! ncheck = sum(icheck)
! if (POP_GlobalSum(ncheck, POP_distrbClinic, errorCode) /= 0) then
! call POP_ErrorSet(errorCode, &
! 'POP_SolversInit: PC and KMU are incompatible')
! return
! endif
!
! deallocate(icheck)
!
! call redistribute_blocks(precondCenter ,distrb_tropic,workC ,distrb_clinic)
! call redistribute_blocks(precondNorth ,distrb_tropic,workN ,distrb_clinic)
! call redistribute_blocks(precondEast ,distrb_tropic,workE ,distrb_clinic)
! call redistribute_blocks(precondSouth ,distrb_tropic,WORKS ,distrb_clinic)
! call redistribute_blocks(precondWest ,distrb_tropic,WORKW ,distrb_clinic)
! call redistribute_blocks(precondNE,distrb_tropic,workNE,distrb_clinic)
! call redistribute_blocks(precondNW,distrb_tropic,workNW,distrb_clinic)
! call redistribute_blocks(precondSE,distrb_tropic,WORKSE,distrb_clinic)
! call redistribute_blocks(precondSW,distrb_tropic,WORKSW,distrb_clinic)
!
! deallocate(WORKS, WORKW, workNW, WORKSE, WORKSW)
!
endif
!-----------------------------------------------------------------------
!
! clean up temporary arrays
!
!-----------------------------------------------------------------------
deallocate(work0, workNorth, workEast, workNE, stat=istat)
if (istat > 0) then
call POP_ErrorSet(errorCode, &
'POP_SolversInit: error deallocating temp mask')
return
endif
!-----------------------------------------------------------------------
!EOC
end subroutine POP_SolversInit
!***********************************************************************
!BOP
! !IROUTINE: POP_SolversDiagonal
! !INTERFACE:
subroutine POP_SolversDiagonal(diagonalCorrection, blockIndx, errorCode) 1
! !DESCRIPTION:
! This routine corrects the center barotropic operator by
! subtracting off the time dependent diagnonal term computed in
! the barotropic driver.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (POP_r8), dimension(:,:), intent(in) :: &
diagonalCorrection ! time dependent diagonal term
integer (POP_i4), intent(in) :: &
blockIndx ! local block index
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! subtract the time dependent diagonal term from the center operator
! coefficient
!
!-----------------------------------------------------------------------
errorCode = POP_Success
centerWgtClinic(:,:,blockIndx) = &
centerWgtClinicIndep(:,:,blockIndx) - &
diagonalCorrection(:,:)
!-----------------------------------------------------------------------
!EOC
end subroutine POP_SolversDiagonal
!***********************************************************************
!BOP
! !IROUTINE: POP_SolversGetDiagnostics
! !INTERFACE:
subroutine POP_SolversGetDiagnostics(iterationCount, residual, & 1
errorCode)
! !DESCRIPTION:
! This routine returns the latest iteration count and residual
! for diagnostic purposes.
!
! !REVISION HISTORY:
! same as module
! !OUTPUT PARAMETERS:
real (POP_r8), intent(out) :: &
residual ! latest residual from solver
integer (POP_i4), intent(out) :: &
iterationCount, &! latest iteration count from solver
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! return the diagnostic quantities
!
!-----------------------------------------------------------------------
errorCode = POP_Success
residual = rmsResidual
iterationCount = numIterations
!-----------------------------------------------------------------------
!EOC
end subroutine POP_SolversGetDiagnostics
!***********************************************************************
!BOP
! !IROUTINE: pcg
! !INTERFACE:
subroutine pcg(X,B,errorCode) 1,18
! !DESCRIPTION:
! This routine uses a preconditioned conjugate-gradient solver to
! solve the equation $Ax=b$. Both the operator $A$ and preconditioner
! are nine-point stencils. If no preconditioner has been supplied,
! a diagonal preconditioner is applied. Convergence is checked
! every {\em ncheck} steps.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(in) :: &
B ! right hand side of linear system
! !INPUT/OUTPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(inout) :: &
X ! on input, an initial guess for the solution
! on output, solution of the linear system
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
i,j,m, &! local iteration counter
nx, ny, &! horizontal extents
numBlocks, &! number of local blocks
iblock ! local block counter
real (POP_r8) :: &
eta0,eta1,rr ! scalar inner product results
real (POP_r8), dimension(size(X,dim=1),size(X,dim=2), &
size(X,dim=3)) :: &
R, &! residual (b-Ax)
S, &! conjugate direction vector
Q,work0,work1 ! various cg intermediate results
!-----------------------------------------------------------------------
!
! compute initial residual and initialize S
!
!-----------------------------------------------------------------------
errorCode = POP_Success
call POP_DistributionGet(POP_distrbTropic, errorCode, &
numLocalBlocks = numBlocks)
nx = size(X,dim=1)
ny = size(X,dim=2)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: error retrieving local block count')
return
endif
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,numBlocks
call btropOperator(S,X,iblock)
do j=1,ny
do i=1,nx
R(i,j,iblock) = B(i,j,iblock) - S(i,j,iblock)
S(i,j,iblock) = 0.0_POP_r8
end do
end do
end do ! block loop
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!
! initialize fields and scalars
!
!-----------------------------------------------------------------------
call POP_HaloUpdate(R, POP_haloTropic, POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: error updating initial residual halo')
return
endif
eta0 =1.0_POP_r8
numIterations = maxIterations
!-----------------------------------------------------------------------
!
! iterate
!
!-----------------------------------------------------------------------
iterationLoop: do m = 1, maxIterations
!-----------------------------------------------------------------------
!
! calculate (PC)r
! diagonal preconditioner if preconditioner not specified
!
!-----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,numBlocks
if (usePreconditioner) then
call preconditioner(work1,R,iblock)
else
do j=1,ny
do i=1,nx
if (btropWgtCenter(i,j,iblock) /= 0.0_POP_r8) then
work1(i,j,iblock) = R(i,j,iblock)/ &
btropWgtCenter(i,j,iblock)
else
work1(i,j,iblock) = 0.0_POP_r8
endif
end do
end do
endif
do j=1,ny
do i=1,nx
work0(i,j,iblock) = R(i,j,iblock)*work1(i,j,iblock)
end do
end do
end do ! block loop
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!
! update conjugate direction vector s
!
!-----------------------------------------------------------------------
if (usePreconditioner) then
call POP_HaloUpdate(work1, POP_haloTropic, &
POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: error updating halo for preconditioner')
return
endif
endif
!*** (r,(PC)r)
eta1 = POP_GlobalSum(work0, POP_distrbTropic, &
POP_gridHorzLocCenter, &
errorCode, mMask = mMaskTropic)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: error in initial dot product')
return
endif
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,numBlocks
do j=1,ny
do i=1,nx
S(i,j,iblock) = work1(i,j,iblock) + &
S(i,j,iblock)*(eta1/eta0)
end do
end do
!-----------------------------------------------------------------------
!
! compute As
!
!-----------------------------------------------------------------------
call btropOperator(Q,S,iblock)
do j=1,ny
do i=1,nx
work0(i,j,iblock) = Q(i,j,iblock)*S(i,j,iblock)
end do
end do
end do ! block loop
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!
! compute next solution and residual
!
!-----------------------------------------------------------------------
call POP_HaloUpdate(Q, POP_haloTropic, POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: error updating Q halo')
return
endif
eta0 = eta1
eta1 = eta0/POP_GlobalSum(work0, POP_distrbTropic, &
POP_gridHorzLocCenter, errorCode, &
mMask = mMaskTropic)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: error in dot product')
return
endif
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,numBlocks
do j=1,ny
do i=1,nx
X(i,j,iblock) = X(i,j,iblock) + eta1*S(i,j,iblock)
R(i,j,iblock) = R(i,j,iblock) - eta1*Q(i,j,iblock)
end do
end do
if (mod(m,convergenceCheckFreq) == 0) then
call btropOperator(R,X,iblock)
do j=1,ny
do i=1,nx
R(i,j,iblock) = B(i,j,iblock) - R(i,j,iblock)
work0(i,j,iblock) = R(i,j,iblock)*R(i,j,iblock)
end do
end do
endif
end do ! block loop
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!
! test for convergence
!
!-----------------------------------------------------------------------
if (mod(m,convergenceCheckFreq) == 0) then
call POP_HaloUpdate(R, POP_haloTropic, POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: error updating residual halo in convrg')
return
endif
rr = POP_GlobalSum(work0, POP_distrbTropic, &
POP_gridHorzLocCenter, &
errorCode, mMask = mMaskTropic) ! (r,r)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: error computing convergence dot prod')
return
endif
if (rr < convergenceCriterion) then
numIterations = m
exit iterationLoop
endif
endif
enddo iterationLoop
rmsResidual = sqrt(rr*residualNorm)
if (numIterations == maxIterations) then
if (convergenceCriterion /= 0.0_POP_r8) then
call POP_ErrorSet(errorCode, &
'POP_SolversPCG: solver not converged')
return
endif
endif
!-----------------------------------------------------------------------
!EOC
end subroutine pcg
!***********************************************************************
!BOP
! !IROUTINE: ChronGear
! !INTERFACE:
subroutine ChronGear(X,B,errorCode) 1,20
! !DESCRIPTION:
! This routine implements the Chronopoulos-Gear conjugate-gradient
! solver with preconditioner for solving the linear system $Ax=b$.
! It is a rearranged conjugate gradient solver that reduces the
! number of inner products per iteration from two to one (not
! counting convergence check). Both the operator $A$ and
! preconditioner are nine-point stencils. If no preconditioner has
! been supplied, a diagonal preconditioner is applied. Convergence
! is checked every {\em ncheck} steps.
!
!
! References:
! Dongarra, J. and V. Eijkhout. LAPACK Working Note 159.
! Finite-choice algorithm optimization in conjugate gradients.
! Tech. Rep. ut-cs-03-502. Computer Science Department.
! University of Tennessee, Knoxville. 2003.
!
! D Azevedo, E.F., V.L. Eijkhout, and C.H. Romine. LAPACK Working
! Note 56. Conjugate gradient algorithms with reduced
! synchronization overhead on distributed memory multiprocessors.
! Tech. Rep. CS-93-185. Computer Science Department.
! University of Tennessee, Knoxville. 1993.
!
! !REVISION HISTORY:
! this routine implemented by Frank Bryan et al., NCAR
! !INPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(in) :: &
B ! right hand side of linear system
! !INPUT/OUTPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(inout) :: &
X ! on input, an initial guess for the solution
! on output, solution of the linear system
! !OUTPUT PARAMETERS:
integer (POP_i4), intent(out) :: &
errorCode ! returned error code
!EOP
!BOC
!-----------------------------------------------------------------------
integer (POP_i4) :: &
i,j,m, &! local iteration counter
nx, ny, &! horizontal extents
numBlocks, &! number of local blocks
iblock ! local block counter
real (POP_r8) :: & ! scalar results
cgAlpha, cgBeta, cgSigma, cgDelta, cgRhoOld, cgRho, rr
real (POP_r8), dimension(size(X,dim=1),size(X,dim=2), &
size(X,dim=3)) :: &
R, &! residual (b-Ax)
S, &! conjugate direction vector
Q,Z,AZ,WORK0, &! various cg intermediate results
A0R ! diagonal preconditioner
real (POP_r8), dimension(size(X,dim=1),size(X,dim=2), &
2,size(X,dim=3)) :: &
WORKN ! WORK array
real (POP_r8), dimension(2) :: &
sumN ! global sum results for multiple arrays
!-----------------------------------------------------------------------
!
! initialize some scalars
!
!-----------------------------------------------------------------------
errorCode = POP_Success
call POP_DistributionGet(POP_distrbTropic, errorCode, &
numLocalBlocks = numBlocks)
nx = size(X,dim=1)
ny = size(X,dim=2)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: error retrieving local block count')
return
endif
cgRho = 1.0_POP_r8
numIterations = maxIterations
!-----------------------------------------------------------------------
!
! compute initial residual and initialize other arrays
!
!-----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,numBlocks
do j=1,ny
do i=1,nx
R (i,j,iblock) = 0.0_POP_r8
S (i,j,iblock) = 0.0_POP_r8
Z (i,j,iblock) = 0.0_POP_r8
Q (i,j,iblock) = 0.0_POP_r8
AZ (i,j,iblock) = 0.0_POP_r8
WORK0(i,j,iblock) = 0.0_POP_r8
WORKN(i,j,1,iblock) = 0.0_POP_r8
WORKN(i,j,2,iblock) = 0.0_POP_r8
end do
end do
if (.not. usePreconditioner) then
!--- diagonal preconditioner if preconditioner not specified
do j=1,ny
do i=1,nx
if (btropWgtCenter(i,j,iblock) /= 0.0_POP_r8) then
A0R(i,j,iblock) = 1.0_POP_r8/btropWgtCenter(i,j,iblock)
else
A0R(i,j,iblock) = 0.0_POP_r8
endif
end do
end do
endif
! use S as a temp here for Ax
call btropOperator(S,X,iblock)
do j=1,ny
do i=1,nx
R(i,j,iblock) = B(i,j,iblock) - S(i,j,iblock) ! b-Ax
end do
end do
end do ! block loop
!$OMP END PARALLEL DO
call POP_HaloUpdate(R, POP_haloTropic, POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: error updating initial residual halo')
return
endif
!-----------------------------------------------------------------------
!
! take one pass of standard algorithm
!
!-----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,numBlocks
!---- calculate (PC)r store in Z
if (usePreconditioner) then
call preconditioner(Z,R,iblock)
else ! use diagonal preconditioner
do j=1,ny
do i=1,nx
Z(i,j,iblock) = R(i,j,iblock)*A0R(i,j,iblock)
end do
end do
endif
!---- Compute intermediate result for dot product
!---- update conjugate direction vector S
do j=1,ny
do i=1,nx
WORKN(i,j,1,iblock) = R(i,j,iblock)*Z(i,j,iblock)
S(i,j,iblock) = Z(i,j,iblock)
end do
end do
!---- compute Q = A * S
call btropOperator(Q,S,iblock)
!---- compute intermediate result for dot product
do j=1,ny
do i=1,nx
WORKN(i,j,2,iblock) = S(i,j,iblock)*Q(i,j,iblock)
end do
end do
end do
!$OMP END PARALLEL DO
call POP_HaloUpdate(Q, POP_haloTropic, POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: error updating Q halo')
return
endif
!---- Form dot products
sumN = POP_GlobalSum(WORKN, POP_distrbTropic, &
POP_gridHorzLocCenter, &
errorCode, mMask = mMaskTropic)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: error in initial dot products')
return
endif
cgRhoOld = sumN(1) !(r,PCr)
cgSigma = sumN(2) !(s,As)
cgAlpha = cgRhoOld/cgSigma
!---- compute first solution and residual
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,numBlocks
do j=1,ny
do i=1,nx
X(i,j,iblock) = X(i,j,iblock) + cgAlpha*S(i,j,iblock)
R(i,j,iblock) = R(i,j,iblock) - cgAlpha*Q(i,j,iblock)
end do
end do
end do
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!
! iterate
!
!-----------------------------------------------------------------------
iterationLoop: do m = 1, maxIterations
!-----------------------------------------------------------------------
!
! calculate (PC)r and A*(Pc)r
!
!-----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock,i,j)
do iblock=1,numBlocks
if (usePreconditioner) then
call preconditioner(Z,R,iblock)
else
do j=1,ny
do i=1,nx
Z(i,j,iblock) = R(i,j,iblock)*A0R(i,j,iblock)
end do
end do
endif
call btropOperator(AZ,Z,iblock)
!--- intermediate results for inner products
do j=1,ny
do i=1,nx
WORKN(i,j,1,iblock) = R(i,j,iblock)*Z(i,j,iblock)
WORKN(i,j,2,iblock) = AZ(i,j,iblock)*Z(i,j,iblock)
end do
end do
end do
!$OMP END PARALLEL DO
call POP_HaloUpdate(AZ, POP_haloTropic, POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: error updating AZ halo')
return
endif
sumN = POP_GlobalSum(WORKN, POP_distrbTropic, &
POP_gridHorzLocCenter, &
errorCode, mMask = mMaskTropic)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: error in dot products')
return
endif
cgRho = sumN(1) ! (r,(PC)r)
cgDelta = sumN(2) ! (A (PC)r,(PC)r)
cgBeta = cgRho/cgRhoOld
cgSigma = cgDelta - (cgBeta**2)*cgSigma
cgAlpha = cgRho/cgSigma
cgRhoOld = cgRho
!-----------------------------------------------------------------------
!
! compute S and Q
! compute next solution and residual
!
!-----------------------------------------------------------------------
!$OMP PARALLEL DO PRIVATE(iblock, i,j)
do iblock=1,numBlocks
do j=1,ny
do i=1,nx
S(i,j,iblock) = Z(i,j,iblock) + cgBeta *S(i,j,iblock)
Q(i,j,iblock) = AZ(i,j,iblock) + cgBeta *Q(i,j,iblock)
X(i,j,iblock) = X(i,j,iblock) + cgAlpha*S(i,j,iblock)
R(i,j,iblock) = R(i,j,iblock) - cgAlpha*Q(i,j,iblock)
end do
end do
!--- recompute residual as b-Ax for convergence check
if (mod(m,convergenceCheckFreq) == 0) then
!--- Reset residual using r = b - Ax
!--- (r,r) for norm of residual
call btropOperator(Z,X,iblock)
do j=1,ny
do i=1,nx
R(i,j,iblock) = B(i,j,iblock) - Z(i,j,iblock)
WORK0(i,j,iblock) = R(i,j,iblock)*R(i,j,iblock)
end do
end do
endif
end do
!$OMP END PARALLEL DO
!-----------------------------------------------------------------------
!
! test for convergence if it is time
!
!-----------------------------------------------------------------------
if (mod(m,convergenceCheckFreq) == 0) then
!--- update ghost cells for next iteration
call POP_HaloUpdate(R, POP_haloTropic, POP_gridHorzLocCenter, &
POP_fieldKindScalar, errorCode)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: error updating residual halo in convrg')
return
endif
!--- residual norm for convergence
rr = POP_GlobalSum(work0, POP_distrbTropic, &! (r,r)
POP_gridHorzLocCenter, &
errorCode, mMask = mMaskTropic)
if (errorCode /= POP_Success) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: error computing convergence dot prod')
return
endif
if (rr < convergenceCriterion) then
numIterations = m
exit iterationLoop
endif
endif
end do iterationLoop
rmsResidual = sqrt(rr*residualNorm)
if (numIterations == maxIterations) then
if (convergenceCriterion /= 0.0_POP_r8) then
call POP_ErrorSet(errorCode, &
'POP_SolversChronGear: solver not converged')
return
endif
endif
!-----------------------------------------------------------------------
!EOC
end subroutine ChronGear
!***********************************************************************
!BOP
! !IROUTINE: preconditioner
! !INTERFACE:
subroutine preconditioner(PX,X,bid) 4
! !DESCRIPTION:
! This function applies a precomputed preconditioner as a nine-point
! stencil operator.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(in) :: &
X ! array to be operated on
integer (POP_i4), intent(in) :: &
bid ! local block address for this block
! !OUTPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(out) :: &
PX ! nine point operator result
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
i,j ! dummy counters
!-----------------------------------------------------------------------
PX(:,:,bid) = 0.0_POP_r8
do j=2,size(X,dim=2)-1
do i=2,size(X,dim=1)-1
PX(i,j,bid) = precondNE (i,j,bid)*X(i+1,j+1,bid) + &
precondNW (i,j,bid)*X(i-1,j+1,bid) + &
precondSE (i,j,bid)*X(i+1,j-1,bid) + &
precondSW (i,j,bid)*X(i-1,j-1,bid) + &
precondNorth (i,j,bid)*X(i ,j+1,bid) + &
precondSouth (i,j,bid)*X(i ,j-1,bid) + &
precondEast (i,j,bid)*X(i+1,j ,bid) + &
precondWest (i,j,bid)*X(i-1,j ,bid) + &
precondCenter(i,j,bid)*X(i ,j ,bid)
end do
end do
!-----------------------------------------------------------------------
!EOC
end subroutine preconditioner
!***********************************************************************
!BOP
! !IROUTINE: btropOperator
! !INTERFACE:
subroutine btropOperator(AX,X,bid) 7
! !DESCRIPTION:
! This routine applies the nine-point stencil operator for the
! barotropic solver. It takes advantage of some 9pt weights being
! shifted versions of others.
!
! !REVISION HISTORY:
! same as module
! !INPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(in) :: &
X ! array to be operated on
integer (POP_i4), intent(in) :: &
bid ! local block address for this block
! !OUTPUT PARAMETERS:
real (POP_r8), dimension(:,:,:), intent(out) :: &
AX ! nine point operator result (Ax)
!EOP
!BOC
!-----------------------------------------------------------------------
!
! local variables
!
!-----------------------------------------------------------------------
integer (POP_i4) :: &
i,j ! dummy counters
!-----------------------------------------------------------------------
AX(:,:,bid) = 0.0_POP_r8
do j=2,size(X,dim=2)-1
do i=2,size(X,dim=1)-1
AX(i,j,bid) = btropWgtCenter(i ,j ,bid)*X(i ,j ,bid) + &
btropWgtNorth (i ,j ,bid)*X(i ,j+1,bid) + &
btropWgtNorth (i ,j-1,bid)*X(i ,j-1,bid) + &
btropWgtEast (i ,j ,bid)*X(i+1,j ,bid) + &
btropWgtEast (i-1,j ,bid)*X(i-1,j ,bid) + &
btropWgtNE (i ,j ,bid)*X(i+1,j+1,bid) + &
btropWgtNE (i ,j-1,bid)*X(i+1,j-1,bid) + &
btropWgtNE (i-1,j ,bid)*X(i-1,j+1,bid) + &
btropWgtNE (i-1,j-1,bid)*X(i-1,j-1,bid)
end do
end do
!-----------------------------------------------------------------------
!EOC
end subroutine btropOperator
!***********************************************************************
end module POP_SolversMod
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||