model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.71 2009/11/29 03:17:05 jmc 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, ioUnit
      CHARACTER*10 sufx
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_NONHYDROSTATIC
      INTEGER kp1
      _RL     wFacKm, wFacKp
      LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
      _RL     tmpVar(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _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 = use3Dsolver .AND. exactConserv
c    &                         .AND. select_rStar.EQ.0
        zeroMeanPnh = .FALSE.
c       zeroMeanPnh = use3Dsolver .AND. select_rStar.NE.0
c       oldFreeSurfTerm = use3Dsolver .AND. select_rStar.EQ.0
c    &                                .AND. .NOT.zeroPsNH
        oldFreeSurfTerm = use3Dsolver .AND. .NOT.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

      IF ( myIter.EQ.1+nIter0 .AND. debugLevel .GE. debLevA ) THEN
        _BEGIN_MASTER( myThid )
        ioUnit = standardMessageUnit
        WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
     &       ' putPmEinXvector =', putPmEinXvector
        CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
#ifdef ALLOW_NONHYDROSTATIC
        WRITE(msgBuf,'(A,2(A,L5))') 'SOLVE_FOR_PRESSURE:',
     &       ' zeroPsNH=', zeroPsNH, ' , zeroMeanPnh=', zeroMeanPnh
        CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
        WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
     &       ' oldFreeSurfTerm =', oldFreeSurfTerm
        CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
#endif
        _END_MASTER( myThid )
      ENDIF

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 ( oldFreeSurfTerm ) 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)

C--   Update or Add free-surface contribution to cg3d_b:
c        IF ( select_rStar.EQ.0 .AND. exactConserv ) THEN
         IF ( select_rStar.EQ.0 .AND. .NOT.oldFreeSurfTerm ) THEN
           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)-etaH(i,j,bi,bj))
     &                              *_rA(i,j,bi,bj)*deepFac2F(ks)
     &                              /deltaTMom/deltaTfreesurf
              ENDIF
            ENDDO
           ENDDO
#ifdef NONLIN_FRSURF
         ELSEIF ( select_rStar.NE.0 ) THEN
           tmpFac = 0.
           DO j=1,sNy
            DO i=1,sNx
              ks = ksurfC(i,j,bi,bj)
              tmpVar(i,j) = freeSurfFac
     &                    *( etaN(i,j,bi,bj) - etaH(i,j,bi,bj) )
     &                    *_rA(i,j,bi,bj)*deepFac2F(ks)
     &                    /deltaTMom/deltaTfreesurf
     &                    *recip_Rcol(i,j,bi,bj)
            ENDDO
           ENDDO
           DO k=1,Nr
            DO j=1,sNy
             DO i=1,sNx
              cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
     &          + tmpVar(i,j)*drF(k)*h0FacC(i,j,k,bi,bj)
             ENDDO
            ENDDO
           ENDDO
#endif /* NONLIN_FRSURF */
         ELSEIF ( 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

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Finish updating cg3d_b: 1) increment in horiz velocity due to new cg2d_x
C                             2) add vertical velocity contribution.
         DO j=1,sNy+1
          DO i=1,sNx+1
           uf(i,j) = -_recip_dxC(i,j,bi,bj)
     &             * implicSurfPress*implicDiv2DFlow
     &             *(cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
           vf(i,j) = -_recip_dyC(i,j,bi,bj)
     &             * implicSurfPress*implicDiv2DFlow
     &             *(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)
     &      vf(i,OB_Jn(i,bi,bj)) = 0.
C Southern boundary
           IF (OB_Js(i,bi,bj).NE.0)
     &      vf(i,OB_Js(i,bi,bj)+1) = 0.
          ENDDO
          DO j=1,sNy+1
C Eastern boundary
           IF (OB_Ie(j,bi,bj).NE.0)
     &      uf(OB_Ie(j,bi,bj),j) = 0.
C Western boundary
           IF (OB_Iw(j,bi,bj).NE.0)
     &      uf(OB_Iw(j,bi,bj)+1,j) = 0.
          ENDDO
         ENDIF
#endif /* ALLOW_OBCS */

C Note: with implicDiv2DFlow < 1, wVel contribution to cg3d_b is similar to
C       uVel,vVel contribution to cg2d_b when exactConserv=T, since wVel is
C       always recomputed from continuity eq (like eta when exactConserv=T)
         k=1
         kp1 = MIN(k+1,Nr)
         wFacKp = implicDiv2DFlow*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 = implicDiv2DFlow*deepFac2F( k )*rhoFacF( k )
          wFacKp = implicDiv2DFlow*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)
     &       cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj) = 0.
