model/src/solve_for_pressure.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.12 2000/03/14 17:47:26 adcroft Exp $

#include "CPP_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SOLVE_FOR_PRESSURE( myThid )
C     /==========================================================\
C     | SUBROUTINE SOLVE_FOR_PRESSURE                            |
C     | o Controls inversion of two and/or three-dimensional     |
C     |   elliptic problems for the pressure field.              |
C     \==========================================================/
      IMPLICIT NONE

C     == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "CG2D.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "CG3D.h"
#include "GW.h"
#endif
#ifdef ALLOW_OBCS
#include "OBCS.h"
#endif

C     == Routine arguments ==
C     myThid - Number of this instance of SOLVE_FOR_PRESSURE
      INTEGER myThid
CEndOfInterface

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)

#ifndef DIVG_IN_DYNAMICS
      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,
     I       myThid)
        ENDDO
       ENDDO
      ENDDO
#endif

#ifdef INCLUDE_CD_CODE
C--   Save previous solution.
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#endif

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
        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)*horiVertRatio*(
     &          -cg2d_x(I,J,bi,bj)
     &          -cg3d_x(I,J,1,bi,bj)
     &        )/deltaTMom/deltaTMom
          cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
     &      +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
     &         -cg2d_x(I,J,bi,bj)
     &         -cg3d_x(I,J,1,bi,bj)
     &       )/deltaTMom/deltaTMom
         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)*horiVertRatio*(
     &          -cg2d_x(I,J,bi,bj)
     &        )/deltaTMom/deltaTMom
         ENDDO
        ENDDO
#endif

#ifdef ALLOW_OBCS
        IF (openBoundaries) 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.
          ENDIF
