_RL psFac, nhFac
psFac = pfFacMom*implicSurfPress
& -deltaTmom*psFac*phiSurfX(i,j)
& -deltaTmom*psFac*phiSurfY(i,j)
& - psFac*phiSurfX(i,j)
& - psFac*phiSurfY(i,j)
_RL phxFac,phyFac, psFac
psFac = pfFacMom*(1. _d 0 - implicSurfPress)