C Southern boundary
            IF (OB_Js(i,bi,bj).NE.0)
     &       cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj) = 0.
           ENDDO
           DO j=1,sNy
C Eastern boundary
            IF (OB_Ie(j,bi,bj).NE.0)
     &       cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj) = 0.
C Western boundary
            IF (OB_Iw(j,bi,bj).NE.0)
     &       cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj) = 0.
           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           myIter, 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 .OR. zeroMeanPnh ) 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 ( zeroPsNH .AND. usingZCoords ) THEN
C-       Z coordinate: assume surface @ level k=1
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             tmpVar(i,j) = phi_nh(i,j,1,bi,bj)
           ENDDO
          ENDDO
         ELSEIF ( zeroPsNH ) THEN
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
             tmpVar(i,j) = phi_nh(i,j,ks,bi,bj)
            ELSE
             tmpVar(i,j) = 0.
            ENDIF
           ENDDO
          ENDDO
#ifdef NONLIN_FRSURF
         ELSE
C        zeroMeanPnh : transfert vertical average of P_NH to EtaN
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             tmpVar(i,j) = 0.
           ENDDO
          ENDDO
          DO k=1,Nr
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             tmpVar(i,j) = tmpVar(i,j)
     &         + phi_nh(i,j,k,bi,bj)*drF(k)*h0FacC(i,j,k,bi,bj)
            ENDDO
           ENDDO
          ENDDO
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             tmpVar(i,j) = tmpVar(i,j)*recip_Rcol(i,j,bi,bj)
           ENDDO
          ENDDO
#endif /* NONLIN_FRSURF */
         ENDIF
         DO k=1,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)
     &                            - tmpVar(i,j)
     &                            )*maskC(i,j,k,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) + tmpVar(i,j) )
          ENDDO
         ENDDO

        ENDDO
       ENDDO
      ENDIF

      ENDIF
#endif /* ALLOW_NONHYDROSTATIC */

#ifdef ALLOW_SHOWFLOPS
      CALL SHOWFLOPS_INSOLVE( myThid)
#endif

      RETURN
      END
