model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.58 2006/12/05 05:25:08 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
      _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL firstResidual,lastResidual
      _RL tmpFac
      _RL sumEmP, tileEmP
      LOGICAL putPmEinXvector
      INTEGER numIters, ks
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_NONHYDROSTATIC
      INTEGER kp1
      _RL     wFacKm, wFacKp
      LOGICAL zeroPsNH
#endif
CEOP

#ifdef TIME_PER_TIMESTEP_SFP
CCE107 common block for per timestep timing
C     !TIMING VARIABLES
C     == Timing variables ==
      REAL*8 utnew, utold, stnew, stold, wtnew, wtold
      COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
#endif
#ifdef USE_PAPI_FLOPS_SFP
CCE107 common block for PAPI summary performance
#include <fpapi.h>
      INTEGER*8 flpops, instr
      INTEGER check
      REAL*4 real_time, proc_time, mflops, ipc
      COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
#else
#ifdef USE_PCL_FLOPS_SFP
CCE107 common block for PCL summary performance
#include <pclh.f>
      INTEGER pcl_counter_list(5), flags, nevents, res, ipcl
      INTEGER*8 i_result(5), descr
      REAL*8 fp_result(5)
      COMMON /pclvars/ i_result, descr, fp_result, pcl_counter_list,
     $     flags, nevents
      INTEGER nmaxevents
      PARAMETER (nmaxevents = 61)
      CHARACTER*22 pcl_counter_name(0:nmaxevents-1)
      COMMON /pclnames/ pcl_counter_name
#endif
#endif

#ifdef ALLOW_NONHYDROSTATIC
c       zeroPsNH = .FALSE.
        zeroPsNH = exactConserv
