model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.69 2009/11/23 16:15:54 mlosch Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: SOLVE_FOR_PRESSURE
C     !INTERFACE:
      SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SOLVE_FOR_PRESSURE
C     | o Controls inversion of two and/or three-dimensional
C     |   elliptic problems for the pressure field.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "FFIELDS.h"
#include "DYNVARS.h"
#include "SOLVE_FOR_PRESSURE.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "SOLVE_FOR_PRESSURE3D.h"
#include "NH_VARS.h"
#endif
#ifdef ALLOW_CD_CODE
#include "CD_CODE_VARS.h"
#endif
#ifdef ALLOW_OBCS
#include "OBCS.h"
#endif

C     === Functions ====
      LOGICAL  DIFFERENT_MULTIPLE
      EXTERNAL DIFFERENT_MULTIPLE

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myTime :: Current time in simulation
C     myIter :: Current iteration number in simulation
C     myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
      _RL myTime
      INTEGER myIter
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER i,j,k,bi,bj
      _RL firstResidual,lastResidual
      _RL tmpFac
      _RL sumEmP, tileEmP(nSx,nSy)
      LOGICAL putPmEinXvector
      INTEGER numIters, ks
      CHARACTER*10 sufx
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_NONHYDROSTATIC
      INTEGER kp1
      _RL     wFacKm, wFacKp
      LOGICAL zeroPsNH
      _RL     uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL     vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
#else
      _RL     cg3d_b(1)
#endif
CEOP

#ifdef ALLOW_NONHYDROSTATIC
        zeroPsNH = .FALSE.
c       zeroPsNH = exactConserv
#else
        cg3d_b(1) = 0.
#endif

C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
C anelastic (always Z-coordinate):
C     1) assume that rhoFacF(1)=1 (and ksurf == 1);
C        (this reduces the number of lines of code to modify)
C     2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
C        (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
C       => 2 factors cancel in elliptic eq. for Phi_s ,
C       but 1rst factor(a) remains in RHS cg2d_b.

C--   Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
C     instead of simply etaN ; This can speed-up the solver convergence in
C     the case where |Global_mean_PmE| is large.
      putPmEinXvector = .FALSE.
c     putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater

C--   Save previous solution & Initialise Vector solution and source term :
      sumEmP = 0.
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
#ifdef ALLOW_CD_CODE
          etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
#endif
          cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
          cg2d_b(i,j,bi,bj) = 0.
         ENDDO
        ENDDO
        IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
         tmpFac = freeSurfFac*mass2rUnit
         IF (exactConserv)
     &        tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
         DO j=1,sNy
          DO i=1,sNx
           cg2d_b(i,j,bi,bj) =
     &       tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
          ENDDO
         ENDDO
        ENDIF
        IF ( putPmEinXvector ) THEN
         tileEmP(bi,bj) = 0.
         DO j=1,sNy
          DO i=1,sNx
            tileEmP(bi,bj) = tileEmP(bi,bj)
     &                     + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
     &                                    *maskH(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDIF
       ENDDO
      ENDDO
      IF ( putPmEinXvector ) THEN
        CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
      ENDIF

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        IF ( putPmEinXvector ) THEN
          tmpFac = 0.
          IF (globalArea.GT.0.) tmpFac =
     &      freeSurfFac*deltaTfreesurf*mass2rUnit*sumEmP/globalArea
          DO j=1,sNy
           DO i=1,sNx
            cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
     &                        - tmpFac*Bo_surf(i,j,bi,bj)
           ENDDO
          ENDDO
        ENDIF
C- RHS: similar to the divergence of the vertically integrated mass transport:
C       del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] }  / deltaT
        DO k=Nr,1,-1
         CALL CALC_DIV_GHAT(
     I                       bi,bj,k,
     U                       cg2d_b, cg3d_b,
     I                       myThid )
        ENDDO
       ENDDO
      ENDDO

