model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.76 2011/05/18 23:41:26 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

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
      INTEGER ks
      INTEGER numIters
      _RL firstResidual,lastResidual
      _RL tmpFac
      _RL sumEmP, tileEmP(nSx,nSy)
      LOGICAL putPmEinXvector
      INTEGER ioUnit
      CHARACTER*10 sufx
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_NONHYDROSTATIC
      LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
#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)
     &                                    *maskInC(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

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_NONHYDROSTATIC
        IF ( oldFreeSurfTerm ) THEN
C--   Add source term arising from w=d/dt (p_s + p_nh)
         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
C--   Add source term arising from w=d/dt (p_s)
        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
C- Note: solver matrix is trivial outside OB region (main diagonal only)
C     => no real need to reset RHS (=cg2d_b) & cg2d_x, except that:
C    a) normalisation is fct of Max(RHS), which can be large ouside OB region
C      (would be different if we were solving for increment of eta/g
C       instead of directly for eta/g).
C       => need to reset RHS to ensure that interior solution does not depend
C       on ouside OB region.
C    b) provide directly the trivial solution cg2d_x == 0 for outside OB region
C      (=> no residual => no effect on solver convergence and interior solution)
        IF (useOBCS) THEN
         DO j=1,sNy
          DO i=1,sNx
           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*maskInC(i,j,bi,bj)
           cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*maskInC(i,j,bi,bj)
         ENDDO
         ENDDO
        ENDIF
#endif /* ALLOW_OBCS */
C-    end bi,bj loops
       ENDDO
      ENDDO

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevD ) 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. debLevD ) 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.

C--   Finish updating cg3d_b: 1) Add EmPmR contribution to top level cg3d_b:
C                             2) Update or Add free-surface contribution
C                             3) increment in horiz velocity due to new cg2d_x
C                             4) add vertical velocity contribution.
       CALL PRE_CG3D(
     I                oldFreeSurfTerm,
     I                cg2d_x,
     U                cg3d_b,
     I                myTime, myIter, myThid )

#ifdef ALLOW_DEBUG
       IF ( debugLevel .GE. debLevD ) 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--   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
        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
        CALL POST_CG3D(
     I                  zeroPsNH, zeroMeanPnh,
     I                  myTime, myIter, myThid )
       ENDIF

      ENDIF
#endif /* ALLOW_NONHYDROSTATIC */

#ifdef ALLOW_SHOWFLOPS
      CALL SHOWFLOPS_INSOLVE( myThid)
#endif

      RETURN
      END
