model/src/solve_for_pressure.F

C $Header: /u/gcmpack/MITgcm/model/src/solve_for_pressure.F,v 1.44 2005/01/13 00:46:33 ce107 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  DIFF_BASE_MULTIPLE
      EXTERNAL DIFF_BASE_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
      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--   Save previous solution & Initialise Vector solution and source term :
      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
       ENDDO
      ENDDO

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