C Southern boundary
          IF (OB_Js(I,bi,bj).NE.0) THEN
           cg2d_b(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.
          ENDIF
C Western boundary
          IF (OB_Iw(J,bi,bj).NE.0) THEN
           cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
          ENDIF
         ENDDO
        ENDIF
#endif
       ENDDO
      ENDDO


C--   Find the surface pressure using a two-dimensional conjugate
C--   gradient solver.
C     see CG2D.h for the interface to this routine.
      CALL CG2D(
     I           cg2d_b,
     U           cg2d_x,
     I           myThid )

      _EXCH_XY_R8(cg2d_x, myThid )

#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)=-gBaro*_recip_dxC(i,j,bi,bj)*
     &         (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
           vf(i,j)=-gBaro*_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 (openBoundaries) 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 )
     &       +(
     &         -wVel(i,j,k+1,bi,bj)
     &        )*_rA(i,j,bi,bj)/deltaTmom
     &       +freeSurfFac*_rA(I,J,bi,bj)*horiVertRatio*(
     &          +cg2d_x(I,J,bi,bj)
     &        )/deltaTMom/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 (openBoundaries) 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

       CALL CG3D( myThid )
       _EXCH_XYZ_R8(cg3d_x, myThid )

      ENDIF
#endif

      RETURN
      END
1	adcroft	1.13	C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.12 2000/03/14 17:47:26 adcroft Exp $
2	cnh	1.1
3	adcroft	1.5	#include "CPP_OPTIONS.h"
4	cnh	1.1
5			CStartOfInterface
6			SUBROUTINE SOLVE_FOR_PRESSURE( myThid )
7			C /==========================================================\
8			C \| SUBROUTINE SOLVE_FOR_PRESSURE \|
9			C \| o Controls inversion of two and/or three-dimensional \|
10			C \| elliptic problems for the pressure field. \|
11			C \==========================================================/
12	adcroft	1.8	IMPLICIT NONE
13	cnh	1.1
14	cnh	1.4	C == Global variables
15			#include "SIZE.h"
16			#include "EEPARAMS.h"
17			#include "PARAMS.h"
18			#include "DYNVARS.h"
19	adcroft	1.12	#include "GRID.h"
20	cnh	1.4	#include "CG2D.h"
21	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
22			#include "CG3D.h"
23			#include "GW.h"
24	adcroft	1.12	#endif
25	adcroft	1.11	#ifdef ALLOW_OBCS
26	adcroft	1.9	#include "OBCS.h"
27	adcroft	1.11	#endif
28	cnh	1.4
29	cnh	1.1	C == Routine arguments ==
30			C myThid - Number of this instance of SOLVE_FOR_PRESSURE
31			INTEGER myThid
32			CEndOfInterface
33	cnh	1.4
34			C Local variables
35	cnh	1.6	INTEGER i,j,k,bi,bj
36	adcroft	1.9	_RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37			_RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
38	adcroft	1.12
39			#ifndef DIVG_IN_DYNAMICS
40			DO bj=myByLo(myThid),myByHi(myThid)
41			DO bi=myBxLo(myThid),myBxHi(myThid)
42			DO K=Nr,1,-1
43			DO j=1,sNy+1
44			DO i=1,sNx+1
45			uf(i,j) = _dyG(i,j,bi,bj)
46			& drF(k)_hFacW(i,j,k,bi,bj)
47			vf(i,j) = _dxG(i,j,bi,bj)
48			& drF(k)_hFacS(i,j,k,bi,bj)
49			ENDDO
50			ENDDO
51			CALL CALC_DIV_GHAT(
52			I bi,bj,1,sNx,1,sNy,K,
53			I uf,vf,
54			I myThid)
55			ENDDO
56			ENDDO
57			ENDDO
58	adcroft	1.9	#endif
59	cnh	1.4
60	cnh	1.7	#ifdef INCLUDE_CD_CODE
61	cnh	1.4	C-- Save previous solution.
62			DO bj=myByLo(myThid),myByHi(myThid)
63			DO bi=myBxLo(myThid),myBxHi(myThid)
64			DO j=1-OLy,sNy+OLy
65			DO i=1-OLx,sNx+OLx
66			cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
67			ENDDO
68			ENDDO
69			ENDDO
70			ENDDO
71			#endif
72	cnh	1.1
73	adcroft	1.12	C-- Add source term arising from w=d/dt (p_s + p_nh)
74			DO bj=myByLo(myThid),myByHi(myThid)
75			DO bi=myBxLo(myThid),myBxHi(myThid)
76	adcroft	1.13	#ifdef ALLOW_NONHYDROSTATIC
77	adcroft	1.12	DO j=1,sNy
78			DO i=1,sNx
79			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
80			& +freeSurfFac_rA(I,J,bi,bj)horiVertRatio*(
81			& -cg2d_x(I,J,bi,bj)
82	adcroft	1.13	& -cg3d_x(I,J,1,bi,bj)
83	adcroft	1.12	& )/deltaTMom/deltaTMom
84	adcroft	1.13	cg3d_b(i,j,1,bi,bj) = cg3d_b(i,j,1,bi,bj)
85			& +freeSurfFac_rA(I,J,bi,bj)horiVertRatio*(
86			& -cg2d_x(I,J,bi,bj)
87			& -cg3d_x(I,J,1,bi,bj)
88			& )/deltaTMom/deltaTMom
89	adcroft	1.12	ENDDO
90			ENDDO
91	adcroft	1.13	#else
92	adcroft	1.12	DO j=1,sNy
93			DO i=1,sNx
94	adcroft	1.13	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
95			& +freeSurfFac_rA(I,J,bi,bj)horiVertRatio*(
96			& -cg2d_x(I,J,bi,bj)
97			& )/deltaTMom/deltaTMom
98	adcroft	1.12	ENDDO
99			ENDDO
100			#endif
101
102			#ifdef ALLOW_OBCS
103			IF (openBoundaries) THEN
104			DO i=1,sNx
105			C Northern boundary
106			IF (OB_Jn(I,bi,bj).NE.0) THEN
107			cg2d_b(I,OB_Jn(I,bi,bj),bi,bj)=0.
108			ENDIF
109			C Southern boundary
110			IF (OB_Js(I,bi,bj).NE.0) THEN
111			cg2d_b(I,OB_Js(I,bi,bj),bi,bj)=0.
112			ENDIF
113			ENDDO
114			DO j=1,sNy
115			C Eastern boundary
116			IF (OB_Ie(J,bi,bj).NE.0) THEN
117			cg2d_b(OB_Ie(J,bi,bj),J,bi,bj)=0.
118			ENDIF
119			C Western boundary
120			IF (OB_Iw(J,bi,bj).NE.0) THEN
121			cg2d_b(OB_Iw(J,bi,bj),J,bi,bj)=0.
122			ENDIF
123			ENDDO
124			ENDIF
125			#endif
126			ENDDO
127			ENDDO
128
129
130	cnh	1.1	C-- Find the surface pressure using a two-dimensional conjugate
131			C-- gradient solver.
132			C see CG2D.h for the interface to this routine.
133			CALL CG2D(
134	cnh	1.6	I cg2d_b,
135			U cg2d_x,
136	cnh	1.1	I myThid )
137
138	adcroft	1.10	_EXCH_XY_R8(cg2d_x, myThid )
139
140	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
141			IF ( nonHydrostatic ) THEN
142
143			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
144			C see CG3D.h for the interface to this routine.
145			DO bj=myByLo(myThid),myByHi(myThid)
146			DO bi=myBxLo(myThid),myBxHi(myThid)
147			DO j=1,sNy+1
148			DO i=1,sNx+1
149			uf(i,j)=-gBaro_recip_dxC(i,j,bi,bj)
150			& (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
151			vf(i,j)=-gBaro_recip_dyC(i,j,bi,bj)
152			& (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
153			ENDDO
154			ENDDO
155
156	adcroft	1.12	#ifdef ALLOW_OBCS
157	adcroft	1.9	IF (openBoundaries) THEN
158			DO i=1,sNx+1
159			C Northern boundary
160			IF (OB_Jn(I,bi,bj).NE.0) THEN
161			vf(I,OB_Jn(I,bi,bj))=0.
162			ENDIF
163			C Southern boundary
164			IF (OB_Js(I,bi,bj).NE.0) THEN
165			vf(I,OB_Js(I,bi,bj)+1)=0.
166			ENDIF
167			ENDDO
168			DO j=1,sNy+1
169			C Eastern boundary
170			IF (OB_Ie(J,bi,bj).NE.0) THEN
171			uf(OB_Ie(J,bi,bj),J)=0.
172			ENDIF
173			C Western boundary
174			IF (OB_Iw(J,bi,bj).NE.0) THEN
175			uf(OB_Iw(J,bi,bj)+1,J)=0.
176			ENDIF
177			ENDDO
178			ENDIF
179	adcroft	1.12	#endif
180	adcroft	1.9
181	adcroft	1.12	K=1
182			DO j=1,sNy
183			DO i=1,sNx
184			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
185			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
186			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
187			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
188			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
189			& +(
190			& -wVel(i,j,k+1,bi,bj)
191			& )*_rA(i,j,bi,bj)/deltaTmom
192			& +freeSurfFac_rA(I,J,bi,bj)horiVertRatio*(
193			& +cg2d_x(I,J,bi,bj)
194			& )/deltaTMom/deltaTMom
195			ENDDO
196			ENDDO
197			DO K=2,Nr-1
198	adcroft	1.9	DO j=1,sNy
199			DO i=1,sNx
200			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
201			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
202			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
203			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
204			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
205	adcroft	1.12	& +( wVel(i,j,k ,bi,bj)
206			& -wVel(i,j,k+1,bi,bj)
207			& )*_rA(i,j,bi,bj)/deltaTmom
208
209	adcroft	1.9	ENDDO
210			ENDDO
211			ENDDO
212	adcroft	1.12	K=Nr
213			DO j=1,sNy
214			DO i=1,sNx
215			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
216			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
217			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
218			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
219			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
220			& +( wVel(i,j,k ,bi,bj)
221			& )*_rA(i,j,bi,bj)/deltaTmom
222
223			ENDDO
224			ENDDO
225
226			#ifdef ALLOW_OBCS
227			IF (openBoundaries) THEN
228			DO K=1,Nr
229			DO i=1,sNx
230			C Northern boundary
231			IF (OB_Jn(I,bi,bj).NE.0) THEN
232			cg3d_b(I,OB_Jn(I,bi,bj),K,bi,bj)=0.
233			ENDIF
234			C Southern boundary
235			IF (OB_Js(I,bi,bj).NE.0) THEN
236			cg3d_b(I,OB_Js(I,bi,bj),K,bi,bj)=0.
237			ENDIF
238			ENDDO
239			DO j=1,sNy
240			C Eastern boundary
241			IF (OB_Ie(J,bi,bj).NE.0) THEN
242			cg3d_b(OB_Ie(J,bi,bj),J,K,bi,bj)=0.
243			ENDIF
244			C Western boundary
245			IF (OB_Iw(J,bi,bj).NE.0) THEN
246			cg3d_b(OB_Iw(J,bi,bj),J,K,bi,bj)=0.
247			ENDIF
248			ENDDO
249			ENDDO
250			ENDIF
251			#endif
252	adcroft	1.9
253			ENDDO ! bi
254			ENDDO ! bj
255
256			CALL CG3D( myThid )
257	adcroft	1.10	_EXCH_XYZ_R8(cg3d_x, myThid )
258	adcroft	1.9
259			ENDIF
260			#endif
261	cnh	1.1
262			RETURN
263			END