1	jmc	1.77	C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.76 2011/05/18 23:41:26 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	cnh	1.4
39	jmc	1.32	C === Functions ====
40	jmc	1.46	LOGICAL DIFFERENT_MULTIPLE
41			EXTERNAL DIFFERENT_MULTIPLE
42	jmc	1.32
43	cnh	1.27	C !INPUT/OUTPUT PARAMETERS:
44	cnh	1.1	C == Routine arguments ==
45	jmc	1.58	C myTime :: Current time in simulation
46			C myIter :: Current iteration number in simulation
47			C myThid :: Thread number for this instance of SOLVE_FOR_PRESSURE
48	jmc	1.28	_RL myTime
49			INTEGER myIter
50	jmc	1.29	INTEGER myThid
51	cnh	1.4
52	cnh	1.27	C !LOCAL VARIABLES:
53	adcroft	1.22	C == Local variables ==
54	cnh	1.6	INTEGER i,j,k,bi,bj
55	jmc	1.73	INTEGER ks
56			INTEGER numIters
57	adcroft	1.22	_RL firstResidual,lastResidual
58	jmc	1.36	_RL tmpFac
59	jmc	1.65	_RL sumEmP, tileEmP(nSx,nSy)
60	jmc	1.47	LOGICAL putPmEinXvector
61	jmc	1.73	INTEGER ioUnit
62	jmc	1.61	CHARACTER*10 sufx
63	adcroft	1.25	CHARACTER*(MAX_LEN_MBUF) msgBuf
64	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
65	jmc	1.71	LOGICAL zeroPsNH, zeroMeanPnh, oldFreeSurfTerm
66	jmc	1.63	#else
67			_RL cg3d_b(1)
68	jmc	1.49	#endif
69	cnh	1.27	CEOP
70	jmc	1.17
71	jmc	1.49	#ifdef ALLOW_NONHYDROSTATIC
72	jmc	1.68	zeroPsNH = .FALSE.
73	jmc	1.71	c zeroPsNH = use3Dsolver .AND. exactConserv
74			c & .AND. select_rStar.EQ.0
75			zeroMeanPnh = .FALSE.
76			c zeroMeanPnh = use3Dsolver .AND. select_rStar.NE.0
77	jmc	1.72	c oldFreeSurfTerm = use3Dsolver .AND. select_rStar.EQ.0
78			c & .AND. .NOT.zeroPsNH
79			oldFreeSurfTerm = use3Dsolver .AND. .NOT.exactConserv
80	jmc	1.63	#else
81			cg3d_b(1) = 0.
82	jmc	1.49	#endif
83
84	jmc	1.58	C deepAtmosphere & useRealFreshWaterFlux: only valid if deepFac2F(ksurf)=1
85			C anelastic (always Z-coordinate):
86			C 1) assume that rhoFacF(1)=1 (and ksurf == 1);
87			C (this reduces the number of lines of code to modify)
88			C 2) (a) 2-D continuity eq. compute div. of mass transport (<- add rhoFac)
89			C (b) gradient of surf.Press in momentum eq. (<- add 1/rhoFac)
90			C => 2 factors cancel in elliptic eq. for Phi_s ,
91			C but 1rst factor(a) remains in RHS cg2d_b.
92
93	jmc	1.47	C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
94			C instead of simply etaN ; This can speed-up the solver convergence in
95			C the case where \|Global_mean_PmE\| is large.
96			putPmEinXvector = .FALSE.
97	jmc	1.64	c putPmEinXvector = useRealFreshWaterFlux.AND.fluidIsWater
98	jmc	1.47
99	jmc	1.71	IF ( myIter.EQ.1+nIter0 .AND. debugLevel .GE. debLevA ) THEN
100			_BEGIN_MASTER( myThid )
101			ioUnit = standardMessageUnit
102			WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
103			& ' putPmEinXvector =', putPmEinXvector
104			CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
105			#ifdef ALLOW_NONHYDROSTATIC
106			WRITE(msgBuf,'(A,2(A,L5))') 'SOLVE_FOR_PRESSURE:',
107			& ' zeroPsNH=', zeroPsNH, ' , zeroMeanPnh=', zeroMeanPnh
108			CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
109			WRITE(msgBuf,'(2A,L5)') 'SOLVE_FOR_PRESSURE:',
110			& ' oldFreeSurfTerm =', oldFreeSurfTerm
111			CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
112			#endif
113			_END_MASTER( myThid )
114			ENDIF
115
116	jmc	1.17	C-- Save previous solution & Initialise Vector solution and source term :
117	jmc	1.47	sumEmP = 0.
118	jmc	1.17	DO bj=myByLo(myThid),myByHi(myThid)
119			DO bi=myBxLo(myThid),myBxHi(myThid)
120			DO j=1-OLy,sNy+OLy
121			DO i=1-OLx,sNx+OLx
122	edhill	1.40	#ifdef ALLOW_CD_CODE
123	jmc	1.17	etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
124	jmc	1.26	#endif
125	jmc	1.18	cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
126	jmc	1.17	cg2d_b(i,j,bi,bj) = 0.
127			ENDDO
128			ENDDO
129	jmc	1.64	IF (useRealFreshWaterFlux.AND.fluidIsWater) THEN
130	jmc	1.62	tmpFac = freeSurfFac*mass2rUnit
131	jmc	1.58	IF (exactConserv)
132	jmc	1.62	& tmpFac = freeSurfFacmass2rUnitimplicDiv2DFlow
133	jmc	1.29	DO j=1,sNy
134			DO i=1,sNx
135	jmc	1.58	cg2d_b(i,j,bi,bj) =
136	jmc	1.29	& tmpFac_rA(i,j,bi,bj)EmPmR(i,j,bi,bj)/deltaTMom
137			ENDDO
138			ENDDO
139			ENDIF
140	jmc	1.47	IF ( putPmEinXvector ) THEN
141	jmc	1.65	tileEmP(bi,bj) = 0.
142	jmc	1.47	DO j=1,sNy
143			DO i=1,sNx
144	jmc	1.67	tileEmP(bi,bj) = tileEmP(bi,bj)
145	jmc	1.65	& + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
146	jmc	1.74	& *maskInC(i,j,bi,bj)
147	jmc	1.47	ENDDO
148			ENDDO
149			ENDIF
150	jmc	1.17	ENDDO
151			ENDDO
152	jmc	1.47	IF ( putPmEinXvector ) THEN
153	jmc	1.65	CALL GLOBAL_SUM_TILE_RL( tileEmP, sumEmP, myThid )
154	jmc	1.47	ENDIF
155	adcroft	1.12
156			DO bj=myByLo(myThid),myByHi(myThid)
157			DO bi=myBxLo(myThid),myBxHi(myThid)
158	jmc	1.47	IF ( putPmEinXvector ) THEN
159			tmpFac = 0.
160	jmc	1.62	IF (globalArea.GT.0.) tmpFac =
161			& freeSurfFacdeltaTfreesurfmass2rUnit*sumEmP/globalArea
162	jmc	1.47	DO j=1,sNy
163			DO i=1,sNx
164			cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
165			& - tmpFac*Bo_surf(i,j,bi,bj)
166			ENDDO
167			ENDDO
168			ENDIF
169	jmc	1.58	C- RHS: similar to the divergence of the vertically integrated mass transport:
170			C del_i { Sum_k [ rhoFac.(dr.hFac).(dy.deepFac).(u*) ] } / deltaT
171	jmc	1.63	DO k=Nr,1,-1
172	adcroft	1.12	CALL CALC_DIV_GHAT(
173	jmc	1.63	I bi,bj,k,
174			U cg2d_b, cg3d_b,
175			I myThid )
176	adcroft	1.12	ENDDO
177			ENDDO
178			ENDDO
179	cnh	1.4
180	adcroft	1.12	DO bj=myByLo(myThid),myByHi(myThid)
181			DO bi=myBxLo(myThid),myBxHi(myThid)
182	adcroft	1.13	#ifdef ALLOW_NONHYDROSTATIC
183	jmc	1.71	IF ( oldFreeSurfTerm ) THEN
184	jmc	1.73	C-- Add source term arising from w=d/dt (p_s + p_nh)
185	jmc	1.28	DO j=1,sNy
186			DO i=1,sNx
187	jmc	1.77	ks = kSurfC(i,j,bi,bj)
188	jmc	1.51	IF ( ks.LE.Nr ) THEN
189			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
190	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
191			& /deltaTMom/deltaTfreesurf
192	jmc	1.28	& *( etaN(i,j,bi,bj)
193	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
194	jmc	1.51	cg3d_b(i,j,ks,bi,bj) = cg3d_b(i,j,ks,bi,bj)
195	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
196			& /deltaTMom/deltaTfreesurf
197	jmc	1.28	& *( etaN(i,j,bi,bj)
198	jmc	1.59	& +phi_nh(i,j,ks,bi,bj)*recip_Bo(i,j,bi,bj) )
199	jmc	1.51	ENDIF
200	jmc	1.28	ENDDO
201	adcroft	1.12	ENDDO
202	jmc	1.28	ELSEIF ( exactConserv ) THEN
203	adcroft	1.13	#else
204	jmc	1.73	C-- Add source term arising from w=d/dt (p_s)
205	jmc	1.26	IF ( exactConserv ) THEN
206	edhill	1.39	#endif /* ALLOW_NONHYDROSTATIC */
207	jmc	1.26	DO j=1,sNy
208			DO i=1,sNx
209	jmc	1.77	ks = kSurfC(i,j,bi,bj)
210	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
211	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
212			& /deltaTMom/deltaTfreesurf
213	jmc	1.26	& * etaH(i,j,bi,bj)
214			ENDDO
215			ENDDO
216			ELSE
217			DO j=1,sNy
218			DO i=1,sNx
219	jmc	1.77	ks = kSurfC(i,j,bi,bj)
220	jmc	1.26	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
221	jmc	1.58	& -freeSurfFac_rA(i,j,bi,bj)deepFac2F(ks)
222			& /deltaTMom/deltaTfreesurf
223	jmc	1.26	& * etaN(i,j,bi,bj)
224			ENDDO
225	adcroft	1.12	ENDDO
226	jmc	1.26	ENDIF
227	adcroft	1.12
228			#ifdef ALLOW_OBCS
229	jmc	1.76	C- Note: solver matrix is trivial outside OB region (main diagonal only)
230			C => no real need to reset RHS (=cg2d_b) & cg2d_x, except that:
231			C a) normalisation is fct of Max(RHS), which can be large ouside OB region
232			C (would be different if we were solving for increment of eta/g
233			C instead of directly for eta/g).
234			C => need to reset RHS to ensure that interior solution does not depend
235			C on ouside OB region.
236			C b) provide directly the trivial solution cg2d_x == 0 for outside OB region
237			C (=> no residual => no effect on solver convergence and interior solution)
238	adcroft	1.14	IF (useOBCS) THEN
239	jmc	1.76	DO j=1,sNy
240			DO i=1,sNx
241			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)*maskInC(i,j,bi,bj)
242			cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)*maskInC(i,j,bi,bj)
243	adcroft	1.12	ENDDO
244			ENDDO
245			ENDIF
246	jmc	1.49	#endif /* ALLOW_OBCS */
247			C- end bi,bj loops
248	adcroft	1.12	ENDDO
249			ENDDO
250
251	edhill	1.42	#ifdef ALLOW_DEBUG
252	jmc	1.77	IF ( debugLevel .GE. debLevD ) THEN
253	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
254			& myThid)
255	adcroft	1.24	ENDIF
256	adcroft	1.23	#endif
257	jmc	1.61	IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
258			WRITE(sufx,'(I10.10)') myIter
259	jmc	1.67	CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
260	jmc	1.61	ENDIF
261	adcroft	1.12
262	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
263			C-- gradient solver.
264	adcroft	1.22	C see CG2D.h for the interface to this routine.
265			firstResidual=0.
266			lastResidual=0.
267	adcroft	1.19	numIters=cg2dMaxIters
268	jmc	1.50	c CALL TIMER_START('CG2D [SOLVE_FOR_PRESSURE]',myThid)
269	heimbach	1.56	#ifdef ALLOW_CG2D_NSA
270			C-- Call the not-self-adjoint version of cg2d
271			CALL CG2D_NSA(
272			U cg2d_b,
273			U cg2d_x,
274			O firstResidual,
275			O lastResidual,
276			U numIters,
277			I myThid )
278			#else /* not ALLOW_CG2D_NSA = default */
279	mlosch	1.69	#ifdef ALLOW_SRCG
280			IF ( useSRCGSolver ) THEN
281			C-- Call the single reduce CG solver
282			CALL CG2D_SR(
283	adcroft	1.22	U cg2d_b,
284	cnh	1.6	U cg2d_x,
285	adcroft	1.22	O firstResidual,
286			O lastResidual,
287	adcroft	1.19	U numIters,
288	cnh	1.1	I myThid )
289	mlosch	1.69	ELSE
290			#else
291			IF (.TRUE.) THEN
292			C-- Call the default CG solver
293			#endif /* ALLOW_SRCG */
294			CALL CG2D(
295			U cg2d_b,
296			U cg2d_x,
297			O firstResidual,
298			O lastResidual,
299			U numIters,
300			I myThid )
301			ENDIF
302	heimbach	1.56	#endif /* ALLOW_CG2D_NSA */
303	jmc	1.67	_EXCH_XY_RL( cg2d_x, myThid )
304	jmc	1.50	c CALL TIMER_STOP ('CG2D [SOLVE_FOR_PRESSURE]',myThid)
305	adcroft	1.23
306	edhill	1.42	#ifdef ALLOW_DEBUG
307	jmc	1.77	IF ( debugLevel .GE. debLevD ) THEN
308	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
309			& myThid)
310	adcroft	1.24	ENDIF
311	adcroft	1.23	#endif
312	cnh	1.1
313	jmc	1.32	C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
314	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
315	jmc	1.45	& ) THEN
316	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
317			_BEGIN_MASTER( myThid )
318			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
319			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
320			WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
321			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
322			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
323			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
324	edhill	1.43	_END_MASTER( myThid )
325	heimbach	1.38	ENDIF
326	jmc	1.32	ENDIF
327	jmc	1.17
328			C-- Transfert the 2D-solution to "etaN" :
329			DO bj=myByLo(myThid),myByHi(myThid)
330			DO bi=myBxLo(myThid),myBxHi(myThid)
331			DO j=1-OLy,sNy+OLy
332			DO i=1-OLx,sNx+OLx
333	jmc	1.18	etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
334	jmc	1.17	ENDDO
335			ENDDO
336			ENDDO
337			ENDDO
338	adcroft	1.10
339	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
340	jmc	1.53	IF ( use3Dsolver ) THEN
341	jmc	1.67	IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
342			WRITE(sufx,'(I10.10)') myIter
343			CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
344			ENDIF
345	adcroft	1.9
346			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
347			C see CG3D.h for the interface to this routine.
348	jmc	1.71
349	jmc	1.73	C-- Finish updating cg3d_b: 1) Add EmPmR contribution to top level cg3d_b:
350			C 2) Update or Add free-surface contribution
351			C 3) increment in horiz velocity due to new cg2d_x
352			C 4) add vertical velocity contribution.
353			CALL PRE_CG3D(
354			I oldFreeSurfTerm,
355			I cg2d_x,
356			U cg3d_b,
357			I myTime, myIter, myThid )
358	adcroft	1.9
359	jmc	1.67	#ifdef ALLOW_DEBUG
360	jmc	1.77	IF ( debugLevel .GE. debLevD ) THEN
361	jmc	1.73	CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
362			& myThid)
363			ENDIF
364	jmc	1.67	#endif
365			IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
366			WRITE(sufx,'(I10.10)') myIter
367	jmc	1.73	CALL WRITE_FLD_XYZ_RL('cg3d_b.',sufx, cg3d_b, myIter,myThid )
368	jmc	1.67	ENDIF
369	jmc	1.49
370	jmc	1.73	firstResidual=0.
371			lastResidual=0.
372			numIters=cg3dMaxIters
373			CALL TIMER_START('CG3D [SOLVE_FOR_PRESSURE]',myThid)
374			CALL CG3D(
375			U cg3d_b,
376			U phi_nh,
377			O firstResidual,
378			O lastResidual,
379			U numIters,
380			I myIter, myThid )
381			_EXCH_XYZ_RL( phi_nh, myThid )
382			CALL TIMER_STOP ('CG3D [SOLVE_FOR_PRESSURE]',myThid)
383
384			IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
385			& ) THEN
386			IF ( debugLevel .GE. debLevA ) THEN
387			_BEGIN_MASTER( myThid )
388			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
389			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
390			WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
391			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
392			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
393			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
394			_END_MASTER( myThid )
395			ENDIF
396			ENDIF
397	jmc	1.49
398	jmc	1.73	C-- Separate the Hydrostatic Surface Pressure adjusment (=> put it in dPhiNH)
399			C from the Non-hydrostatic pressure (since cg3d_x contains both contribution)
400			IF ( nonHydrostatic .AND. exactConserv ) THEN
401			IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
402			WRITE(sufx,'(I10.10)') myIter
403			CALL WRITE_FLD_XYZ_RL('cg3d_x.',sufx, phi_nh, myIter,myThid )
404			ENDIF
405			CALL POST_CG3D(
406			I zeroPsNH, zeroMeanPnh,
407			I myTime, myIter, myThid )
408			ENDIF
409	jmc	1.49
410			ENDIF
411			#endif /* ALLOW_NONHYDROSTATIC */
412	cnh	1.1
413	heimbach	1.60	#ifdef ALLOW_SHOWFLOPS
414			CALL SHOWFLOPS_INSOLVE( myThid)
415	ce107	1.52	#endif
416	heimbach	1.60
417	cnh	1.1	RETURN
418			END