model/src/solve_for_pressure.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.10 1999/05/24 14:15:15 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 "CG2D.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "CG3D.h"
#include "GW.h"
#ifdef ALLOW_OBCS
#include "OBCS.h"
#endif
#include "GRID.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
#ifdef ALLOW_NONHYDROSTATIC
      _RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
#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--   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

         IF (openBoundaries) THEN
          DO i=1,sNx+1
C Northern boundary
          IF (OB_Jn(I,bi,bj).NE.0) THEN
           uf(I,OB_Jn(I,bi,bj))=0.
           vf(I,OB_Jn(I,bi,bj))=0.
          ENDIF
C Southern boundary
          IF (OB_Js(I,bi,bj).NE.0) THEN
           uf(I,OB_Js(I,bi,bj))=0.
           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.
           vf(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.
           vf(OB_Iw(J,bi,bj),J)=0.
          ENDIF
          ENDDO
         ENDIF

         DO K=1,Nr
          DO j=1,sNy
           DO i=1,sNx
c           cg3d_x(i,j,k,bi,bj) = 0.
            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 )
           ENDDO
          ENDDO
         ENDDO

        ENDDO ! bi
       ENDDO ! bj

       CALL CG3D( myThid )

       _EXCH_XYZ_R8(cg3d_x, myThid )

      ENDIF
#endif

      RETURN
      END
1	adcroft	1.11	C $Header: /u/gcmpack/models/MITgcmUV/model/src/solve_for_pressure.F,v 1.10 1999/05/24 14:15:15 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			#include "CG2D.h"
20	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
21			#include "CG3D.h"
22			#include "GW.h"
23	adcroft	1.11	#ifdef ALLOW_OBCS
24	adcroft	1.9	#include "OBCS.h"
25	adcroft	1.11	#endif
26	adcroft	1.9	#include "GRID.h"
27			#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	#ifdef ALLOW_NONHYDROSTATIC
37			_RS uf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
38			_RS vf(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
39			#endif
40	cnh	1.4
41	cnh	1.7	#ifdef INCLUDE_CD_CODE
42	cnh	1.4	C-- Save previous solution.
43			DO bj=myByLo(myThid),myByHi(myThid)
44			DO bi=myBxLo(myThid),myBxHi(myThid)
45			DO j=1-OLy,sNy+OLy
46			DO i=1-OLx,sNx+OLx
47			cg2d_xNM1(i,j,bi,bj) = cg2d_x(i,j,bi,bj)
48			ENDDO
49			ENDDO
50			ENDDO
51			ENDDO
52			#endif
53	cnh	1.1
54			C-- Find the surface pressure using a two-dimensional conjugate
55			C-- gradient solver.
56			C see CG2D.h for the interface to this routine.
57			CALL CG2D(
58	cnh	1.6	I cg2d_b,
59			U cg2d_x,
60	cnh	1.1	I myThid )
61
62	adcroft	1.10	_EXCH_XY_R8(cg2d_x, myThid )
63
64	adcroft	1.9	#ifdef ALLOW_NONHYDROSTATIC
65			IF ( nonHydrostatic ) THEN
66
67			C-- Solve for a three-dimensional pressure term (NH or IGW or both ).
68			C see CG3D.h for the interface to this routine.
69			DO bj=myByLo(myThid),myByHi(myThid)
70			DO bi=myBxLo(myThid),myBxHi(myThid)
71			DO j=1,sNy+1
72			DO i=1,sNx+1
73			uf(i,j)=-gBaro_recip_dxC(i,j,bi,bj)
74			& (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))
75			vf(i,j)=-gBaro_recip_dyC(i,j,bi,bj)
76			& (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))
77			ENDDO
78			ENDDO
79
80			IF (openBoundaries) THEN
81			DO i=1,sNx+1
82			C Northern boundary
83			IF (OB_Jn(I,bi,bj).NE.0) THEN
84			uf(I,OB_Jn(I,bi,bj))=0.
85			vf(I,OB_Jn(I,bi,bj))=0.
86			ENDIF
87			C Southern boundary
88			IF (OB_Js(I,bi,bj).NE.0) THEN
89			uf(I,OB_Js(I,bi,bj))=0.
90			vf(I,OB_Js(I,bi,bj)+1)=0.
91			ENDIF
92			ENDDO
93			DO j=1,sNy+1
94			C Eastern boundary
95			IF (OB_Ie(J,bi,bj).NE.0) THEN
96			uf(OB_Ie(J,bi,bj),J)=0.
97			vf(OB_Ie(J,bi,bj),J)=0.
98			ENDIF
99			C Western boundary
100			IF (OB_Iw(J,bi,bj).NE.0) THEN
101			uf(OB_Iw(J,bi,bj)+1,J)=0.
102			vf(OB_Iw(J,bi,bj),J)=0.
103			ENDIF
104			ENDDO
105			ENDIF
106
107			DO K=1,Nr
108			DO j=1,sNy
109			DO i=1,sNx
110			c cg3d_x(i,j,k,bi,bj) = 0.
111			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
112			& +dRF(K)dYG(i+1,j,bi,bj)hFacW(i+1,j,k,bi,bj)*uf(i+1,j)
113			& -dRF(K)dYG( i ,j,bi,bj)hFacW( i ,j,k,bi,bj)*uf( i ,j)
114			& +dRF(K)dXG(i,j+1,bi,bj)hFacS(i,j+1,k,bi,bj)*vf(i,j+1)
115			& -dRF(K)dXG(i, j ,bi,bj)hFacS(i, j ,k,bi,bj)*vf(i, j )
116			ENDDO
117			ENDDO
118			ENDDO
119
120			ENDDO ! bi
121			ENDDO ! bj
122
123			CALL CG3D( myThid )
124	adcroft	1.10
125			_EXCH_XYZ_R8(cg3d_x, myThid )
126	adcroft	1.9
127			ENDIF
128			#endif
129	cnh	1.1
130			RETURN
131			END