C--   Add source term arising from w=d/dt (p_s + p_nh)
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_NONHYDROSTATIC
C--   Add EmPmR contribution to top level cg3d_b:
C      (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
        IF ( use3Dsolver .AND.
     &       useRealFreshWaterFlux.AND.fluidIsWater ) THEN
          tmpFac = freeSurfFac*mass2rUnit
          IF (exactConserv)
     &     tmpFac = freeSurfFac*mass2rUnit*implicDiv2DFlow
          ks = 1
          IF ( usingPCoords ) ks = Nr
          DO j=1,sNy
           DO i=1,sNx
            cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
     &        + tmpFac*_rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)/deltaTMom
           ENDDO
          ENDDO
        ENDIF
        IF ( use3Dsolver .AND. zeroPsNH ) THEN
         DO j=1,sNy
          DO i=1,sNx
           ks = ksurfC(i,j,bi,bj)
           IF ( ks.LE.Nr ) THEN
            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
     &         /deltaTMom/deltaTfreesurf
     &         * etaH(i,j,bi,bj)
            cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
     &         /deltaTMom/deltaTfreesurf
     &         * etaH(i,j,bi,bj)
           ENDIF
          ENDDO
         ENDDO
        ELSEIF ( use3Dsolver ) THEN
         DO j=1,sNy
          DO i=1,sNx
           ks = ksurfC(i,j,bi,bj)
           IF ( ks.LE.Nr ) THEN
            cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
     &         /deltaTMom/deltaTfreesurf
     &         *( etaN(i,j,bi,bj)
     &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
            cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
     &         /deltaTMom/deltaTfreesurf
     &         *( etaN(i,j,bi,bj)
     &           +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
           ENDIF
          ENDDO
         ENDDO
        ELSEIF ( exactConserv ) THEN
#else
        IF ( exactConserv ) THEN
#endif /* ALLOW_NONHYDROSTATIC */
         DO j=1,sNy
          DO i=1,sNx
           ks = ksurfC(i,j,bi,bj)
           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
     &         /deltaTMom/deltaTfreesurf
     &         * etaH(i,j,bi,bj)
          ENDDO
         ENDDO
        ELSE
         DO j=1,sNy
          DO i=1,sNx
           ks = ksurfC(i,j,bi,bj)
           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)*deepFac2F(ks)
     &         /deltaTMom/deltaTfreesurf
     &         * etaN(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDIF

#ifdef ALLOW_OBCS
        IF (useOBCS) THEN
         DO i=1,sNx
C Northern boundary
          IF (OB_Jn(i,bi,bj).NE.0) THEN
           cg2d_b(i,OB_Jn(i,bi,bj),bi,bj)=0.
           cg2d_x(i,OB_Jn(i,bi,bj),bi,bj)=0.
          ENDIF
C Southern boundary
          IF (OB_Js(i,bi,bj).NE.0) THEN
           cg2d_b(i,OB_Js(i,bi,bj),bi,bj)=0.
           cg2d_x(i,OB_Js(i,bi,bj),bi,bj)=0.
          ENDIF
         ENDDO
         DO j=1,sNy
C Eastern boundary
          IF (OB_Ie(j,bi,bj).NE.0) THEN
           cg2d_b(OB_Ie(j,bi,bj),j,bi,bj)=0.
           cg2d_x(OB_Ie(j,bi,bj),j,bi,bj)=0.
          ENDIF
C Western boundary
          IF (OB_Iw(j,bi,bj).NE.0) THEN
           cg2d_b(OB_Iw(j,bi,bj),j,bi,bj)=0.
           cg2d_x(OB_Iw(j,bi,bj),j,bi,bj)=0.
          ENDIF
         ENDDO
        ENDIF
#endif /* ALLOW_OBCS */
C-    end bi,bj loops
       ENDDO
      ENDDO

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevB ) THEN
       CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
     &                        myThid)
      ENDIF
#endif
      IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
       WRITE(sufx,'(I10.10)') myIter
       CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
      ENDIF

C--   Find the surface pressure using a two-dimensional conjugate
C--   gradient solver.
C     see CG2D.h for the interface to this routine.
      firstResidual=0.
      lastResidual=0.
      numIters=cg2dMaxIters
c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
#ifdef ALLOW_CG2D_NSA
C--   Call the not-self-adjoint version of cg2d
      CALL CG2D_NSA(
     U           cg2d_b,
     U           cg2d_x,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
#else /* not ALLOW_CG2D_NSA = default */
#ifdef ALLOW_SRCG
      IF ( useSRCGSolver ) THEN
C--   Call the single reduce CG solver
       CALL CG2D_SR(
     U           cg2d_b,
     U           cg2d_x,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
      ELSE
#else
      IF (.TRUE.) THEN
C--   Call the default CG solver
#endif /* ALLOW_SRCG */
       CALL CG2D(
     U           cg2d_b,
     U           cg2d_x,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
      ENDIF
#endif /* ALLOW_CG2D_NSA */
      _EXCH_XY_RL( cg2d_x, myThid )
c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevB ) THEN
       CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
     &                        myThid)
      ENDIF
#endif

C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
      IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
     &   ) THEN
       IF ( debugLevel .GE. debLevA ) THEN
        _BEGIN_MASTER( myThid )
        WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        _END_MASTER( myThid )
       ENDIF
      ENDIF

C--   Transfert the 2D-solution to "etaN" :
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_NONHYDROSTATIC
      IF ( use3Dsolver ) THEN
       IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
        WRITE(sufx,'(I10.10)') myIter
        CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
       ENDIF

