model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm_contrib/nesting_sannino/code_nest_merged/solve_for_pressure.F,v 1.3 2010/11/28 01:56:55 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
      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
        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
#ifndef ALLOW_NEST_CHILD
         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 /* ALLOW_NEST_CHILD */
        ENDIF
#endif /* ALLOW_OBCS */
C-    end bi,bj loops
       ENDDO
      ENDDO

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevB ) THEN
       CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
     &                        myThid)
      ENDIF
#endif
      IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
       WRITE(sufx,'(I10.10)') myIter
       CALL WRITE_FLD_XY_RL( 'cg2d_b.', sufx, cg2d_b, myIter, myThid )
      ENDIF

C--   Find the surface pressure using a two-dimensional conjugate
C--   gradient solver.
C     see CG2D.h for the interface to this routine.
      firstResidual=0.
      lastResidual=0.
      numIters=cg2dMaxIters
c     CALL TIMER_START('CG2D   [SOLVE_FOR_PRESSURE]',myThid)
#ifdef ALLOW_CG2D_NSA
C--   Call the not-self-adjoint version of cg2d
      CALL CG2D_NSA(
     U           cg2d_b,
     U           cg2d_x,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
#else /* not ALLOW_CG2D_NSA = default */
#ifdef ALLOW_SRCG
      IF ( useSRCGSolver ) THEN
C--   Call the single reduce CG solver
       CALL CG2D_SR(
     U           cg2d_b,
     U           cg2d_x,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
      ELSE
#else
      IF (.TRUE.) THEN
C--   Call the default CG solver
#endif /* ALLOW_SRCG */
       CALL CG2D(
     U           cg2d_b,
     U           cg2d_x,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
      ENDIF
#endif /* ALLOW_CG2D_NSA */
      _EXCH_XY_RL( cg2d_x, myThid )
c     CALL TIMER_STOP ('CG2D   [SOLVE_FOR_PRESSURE]',myThid)

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevB ) THEN
       CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
     &                        myThid)
      ENDIF
#endif

C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
      IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
     &   ) THEN
       IF ( debugLevel .GE. debLevA ) THEN
        _BEGIN_MASTER( myThid )
        WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
        CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
        _END_MASTER( myThid )
       ENDIF
      ENDIF

C--   Transfert the 2D-solution to "etaN" :
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_NONHYDROSTATIC
      IF ( use3Dsolver ) THEN
       IF ( DIFFERENT_MULTIPLE(diagFreq, myTime, deltaTClock) ) THEN
        WRITE(sufx,'(I10.10)') myIter
        CALL WRITE_FLD_XY_RL( 'cg2d_x.',sufx, cg2d_x, myIter, myThid )
       ENDIF

C--   Solve for a three-dimensional pressure term (NH or IGW or both ).
C     see CG3D.h for the interface to this routine.

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. debLevB ) THEN
        CALL DEBUG_STATS_RL(Nr,cg3d_b,'cg3d_b (SOLVE_FOR_PRESSURE)',
     &                         myThid)
       ENDIF
#endif
       IF ( DIFFERENT_MULTIPLE( diagFreq, myTime, deltaTClock) ) THEN
        WRITE(sufx,'(I10.10)') myIter
        CALL WRITE_FLD_XYZ_RL('cg3d_b.',sufx, cg3d_b, myIter,myThid )
       ENDIF

       firstResidual=0.
       lastResidual=0.
       numIters=cg3dMaxIters
       CALL TIMER_START('CG3D   [SOLVE_FOR_PRESSURE]',myThid)
       CALL CG3D(
     U            cg3d_b,
     U            phi_nh,
     O            firstResidual,
     O            lastResidual,
     U            numIters,
     I            myIter, myThid )
       _EXCH_XYZ_RL( phi_nh, myThid )
       CALL TIMER_STOP ('CG3D   [SOLVE_FOR_PRESSURE]',myThid)

       IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
     &    ) THEN
        IF ( debugLevel .GE. debLevA ) THEN
         _BEGIN_MASTER( myThid )
         WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
         WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
         WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
         CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
         _END_MASTER( myThid )
        ENDIF
       ENDIF

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