model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.65 2009/04/22 01:30:00 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
      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
c       zeroPsNH = .FALSE.
        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
        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_RS( '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 */
      CALL CG2D(
     U           cg2d_b,
     U           cg2d_x,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
#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

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

      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
       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.66	C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.65 2009/04/22 01:30:00 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.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			c zeroPsNH = .FALSE.
78			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.65	tileEmP(bi,bj) = tileEmP(bi,bj)
127			& + 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.53	IF ( use3Dsolver .AND. zeroPsNH ) THEN
167	jmc	1.49	DO j=1,sNy
168			DO i=1,sNx
169	jmc	1.51	ks = ksurfC(i,j,bi,bj)
170			IF ( ks.LE.Nr ) THEN
171			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
172	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
173			& /deltaTMom/deltaTfreesurf
174	jmc	1.49	& * etaH(i,j,bi,bj)
175	jmc	1.51	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
176	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
177			& /deltaTMom/deltaTfreesurf
178	jmc	1.49	& * etaH(i,j,bi,bj)
179	jmc	1.51	ENDIF
180	jmc	1.49	ENDDO
181			ENDDO
182	jmc	1.53	ELSEIF ( use3Dsolver ) THEN
183	jmc	1.28	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.28	& *( etaN(i,j,bi,bj)
191	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
192	jmc	1.51	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
193	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
194			& /deltaTMom/deltaTfreesurf
195	jmc	1.28	& *( etaN(i,j,bi,bj)
196	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
197	jmc	1.51	ENDIF
198	jmc	1.28	ENDDO
199	adcroft	1.12	ENDDO
200	jmc	1.28	ELSEIF ( exactConserv ) THEN
201	adcroft	1.13	#else
202	jmc	1.26	IF ( exactConserv ) THEN
203	edhill	1.39	#endif /* ALLOW_NONHYDROSTATIC */
204	jmc	1.26	DO j=1,sNy
205			DO i=1,sNx
206	jmc	1.58	ks = ksurfC(i,j,bi,bj)
207	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
208	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
209			& /deltaTMom/deltaTfreesurf
210	jmc	1.26	& * etaH(i,j,bi,bj)
211			ENDDO
212			ENDDO
213			ELSE
214			DO j=1,sNy
215			DO i=1,sNx
216	jmc	1.58	ks = ksurfC(i,j,bi,bj)
217	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
218	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
219			& /deltaTMom/deltaTfreesurf
220	jmc	1.26	& * etaN(i,j,bi,bj)
221			ENDDO
222	adcroft	1.12	ENDDO
223	jmc	1.26	ENDIF
224	adcroft	1.12
225			#ifdef ALLOW_OBCS
226	adcroft	1.14	IF (useOBCS) THEN
227	adcroft	1.12	DO i=1,sNx
228			C Northern boundary
229	jmc	1.63	IF (OB_Jn(i,bi,bj).NE.0) THEN
230			cg2d_b(i,OB_Jn(i,bi,bj),bi,bj)=0.
231			cg2d_x(i,OB_Jn(i,bi,bj),bi,bj)=0.
232	adcroft	1.12	ENDIF
233			C Southern boundary
234	jmc	1.63	IF (OB_Js(i,bi,bj).NE.0) THEN
235			cg2d_b(i,OB_Js(i,bi,bj),bi,bj)=0.
236			cg2d_x(i,OB_Js(i,bi,bj),bi,bj)=0.
237	adcroft	1.12	ENDIF
238			ENDDO
239			DO j=1,sNy
240			C Eastern boundary
241	jmc	1.63	IF (OB_Ie(j,bi,bj).NE.0) THEN
242			cg2d_b(OB_Ie(j,bi,bj),j,bi,bj)=0.
243			cg2d_x(OB_Ie(j,bi,bj),j,bi,bj)=0.
244	adcroft	1.12	ENDIF
245			C Western boundary
246	jmc	1.63	IF (OB_Iw(j,bi,bj).NE.0) THEN
247			cg2d_b(OB_Iw(j,bi,bj),j,bi,bj)=0.
248			cg2d_x(OB_Iw(j,bi,bj),j,bi,bj)=0.
249	adcroft	1.12	ENDIF
250			ENDDO
251			ENDIF
252	jmc	1.49	#endif /* ALLOW_OBCS */
253			C- end bi,bj loops
254	adcroft	1.12	ENDDO
255			ENDDO
256
257	edhill	1.42	#ifdef ALLOW_DEBUG
258	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
259	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
260			& myThid)
261	adcroft	1.24	ENDIF
262	adcroft	1.23	#endif
263	jmc	1.61	IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
264			WRITE(sufx,'(I10.10)') myIter
265			CALL WRITE_FLD_XY_RS( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
266			ENDIF
267	adcroft	1.12
268	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
269			C-- gradient solver.
270	adcroft	1.22	C see CG2D.h for the interface to this routine.
271			firstResidual=0.
272			lastResidual=0.
273	adcroft	1.19	numIters=cg2dMaxIters
274	jmc	1.50	c CALL TIMER_START('CG2D [SOLVE_FOR_PRESSURE]',myThid)
275	heimbach	1.56	#ifdef ALLOW_CG2D_NSA
276			C-- Call the not-self-adjoint version of cg2d
277			CALL CG2D_NSA(
278			U cg2d_b,
279			U cg2d_x,
280			O firstResidual,
281			O lastResidual,
282			U numIters,
283			I myThid )
284			#else /* not ALLOW_CG2D_NSA = default */
285	cnh	1.1	CALL CG2D(
286	adcroft	1.22	U cg2d_b,
287	cnh	1.6	U cg2d_x,
288	adcroft	1.22	O firstResidual,
289			O lastResidual,
290	adcroft	1.19	U numIters,
291	cnh	1.1	I myThid )
292	heimbach	1.56	#endif /* ALLOW_CG2D_NSA */
293	jmc	1.66	_EXCH_XY_RL(cg2d_x, myThid )
294	jmc	1.50	c CALL TIMER_STOP ('CG2D [SOLVE_FOR_PRESSURE]',myThid)
295	adcroft	1.23
296	edhill	1.42	#ifdef ALLOW_DEBUG
297	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
298	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
299			& myThid)
300	adcroft	1.24	ENDIF
301	adcroft	1.23	#endif
302	cnh	1.1
303	jmc	1.32	C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
304	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
305	jmc	1.45	& ) THEN
306	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
307			_BEGIN_MASTER( myThid )
308			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
309			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
310			WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
311			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
312			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
313			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
314	edhill	1.43	_END_MASTER( myThid )
315	heimbach	1.38	ENDIF
316	jmc	1.32	ENDIF
317	jmc	1.17
318			C-- Transfert the 2D-solution to "etaN" :
319			DO bj=myByLo(myThid),myByHi(myThid)
320			DO bi=myBxLo(myThid),myBxHi(myThid)
321			DO j=1-OLy,sNy+OLy
322			DO i=1-OLx,sNx+OLx
323	jmc	1.18	etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
324	jmc	1.17	ENDDO
325			ENDDO
326			ENDDO
327			ENDDO
328	adcroft	1.10
329	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
330	jmc	1.53	IF ( use3Dsolver ) THEN
331	adcroft	1.9
332			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
333			C see CG3D.h for the interface to this routine.
334			DO bj=myByLo(myThid),myByHi(myThid)
335			DO bi=myBxLo(myThid),myBxHi(myThid)
336			DO j=1,sNy+1
337			DO i=1,sNx+1
338	jmc	1.18	uf(i,j)=-_recip_dxC(i,j,bi,bj)*
339	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
340	jmc	1.18	vf(i,j)=-_recip_dyC(i,j,bi,bj)*
341	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
342			ENDDO
343			ENDDO
344
345	adcroft	1.12	#ifdef ALLOW_OBCS
346	adcroft	1.14	IF (useOBCS) THEN
347	adcroft	1.9	DO i=1,sNx+1
348			C Northern boundary
349	jmc	1.63	IF (OB_Jn(i,bi,bj).NE.0) THEN
350			vf(i,OB_Jn(i,bi,bj))=0.
351	adcroft	1.9	ENDIF
352			C Southern boundary
353	jmc	1.63	IF (OB_Js(i,bi,bj).NE.0) THEN
354			vf(i,OB_Js(i,bi,bj)+1)=0.
355	adcroft	1.9	ENDIF
356			ENDDO
357			DO j=1,sNy+1
358			C Eastern boundary
359	jmc	1.63	IF (OB_Ie(j,bi,bj).NE.0) THEN
360			uf(OB_Ie(j,bi,bj),j)=0.
361	adcroft	1.9	ENDIF
362			C Western boundary
363	jmc	1.63	IF (OB_Iw(j,bi,bj).NE.0) THEN
364			uf(OB_Iw(j,bi,bj)+1,J)=0.
365	adcroft	1.9	ENDIF
366			ENDDO
367			ENDIF
368	jmc	1.49	#endif /* ALLOW_OBCS */
369	adcroft	1.9
370	jmc	1.58	IF ( usingZCoords ) THEN
371	jmc	1.51	C- Z coordinate: assume surface @ level k=1
372	jmc	1.58	tmpFac = freeSurfFac*deepFac2F(1)
373	jmc	1.51	ELSE
374			C- Other than Z coordinate: no assumption on surface level index
375	jmc	1.58	tmpFac = 0.
376	jmc	1.51	DO j=1,sNy
377			DO i=1,sNx
378			ks = ksurfC(i,j,bi,bj)
379			IF ( ks.LE.Nr ) THEN
380			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
381			& +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
382	jmc	1.58	& _rA(i,j,bi,bj)deepFac2F(ks)/deltaTmom
383	jmc	1.51	ENDIF
384			ENDDO
385			ENDDO
386			ENDIF
387	jmc	1.63	k=1
388	jmc	1.51	kp1 = MIN(k+1,Nr)
389	jmc	1.58	wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
390			IF (k.GE.Nr) wFacKp = 0.
391	adcroft	1.12	DO j=1,sNy
392			DO i=1,sNx
393	jmc	1.51	cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
394	jmc	1.63	& +drF(k)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
395			& -drF(k)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
396			& +drF(k)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
397			& -drF(k)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
398	jmc	1.51	& +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
399	jmc	1.58	& -wVel(i,j,kp1,bi,bj)*wFacKp
400	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
401			ENDDO
402			ENDDO
403	jmc	1.63	DO k=2,Nr
404	jmc	1.51	kp1 = MIN(k+1,Nr)
405	jmc	1.58	C- deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
406			C both appear in wVel term, but at 2 different levels
407			wFacKm = deepFac2F( k )*rhoFacF( k )
408			wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
409			IF (k.GE.Nr) wFacKp = 0.
410	adcroft	1.9	DO j=1,sNy
411			DO i=1,sNx
412			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
413	jmc	1.63	& +drF(k)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
414			& -drF(k)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
415			& +drF(k)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
416			& -drF(k)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
417	jmc	1.58	& +( wVel(i,j, k ,bi,bj)wFacKmmaskC(i,j,k-1,bi,bj)
418			& -wVel(i,j,kp1,bi,bj)*wFacKp
419	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
420
421	adcroft	1.9	ENDDO
422			ENDDO
423			ENDDO
424	adcroft	1.12
425			#ifdef ALLOW_OBCS
426	adcroft	1.14	IF (useOBCS) THEN
427	jmc	1.63	DO k=1,Nr
428	adcroft	1.12	DO i=1,sNx
429			C Northern boundary
430	jmc	1.63	IF (OB_Jn(i,bi,bj).NE.0) THEN
431			cg3d_b(i,OB_Jn(i,bi,bj),k,bi,bj)=0.
432	adcroft	1.12	ENDIF
433			C Southern boundary
434	jmc	1.63	IF (OB_Js(i,bi,bj).NE.0) THEN
435			cg3d_b(i,OB_Js(i,bi,bj),k,bi,bj)=0.
436	adcroft	1.12	ENDIF
437			ENDDO
438			DO j=1,sNy
439			C Eastern boundary
440	jmc	1.63	IF (OB_Ie(j,bi,bj).NE.0) THEN
441			cg3d_b(OB_Ie(j,bi,bj),j,k,bi,bj)=0.
442	adcroft	1.12	ENDIF
443			C Western boundary
444	jmc	1.63	IF (OB_Iw(j,bi,bj).NE.0) THEN
445			cg3d_b(OB_Iw(j,bi,bj),j,k,bi,bj)=0.
446	adcroft	1.12	ENDIF
447			ENDDO
448			ENDDO
449			ENDIF
450	jmc	1.49	#endif /* ALLOW_OBCS */
451			C- end bi,bj loops
452			ENDDO
453			ENDDO
454	adcroft	1.9
455	adcroft	1.25	firstResidual=0.
456			lastResidual=0.
457	jmc	1.49	numIters=cg3dMaxIters
458	jmc	1.50	CALL TIMER_START('CG3D [SOLVE_FOR_PRESSURE]',myThid)
459	adcroft	1.25	CALL CG3D(
460			U cg3d_b,
461			U phi_nh,
462			O firstResidual,
463			O lastResidual,
464			U numIters,
465			I myThid )
466	jmc	1.66	_EXCH_XYZ_RL(phi_nh, myThid )
467	jmc	1.50	CALL TIMER_STOP ('CG3D [SOLVE_FOR_PRESSURE]',myThid)
468	adcroft	1.25
469	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
470	jmc	1.45	& ) THEN
471	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
472			_BEGIN_MASTER( myThid )
473			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
474			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
475			WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
476			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
477			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
478			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
479	edhill	1.43	_END_MASTER( myThid )
480	heimbach	1.38	ENDIF
481	mlosch	1.37	ENDIF
482	adcroft	1.9
483	jmc	1.49	C-- Update surface pressure (account for NH-p @ surface level) and NH pressure:
484			IF ( zeroPsNH ) THEN
485			DO bj=myByLo(myThid),myByHi(myThid)
486			DO bi=myBxLo(myThid),myBxHi(myThid)
487
488			IF ( usingZCoords ) THEN
489			C- Z coordinate: assume surface @ level k=1
490			DO k=2,Nr
491			DO j=1-OLy,sNy+OLy
492			DO i=1-OLx,sNx+OLx
493			phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
494			& - phi_nh(i,j,1,bi,bj)
495			ENDDO
496			ENDDO
497			ENDDO
498			DO j=1-OLy,sNy+OLy
499			DO i=1-OLx,sNx+OLx
500			etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
501			& *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
502			phi_nh(i,j,1,bi,bj) = 0.
503			ENDDO
504			ENDDO
505			ELSE
506			C- Other than Z coordinate: no assumption on surface level index
507			DO j=1-OLy,sNy+OLy
508			DO i=1-OLx,sNx+OLx
509			ks = ksurfC(i,j,bi,bj)
510			IF ( ks.LE.Nr ) THEN
511			etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
512			& *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
513			DO k=Nr,1,-1
514			phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
515			& - phi_nh(i,j,ks,bi,bj)
516			ENDDO
517			ENDIF
518			ENDDO
519			ENDDO
520			ENDIF
521
522			ENDDO
523			ENDDO
524	adcroft	1.9	ENDIF
525	jmc	1.49
526			ENDIF
527			#endif /* ALLOW_NONHYDROSTATIC */
528	cnh	1.1
529	heimbach	1.60	#ifdef ALLOW_SHOWFLOPS
530			CALL SHOWFLOPS_INSOLVE( myThid)
531	ce107	1.52	#endif
532	heimbach	1.60
533	cnh	1.1	RETURN
534			END