C--   Solve for a three-dimensional pressure term (NH or IGW or both ).
C     see CG3D.h for the interface to this routine.
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1,sNy+1
          DO i=1,sNx+1
           uf(i,j)=-_recip_dxC(i,j,bi,bj)*
     &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
           vf(i,j)=-_recip_dyC(i,j,bi,bj)*
     &         (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
          ENDDO
         ENDDO

#ifdef ALLOW_OBCS
         IF (useOBCS) THEN
          DO i=1,sNx+1
C Northern boundary
          IF (OB_Jn(i,bi,bj).NE.0) THEN
           vf(i,OB_Jn(i,bi,bj))=0.
          ENDIF
C Southern boundary
          IF (OB_Js(i,bi,bj).NE.0) THEN
           vf(i,OB_Js(i,bi,bj)+1)=0.
          ENDIF
          ENDDO
          DO j=1,sNy+1
C Eastern boundary
          IF (OB_Ie(j,bi,bj).NE.0) THEN
           uf(OB_Ie(j,bi,bj),j)=0.
          ENDIF
C Western boundary
          IF (OB_Iw(j,bi,bj).NE.0) THEN
           uf(OB_Iw(j,bi,bj)+1,J)=0.
          ENDIF
          ENDDO
         ENDIF
#endif /* ALLOW_OBCS */

         IF ( usingZCoords ) THEN
C-       Z coordinate: assume surface @ level k=1
           tmpFac = freeSurfFac*deepFac2F(1)
         ELSE
C-       Other than Z coordinate: no assumption on surface level index
           tmpFac = 0.
           DO j=1,sNy
            DO i=1,sNx
              ks = ksurfC(i,j,bi,bj)
              IF ( ks.LE.Nr ) THEN
               cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
     &              +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
     &                  *_rA(i,j,bi,bj)*deepFac2F(ks)/deltaTmom
              ENDIF
            ENDDO
           ENDDO
         ENDIF
         k=1
         kp1 = MIN(k+1,Nr)
         wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
         IF (k.GE.Nr) wFacKp = 0.
         DO j=1,sNy
          DO i=1,sNx
            cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
     &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
     &       -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
     &       +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
     &       -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
     &       +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
     &         -wVel(i,j,kp1,bi,bj)*wFacKp
     &        )*_rA(i,j,bi,bj)/deltaTmom
          ENDDO
         ENDDO
         DO k=2,Nr
          kp1 = MIN(k+1,Nr)
C-       deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
C        both appear in wVel term, but at 2 different levels
          wFacKm = deepFac2F( k )*rhoFacF( k )
          wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
          IF (k.GE.Nr) wFacKp = 0.
          DO j=1,sNy
           DO i=1,sNx
            cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
     &       +drF(k)*dyG(i+1,j,bi,bj)*_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
     &       -drF(k)*dyG( i ,j,bi,bj)*_hFacW( i ,j,k,bi,bj)*uf( i ,j)
     &       +drF(k)*dxG(i,j+1,bi,bj)*_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
     &       -drF(k)*dxG(i, j ,bi,bj)*_hFacS(i, j ,k,bi,bj)*vf(i, j )
     &       +( wVel(i,j, k ,bi,bj)*wFacKm*maskC(i,j,k-1,bi,bj)
     &         -wVel(i,j,kp1,bi,bj)*wFacKp
     &        )*_rA(i,j,bi,bj)/deltaTmom

           ENDDO
          ENDDO
         ENDDO

#ifdef ALLOW_OBCS
         IF (useOBCS) THEN
          DO k=1,Nr
          DO i=1,sNx
C Northern boundary
          IF (OB_Jn(i,bi,bj).NE.0) THEN
           cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj)=0.
          ENDIF
C Southern boundary
          IF (OB_Js(i,bi,bj).NE.0) THEN
           cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj)=0.
          ENDIF
          ENDDO
          DO j=1,sNy
C Eastern boundary
          IF (OB_Ie(j,bi,bj).NE.0) THEN
           cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj)=0.
          ENDIF
C Western boundary
          IF (OB_Iw(j,bi,bj).NE.0) THEN
           cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj)=0.
          ENDIF
          ENDDO
          ENDDO
         ENDIF
#endif /* ALLOW_OBCS */
C-    end bi,bj loops
        ENDDO
       ENDDO

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevB ) THEN
       CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
     &                        myThid)
      ENDIF
#endif
      IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
       WRITE(sufx,'(I10.10)') myIter
       CALL WRITE_FLD_XYZ_RL( 'cg3d_b.',sufx, cg3d_b, myIter, myThid )
      ENDIF

      firstResidual=0.
      lastResidual=0.
      numIters=cg3dMaxIters
      CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
      CALL CG3D(
     U           cg3d_b,
     U           phi_nh,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
      _EXCH_XYZ_RL( phi_nh, myThid )
      CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)

      IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
     &   ) THEN
       IF ( debugLevel .GE. debLevA ) THEN
        _BEGIN_MASTER( myThid )
        WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        _END_MASTER( myThid )
       ENDIF
      ENDIF

C--   Update surface pressure (account for NH-p @ surface level) and NH pressure:
      IF ( zeroPsNH ) THEN
       IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
        WRITE(sufx,'(I10.10)') myIter
        CALL WRITE_FLD_XYZ_RL( 'cg3d_x.',sufx,phi_nh, myIter, myThid )
       ENDIF
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

         IF ( usingZCoords ) THEN
C-       Z coordinate: assume surface @ level k=1
          DO k=2,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
     &                           - phi_nh(i,j,1,bi,bj)
            ENDDO
           ENDDO
          ENDDO
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
     &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
             phi_nh(i,j,1,bi,bj) = 0.
           ENDDO
          ENDDO
         ELSE
C-       Other than Z coordinate: no assumption on surface level index
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            ks = ksurfC(i,j,bi,bj)
            IF ( ks.LE.Nr ) THEN
             etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
     &                       *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
             DO k=Nr,1,-1
              phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
     &                            - phi_nh(i,j,ks,bi,bj)
             ENDDO
            ENDIF
           ENDDO
          ENDDO
         ENDIF

        ENDDO
       ENDDO
      ENDIF

      ENDIF
