model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.46 2005/05/15 03:02: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 "DYNVARS.h"
#ifdef ALLOW_CD_CODE
#include "CD_CODE_VARS.h"
#endif
#include "GRID.h"
#include "SURFACE.h"
#include "FFIELDS.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "SOLVE_FOR_PRESSURE3D.h"
#include "GW.h"
#endif
#ifdef ALLOW_OBCS
#include "OBCS.h"
#endif
#include "SOLVE_FOR_PRESSURE.h"

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
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

#ifdef TIME_PER_TIMESTEP
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

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
        DO K=Nr,1,-1
         DO j=1,sNy+1
          DO i=1,sNx+1
           uf(i,j) = _dyG(i,j,bi,bj)
     &      *drF(k)*_hFacW(i,j,k,bi,bj)
           vf(i,j) = _dxG(i,j,bi,bj)
     &      *drF(k)*_hFacS(i,j,k,bi,bj)
          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 ( nonHydrostatic ) THEN
         DO j=1,sNy
          DO i=1,sNx
           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
     &         *( etaN(i,j,bi,bj)
     &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
           cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
     &         *( etaN(i,j,bi,bj)
     &           +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
          ENDDO
         ENDDO
        ELSEIF ( exactConserv ) THEN
#else
        IF ( exactConserv ) THEN
#endif /* ALLOW_NONHYDROSTATIC */
         DO j=1,sNy
          DO i=1,sNx
           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
     &         * etaH(i,j,bi,bj)
          ENDDO
         ENDDO
        ELSE
         DO j=1,sNy
          DO i=1,sNx
           cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &       -freeSurfFac*_rA(i,j,bi,bj)/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
       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
      CALL CG2D(
     U           cg2d_b,
     U           cg2d_x,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
      _EXCH_XY_R8(cg2d_x, 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 ( nonHydrostatic ) 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

         K=1
         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 )
     &       +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
     &          -wVel(i,j,k+1,bi,bj)
     &        )*_rA(i,j,bi,bj)/deltaTmom
          ENDDO
         ENDDO
         DO K=2,Nr-1
          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)
     &         -wVel(i,j,k+1,bi,bj)
     &        )*_rA(i,j,bi,bj)/deltaTmom

           ENDDO
          ENDDO
         ENDDO
         K=Nr
         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)
     &        )*_rA(i,j,bi,bj)/deltaTmom

          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

        ENDDO ! bi
       ENDDO ! bj

      firstResidual=0.
      lastResidual=0.
      numIters=cg2dMaxIters
      CALL CG3D(
     U           cg3d_b,
     U           phi_nh,
     O           firstResidual,
     O           lastResidual,
     U           numIters,
     I           myThid )
      _EXCH_XYZ_R8(phi_nh, 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

      ENDIF
#endif

#ifdef TIME_PER_TIMESTEP
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

      RETURN
      END

#ifdef TIME_PER_TIMESTEP
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
1	jmc	1.47	C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.46 2005/05/15 03:02: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			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			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			#include "DYNVARS.h"
27	edhill	1.41	#ifdef ALLOW_CD_CODE
28			#include "CD_CODE_VARS.h"
29			#endif
30	adcroft	1.12	#include "GRID.h"
31	jmc	1.17	#include "SURFACE.h"
32	jmc	1.28	#include "FFIELDS.h"
33	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
34	adcroft	1.25	#include "SOLVE_FOR_PRESSURE3D.h"
35	adcroft	1.9	#include "GW.h"
36	adcroft	1.12	#endif
37	adcroft	1.11	#ifdef ALLOW_OBCS
38	adcroft	1.9	#include "OBCS.h"
39	adcroft	1.11	#endif
40	adcroft	1.22	#include "SOLVE_FOR_PRESSURE.h"
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.28	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			_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	adcroft	1.19	INTEGER numIters
65	adcroft	1.25	CHARACTER*(MAX_LEN_MBUF) msgBuf
66	cnh	1.27	CEOP
67	jmc	1.17
68	ce107	1.44	#ifdef TIME_PER_TIMESTEP
69			CCE107 common block for per timestep timing
70			C !TIMING VARIABLES
71			C == Timing variables ==
72			REAL*8 utnew, utold, stnew, stold, wtnew, wtold
73			COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
74			#endif
75
76	jmc	1.47	C-- Initialise the Vector solution with etaN + deltaT*Global_mean_PmE
77			C instead of simply etaN ; This can speed-up the solver convergence in
78			C the case where \|Global_mean_PmE\| is large.
79			putPmEinXvector = .FALSE.
80			c putPmEinXvector = useRealFreshWaterFlux
81
82	jmc	1.17	C-- Save previous solution & Initialise Vector solution and source term :
83	jmc	1.47	sumEmP = 0.
84	jmc	1.17	DO bj=myByLo(myThid),myByHi(myThid)
85			DO bi=myBxLo(myThid),myBxHi(myThid)
86			DO j=1-OLy,sNy+OLy
87			DO i=1-OLx,sNx+OLx
88	edhill	1.40	#ifdef ALLOW_CD_CODE
89	jmc	1.17	etaNm1(i,j,bi,bj) = etaN(i,j,bi,bj)
90	jmc	1.26	#endif
91	jmc	1.18	cg2d_x(i,j,bi,bj) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
92	jmc	1.17	cg2d_b(i,j,bi,bj) = 0.
93			ENDDO
94			ENDDO
95	jmc	1.29	IF (useRealFreshWaterFlux) THEN
96	jmc	1.36	tmpFac = freeSurfFac*convertEmP2rUnit
97	mlosch	1.35	IF (exactConserv)
98	jmc	1.36	& tmpFac = freeSurfFacconvertEmP2rUnitimplicDiv2DFlow
99	jmc	1.29	DO j=1,sNy
100			DO i=1,sNx
101			cg2d_b(i,j,bi,bj) =
102			& tmpFac_rA(i,j,bi,bj)EmPmR(i,j,bi,bj)/deltaTMom
103			ENDDO
104			ENDDO
105			ENDIF
106	jmc	1.47	IF ( putPmEinXvector ) THEN
107			tileEmP = 0.
108			DO j=1,sNy
109			DO i=1,sNx
110			tileEmP = tileEmP + rA(i,j,bi,bj)*EmPmR(i,j,bi,bj)
111			& *maskH(i,j,bi,bj)
112			ENDDO
113			ENDDO
114			sumEmP = sumEmP + tileEmP
115			ENDIF
116	jmc	1.17	ENDDO
117			ENDDO
118	jmc	1.47	IF ( putPmEinXvector ) THEN
119			_GLOBAL_SUM_R8( sumEmP, myThid )
120			ENDIF
121	adcroft	1.12
122			DO bj=myByLo(myThid),myByHi(myThid)
123			DO bi=myBxLo(myThid),myBxHi(myThid)
124	jmc	1.47	IF ( putPmEinXvector ) THEN
125			tmpFac = 0.
126			IF (globalArea.GT.0.) tmpFac = freeSurfFac*deltaTfreesurf
127			& convertEmP2rUnitsumEmP/globalArea
128			DO j=1,sNy
129			DO i=1,sNx
130			cg2d_x(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
131			& - tmpFac*Bo_surf(i,j,bi,bj)
132			ENDDO
133			ENDDO
134			ENDIF
135	adcroft	1.12	DO K=Nr,1,-1
136			DO j=1,sNy+1
137			DO i=1,sNx+1
138			uf(i,j) = _dyG(i,j,bi,bj)
139			& drF(k)_hFacW(i,j,k,bi,bj)
140			vf(i,j) = _dxG(i,j,bi,bj)
141			& drF(k)_hFacS(i,j,k,bi,bj)
142			ENDDO
143			ENDDO
144			CALL CALC_DIV_GHAT(
145			I bi,bj,1,sNx,1,sNy,K,
146			I uf,vf,
147	jmc	1.17	U cg2d_b,
148	adcroft	1.12	I myThid)
149			ENDDO
150			ENDDO
151			ENDDO
152	cnh	1.4
153	adcroft	1.12	C-- Add source term arising from w=d/dt (p_s + p_nh)
154			DO bj=myByLo(myThid),myByHi(myThid)
155			DO bi=myBxLo(myThid),myBxHi(myThid)
156	adcroft	1.13	#ifdef ALLOW_NONHYDROSTATIC
157	jmc	1.28	IF ( nonHydrostatic ) THEN
158			DO j=1,sNy
159			DO i=1,sNx
160			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
161	adcroft	1.33	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
162	jmc	1.28	& *( etaN(i,j,bi,bj)
163			& +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
164			cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
165	adcroft	1.33	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
166	jmc	1.28	& *( etaN(i,j,bi,bj)
167			& +phi_nh(i,j,1,bi,bj)*horiVertRatio/gravity )
168			ENDDO
169	adcroft	1.12	ENDDO
170	jmc	1.28	ELSEIF ( exactConserv ) THEN
171	adcroft	1.13	#else
172	jmc	1.26	IF ( exactConserv ) THEN
173	edhill	1.39	#endif /* ALLOW_NONHYDROSTATIC */
174	jmc	1.26	DO j=1,sNy
175			DO i=1,sNx
176			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
177	adcroft	1.33	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
178	jmc	1.26	& * etaH(i,j,bi,bj)
179			ENDDO
180			ENDDO
181			ELSE
182			DO j=1,sNy
183			DO i=1,sNx
184			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
185	adcroft	1.33	& -freeSurfFac*_rA(i,j,bi,bj)/deltaTMom/deltaTfreesurf
186	jmc	1.26	& * etaN(i,j,bi,bj)
187			ENDDO
188	adcroft	1.12	ENDDO
189	jmc	1.26	ENDIF
190	adcroft	1.12
191			#ifdef ALLOW_OBCS
192	adcroft	1.14	IF (useOBCS) THEN
193	adcroft	1.12	DO i=1,sNx
194			C Northern boundary
195			IF (OB_Jn(I,bi,bj).NE.0) THEN
196			cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
197	jmc	1.31	cg2d_x(I,OB_Jn(I,bi,bj),bi,bj)=0.
198	adcroft	1.12	ENDIF
199			C Southern boundary
200			IF (OB_Js(I,bi,bj).NE.0) THEN
201			cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
202	jmc	1.31	cg2d_x(I,OB_Js(I,bi,bj),bi,bj)=0.
203	adcroft	1.12	ENDIF
204			ENDDO
205			DO j=1,sNy
206			C Eastern boundary
207			IF (OB_Ie(J,bi,bj).NE.0) THEN
208			cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
209	jmc	1.31	cg2d_x(OB_Ie(J,bi,bj),J,bi,bj)=0.
210	adcroft	1.12	ENDIF
211			C Western boundary
212			IF (OB_Iw(J,bi,bj).NE.0) THEN
213			cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
214	jmc	1.31	cg2d_x(OB_Iw(J,bi,bj),J,bi,bj)=0.
215	adcroft	1.12	ENDIF
216			ENDDO
217			ENDIF
218			#endif
219			ENDDO
220			ENDDO
221
222	edhill	1.42	#ifdef ALLOW_DEBUG
223	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
224	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_b,'cg2d_b (SOLVE_FOR_PRESSURE)',
225			& myThid)
226	adcroft	1.24	ENDIF
227	adcroft	1.23	#endif
228	adcroft	1.12
229	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
230			C-- gradient solver.
231	adcroft	1.22	C see CG2D.h for the interface to this routine.
232			firstResidual=0.
233			lastResidual=0.
234	adcroft	1.19	numIters=cg2dMaxIters
235	cnh	1.1	CALL CG2D(
236	adcroft	1.22	U cg2d_b,
237	cnh	1.6	U cg2d_x,
238	adcroft	1.22	O firstResidual,
239			O lastResidual,
240	adcroft	1.19	U numIters,
241	cnh	1.1	I myThid )
242	adcroft	1.19	_EXCH_XY_R8(cg2d_x, myThid )
243	adcroft	1.23
244	edhill	1.42	#ifdef ALLOW_DEBUG
245	heimbach	1.38	IF ( debugLevel .GE. debLevB ) THEN
246	adcroft	1.23	CALL DEBUG_STATS_RL(1,cg2d_x,'cg2d_x (SOLVE_FOR_PRESSURE)',
247			& myThid)
248	adcroft	1.24	ENDIF
249	adcroft	1.23	#endif
250	cnh	1.1
251	jmc	1.32	C- dump CG2D output at monitorFreq (to reduce size of STD-OUTPUT files) :
252	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
253	jmc	1.45	& ) THEN
254	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
255			_BEGIN_MASTER( myThid )
256			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_init_res =',firstResidual
257			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
258			WRITE(msgBuf,'(A34,I6)') 'cg2d_iters =',numIters
259			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
260			WRITE(msgBuf,'(A34,1PE24.14)') 'cg2d_res =',lastResidual
261			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
262	edhill	1.43	_END_MASTER( myThid )
263	heimbach	1.38	ENDIF
264	jmc	1.32	ENDIF
265	jmc	1.17
266			C-- Transfert the 2D-solution to "etaN" :
267			DO bj=myByLo(myThid),myByHi(myThid)
268			DO bi=myBxLo(myThid),myBxHi(myThid)
269			DO j=1-OLy,sNy+OLy
270			DO i=1-OLx,sNx+OLx
271	jmc	1.18	etaN(i,j,bi,bj) = recip_Bo(i,j,bi,bj)*cg2d_x(i,j,bi,bj)
272	jmc	1.17	ENDDO
273			ENDDO
274			ENDDO
275			ENDDO
276	adcroft	1.10
277	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
278			IF ( nonHydrostatic ) THEN
279
280			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
281			C see CG3D.h for the interface to this routine.
282			DO bj=myByLo(myThid),myByHi(myThid)
283			DO bi=myBxLo(myThid),myBxHi(myThid)
284			DO j=1,sNy+1
285			DO i=1,sNx+1
286	jmc	1.18	uf(i,j)=-_recip_dxC(i,j,bi,bj)*
287	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
288	jmc	1.18	vf(i,j)=-_recip_dyC(i,j,bi,bj)*
289	adcroft	1.9	& (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
290			ENDDO
291			ENDDO
292
293	adcroft	1.12	#ifdef ALLOW_OBCS
294	adcroft	1.14	IF (useOBCS) THEN
295	adcroft	1.9	DO i=1,sNx+1
296			C Northern boundary
297			IF (OB_Jn(I,bi,bj).NE.0) THEN
298			vf(I,OB_Jn(I,bi,bj))=0.
299			ENDIF
300			C Southern boundary
301			IF (OB_Js(I,bi,bj).NE.0) THEN
302			vf(I,OB_Js(I,bi,bj)+1)=0.
303			ENDIF
304			ENDDO
305			DO j=1,sNy+1
306			C Eastern boundary
307			IF (OB_Ie(J,bi,bj).NE.0) THEN
308			uf(OB_Ie(J,bi,bj),J)=0.
309			ENDIF
310			C Western boundary
311			IF (OB_Iw(J,bi,bj).NE.0) THEN
312			uf(OB_Iw(J,bi,bj)+1,J)=0.
313			ENDIF
314			ENDDO
315			ENDIF
316	adcroft	1.12	#endif
317	adcroft	1.9
318	adcroft	1.12	K=1
319			DO j=1,sNy
320			DO i=1,sNx
321			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
322			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
323			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
324			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
325			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
326	jmc	1.18	& +( freeSurfFac*etaN(i,j,bi,bj)/deltaTMom
327			& -wVel(i,j,k+1,bi,bj)
328	adcroft	1.12	& )*_rA(i,j,bi,bj)/deltaTmom
329			ENDDO
330			ENDDO
331			DO K=2,Nr-1
332	adcroft	1.9	DO j=1,sNy
333			DO i=1,sNx
334			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
335			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
336			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
337			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
338			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
339	adcroft	1.12	& +( wVel(i,j,k ,bi,bj)
340			& -wVel(i,j,k+1,bi,bj)
341			& )*_rA(i,j,bi,bj)/deltaTmom
342
343	adcroft	1.9	ENDDO
344			ENDDO
345			ENDDO
346	adcroft	1.12	K=Nr
347			DO j=1,sNy
348			DO i=1,sNx
349			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
350			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
351			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
352			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
353			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
354			& +( wVel(i,j,k ,bi,bj)
355			& )*_rA(i,j,bi,bj)/deltaTmom
356
357			ENDDO
358			ENDDO
359
360			#ifdef ALLOW_OBCS
361	adcroft	1.14	IF (useOBCS) THEN
362	adcroft	1.12	DO K=1,Nr
363			DO i=1,sNx
364			C Northern boundary
365			IF (OB_Jn(I,bi,bj).NE.0) THEN
366			cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
367			ENDIF
368			C Southern boundary
369			IF (OB_Js(I,bi,bj).NE.0) THEN
370			cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
371			ENDIF
372			ENDDO
373			DO j=1,sNy
374			C Eastern boundary
375			IF (OB_Ie(J,bi,bj).NE.0) THEN
376			cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
377			ENDIF
378			C Western boundary
379			IF (OB_Iw(J,bi,bj).NE.0) THEN
380			cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
381			ENDIF
382			ENDDO
383			ENDDO
384			ENDIF
385			#endif
386	adcroft	1.9
387			ENDDO ! bi
388			ENDDO ! bj
389
390	adcroft	1.25	firstResidual=0.
391			lastResidual=0.
392			numIters=cg2dMaxIters
393			CALL CG3D(
394			U cg3d_b,
395			U phi_nh,
396			O firstResidual,
397			O lastResidual,
398			U numIters,
399			I myThid )
400			_EXCH_XYZ_R8(phi_nh, myThid )
401
402	jmc	1.46	IF ( DIFFERENT_MULTIPLE(monitorFreq,myTime,deltaTClock)
403	jmc	1.45	& ) THEN
404	heimbach	1.38	IF ( debugLevel .GE. debLevA ) THEN
405			_BEGIN_MASTER( myThid )
406			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_init_res =',firstResidual
407			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
408			WRITE(msgBuf,'(A34,I6)') 'cg3d_iters =',numIters
409			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
410			WRITE(msgBuf,'(A34,1PE24.14)') 'cg3d_res =',lastResidual
411			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
412	edhill	1.43	_END_MASTER( myThid )
413	heimbach	1.38	ENDIF
414	mlosch	1.37	ENDIF
415	adcroft	1.9
416			ENDIF
417			#endif
418	cnh	1.1
419	ce107	1.44	#ifdef TIME_PER_TIMESTEP
420			CCE107 Time per timestep information
421			_BEGIN_MASTER( myThid )
422			CALL TIMER_GET_TIME( utnew, stnew, wtnew )
423			C Only output timing information after the 1st timestep
424			IF ( wtold .NE. 0.0D0 ) THEN
425			WRITE(msgBuf,'(A34,3F10.6)')
426			$ 'User, system and wallclock time:', utnew - utold,
427			$ stnew - stold, wtnew - wtold
428			CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1)
429			ENDIF
430			utold = utnew
431			stold = stnew
432			wtold = wtnew
433			_END_MASTER( myThid )
434			#endif
435
436	cnh	1.1	RETURN
437			END
438	ce107	1.44
439			#ifdef TIME_PER_TIMESTEP
440			CCE107 Initialization of common block for per timestep timing
441			BLOCK DATA settimers
442			C !TIMING VARIABLES
443			C == Timing variables ==
444			REAL*8 utnew, utold, stnew, stold, wtnew, wtold
445			COMMON /timevars/ utnew, utold, stnew, stold, wtnew, wtold
446			DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/
447			END
448			#endif