C $Header: /u/gcmpack/MITgcm/model/src/post_cg3d.F,v 1.8 2017/03/24 15:41:19 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: POST_CG3D
C     !INTERFACE:
      SUBROUTINE POST_CG3D(
     I                      zeroPsNH, zeroMeanPnh,
     I                      myTime, myIter, myThid )

C     !DESCRIPTION:
C     Called from SOLVE_FOR_PRESSURE, after 3-D solver (cg3d):
C     Finish computation of Non-hydrostatic pressure from 3-D solver solution

C     !USES:
      IMPLICIT NONE
C     == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
c#include "FFIELDS.h"
#include "DYNVARS.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "NH_VARS.h"
#endif

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

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     zeroPsNH    :: account for Hyd.component of cg3d_x by updating NH & Surf.Press
C     zeroMeanPnh :: account for Hyd.component of cg3d_x by updating NH & Surf.Press
C     myTime      :: Current time in simulation
C     myIter      :: Current iteration number in simulation
C     myThid      :: My Thread Id. number
      LOGICAL zeroPsNH, zeroMeanPnh
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

#ifdef ALLOW_NONHYDROSTATIC
C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER i,j,k,bi,bj
      INTEGER ks
c     CHARACTER*(MAX_LEN_MBUF) msgBuf
      _RL     locGamma
CEOP

C--   Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
C     from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
      IF ( nonHydrostatic .AND. exactConserv ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

         IF ( selectNHfreeSurf.GE.1 ) THEN
          DO j=1,sNy
           DO i=1,sNx
            locGamma = drC(1)*recip_Bo(i,j,bi,bj)
     &               /( deltaTMom*deltaTFreeSurf
     &                 *implicitNHPress*implicDiv2DFlow )
            ks = 1
c           ks = kSurfC(i,j,bi,bj)
c           IF ( ks.LE.Nr ) THEN
             dPhiNH(i,j,bi,bj) = ( phi_nh(i,j,ks,bi,bj)
     &           + locGamma*Bo_surf(i,j,bi,bj)
     &                     *implicDiv2DFlow*deltaTFreeSurf
c    &                     *( wVel(i,j,ks,bi,bj) - wSurfP2d(i,j) )
     &                     *( wVel(i,j,ks,bi,bj) - dPhiNH(i,j,bi,bj) )
     &                           )/(1. _d 0 + locGamma )
c           ENDIF
           ENDDO
          ENDDO
         ELSEIF ( uniformFreeSurfLev ) THEN
C-       Z coordinate: assume surface @ level k=1
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             dPhiNH(i,j,bi,bj) = phi_nh(i,j,1,bi,bj)
           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
             dPhiNH(i,j,bi,bj) = phi_nh(i,j,ks,bi,bj)
            ELSE
             dPhiNH(i,j,bi,bj) = 0.
            ENDIF
           ENDDO
          ENDDO
         ENDIF

        ENDDO
       ENDDO
       IF ( selectNHfreeSurf.GE.1 .AND.
     &  ( implicitNHPress.LT.oneRL .OR. selectP_inEOS_Zc.EQ.3 ) ) THEN
         CALL EXCH_XY_RL( dPhiNH, 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
        CALL WRITE_FLD_XYZ_RL( 'cg3d_x','I10', phi_nh, myIter, myThid )
       ENDIF
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

         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)
     &                            - dPhiNH(i,j,bi,bj)
     &                            )*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) = etaN(i,j,bi,bj)
     &                      + recip_Bo(i,j,bi,bj)*dPhiNH(i,j,bi,bj)
            dPhiNH(i,j,bi,bj) = 0. _d 0
          ENDDO
         ENDDO

        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END