1	jmc	1.72	C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.71 2009/11/29 03:17:05 jmc 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.71	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.71	INTEGER numIters, ks, ioUnit
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.71	LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
69			_RL tmpVar(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70	jmc	1.63	_RL uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
71			_RL vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
72			#else
73			_RL cg3d_b(1)
74	jmc	1.49	#endif
75	cnh	1.27	CEOP
76	jmc	1.17
77	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
78	jmc	1.68	zeroPsNH = .FALSE.
79	jmc	1.71	c zeroPsNH = use3Dsolver .AND. exactConserv
80			c & .AND. select_rStar.EQ.0
81			zeroMeanPnh = .FALSE.
82			c zeroMeanPnh = use3Dsolver .AND. select_rStar.NE.0
83	jmc	1.72	c oldFreeSurfTerm = use3Dsolver .AND. select_rStar.EQ.0
84			c & .AND. .NOT.zeroPsNH
85			oldFreeSurfTerm = use3Dsolver .AND. .NOT.exactConserv
86	jmc	1.63	#else
87			cg3d_b(1) = 0.
88	jmc	1.49	#endif
89
90	jmc	1.58	C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
91			C anelastic (always Z-coordinate):
92			C 1) assume that rhoFacF(1)=1 (and ksurf == 1);
93			C (this reduces the number of lines of code to modify)
94			C 2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
95			C (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
96			C => 2 factors cancel in elliptic eq. for Phi_s ,
97			C but 1rst factor(a) remains in RHS cg2d_b.
98
99	jmc	1.47	C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
100			C instead of simply etaN ; This can speed-up the solver convergence in
101			C the case where \|Global_mean_PmE\| is large.
102			putPmEinXvector = .FALSE.
103	jmc	1.64	c putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
104	jmc	1.47
105	jmc	1.71	IF ( myIter.EQ.1+nIter0 .AND. debugLevel .GE. debLevA ) THEN
106			_BEGIN_MASTER( myThid )
107			ioUnit = standardMessageUnit
108			WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
109			& ' putPmEinXvector =', putPmEinXvector
110			CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
111			#ifdef ALLOW_NONHYDROSTATIC
112			WRITE(msgBuf,'(A,2(A,L5))') 'SOLVE_FOR_PRESSURE:',
113			& ' zeroPsNH=', zeroPsNH, ' , zeroMeanPnh=', zeroMeanPnh
114			CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
115			WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
116			& ' oldFreeSurfTerm =', oldFreeSurfTerm
117			CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
118			#endif
119			_END_MASTER( myThid )
120			ENDIF
121
122	jmc	1.17	C-- Save previous solution & Initialise Vector solution and source term :
123	jmc	1.47	sumEmP = 0.
124	jmc	1.17	DO bj=myByLo(myThid),myByHi(myThid)
125			DO bi=myBxLo(myThid),myBxHi(myThid)
126			DO j=1-OLy,sNy+OLy
127			DO i=1-OLx,sNx+OLx
128	edhill	1.40	#ifdef ALLOW_CD_CODE
129	jmc	1.17	etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
130	jmc	1.26	#endif
131	jmc	1.18	cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
132	jmc	1.17	cg2d_b(i,j,bi,bj) = 0.
133			ENDDO
134			ENDDO
135	jmc	1.64	IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
136	jmc	1.62	tmpFac = freeSurfFac*mass2rUnit
137	jmc	1.58	IF (exactConserv)
138	jmc	1.62	& tmpFac = freeSurfFacmass2rUnitimplicDiv2DFlow
139	jmc	1.29	DO j=1,sNy
140			DO i=1,sNx
141	jmc	1.58	cg2d_b(i,j,bi,bj) =
142	jmc	1.29	& tmpFac_rA(i,j,bi,bj)EmPmR(i,j,bi,bj)/deltaTMom
143			ENDDO
144			ENDDO
145			ENDIF
146	jmc	1.47	IF ( putPmEinXvector ) THEN
147	jmc	1.65	tileEmP(bi,bj) = 0.
148	jmc	1.47	DO j=1,sNy
149			DO i=1,sNx
150	jmc	1.67	tileEmP(bi,bj) = tileEmP(bi,bj)
151	jmc	1.65	& + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
152			& *maskH(i,j,bi,bj)
153	jmc	1.47	ENDDO
154			ENDDO
155			ENDIF
156	jmc	1.17	ENDDO
157			ENDDO
158	jmc	1.47	IF ( putPmEinXvector ) THEN
159	jmc	1.65	CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
160	jmc	1.47	ENDIF
161	adcroft	1.12
162			DO bj=myByLo(myThid),myByHi(myThid)
163			DO bi=myBxLo(myThid),myBxHi(myThid)
164	jmc	1.47	IF ( putPmEinXvector ) THEN
165			tmpFac = 0.
166	jmc	1.62	IF (globalArea.GT.0.) tmpFac =
167			& freeSurfFacdeltaTfreesurfmass2rUnit*sumEmP/globalArea
168	jmc	1.47	DO j=1,sNy
169			DO i=1,sNx
170			cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
171			& - tmpFac*Bo_surf(i,j,bi,bj)
172			ENDDO
173			ENDDO
174			ENDIF
175	jmc	1.58	C- RHS: similar to the divergence of the vertically integrated mass transport:
176			C del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] } / deltaT
177	jmc	1.63	DO k=Nr,1,-1
178	adcroft	1.12	CALL CALC_DIV_GHAT(
179	jmc	1.63	I bi,bj,k,
180			U cg2d_b, cg3d_b,
181			I myThid )
182	adcroft	1.12	ENDDO
183			ENDDO
184			ENDDO
185	cnh	1.4
186	adcroft	1.12	C-- Add source term arising from w=d/dt (p_s + p_nh)
187			DO bj=myByLo(myThid),myByHi(myThid)
188			DO bi=myBxLo(myThid),myBxHi(myThid)
189	adcroft	1.13	#ifdef ALLOW_NONHYDROSTATIC
190	jmc	1.70	C-- Add EmPmR contribution to top level cg3d_b:
191			C (has been done for cg2d_b ; and addMass was added by CALC_DIV_GHAT)
192			IF ( use3Dsolver .AND.
193			& useRealFreshWaterFlux.AND.fluidIsWater ) THEN
194			tmpFac = freeSurfFac*mass2rUnit
195			IF (exactConserv)
196			& tmpFac = freeSurfFacmass2rUnitimplicDiv2DFlow
197			ks = 1
198			IF ( usingPCoords ) ks = Nr
199			DO j=1,sNy
200			DO i=1,sNx
201			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
202			& + tmpFac_rA(i,j,bi,bj)EmPmR(i,j,bi,bj)/deltaTMom
203			ENDDO
204			ENDDO
205			ENDIF
206	jmc	1.71
207			IF ( oldFreeSurfTerm ) THEN
208	jmc	1.28	DO j=1,sNy
209			DO i=1,sNx
210	jmc	1.51	ks = ksurfC(i,j,bi,bj)
211			IF ( ks.LE.Nr ) THEN
212			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
213	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
214			& /deltaTMom/deltaTfreesurf
215	jmc	1.28	& *( etaN(i,j,bi,bj)
216	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
217	jmc	1.51	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
218	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
219			& /deltaTMom/deltaTfreesurf
220	jmc	1.28	& *( etaN(i,j,bi,bj)
221	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
222	jmc	1.51	ENDIF
223	jmc	1.28	ENDDO
224	adcroft	1.12	ENDDO
225	jmc	1.28	ELSEIF ( exactConserv ) THEN
226	adcroft	1.13	#else
227	jmc	1.26	IF ( exactConserv ) THEN
228	edhill	1.39	#endif /* ALLOW_NONHYDROSTATIC */
229	jmc	1.26	DO j=1,sNy
230			DO i=1,sNx
231	jmc	1.58	ks = ksurfC(i,j,bi,bj)
232	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
233	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
234			& /deltaTMom/deltaTfreesurf
235	jmc	1.26	& * etaH(i,j,bi,bj)
236			ENDDO
237			ENDDO
238			ELSE
239			DO j=1,sNy
240			DO i=1,sNx
241	jmc	1.58	ks = ksurfC(i,j,bi,bj)
242	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
243	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
244			& /deltaTMom/deltaTfreesurf
245	jmc	1.26	& * etaN(i,j,bi,bj)
246			ENDDO
247	adcroft	1.12	ENDDO
248	jmc	1.26	ENDIF
249	adcroft	1.12
250			#ifdef ALLOW_OBCS
251	adcroft	1.14	IF (useOBCS) THEN
252	adcroft	1.12	DO i=1,sNx
253			C Northern boundary
254	jmc	1.63	IF (OB_Jn(i,bi,bj).NE.0) THEN
255			cg2d_b(i,OB_Jn(i,bi,bj),bi,bj)=0.
256			cg2d_x(i,OB_Jn(i,bi,bj),bi,bj)=0.
257	adcroft	1.12	ENDIF
258			C Southern boundary
259	jmc	1.63	IF (OB_Js(i,bi,bj).NE.0) THEN
260			cg2d_b(i,OB_Js(i,bi,bj),bi,bj)=0.
261			cg2d_x(i,OB_Js(i,bi,bj),bi,bj)=0.
262	adcroft	1.12	ENDIF
263			ENDDO
264			DO j=1,sNy
265			C Eastern boundary
266	jmc	1.63	IF (OB_Ie(j,bi,bj).NE.0) THEN
267			cg2d_b(OB_Ie(j,bi,bj),j,bi,bj)=0.
268			cg2d_x(OB_Ie(j,bi,bj),j,bi,bj)=0.
269	adcroft	1.12	ENDIF
270			C Western boundary
271	jmc	1.63	IF (OB_Iw(j,bi,bj).NE.0) THEN
272			cg2d_b(OB_Iw(j,bi,bj),j,bi,bj)=0.
273			cg2d_x(OB_Iw(j,bi,bj),j,bi,bj)=0.
274	adcroft	1.12	ENDIF
275			ENDDO
276			ENDIF
277	jmc	1.49	#endif /* ALLOW_OBCS */
278			C- end bi,bj loops
279	adcroft	1.12	ENDDO
280			ENDDO
281
282	edhill	1.42	#ifdef ALLOW_DEBUG
283	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
284	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
285			& myThid)
286	adcroft	1.24	ENDIF
287	adcroft	1.23	#endif
288	jmc	1.61	IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
289			WRITE(sufx,'(I10.10)') myIter
290	jmc	1.67	CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
291	jmc	1.61	ENDIF
292	adcroft	1.12
293	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
294			C-- gradient solver.
295	adcroft	1.22	C see CG2D.h for the interface to this routine.
296			firstResidual=0.
297			lastResidual=0.
298	adcroft	1.19	numIters=cg2dMaxIters
299	jmc	1.50	c CALL TIMER_START('CG2D [SOLVE_FOR_PRESSURE]',myThid)
300	heimbach	1.56	#ifdef ALLOW_CG2D_NSA
301			C-- Call the not-self-adjoint version of cg2d
302			CALL CG2D_NSA(
303			U cg2d_b,
304			U cg2d_x,
305			O firstResidual,
306			O lastResidual,
307			U numIters,
308			I myThid )
309			#else /* not ALLOW_CG2D_NSA = default */
310	mlosch	1.69	#ifdef ALLOW_SRCG
311			IF ( useSRCGSolver ) THEN
312			C-- Call the single reduce CG solver
313			CALL CG2D_SR(
314	adcroft	1.22	U cg2d_b,
315	cnh	1.6	U cg2d_x,
316	adcroft	1.22	O firstResidual,
317			O lastResidual,
318	adcroft	1.19	U numIters,
319	cnh	1.1	I myThid )
320	mlosch	1.69	ELSE
321			#else
322			IF (.TRUE.) THEN
323			C-- Call the default CG solver
324			#endif /* ALLOW_SRCG */
325			CALL CG2D(
326			U cg2d_b,
327			U cg2d_x,
328			O firstResidual,
329			O lastResidual,
330			U numIters,
331			I myThid )
332			ENDIF
333	heimbach	1.56	#endif /* ALLOW_CG2D_NSA */
334	jmc	1.67	_EXCH_XY_RL( cg2d_x, myThid )
335	jmc	1.50	c CALL TIMER_STOP ('CG2D [SOLVE_FOR_PRESSURE]',myThid)
336	adcroft	1.23
337	edhill	1.42	#ifdef ALLOW_DEBUG
338	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
339	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
340			& myThid)
341	adcroft	1.24	ENDIF
342	adcroft	1.23	#endif
343	cnh	1.1
344	jmc	1.32	C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
345	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
346	jmc	1.45	& ) THEN
347	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
348			_BEGIN_MASTER( myThid )
349			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
350			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
351			WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
352			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
353			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
354			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
355	edhill	1.43	_END_MASTER( myThid )
356	heimbach	1.38	ENDIF
357	jmc	1.32	ENDIF
358	jmc	1.17
359			C-- Transfert the 2D-solution to "etaN" :
360			DO bj=myByLo(myThid),myByHi(myThid)
361			DO bi=myBxLo(myThid),myBxHi(myThid)
362			DO j=1-OLy,sNy+OLy
363			DO i=1-OLx,sNx+OLx
364	jmc	1.18	etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
365	jmc	1.17	ENDDO
366			ENDDO
367			ENDDO
368			ENDDO
369	adcroft	1.10
370	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
371	jmc	1.53	IF ( use3Dsolver ) THEN
372	jmc	1.67	IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
373			WRITE(sufx,'(I10.10)') myIter
374			CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
375			ENDIF
376	adcroft	1.9
377			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
378			C see CG3D.h for the interface to this routine.
379			DO bj=myByLo(myThid),myByHi(myThid)
380			DO bi=myBxLo(myThid),myBxHi(myThid)
381	jmc	1.71
382			C-- Update or Add free-surface contribution to cg3d_b:
383			c IF ( select_rStar.EQ.0 .AND. exactConserv ) THEN
384			IF ( select_rStar.EQ.0 .AND. .NOT.oldFreeSurfTerm ) THEN
385			tmpFac = 0.
386			DO j=1,sNy
387			DO i=1,sNx
388			ks = ksurfC(i,j,bi,bj)
389			IF ( ks.LE.Nr ) THEN
390			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
391			& +freeSurfFac*(etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
392			& _rA(i,j,bi,bj)deepFac2F(ks)
393			& /deltaTMom/deltaTfreesurf
394			ENDIF
395			ENDDO
396			ENDDO
397			#ifdef NONLIN_FRSURF
398			ELSEIF ( select_rStar.NE.0 ) THEN
399			tmpFac = 0.
400			DO j=1,sNy
401			DO i=1,sNx
402			ks = ksurfC(i,j,bi,bj)
403			tmpVar(i,j) = freeSurfFac
404			& *( etaN(i,j,bi,bj) - etaH(i,j,bi,bj) )
405			& _rA(i,j,bi,bj)deepFac2F(ks)
406			& /deltaTMom/deltaTfreesurf
407			& *recip_Rcol(i,j,bi,bj)
408			ENDDO
409			ENDDO
410			DO k=1,Nr
411			DO j=1,sNy
412			DO i=1,sNx
413			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
414			& + tmpVar(i,j)drF(k)h0FacC(i,j,k,bi,bj)
415			ENDDO
416			ENDDO
417			ENDDO
418			#endif /* NONLIN_FRSURF */
419			ELSEIF ( usingZCoords ) THEN
420			C- Z coordinate: assume surface @ level k=1
421			tmpFac = freeSurfFac*deepFac2F(1)
422			ELSE
423			C- Other than Z coordinate: no assumption on surface level index
424			tmpFac = 0.
425			DO j=1,sNy
426			DO i=1,sNx
427			ks = ksurfC(i,j,bi,bj)
428			IF ( ks.LE.Nr ) THEN
429			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
430			& +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
431			& _rA(i,j,bi,bj)deepFac2F(ks)/deltaTmom
432			ENDIF
433			ENDDO
434			ENDDO
435			ENDIF
436
437			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
438
439			C-- Finish updating cg3d_b: 1) increment in horiz velocity due to new cg2d_x
440			C 2) add vertical velocity contribution.
441	adcroft	1.9	DO j=1,sNy+1
442			DO i=1,sNx+1
443	jmc	1.71	uf(i,j) = -_recip_dxC(i,j,bi,bj)
444			& * implicSurfPress*implicDiv2DFlow
445			& *(cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
446			vf(i,j) = -_recip_dyC(i,j,bi,bj)
447			& * implicSurfPress*implicDiv2DFlow
448			& *(cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
449	adcroft	1.9	ENDDO
450			ENDDO
451
452	adcroft	1.12	#ifdef ALLOW_OBCS
453	adcroft	1.14	IF (useOBCS) THEN
454	adcroft	1.9	DO i=1,sNx+1
455			C Northern boundary
456	jmc	1.71	IF (OB_Jn(i,bi,bj).NE.0)
457			& vf(i,OB_Jn(i,bi,bj)) = 0.
458	adcroft	1.9	C Southern boundary
459	jmc	1.71	IF (OB_Js(i,bi,bj).NE.0)
460			& vf(i,OB_Js(i,bi,bj)+1) = 0.
461	adcroft	1.9	ENDDO
462			DO j=1,sNy+1
463			C Eastern boundary
464	jmc	1.71	IF (OB_Ie(j,bi,bj).NE.0)
465			& uf(OB_Ie(j,bi,bj),j) = 0.
466	adcroft	1.9	C Western boundary
467	jmc	1.71	IF (OB_Iw(j,bi,bj).NE.0)
468			& uf(OB_Iw(j,bi,bj)+1,j) = 0.
469	adcroft	1.9	ENDDO
470			ENDIF
471	jmc	1.49	#endif /* ALLOW_OBCS */
472	adcroft	1.9
473	jmc	1.71	C Note: with implicDiv2DFlow < 1, wVel contribution to cg3d_b is similar to
474			C uVel,vVel contribution to cg2d_b when exactConserv=T, since wVel is
475			C always recomputed from continuity eq (like eta when exactConserv=T)
476	jmc	1.63	k=1
477	jmc	1.51	kp1 = MIN(k+1,Nr)
478	jmc	1.71	wFacKp = implicDiv2DFlowdeepFac2F(kp1)rhoFacF(kp1)
479	jmc	1.58	IF (k.GE.Nr) wFacKp = 0.
480	adcroft	1.12	DO j=1,sNy
481			DO i=1,sNx
482	jmc	1.51	cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
483	jmc	1.63	& +drF(k)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
484			& -drF(k)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
485			& +drF(k)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
486			& -drF(k)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
487	jmc	1.51	& +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
488	jmc	1.58	& -wVel(i,j,kp1,bi,bj)*wFacKp
489	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
490			ENDDO
491			ENDDO
492	jmc	1.63	DO k=2,Nr
493	jmc	1.51	kp1 = MIN(k+1,Nr)
494	jmc	1.58	C- deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
495			C both appear in wVel term, but at 2 different levels
496	jmc	1.71	wFacKm = implicDiv2DFlowdeepFac2F( k )rhoFacF( k )
497			wFacKp = implicDiv2DFlowdeepFac2F(kp1)rhoFacF(kp1)
498	jmc	1.58	IF (k.GE.Nr) wFacKp = 0.
499	adcroft	1.9	DO j=1,sNy
500			DO i=1,sNx
501			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
502	jmc	1.63	& +drF(k)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
503			& -drF(k)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
504			& +drF(k)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
505			& -drF(k)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
506	jmc	1.58	& +( wVel(i,j, k ,bi,bj)wFacKmmaskC(i,j,k-1,bi,bj)
507			& -wVel(i,j,kp1,bi,bj)*wFacKp
508	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
509
510	adcroft	1.9	ENDDO
511			ENDDO
512			ENDDO
513	adcroft	1.12
514			#ifdef ALLOW_OBCS
515	adcroft	1.14	IF (useOBCS) THEN
516	jmc	1.63	DO k=1,Nr
517	jmc	1.71	DO i=1,sNx
518	adcroft	1.12	C Northern boundary
519	jmc	1.71	IF (OB_Jn(i,bi,bj).NE.0)
520			& cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj) = 0.
521	adcroft	1.12	C Southern boundary
522	jmc	1.71	IF (OB_Js(i,bi,bj).NE.0)
523			& cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj) = 0.
524			ENDDO
525			DO j=1,sNy
526	adcroft	1.12	C Eastern boundary
527	jmc	1.71	IF (OB_Ie(j,bi,bj).NE.0)
528			& cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj) = 0.
529	adcroft	1.12	C Western boundary
530	jmc	1.71	IF (OB_Iw(j,bi,bj).NE.0)
531			& cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj) = 0.
532			ENDDO
533	adcroft	1.12	ENDDO
534			ENDIF
535	jmc	1.49	#endif /* ALLOW_OBCS */
536			C- end bi,bj loops
537			ENDDO
538			ENDDO
539	adcroft	1.9
540	jmc	1.67	#ifdef ALLOW_DEBUG
541			IF ( debugLevel .GE. debLevB ) THEN
542			CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
543			& myThid)
544			ENDIF
545			#endif
546			IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
547			WRITE(sufx,'(I10.10)') myIter
548			CALL WRITE_FLD_XYZ_RL( 'cg3d_b.',sufx, cg3d_b, myIter, myThid )
549			ENDIF
550
551	adcroft	1.25	firstResidual=0.
552			lastResidual=0.
553	jmc	1.49	numIters=cg3dMaxIters
554	jmc	1.50	CALL TIMER_START('CG3D [SOLVE_FOR_PRESSURE]',myThid)
555	adcroft	1.25	CALL CG3D(
556			U cg3d_b,
557			U phi_nh,
558			O firstResidual,
559			O lastResidual,
560			U numIters,
561	jmc	1.71	I myIter, myThid )
562	jmc	1.67	_EXCH_XYZ_RL( phi_nh, myThid )
563	jmc	1.50	CALL TIMER_STOP ('CG3D [SOLVE_FOR_PRESSURE]',myThid)
564	adcroft	1.25
565	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
566	jmc	1.45	& ) THEN
567	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
568			_BEGIN_MASTER( myThid )
569			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
570			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
571			WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
572			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
573			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
574			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
575	edhill	1.43	_END_MASTER( myThid )
576	heimbach	1.38	ENDIF
577	mlosch	1.37	ENDIF
578	adcroft	1.9
579	jmc	1.49	C-- Update surface pressure (account for NH-p @ surface level) and NH pressure:
580	jmc	1.71	IF ( zeroPsNH .OR. zeroMeanPnh ) THEN
581	jmc	1.67	IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
582			WRITE(sufx,'(I10.10)') myIter
583	jmc	1.71	CALL WRITE_FLD_XYZ_RL( 'cg3d_x.',sufx, phi_nh, myIter, myThid )
584	jmc	1.67	ENDIF
585	jmc	1.49	DO bj=myByLo(myThid),myByHi(myThid)
586			DO bi=myBxLo(myThid),myBxHi(myThid)
587
588	jmc	1.71	IF ( zeroPsNH .AND. usingZCoords ) THEN
589	jmc	1.49	C- Z coordinate: assume surface @ level k=1
590	jmc	1.71	DO j=1-OLy,sNy+OLy
591			DO i=1-OLx,sNx+OLx
592			tmpVar(i,j) = phi_nh(i,j,1,bi,bj)
593			ENDDO
594			ENDDO
595			ELSEIF ( zeroPsNH ) THEN
596			C- Other than Z coordinate: no assumption on surface level index
597			DO j=1-OLy,sNy+OLy
598			DO i=1-OLx,sNx+OLx
599			ks = ksurfC(i,j,bi,bj)
600			IF ( ks.LE.Nr ) THEN
601			tmpVar(i,j) = phi_nh(i,j,ks,bi,bj)
602			ELSE
603			tmpVar(i,j) = 0.
604			ENDIF
605			ENDDO
606			ENDDO
607			#ifdef NONLIN_FRSURF
608			ELSE
609			C zeroMeanPnh : transfert vertical average of P_NH to EtaN
610			DO j=1-OLy,sNy+OLy
611			DO i=1-OLx,sNx+OLx
612			tmpVar(i,j) = 0.
613			ENDDO
614			ENDDO
615			DO k=1,Nr
616	jmc	1.49	DO j=1-OLy,sNy+OLy
617			DO i=1-OLx,sNx+OLx
618	jmc	1.71	tmpVar(i,j) = tmpVar(i,j)
619			& + phi_nh(i,j,k,bi,bj)drF(k)h0FacC(i,j,k,bi,bj)
620	jmc	1.49	ENDDO
621			ENDDO
622			ENDDO
623			DO j=1-OLy,sNy+OLy
624			DO i=1-OLx,sNx+OLx
625	jmc	1.71	tmpVar(i,j) = tmpVar(i,j)*recip_Rcol(i,j,bi,bj)
626	jmc	1.49	ENDDO
627			ENDDO
628	jmc	1.71	#endif /* NONLIN_FRSURF */
629			ENDIF
630			DO k=1,Nr
631	jmc	1.49	DO j=1-OLy,sNy+OLy
632			DO i=1-OLx,sNx+OLx
633	jmc	1.71	phi_nh(i,j,k,bi,bj) = ( phi_nh(i,j,k,bi,bj)
634			& - tmpVar(i,j)
635			& )*maskC(i,j,k,bi,bj)
636	jmc	1.49	ENDDO
637			ENDDO
638	jmc	1.71	ENDDO
639			DO j=1-OLy,sNy+OLy
640			DO i=1-OLx,sNx+OLx
641			etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
642			& *( cg2d_x(i,j,bi,bj) + tmpVar(i,j) )
643			ENDDO
644			ENDDO
645	jmc	1.49
646			ENDDO
647			ENDDO
648	adcroft	1.9	ENDIF
649	jmc	1.49
650			ENDIF
651			#endif /* ALLOW_NONHYDROSTATIC */
652	cnh	1.1
653	heimbach	1.60	#ifdef ALLOW_SHOWFLOPS
654			CALL SHOWFLOPS_INSOLVE( myThid)
655	ce107	1.52	#endif
656	heimbach	1.60
657	cnh	1.1	RETURN
658			END