#endif /* ALLOW_NONHYDROSTATIC */

#ifdef ALLOW_SHOWFLOPS
      CALL SHOWFLOPS_INSOLVE( myThid)
#endif

      RETURN
      END
1	jmc	1.70	C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.69 2009/11/23 16:15:54 mlosch Exp $
2	heimbach	1.21	C $Name: $
3	cnh	1.1
4	edhill	1.39	#include "PACKAGES_CONFIG.h"
5	adcroft	1.5	#include "CPP_OPTIONS.h"
6	cnh	1.1
7	cnh	1.27	CBOP
8			C !ROUTINE: SOLVE_FOR_PRESSURE
9			C !INTERFACE:
10	jmc	1.29	SUBROUTINE SOLVE_FOR_PRESSURE(myTime, myIter, myThid)
11	cnh	1.27
12			C !DESCRIPTION: \bv
13			C ==========================================================
14	jmc	1.58	C \| SUBROUTINE SOLVE_FOR_PRESSURE
15			C \| o Controls inversion of two and/or three-dimensional
16			C \| elliptic problems for the pressure field.
17	cnh	1.27	C ==========================================================
18			C \ev
19
20			C !USES:
21	adcroft	1.8	IMPLICIT NONE
22	cnh	1.4	C == Global variables
23			#include "SIZE.h"
24			#include "EEPARAMS.h"
25			#include "PARAMS.h"
26	adcroft	1.12	#include "GRID.h"
27	jmc	1.17	#include "SURFACE.h"
28	jmc	1.28	#include "FFIELDS.h"
29	jmc	1.48	#include "DYNVARS.h"
30			#include "SOLVE_FOR_PRESSURE.h"
31	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
32	adcroft	1.25	#include "SOLVE_FOR_PRESSURE3D.h"
33	jmc	1.48	#include "NH_VARS.h"
34			#endif
35			#ifdef ALLOW_CD_CODE
36			#include "CD_CODE_VARS.h"
37	adcroft	1.12	#endif
38	adcroft	1.11	#ifdef ALLOW_OBCS
39	adcroft	1.9	#include "OBCS.h"
40	adcroft	1.11	#endif
41	cnh	1.4
42	jmc	1.32	C === Functions ====
43	jmc	1.46	LOGICAL DIFFERENT_MULTIPLE
44			EXTERNAL DIFFERENT_MULTIPLE
45	jmc	1.32
46	cnh	1.27	C !INPUT/OUTPUT PARAMETERS:
47	cnh	1.1	C == Routine arguments ==
48	jmc	1.58	C myTime :: Current time in simulation
49			C myIter :: Current iteration number in simulation
50			C myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
51	jmc	1.28	_RL myTime
52			INTEGER myIter
53	jmc	1.29	INTEGER myThid
54	cnh	1.4
55	cnh	1.27	C !LOCAL VARIABLES:
56	adcroft	1.22	C == Local variables ==
57	cnh	1.6	INTEGER i,j,k,bi,bj
58	adcroft	1.22	_RL firstResidual,lastResidual
59	jmc	1.36	_RL tmpFac
60	jmc	1.65	_RL sumEmP, tileEmP(nSx,nSy)
61	jmc	1.47	LOGICAL putPmEinXvector
62	jmc	1.58	INTEGER numIters, ks
63	jmc	1.61	CHARACTER*10 sufx
64	adcroft	1.25	CHARACTER*(MAX_LEN_MBUF) msgBuf
65	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
66	jmc	1.58	INTEGER kp1
67			_RL wFacKm, wFacKp
68	jmc	1.49	LOGICAL zeroPsNH
69	jmc	1.63	_RL uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70			_RL vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71			#else
72			_RL cg3d_b(1)
73	jmc	1.49	#endif
74	cnh	1.27	CEOP
75	jmc	1.17
76	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
77	jmc	1.68	zeroPsNH = .FALSE.
78			c zeroPsNH = exactConserv
79	jmc	1.63	#else
80			cg3d_b(1) = 0.
81	jmc	1.49	#endif
82
83	jmc	1.58	C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
84			C anelastic (always Z-coordinate):
85			C 1) assume that rhoFacF(1)=1 (and ksurf == 1);
86			C (this reduces the number of lines of code to modify)
87			C 2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
88			C (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
89			C => 2 factors cancel in elliptic eq. for Phi_s ,
90			C but 1rst factor(a) remains in RHS cg2d_b.
91
92	jmc	1.47	C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
93			C instead of simply etaN ; This can speed-up the solver convergence in
94			C the case where \|Global_mean_PmE\| is large.
95			putPmEinXvector = .FALSE.
96	jmc	1.64	c putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
97	jmc	1.47
98	jmc	1.17	C-- Save previous solution & Initialise Vector solution and source term :
99	jmc	1.47	sumEmP = 0.
100	jmc	1.17	DO bj=myByLo(myThid),myByHi(myThid)
101			DO bi=myBxLo(myThid),myBxHi(myThid)
102			DO j=1-OLy,sNy+OLy
103			DO i=1-OLx,sNx+OLx
104	edhill	1.40	#ifdef ALLOW_CD_CODE
105	jmc	1.17	etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
106	jmc	1.26	#endif
107	jmc	1.18	cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
108	jmc	1.17	cg2d_b(i,j,bi,bj) = 0.
109			ENDDO
110			ENDDO
111	jmc	1.64	IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
112	jmc	1.62	tmpFac = freeSurfFac*mass2rUnit
113	jmc	1.58	IF (exactConserv)
114	jmc	1.62	& tmpFac = freeSurfFacmass2rUnitimplicDiv2DFlow
115	jmc	1.29	DO j=1,sNy
116			DO i=1,sNx
117	jmc	1.58	cg2d_b(i,j,bi,bj) =
118	jmc	1.29	& tmpFac_rA(i,j,bi,bj)EmPmR(i,j,bi,bj)/deltaTMom
119			ENDDO
120			ENDDO
121			ENDIF
122	jmc	1.47	IF ( putPmEinXvector ) THEN
123	jmc	1.65	tileEmP(bi,bj) = 0.
124	jmc	1.47	DO j=1,sNy
125			DO i=1,sNx
126	jmc	1.67	tileEmP(bi,bj) = tileEmP(bi,bj)
127	jmc	1.65	& + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
128			& *maskH(i,j,bi,bj)
129	jmc	1.47	ENDDO
130			ENDDO
131			ENDIF
132	jmc	1.17	ENDDO
133			ENDDO
134	jmc	1.47	IF ( putPmEinXvector ) THEN
135	jmc	1.65	CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
136	jmc	1.47	ENDIF
137	adcroft	1.12
138			DO bj=myByLo(myThid),myByHi(myThid)
139			DO bi=myBxLo(myThid),myBxHi(myThid)
140	jmc	1.47	IF ( putPmEinXvector ) THEN
141			tmpFac = 0.
142	jmc	1.62	IF (globalArea.GT.0.) tmpFac =
143			& freeSurfFacdeltaTfreesurfmass2rUnit*sumEmP/globalArea
144	jmc	1.47	DO j=1,sNy
145			DO i=1,sNx
146			cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
147			& - tmpFac*Bo_surf(i,j,bi,bj)
148			ENDDO
149			ENDDO
150			ENDIF
151	jmc	1.58	C- RHS: similar to the divergence of the vertically integrated mass transport:
152			C del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] } / deltaT
153	jmc	1.63	DO k=Nr,1,-1
154	adcroft	1.12	CALL CALC_DIV_GHAT(
155	jmc	1.63	I bi,bj,k,
156			U cg2d_b, cg3d_b,
157			I myThid )
158	adcroft	1.12	ENDDO
159			ENDDO
160			ENDDO
161	cnh	1.4
162	adcroft	1.12	C-- Add source term arising from w=d/dt (p_s + p_nh)
163			DO bj=myByLo(myThid),myByHi(myThid)
164			DO bi=myBxLo(myThid),myBxHi(myThid)
165	adcroft	1.13	#ifdef ALLOW_NONHYDROSTATIC
166	jmc	1.70	C-- Add EmPmR contribution to top level cg3d_b:
167			C (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
168			IF ( use3Dsolver .AND.
169			& useRealFreshWaterFlux.AND.fluidIsWater ) THEN
170			tmpFac = freeSurfFac*mass2rUnit
171			IF (exactConserv)
172			& tmpFac = freeSurfFacmass2rUnitimplicDiv2DFlow
173			ks = 1
174			IF ( usingPCoords ) ks = Nr
175			DO j=1,sNy
176			DO i=1,sNx
177			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
178			& + tmpFac_rA(i,j,bi,bj)EmPmR(i,j,bi,bj)/deltaTMom
179			ENDDO
180			ENDDO
181			ENDIF
182	jmc	1.53	IF ( use3Dsolver .AND. zeroPsNH ) THEN
183	jmc	1.49	DO j=1,sNy
184			DO i=1,sNx
185	jmc	1.51	ks = ksurfC(i,j,bi,bj)
186			IF ( ks.LE.Nr ) THEN
187			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
188	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
189			& /deltaTMom/deltaTfreesurf
190	jmc	1.49	& * etaH(i,j,bi,bj)
191	jmc	1.51	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
192	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
193			& /deltaTMom/deltaTfreesurf
194	jmc	1.49	& * etaH(i,j,bi,bj)
195	jmc	1.51	ENDIF
196	jmc	1.49	ENDDO
197			ENDDO
198	jmc	1.53	ELSEIF ( use3Dsolver ) THEN
199	jmc	1.28	DO j=1,sNy
200			DO i=1,sNx
201	jmc	1.51	ks = ksurfC(i,j,bi,bj)
202			IF ( ks.LE.Nr ) THEN
203			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
204	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
205			& /deltaTMom/deltaTfreesurf
206	jmc	1.28	& *( etaN(i,j,bi,bj)
207	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
208	jmc	1.51	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
209	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
210			& /deltaTMom/deltaTfreesurf
211	jmc	1.28	& *( etaN(i,j,bi,bj)
212	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
213	jmc	1.51	ENDIF
214	jmc	1.28	ENDDO
215	adcroft	1.12	ENDDO
216	jmc	1.28	ELSEIF ( exactConserv ) THEN
217	adcroft	1.13	#else
218	jmc	1.26	IF ( exactConserv ) THEN
219	edhill	1.39	#endif /* ALLOW_NONHYDROSTATIC */
220	jmc	1.26	DO j=1,sNy
221			DO i=1,sNx
222	jmc	1.58	ks = ksurfC(i,j,bi,bj)
223	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
224	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
225			& /deltaTMom/deltaTfreesurf
226	jmc	1.26	& * etaH(i,j,bi,bj)
227			ENDDO
228			ENDDO
229			ELSE
230			DO j=1,sNy
231			DO i=1,sNx
232	jmc	1.58	ks = ksurfC(i,j,bi,bj)
233	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
234	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
235			& /deltaTMom/deltaTfreesurf
236	jmc	1.26	& * etaN(i,j,bi,bj)
237			ENDDO
238	adcroft	1.12	ENDDO
239	jmc	1.26	ENDIF
240	adcroft	1.12
241			#ifdef ALLOW_OBCS
242	adcroft	1.14	IF (useOBCS) THEN
243	adcroft	1.12	DO i=1,sNx
244			C Northern boundary
245	jmc	1.63	IF (OB_Jn(i,bi,bj).NE.0) THEN
246			cg2d_b(i,OB_Jn(i,bi,bj),bi,bj)=0.
247			cg2d_x(i,OB_Jn(i,bi,bj),bi,bj)=0.
248	adcroft	1.12	ENDIF
249			C Southern boundary
250	jmc	1.63	IF (OB_Js(i,bi,bj).NE.0) THEN
251			cg2d_b(i,OB_Js(i,bi,bj),bi,bj)=0.
252			cg2d_x(i,OB_Js(i,bi,bj),bi,bj)=0.
253	adcroft	1.12	ENDIF
254			ENDDO
255			DO j=1,sNy
256			C Eastern boundary
257	jmc	1.63	IF (OB_Ie(j,bi,bj).NE.0) THEN
258			cg2d_b(OB_Ie(j,bi,bj),j,bi,bj)=0.
259			cg2d_x(OB_Ie(j,bi,bj),j,bi,bj)=0.
260	adcroft	1.12	ENDIF
261			C Western boundary
262	jmc	1.63	IF (OB_Iw(j,bi,bj).NE.0) THEN
263			cg2d_b(OB_Iw(j,bi,bj),j,bi,bj)=0.
264			cg2d_x(OB_Iw(j,bi,bj),j,bi,bj)=0.
265	adcroft	1.12	ENDIF
266			ENDDO
267			ENDIF
268	jmc	1.49	#endif /* ALLOW_OBCS */
269			C- end bi,bj loops
270	adcroft	1.12	ENDDO
271			ENDDO
272
273	edhill	1.42	#ifdef ALLOW_DEBUG
274	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
275	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
276			& myThid)
277	adcroft	1.24	ENDIF
278	adcroft	1.23	#endif
279	jmc	1.61	IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
280			WRITE(sufx,'(I10.10)') myIter
281	jmc	1.67	CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
282	jmc	1.61	ENDIF
283	adcroft	1.12
284	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
285			C-- gradient solver.
286	adcroft	1.22	C see CG2D.h for the interface to this routine.
287			firstResidual=0.
288			lastResidual=0.
289	adcroft	1.19	numIters=cg2dMaxIters
290	jmc	1.50	c CALL TIMER_START('CG2D [SOLVE_FOR_PRESSURE]',myThid)
291	heimbach	1.56	#ifdef ALLOW_CG2D_NSA
292			C-- Call the not-self-adjoint version of cg2d
293			CALL CG2D_NSA(
294			U cg2d_b,
295			U cg2d_x,
296			O firstResidual,
297			O lastResidual,
298			U numIters,
299			I myThid )
300			#else /* not ALLOW_CG2D_NSA = default */
301	mlosch	1.69	#ifdef ALLOW_SRCG
302			IF ( useSRCGSolver ) THEN
303			C-- Call the single reduce CG solver
304			CALL CG2D_SR(
305	adcroft	1.22	U cg2d_b,
306	cnh	1.6	U cg2d_x,
307	adcroft	1.22	O firstResidual,
308			O lastResidual,
309	adcroft	1.19	U numIters,
310	cnh	1.1	I myThid )
311	mlosch	1.69	ELSE
312			#else
313			IF (.TRUE.) THEN
314			C-- Call the default CG solver
315			#endif /* ALLOW_SRCG */
316			CALL CG2D(
317			U cg2d_b,
318			U cg2d_x,
319			O firstResidual,
320			O lastResidual,
321			U numIters,
322			I myThid )
323			ENDIF
324	heimbach	1.56	#endif /* ALLOW_CG2D_NSA */
325	jmc	1.67	_EXCH_XY_RL( cg2d_x, myThid )
326	jmc	1.50	c CALL TIMER_STOP ('CG2D [SOLVE_FOR_PRESSURE]',myThid)
327	adcroft	1.23
328	edhill	1.42	#ifdef ALLOW_DEBUG
329	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
330	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
331			& myThid)
332	adcroft	1.24	ENDIF
333	adcroft	1.23	#endif
334	cnh	1.1
335	jmc	1.32	C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
336	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
337	jmc	1.45	& ) THEN
338	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
339			_BEGIN_MASTER( myThid )
340			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
341			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
342			WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
343			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
344			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
345			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
346	edhill	1.43	_END_MASTER( myThid )
347	heimbach	1.38	ENDIF
348	jmc	1.32	ENDIF
349	jmc	1.17
350			C-- Transfert the 2D-solution to "etaN" :
351			DO bj=myByLo(myThid),myByHi(myThid)
352			DO bi=myBxLo(myThid),myBxHi(myThid)
353			DO j=1-OLy,sNy+OLy
354			DO i=1-OLx,sNx+OLx
355	jmc	1.18	etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
356	jmc	1.17	ENDDO
357			ENDDO
358			ENDDO
359			ENDDO
360	adcroft	1.10
361	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
362	jmc	1.53	IF ( use3Dsolver ) THEN
363	jmc	1.67	IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
364			WRITE(sufx,'(I10.10)') myIter
365			CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
366			ENDIF
367	adcroft	1.9
368			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
369			C see CG3D.h for the interface to this routine.
370			DO bj=myByLo(myThid),myByHi(myThid)
371			DO bi=myBxLo(myThid),myBxHi(myThid)
372			DO j=1,sNy+1
373			DO i=1,sNx+1
374	jmc	1.18	uf(i,j)=-_recip_dxC(i,j,bi,bj)*
375	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
376	jmc	1.18	vf(i,j)=-_recip_dyC(i,j,bi,bj)*
377	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
378			ENDDO
379			ENDDO
380
381	adcroft	1.12	#ifdef ALLOW_OBCS
382	adcroft	1.14	IF (useOBCS) THEN
383	adcroft	1.9	DO i=1,sNx+1
384			C Northern boundary
385	jmc	1.63	IF (OB_Jn(i,bi,bj).NE.0) THEN
386			vf(i,OB_Jn(i,bi,bj))=0.
387	adcroft	1.9	ENDIF
388			C Southern boundary
389	jmc	1.63	IF (OB_Js(i,bi,bj).NE.0) THEN
390			vf(i,OB_Js(i,bi,bj)+1)=0.
391	adcroft	1.9	ENDIF
392			ENDDO
393			DO j=1,sNy+1
394			C Eastern boundary
395	jmc	1.63	IF (OB_Ie(j,bi,bj).NE.0) THEN
396			uf(OB_Ie(j,bi,bj),j)=0.
397	adcroft	1.9	ENDIF
398			C Western boundary
399	jmc	1.63	IF (OB_Iw(j,bi,bj).NE.0) THEN
400			uf(OB_Iw(j,bi,bj)+1,J)=0.
401	adcroft	1.9	ENDIF
402			ENDDO
403			ENDIF
404	jmc	1.49	#endif /* ALLOW_OBCS */
405	adcroft	1.9
406	jmc	1.58	IF ( usingZCoords ) THEN
407	jmc	1.51	C- Z coordinate: assume surface @ level k=1
408	jmc	1.58	tmpFac = freeSurfFac*deepFac2F(1)
409	jmc	1.51	ELSE
410			C- Other than Z coordinate: no assumption on surface level index
411	jmc	1.58	tmpFac = 0.
412	jmc	1.51	DO j=1,sNy
413			DO i=1,sNx
414			ks = ksurfC(i,j,bi,bj)
415			IF ( ks.LE.Nr ) THEN
416			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
417			& +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
418	jmc	1.58	& _rA(i,j,bi,bj)deepFac2F(ks)/deltaTmom
419	jmc	1.51	ENDIF
420			ENDDO
421			ENDDO
422			ENDIF
423	jmc	1.63	k=1
424	jmc	1.51	kp1 = MIN(k+1,Nr)
425	jmc	1.58	wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
426			IF (k.GE.Nr) wFacKp = 0.
427	adcroft	1.12	DO j=1,sNy
428			DO i=1,sNx
429	jmc	1.51	cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
430	jmc	1.63	& +drF(k)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
431			& -drF(k)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
432			& +drF(k)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
433			& -drF(k)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
434	jmc	1.51	& +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
435	jmc	1.58	& -wVel(i,j,kp1,bi,bj)*wFacKp
436	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
437			ENDDO
438			ENDDO
439	jmc	1.63	DO k=2,Nr
440	jmc	1.51	kp1 = MIN(k+1,Nr)
441	jmc	1.58	C- deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
442			C both appear in wVel term, but at 2 different levels
443			wFacKm = deepFac2F( k )*rhoFacF( k )
444			wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
445			IF (k.GE.Nr) wFacKp = 0.
446	adcroft	1.9	DO j=1,sNy
447			DO i=1,sNx
448			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
449	jmc	1.63	& +drF(k)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
450			& -drF(k)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
451			& +drF(k)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
452			& -drF(k)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
453	jmc	1.58	& +( wVel(i,j, k ,bi,bj)wFacKmmaskC(i,j,k-1,bi,bj)
454			& -wVel(i,j,kp1,bi,bj)*wFacKp
455	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
456
457	adcroft	1.9	ENDDO
458			ENDDO
459			ENDDO
460	adcroft	1.12
461			#ifdef ALLOW_OBCS
462	adcroft	1.14	IF (useOBCS) THEN
463	jmc	1.63	DO k=1,Nr
464	adcroft	1.12	DO i=1,sNx
465			C Northern boundary
466	jmc	1.63	IF (OB_Jn(i,bi,bj).NE.0) THEN
467			cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj)=0.
468	adcroft	1.12	ENDIF
469			C Southern boundary
470	jmc	1.63	IF (OB_Js(i,bi,bj).NE.0) THEN
471			cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj)=0.
472	adcroft	1.12	ENDIF
473			ENDDO
474			DO j=1,sNy
475			C Eastern boundary
476	jmc	1.63	IF (OB_Ie(j,bi,bj).NE.0) THEN
477			cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj)=0.
478	adcroft	1.12	ENDIF
479			C Western boundary
480	jmc	1.63	IF (OB_Iw(j,bi,bj).NE.0) THEN
481			cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj)=0.
482	adcroft	1.12	ENDIF
483			ENDDO
484			ENDDO
485			ENDIF
486	jmc	1.49	#endif /* ALLOW_OBCS */
487			C- end bi,bj loops
488			ENDDO
489			ENDDO
490	adcroft	1.9
491	jmc	1.67	#ifdef ALLOW_DEBUG
492			IF ( debugLevel .GE. debLevB ) THEN
493			CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
494			& myThid)
495			ENDIF
496			#endif
497			IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
498			WRITE(sufx,'(I10.10)') myIter
499			CALL WRITE_FLD_XYZ_RL( 'cg3d_b.',sufx, cg3d_b, myIter, myThid )
500			ENDIF
501
502	adcroft	1.25	firstResidual=0.
503			lastResidual=0.
504	jmc	1.49	numIters=cg3dMaxIters
505	jmc	1.50	CALL TIMER_START('CG3D [SOLVE_FOR_PRESSURE]',myThid)
506	adcroft	1.25	CALL CG3D(
507			U cg3d_b,
508			U phi_nh,
509			O firstResidual,
510			O lastResidual,
511			U numIters,
512			I myThid )
513	jmc	1.67	_EXCH_XYZ_RL( phi_nh, myThid )
514	jmc	1.50	CALL TIMER_STOP ('CG3D [SOLVE_FOR_PRESSURE]',myThid)
515	adcroft	1.25
516	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
517	jmc	1.45	& ) THEN
518	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
519			_BEGIN_MASTER( myThid )
520			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
521			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
522			WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
523			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
524			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
525			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
526	edhill	1.43	_END_MASTER( myThid )
527	heimbach	1.38	ENDIF
528	mlosch	1.37	ENDIF
529	adcroft	1.9
530	jmc	1.49	C-- Update surface pressure (account for NH-p @ surface level) and NH pressure:
531			IF ( zeroPsNH ) THEN
532	jmc	1.67	IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
533			WRITE(sufx,'(I10.10)') myIter
534			CALL WRITE_FLD_XYZ_RL( 'cg3d_x.',sufx,phi_nh, myIter, myThid )
535			ENDIF
536	jmc	1.49	DO bj=myByLo(myThid),myByHi(myThid)
537			DO bi=myBxLo(myThid),myBxHi(myThid)
538
539			IF ( usingZCoords ) THEN
540			C- Z coordinate: assume surface @ level k=1
541			DO k=2,Nr
542			DO j=1-OLy,sNy+OLy
543			DO i=1-OLx,sNx+OLx
544			phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
545			& - phi_nh(i,j,1,bi,bj)
546			ENDDO
547			ENDDO
548			ENDDO
549			DO j=1-OLy,sNy+OLy
550			DO i=1-OLx,sNx+OLx
551			etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
552			& *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
553			phi_nh(i,j,1,bi,bj) = 0.
554			ENDDO
555			ENDDO
556			ELSE
557			C- Other than Z coordinate: no assumption on surface level index
558			DO j=1-OLy,sNy+OLy
559			DO i=1-OLx,sNx+OLx
560			ks = ksurfC(i,j,bi,bj)
561			IF ( ks.LE.Nr ) THEN
562			etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
563			& *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
564			DO k=Nr,1,-1
565			phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
566			& - phi_nh(i,j,ks,bi,bj)
567			ENDDO
568			ENDIF
569			ENDDO
570			ENDDO
571			ENDIF
572
573			ENDDO
574			ENDDO
575	adcroft	1.9	ENDIF
576	jmc	1.49
577			ENDIF
578			#endif /* ALLOW_NONHYDROSTATIC */
579	cnh	1.1
580	heimbach	1.60	#ifdef ALLOW_SHOWFLOPS
581			CALL SHOWFLOPS_INSOLVE( myThid)
582	ce107	1.52	#endif
583	heimbach	1.60
584	cnh	1.1	RETURN
585			END