#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

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) THEN
         tmpFac = freeSurfFac*convertEmP2rUnit
         IF (exactConserv)
     &        tmpFac = freeSurfFac*convertEmP2rUnit*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 = 0.
         DO j=1,sNy
          DO i=1,sNx
            tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
     &                                       *maskH(i,j,bi,bj)
          ENDDO
         ENDDO
         sumEmP = sumEmP + tileEmP
        ENDIF
       ENDDO
      ENDDO
      IF ( putPmEinXvector ) THEN
        _GLOBAL_SUM_R8( 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
     &                          *convertEmP2rUnit*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
         DO j=1,sNy+1
          DO i=1,sNx+1
           uf(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
     &               *drF(k)*_hFacW(i,j,k,bi,bj)*rhoFacC(k)
           vf(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
     &               *drF(k)*_hFacS(i,j,k,bi,bj)*rhoFacC(k)
          ENDDO
         ENDDO
         CALL CALC_DIV_GHAT(
     I       bi,bj,1,sNx,1,sNy,K,
     I       uf,vf,
     U       cg2d_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

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_R8(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_R8(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 TIME_PER_TIMESTEP_SFP
CCE107 Time per timestep information
      _BEGIN_MASTER( myThid )
      CALL TIMER_GET_TIME( utnew, stnew, wtnew )
C Only output timing information after the 1st timestep
      IF ( wtold .NE. 0.0D0 ) THEN
        WRITE(msgBuf,'(A34,3F10.6)')
     $        'User, system and wallclock time:', utnew - utold,
     $        stnew - stold, wtnew - wtold
         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      ENDIF
      utold = utnew
      stold = stnew
      wtold = wtnew
      _END_MASTER( myThid )
#endif
#ifdef USE_PAPI_FLOPS_SFP
CCE107 PAPI summary performance
      _BEGIN_MASTER( myThid )
#ifdef USE_FLIPS
      call PAPIF_flips(real_time, proc_time, flpops, mflops, check)
#else
      call PAPIF_flops(real_time, proc_time, flpops, mflops, check)
#endif
      WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
     $     'Mflop/s during this timestep:', mflops, ' ', mflops
     $     *proc_time/(real_time + 1E-36)
      CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
#ifdef PAPI_VERSION
      call PAPIF_ipc(real_time, proc_time, instr, ipc, check)
      WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
     $     'IPC during this timestep:', ipc, ' ', ipc*proc_time
     $     /(real_time + 1E-36)
      CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
#endif
      _END_MASTER( myThid )
#else
#ifdef USE_PCL_FLOPS_SFP
CCE107 PCL summary performance
      _BEGIN_MASTER( myThid )
      PCLstop(descr, i_result, fp_result, nevents)
      do ipcl = 1, nevents
         WRITE(msgBuf,'(A22,A26,F10.6)'),
     $        pcl_counter_name(pcl_counter_list(ipcl)),
     $        'during this timestep:', fp_results(ipcl)
         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
      enddo
      PCLstart(descr, pcl_counter_list, nevents, flags)
      _END_MASTER( myThid )
#endif
#endif
      RETURN
      END

#ifdef TIME_PER_TIMESTEP_SFP
CCE107 Initialization of common block for per timestep timing
      BLOCK DATA settimers
C     !TIMING VARIABLES
C     == Timing variables ==
      REAL*8 utnew, utold, stnew, stold, wtnew, wtold
      COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
      DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
      END
#endif
#ifdef USE_PAPI_FLOPS_SFP
CCE107 Initialization of common block for PAPI summary performance
      BLOCK DATA setpapis
      INTEGER*8 flpops, instr
      REAL real_time, proc_time, mflops, ipc
      COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
      DATA flpops, instr, real_time, proc_time, mflops, ipc /2*0,4*0.E0/
      END
#endif
1	jmc	1.59	C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.58 2006/12/05 05:25:08 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.9	_RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59			_RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60	adcroft	1.22	_RL firstResidual,lastResidual
61	jmc	1.36	_RL tmpFac
62	jmc	1.47	_RL sumEmP, tileEmP
63			LOGICAL putPmEinXvector
64	jmc	1.58	INTEGER numIters, ks
65	adcroft	1.25	CHARACTER*(MAX_LEN_MBUF) msgBuf
66	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
67	jmc	1.58	INTEGER kp1
68			_RL wFacKm, wFacKp
69	jmc	1.49	LOGICAL zeroPsNH
70			#endif
71	cnh	1.27	CEOP
72	jmc	1.17
73	ce107	1.52	#ifdef TIME_PER_TIMESTEP_SFP
74	ce107	1.44	CCE107 common block for per timestep timing
75			C !TIMING VARIABLES
76			C == Timing variables ==
77			REAL*8 utnew, utold, stnew, stold, wtnew, wtold
78			COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
79			#endif
80	ce107	1.52	#ifdef USE_PAPI_FLOPS_SFP
81			CCE107 common block for PAPI summary performance
82			#include <fpapi.h>
83	ce107	1.54	INTEGER*8 flpops, instr
84	ce107	1.52	INTEGER check
85	ce107	1.54	REAL*4 real_time, proc_time, mflops, ipc
86			COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
87			#else
88			#ifdef USE_PCL_FLOPS_SFP
89			CCE107 common block for PCL summary performance
90			#include <pclh.f>
91			INTEGER pcl_counter_list(5), flags, nevents, res, ipcl
92			INTEGER*8 i_result(5), descr
93			REAL*8 fp_result(5)
94			COMMON /pclvars/ i_result, descr, fp_result, pcl_counter_list,
95			$ flags, nevents
96			INTEGER nmaxevents
97			PARAMETER (nmaxevents = 61)
98			CHARACTER*22 pcl_counter_name(0:nmaxevents-1)
99			COMMON /pclnames/ pcl_counter_name
100			#endif
101	ce107	1.52	#endif
102	ce107	1.44
103	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
104			c zeroPsNH = .FALSE.
105			zeroPsNH = exactConserv
106			#endif
107
108	jmc	1.58	C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
109			C anelastic (always Z-coordinate):
110			C 1) assume that rhoFacF(1)=1 (and ksurf == 1);
111			C (this reduces the number of lines of code to modify)
112			C 2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
113			C (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
114			C => 2 factors cancel in elliptic eq. for Phi_s ,
115			C but 1rst factor(a) remains in RHS cg2d_b.
116
117	jmc	1.47	C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
118			C instead of simply etaN ; This can speed-up the solver convergence in
119			C the case where \|Global_mean_PmE\| is large.
120			putPmEinXvector = .FALSE.
121			c putPmEinXvector = useRealFreshWaterFlux
122
123	jmc	1.17	C-- Save previous solution & Initialise Vector solution and source term :
124	jmc	1.47	sumEmP = 0.
125	jmc	1.17	DO bj=myByLo(myThid),myByHi(myThid)
126			DO bi=myBxLo(myThid),myBxHi(myThid)
127			DO j=1-OLy,sNy+OLy
128			DO i=1-OLx,sNx+OLx
129	edhill	1.40	#ifdef ALLOW_CD_CODE
130	jmc	1.17	etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
131	jmc	1.26	#endif
132	jmc	1.18	cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
133	jmc	1.17	cg2d_b(i,j,bi,bj) = 0.
134			ENDDO
135			ENDDO
136	jmc	1.29	IF (useRealFreshWaterFlux) THEN
137	jmc	1.36	tmpFac = freeSurfFac*convertEmP2rUnit
138	jmc	1.58	IF (exactConserv)
139	jmc	1.36	& tmpFac = freeSurfFacconvertEmP2rUnitimplicDiv2DFlow
140	jmc	1.29	DO j=1,sNy
141			DO i=1,sNx
142	jmc	1.58	cg2d_b(i,j,bi,bj) =
143	jmc	1.29	& tmpFac_rA(i,j,bi,bj)EmPmR(i,j,bi,bj)/deltaTMom
144			ENDDO
145			ENDDO
146			ENDIF
147	jmc	1.47	IF ( putPmEinXvector ) THEN
148			tileEmP = 0.
149			DO j=1,sNy
150			DO i=1,sNx
151			tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
152			& *maskH(i,j,bi,bj)
153			ENDDO
154			ENDDO
155			sumEmP = sumEmP + tileEmP
156			ENDIF
157	jmc	1.17	ENDDO
158			ENDDO
159	jmc	1.47	IF ( putPmEinXvector ) THEN
160			_GLOBAL_SUM_R8( sumEmP, myThid )
161			ENDIF
162	adcroft	1.12
163			DO bj=myByLo(myThid),myByHi(myThid)
164			DO bi=myBxLo(myThid),myBxHi(myThid)
165	jmc	1.47	IF ( putPmEinXvector ) THEN
166			tmpFac = 0.
167			IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf
168			& convertEmP2rUnitsumEmP/globalArea
169			DO j=1,sNy
170			DO i=1,sNx
171			cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
172			& - tmpFac*Bo_surf(i,j,bi,bj)
173			ENDDO
174			ENDDO
175			ENDIF
176	jmc	1.58	C- RHS: similar to the divergence of the vertically integrated mass transport:
177			C del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] } / deltaT
178	adcroft	1.12	DO K=Nr,1,-1
179			DO j=1,sNy+1
180			DO i=1,sNx+1
181	jmc	1.58	uf(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
182			& drF(k)_hFacW(i,j,k,bi,bj)*rhoFacC(k)
183			vf(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
184			& drF(k)_hFacS(i,j,k,bi,bj)*rhoFacC(k)
185	adcroft	1.12	ENDDO
186			ENDDO
187			CALL CALC_DIV_GHAT(
188			I bi,bj,1,sNx,1,sNy,K,
189			I uf,vf,
190	jmc	1.17	U cg2d_b,
191	adcroft	1.12	I myThid)
192			ENDDO
193			ENDDO
194			ENDDO
195	cnh	1.4
196	adcroft	1.12	C-- Add source term arising from w=d/dt (p_s + p_nh)
197			DO bj=myByLo(myThid),myByHi(myThid)
198			DO bi=myBxLo(myThid),myBxHi(myThid)
199	adcroft	1.13	#ifdef ALLOW_NONHYDROSTATIC
200	jmc	1.53	IF ( use3Dsolver .AND. zeroPsNH ) THEN
201	jmc	1.49	DO j=1,sNy
202			DO i=1,sNx
203	jmc	1.51	ks = ksurfC(i,j,bi,bj)
204			IF ( ks.LE.Nr ) THEN
205			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
206	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
207			& /deltaTMom/deltaTfreesurf
208	jmc	1.49	& * etaH(i,j,bi,bj)
209	jmc	1.51	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
210	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
211			& /deltaTMom/deltaTfreesurf
212	jmc	1.49	& * etaH(i,j,bi,bj)
213	jmc	1.51	ENDIF
214	jmc	1.49	ENDDO
215			ENDDO
216	jmc	1.53	ELSEIF ( use3Dsolver ) THEN
217	jmc	1.28	DO j=1,sNy
218			DO i=1,sNx
219	jmc	1.51	ks = ksurfC(i,j,bi,bj)
220			IF ( ks.LE.Nr ) THEN
221			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
222	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
223			& /deltaTMom/deltaTfreesurf
224	jmc	1.28	& *( etaN(i,j,bi,bj)
225	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
226	jmc	1.51	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
227	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
228			& /deltaTMom/deltaTfreesurf
229	jmc	1.28	& *( etaN(i,j,bi,bj)
230	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
231	jmc	1.51	ENDIF
232	jmc	1.28	ENDDO
233	adcroft	1.12	ENDDO
234	jmc	1.28	ELSEIF ( exactConserv ) THEN
235	adcroft	1.13	#else
236	jmc	1.26	IF ( exactConserv ) THEN
237	edhill	1.39	#endif /* ALLOW_NONHYDROSTATIC */
238	jmc	1.26	DO j=1,sNy
239			DO i=1,sNx
240	jmc	1.58	ks = ksurfC(i,j,bi,bj)
241	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
242	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
243			& /deltaTMom/deltaTfreesurf
244	jmc	1.26	& * etaH(i,j,bi,bj)
245			ENDDO
246			ENDDO
247			ELSE
248			DO j=1,sNy
249			DO i=1,sNx
250	jmc	1.58	ks = ksurfC(i,j,bi,bj)
251	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
252	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
253			& /deltaTMom/deltaTfreesurf
254	jmc	1.26	& * etaN(i,j,bi,bj)
255			ENDDO
256	adcroft	1.12	ENDDO
257	jmc	1.26	ENDIF
258	adcroft	1.12
259			#ifdef ALLOW_OBCS
260	adcroft	1.14	IF (useOBCS) THEN
261	adcroft	1.12	DO i=1,sNx
262			C Northern boundary
263			IF (OB_Jn(I,bi,bj).NE.0) THEN
264			cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
265	jmc	1.31	cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
266	adcroft	1.12	ENDIF
267			C Southern boundary
268			IF (OB_Js(I,bi,bj).NE.0) THEN
269			cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
270	jmc	1.31	cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
271	adcroft	1.12	ENDIF
272			ENDDO
273			DO j=1,sNy
274			C Eastern boundary
275			IF (OB_Ie(J,bi,bj).NE.0) THEN
276			cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
277	jmc	1.31	cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
278	adcroft	1.12	ENDIF
279			C Western boundary
280			IF (OB_Iw(J,bi,bj).NE.0) THEN
281			cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
282	jmc	1.31	cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
283	adcroft	1.12	ENDIF
284			ENDDO
285			ENDIF
286	jmc	1.49	#endif /* ALLOW_OBCS */
287			C- end bi,bj loops
288	adcroft	1.12	ENDDO
289			ENDDO
290
291	edhill	1.42	#ifdef ALLOW_DEBUG
292	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
293	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
294			& myThid)
295	adcroft	1.24	ENDIF
296	adcroft	1.23	#endif
297	adcroft	1.12
298	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
299			C-- gradient solver.
300	adcroft	1.22	C see CG2D.h for the interface to this routine.
301			firstResidual=0.
302			lastResidual=0.
303	adcroft	1.19	numIters=cg2dMaxIters
304	jmc	1.50	c CALL TIMER_START('CG2D [SOLVE_FOR_PRESSURE]',myThid)
305	heimbach	1.56	#ifdef ALLOW_CG2D_NSA
306			C-- Call the not-self-adjoint version of cg2d
307			CALL CG2D_NSA(
308			U cg2d_b,
309			U cg2d_x,
310			O firstResidual,
311			O lastResidual,
312			U numIters,
313			I myThid )
314			#else /* not ALLOW_CG2D_NSA = default */
315	cnh	1.1	CALL CG2D(
316	adcroft	1.22	U cg2d_b,
317	cnh	1.6	U cg2d_x,
318	adcroft	1.22	O firstResidual,
319			O lastResidual,
320	adcroft	1.19	U numIters,
321	cnh	1.1	I myThid )
322	heimbach	1.56	#endif /* ALLOW_CG2D_NSA */
323	adcroft	1.19	_EXCH_XY_R8(cg2d_x, myThid )
324	jmc	1.50	c CALL TIMER_STOP ('CG2D [SOLVE_FOR_PRESSURE]',myThid)
325	adcroft	1.23
326	edhill	1.42	#ifdef ALLOW_DEBUG
327	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
328	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
329			& myThid)
330	adcroft	1.24	ENDIF
331	adcroft	1.23	#endif
332	cnh	1.1
333	jmc	1.32	C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
334	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
335	jmc	1.45	& ) THEN
336	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
337			_BEGIN_MASTER( myThid )
338			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
339			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
340			WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
341			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
342			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
343			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
344	edhill	1.43	_END_MASTER( myThid )
345	heimbach	1.38	ENDIF
346	jmc	1.32	ENDIF
347	jmc	1.17
348			C-- Transfert the 2D-solution to "etaN" :
349			DO bj=myByLo(myThid),myByHi(myThid)
350			DO bi=myBxLo(myThid),myBxHi(myThid)
351			DO j=1-OLy,sNy+OLy
352			DO i=1-OLx,sNx+OLx
353	jmc	1.18	etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
354	jmc	1.17	ENDDO
355			ENDDO
356			ENDDO
357			ENDDO
358	adcroft	1.10
359	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
360	jmc	1.53	IF ( use3Dsolver ) THEN
361	adcroft	1.9
362			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
363			C see CG3D.h for the interface to this routine.
364			DO bj=myByLo(myThid),myByHi(myThid)
365			DO bi=myBxLo(myThid),myBxHi(myThid)
366			DO j=1,sNy+1
367			DO i=1,sNx+1
368	jmc	1.18	uf(i,j)=-_recip_dxC(i,j,bi,bj)*
369	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
370	jmc	1.18	vf(i,j)=-_recip_dyC(i,j,bi,bj)*
371	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
372			ENDDO
373			ENDDO
374
375	adcroft	1.12	#ifdef ALLOW_OBCS
376	adcroft	1.14	IF (useOBCS) THEN
377	adcroft	1.9	DO i=1,sNx+1
378			C Northern boundary
379			IF (OB_Jn(I,bi,bj).NE.0) THEN
380			vf(I,OB_Jn(I,bi,bj))=0.
381			ENDIF
382			C Southern boundary
383			IF (OB_Js(I,bi,bj).NE.0) THEN
384			vf(I,OB_Js(I,bi,bj)+1)=0.
385			ENDIF
386			ENDDO
387			DO j=1,sNy+1
388			C Eastern boundary
389			IF (OB_Ie(J,bi,bj).NE.0) THEN
390			uf(OB_Ie(J,bi,bj),J)=0.
391			ENDIF
392			C Western boundary
393			IF (OB_Iw(J,bi,bj).NE.0) THEN
394			uf(OB_Iw(J,bi,bj)+1,J)=0.
395			ENDIF
396			ENDDO
397			ENDIF
398	jmc	1.49	#endif /* ALLOW_OBCS */
399	adcroft	1.9
400	jmc	1.58	IF ( usingZCoords ) THEN
401	jmc	1.51	C- Z coordinate: assume surface @ level k=1
402	jmc	1.58	tmpFac = freeSurfFac*deepFac2F(1)
403	jmc	1.51	ELSE
404			C- Other than Z coordinate: no assumption on surface level index
405	jmc	1.58	tmpFac = 0.
406	jmc	1.51	DO j=1,sNy
407			DO i=1,sNx
408			ks = ksurfC(i,j,bi,bj)
409			IF ( ks.LE.Nr ) THEN
410			cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
411			& +freeSurfFac*etaN(i,j,bi,bj)/deltaTfreesurf
412	jmc	1.58	& _rA(i,j,bi,bj)deepFac2F(ks)/deltaTmom
413	jmc	1.51	ENDIF
414			ENDDO
415			ENDDO
416			ENDIF
417	adcroft	1.12	K=1
418	jmc	1.51	kp1 = MIN(k+1,Nr)
419	jmc	1.58	wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
420			IF (k.GE.Nr) wFacKp = 0.
421	adcroft	1.12	DO j=1,sNy
422			DO i=1,sNx
423	jmc	1.51	cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
424	heimbach	1.56	& +drF(K)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
425			& -drF(K)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
426			& +drF(K)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
427			& -drF(K)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
428	jmc	1.51	& +( tmpFac*etaN(i,j,bi,bj)/deltaTfreesurf
429	jmc	1.58	& -wVel(i,j,kp1,bi,bj)*wFacKp
430	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
431			ENDDO
432			ENDDO
433	jmc	1.51	DO K=2,Nr
434			kp1 = MIN(k+1,Nr)
435	jmc	1.58	C- deepFac & rhoFac cancel with the ones in uf[=del_i(Phi)/dx],vf ;
436			C both appear in wVel term, but at 2 different levels
437			wFacKm = deepFac2F( k )*rhoFacF( k )
438			wFacKp = deepFac2F(kp1)*rhoFacF(kp1)
439			IF (k.GE.Nr) wFacKp = 0.
440	adcroft	1.9	DO j=1,sNy
441			DO i=1,sNx
442			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
443	heimbach	1.56	& +drF(K)dyG(i+1,j,bi,bj)_hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
444			& -drF(K)dyG( i ,j,bi,bj)_hFacW( i ,j,k,bi,bj)*uf( i ,j)
445			& +drF(K)dxG(i,j+1,bi,bj)_hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
446			& -drF(K)dxG(i, j ,bi,bj)_hFacS(i, j ,k,bi,bj)*vf(i, j )
447	jmc	1.58	& +( wVel(i,j, k ,bi,bj)wFacKmmaskC(i,j,k-1,bi,bj)
448			& -wVel(i,j,kp1,bi,bj)*wFacKp
449	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
450
451	adcroft	1.9	ENDDO
452			ENDDO
453			ENDDO
454	adcroft	1.12
455			#ifdef ALLOW_OBCS
456	adcroft	1.14	IF (useOBCS) THEN
457	adcroft	1.12	DO K=1,Nr
458			DO i=1,sNx
459			C Northern boundary
460			IF (OB_Jn(I,bi,bj).NE.0) THEN
461			cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
462			ENDIF
463			C Southern boundary
464			IF (OB_Js(I,bi,bj).NE.0) THEN
465			cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
466			ENDIF
467			ENDDO
468			DO j=1,sNy
469			C Eastern boundary
470			IF (OB_Ie(J,bi,bj).NE.0) THEN
471			cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
472			ENDIF
473			C Western boundary
474			IF (OB_Iw(J,bi,bj).NE.0) THEN
475			cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
476			ENDIF
477			ENDDO
478			ENDDO
479			ENDIF
480	jmc	1.49	#endif /* ALLOW_OBCS */
481			C- end bi,bj loops
482			ENDDO
483			ENDDO
484	adcroft	1.9
485	adcroft	1.25	firstResidual=0.
486			lastResidual=0.
487	jmc	1.49	numIters=cg3dMaxIters
488	jmc	1.50	CALL TIMER_START('CG3D [SOLVE_FOR_PRESSURE]',myThid)
489	adcroft	1.25	CALL CG3D(
490			U cg3d_b,
491			U phi_nh,
492			O firstResidual,
493			O lastResidual,
494			U numIters,
495			I myThid )
496			_EXCH_XYZ_R8(phi_nh, myThid )
497	jmc	1.50	CALL TIMER_STOP ('CG3D [SOLVE_FOR_PRESSURE]',myThid)
498	adcroft	1.25
499	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
500	jmc	1.45	& ) THEN
501	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
502			_BEGIN_MASTER( myThid )
503			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
504			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
505			WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
506			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
507			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
508			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
509	edhill	1.43	_END_MASTER( myThid )
510	heimbach	1.38	ENDIF
511	mlosch	1.37	ENDIF
512	adcroft	1.9
513	jmc	1.49	C-- Update surface pressure (account for NH-p @ surface level) and NH pressure:
514			IF ( zeroPsNH ) THEN
515			DO bj=myByLo(myThid),myByHi(myThid)
516			DO bi=myBxLo(myThid),myBxHi(myThid)
517
518			IF ( usingZCoords ) THEN
519			C- Z coordinate: assume surface @ level k=1
520			DO k=2,Nr
521			DO j=1-OLy,sNy+OLy
522			DO i=1-OLx,sNx+OLx
523			phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
524			& - phi_nh(i,j,1,bi,bj)
525			ENDDO
526			ENDDO
527			ENDDO
528			DO j=1-OLy,sNy+OLy
529			DO i=1-OLx,sNx+OLx
530			etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
531			& *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,1,bi,bj))
532			phi_nh(i,j,1,bi,bj) = 0.
533			ENDDO
534			ENDDO
535			ELSE
536			C- Other than Z coordinate: no assumption on surface level index
537			DO j=1-OLy,sNy+OLy
538			DO i=1-OLx,sNx+OLx
539			ks = ksurfC(i,j,bi,bj)
540			IF ( ks.LE.Nr ) THEN
541			etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)
542			& *(cg2d_x(i,j,bi,bj) + phi_nh(i,j,ks,bi,bj))
543			DO k=Nr,1,-1
544			phi_nh(i,j,k,bi,bj) = phi_nh(i,j,k,bi,bj)
545			& - phi_nh(i,j,ks,bi,bj)
546			ENDDO
547			ENDIF
548			ENDDO
549			ENDDO
550			ENDIF
551
552			ENDDO
553			ENDDO
554	adcroft	1.9	ENDIF
555	jmc	1.49
556			ENDIF
557			#endif /* ALLOW_NONHYDROSTATIC */
558	cnh	1.1
559	ce107	1.52	#ifdef TIME_PER_TIMESTEP_SFP
560	ce107	1.44	CCE107 Time per timestep information
561			_BEGIN_MASTER( myThid )
562			CALL TIMER_GET_TIME( utnew, stnew, wtnew )
563			C Only output timing information after the 1st timestep
564			IF ( wtold .NE. 0.0D0 ) THEN
565			WRITE(msgBuf,'(A34,3F10.6)')
566			$ 'User, system and wallclock time:', utnew - utold,
567			$ stnew - stold, wtnew - wtold
568			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
569			ENDIF
570			utold = utnew
571			stold = stnew
572			wtold = wtnew
573			_END_MASTER( myThid )
574			#endif
575	ce107	1.52	#ifdef USE_PAPI_FLOPS_SFP
576			CCE107 PAPI summary performance
577			_BEGIN_MASTER( myThid )
578	ce107	1.54	#ifdef USE_FLIPS
579			call PAPIF_flips(real_time, proc_time, flpops, mflops, check)
580			#else
581	ce107	1.52	call PAPIF_flops(real_time, proc_time, flpops, mflops, check)
582	ce107	1.54	#endif
583	ce107	1.55	WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
584			$ 'Mflop/s during this timestep:', mflops, ' ', mflops
585			$ *proc_time/(real_time + 1E-36)
586	ce107	1.52	CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
587	ce107	1.54	#ifdef PAPI_VERSION
588			call PAPIF_ipc(real_time, proc_time, instr, ipc, check)
589	ce107	1.55	WRITE(msgBuf,'(A34,F10.6,A,F10.6)')
590			$ 'IPC during this timestep:', ipc, ' ', ipc*proc_time
591	ce107	1.57	$ /(real_time + 1E-36)
592	ce107	1.54	CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
593			#endif
594	ce107	1.52	_END_MASTER( myThid )
595	ce107	1.54	#else
596			#ifdef USE_PCL_FLOPS_SFP
597			CCE107 PCL summary performance
598			_BEGIN_MASTER( myThid )
599			PCLstop(descr, i_result, fp_result, nevents)
600			do ipcl = 1, nevents
601			WRITE(msgBuf,'(A22,A26,F10.6)'),
602			$ pcl_counter_name(pcl_counter_list(ipcl)),
603			$ 'during this timestep:', fp_results(ipcl)
604			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
605			enddo
606			PCLstart(descr, pcl_counter_list, nevents, flags)
607			_END_MASTER( myThid )
608			#endif
609	ce107	1.52	#endif
610	cnh	1.1	RETURN
611			END
612	ce107	1.44
613	ce107	1.52	#ifdef TIME_PER_TIMESTEP_SFP
614	ce107	1.44	CCE107 Initialization of common block for per timestep timing
615			BLOCK DATA settimers
616			C !TIMING VARIABLES
617			C == Timing variables ==
618			REAL*8 utnew, utold, stnew, stold, wtnew, wtold
619			COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
620			DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
621			END
622			#endif
623	ce107	1.52	#ifdef USE_PAPI_FLOPS_SFP
624			CCE107 Initialization of common block for PAPI summary performance
625			BLOCK DATA setpapis
626	ce107	1.54	INTEGER*8 flpops, instr
627			REAL real_time, proc_time, mflops, ipc
628			COMMON /papivars/ flpops, instr, real_time, proc_time, mflops, ipc
629			DATA flpops, instr, real_time, proc_time, mflops, ipc /20,40.E0/
630	ce107	1.52	END
631			#endif