model/src/solve_for_pressure.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.14 2001/02/02 21:04:48 adcroft Exp $
C $Name: $

#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 (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.
          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 (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 )
     &       +(
     &         -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 (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

       CALL CG3D( myThid )
       _EXCH_XYZ_R8(cg3d_x, myThid )

      ENDIF
#endif

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