model/src/calc_div_ghat.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_div_ghat.F,v 1.12 2001/02/20 15:06:21 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

C     /==========================================================\
C     | S/R CALC_DIV_GHAT                                        |
C     | o Form the right hand-side of the surface pressure eqn.  |
C     |==========================================================|
C     \==========================================================/
      SUBROUTINE CALC_DIV_GHAT( 
     I        bi,bj,iMin,iMax,jMin,jMax,K,
     I        xA,yA,
     U        cg2d_b,
     I        myThid)

      IMPLICIT NONE

C     == Global variables ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_NONHYDROSTATIC
#include "CG3D.h"
#endif

C     == Routine arguments ==
C     bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
C                                      results will be set.
C     k                              - Index of layer.
C     xA, yA                         - Cell face areas
C     cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector 
C     myThid - Instance number for this call of CALC_DIV_GHAT
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER K
      _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER myThid

C     == Local variables ==
      INTEGER i,j
      _RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)

C--   Pressure equation source term
C     Term is the vertical integral of the divergence of the 
C     time tendency terms along with a relaxation term that
C     pulls div(U) + dh/dt back toward zero.

      IF (implicDiv2Dflow.EQ.1.) then
C     Fully Implicit treatment of the Barotropic Flow Divergence
        DO j=1,sNy
         DO i=1,sNx+1
          pf(i,j) = xA(i,j)*gUNm1(i,j,k,bi,bj) / deltaTmom
         ENDDO
        ENDDO
      ELSE
C     Explicit+Implicit part of the Barotropic Flow Divergence
C      => Filtering of uVel,vVel is necessary
#ifdef ALLOW_ZONAL_FILT
        IF (zonal_filt_lat.LT.90.) THEN 
          CALL ZONAL_FILTER(
     &      uVel, hFacW, 1-1, sNy+1, k, k, bi, bj, 1, myThid)
          CALL ZONAL_FILTER(
     &      vVel, hFacS, 1-1, sNy+1, k, k, bi, bj, 2, myThid)
        ENDIF
#endif
        DO j=1,sNy
         DO i=1,sNx+1
          pf(i,j) = ( implicDiv2Dflow * gUNm1(i,j,k,bi,bj)
     &          + (1.-implicDiv2Dflow) * uVel(i,j,k,bi,bj)
     &               ) * xA(i,j) / deltaTmom
         ENDDO
        ENDDO
      ENDIF
      DO j=1,sNy
       DO i=1,sNx
        cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
     &   pf(i+1,j)-pf(i,j)
       ENDDO
      ENDDO

#ifdef ALLOW_NONHYDROSTATIC
      IF (nonHydrostatic) THEN
       DO j=1,sNy
        DO i=1,sNx
         cg3d_b(i,j,k,bi,bj) = pf(i+1,j)-pf(i,j)
        ENDDO
       ENDDO
      ENDIF
#endif

      IF (implicDiv2Dflow.EQ.1.) then
C     Fully Implicit treatment of the Barotropic Flow Divergence
        DO j=1,sNy+1
         DO i=1,sNx
          pf(i,j) = yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom
         ENDDO
        ENDDO
      ELSE
C     Explicit+Implicit part of the Barotropic Flow Divergence
        DO j=1,sNy+1
         DO i=1,sNx
          pf(i,j) = ( implicDiv2Dflow * gVNm1(i,j,k,bi,bj)
     &          + (1.-implicDiv2Dflow) * vVel(i,j,k,bi,bj)
     &               ) * yA(i,j) / deltaTmom
         ENDDO
        ENDDO
      ENDIF
      DO j=1,sNy
       DO i=1,sNx
        cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
     &   pf(i,j+1)-pf(i,j)
       ENDDO
      ENDDO

#ifdef ALLOW_NONHYDROSTATIC
      IF (nonHydrostatic) THEN
       DO j=1,sNy
        DO i=1,sNx
         cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) +
     &    pf(i,j+1)-pf(i,j)
        ENDDO
       ENDDO
      ENDIF
#endif

      RETURN
      END
1	jmc	1.13	C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_div_ghat.F,v 1.12 2001/02/20 15:06:21 jmc Exp $
2	jmc	1.12	C $Name: $
3	cnh	1.1
4	cnh	1.6	#include "CPP_OPTIONS.h"
5	cnh	1.1
6			C /==========================================================\
7			C \| S/R CALC_DIV_GHAT \|
8			C \| o Form the right hand-side of the surface pressure eqn. \|
9	adcroft	1.9	C \|==========================================================\|
10	cnh	1.1	C \==========================================================/
11			SUBROUTINE CALC_DIV_GHAT(
12	jmc	1.13	I bi,bj,iMin,iMax,jMin,jMax,K,
13	cnh	1.1	I xA,yA,
14	jmc	1.13	U cg2d_b,
15	cnh	1.1	I myThid)
16
17			IMPLICIT NONE
18
19			C == Global variables ==
20			#include "SIZE.h"
21			#include "DYNVARS.h"
22			#include "FFIELDS.h"
23			#include "EEPARAMS.h"
24			#include "PARAMS.h"
25			#include "GRID.h"
26	adcroft	1.8	#ifdef ALLOW_NONHYDROSTATIC
27			#include "CG3D.h"
28			#endif
29	cnh	1.1
30			C == Routine arguments ==
31			C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
32			C results will be set.
33	adcroft	1.9	C k - Index of layer.
34			C xA, yA - Cell face areas
35	jmc	1.13	C cg2d_b - Conjugate Gradient 2-D solver : Right-hand side vector
36			C myThid - Instance number for this call of CALC_DIV_GHAT
37	cnh	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
38			INTEGER K
39			_RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40			_RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41	jmc	1.13	_RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
42	cnh	1.1	INTEGER myThid
43
44			C == Local variables ==
45			INTEGER i,j
46			_RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47
48			C-- Pressure equation source term
49			C Term is the vertical integral of the divergence of the
50			C time tendency terms along with a relaxation term that
51			C pulls div(U) + dh/dt back toward zero.
52
53	jmc	1.12	IF (implicDiv2Dflow.EQ.1.) then
54			C Fully Implicit treatment of the Barotropic Flow Divergence
55			DO j=1,sNy
56			DO i=1,sNx+1
57			pf(i,j) = xA(i,j)*gUNm1(i,j,k,bi,bj) / deltaTmom
58			ENDDO
59			ENDDO
60			ELSE
61			C Explicit+Implicit part of the Barotropic Flow Divergence
62			C => Filtering of uVel,vVel is necessary
63			#ifdef ALLOW_ZONAL_FILT
64			IF (zonal_filt_lat.LT.90.) THEN
65			CALL ZONAL_FILTER(
66			& uVel, hFacW, 1-1, sNy+1, k, k, bi, bj, 1, myThid)
67			CALL ZONAL_FILTER(
68			& vVel, hFacS, 1-1, sNy+1, k, k, bi, bj, 2, myThid)
69			ENDIF
70			#endif
71			DO j=1,sNy
72			DO i=1,sNx+1
73			pf(i,j) = ( implicDiv2Dflow * gUNm1(i,j,k,bi,bj)
74			& + (1.-implicDiv2Dflow) * uVel(i,j,k,bi,bj)
75			& ) * xA(i,j) / deltaTmom
76			ENDDO
77			ENDDO
78			ENDIF
79	cnh	1.1	DO j=1,sNy
80			DO i=1,sNx
81	adcroft	1.8	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
82	cnh	1.1	& pf(i+1,j)-pf(i,j)
83			ENDDO
84			ENDDO
85
86	adcroft	1.8	#ifdef ALLOW_NONHYDROSTATIC
87			IF (nonHydrostatic) THEN
88			DO j=1,sNy
89			DO i=1,sNx
90			cg3d_b(i,j,k,bi,bj) = pf(i+1,j)-pf(i,j)
91			ENDDO
92			ENDDO
93			ENDIF
94			#endif
95
96	jmc	1.12	IF (implicDiv2Dflow.EQ.1.) then
97			C Fully Implicit treatment of the Barotropic Flow Divergence
98			DO j=1,sNy+1
99			DO i=1,sNx
100			pf(i,j) = yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom
101			ENDDO
102			ENDDO
103			ELSE
104			C Explicit+Implicit part of the Barotropic Flow Divergence
105			DO j=1,sNy+1
106			DO i=1,sNx
107			pf(i,j) = ( implicDiv2Dflow * gVNm1(i,j,k,bi,bj)
108			& + (1.-implicDiv2Dflow) * vVel(i,j,k,bi,bj)
109			& ) * yA(i,j) / deltaTmom
110			ENDDO
111			ENDDO
112			ENDIF
113	cnh	1.1	DO j=1,sNy
114			DO i=1,sNx
115	adcroft	1.8	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
116	cnh	1.1	& pf(i,j+1)-pf(i,j)
117			ENDDO
118			ENDDO
119	cnh	1.4
120	adcroft	1.8	#ifdef ALLOW_NONHYDROSTATIC
121			IF (nonHydrostatic) THEN
122			DO j=1,sNy
123			DO i=1,sNx
124			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) +
125			& pf(i,j+1)-pf(i,j)
126			ENDDO
127			ENDDO
128			ENDIF
129			#endif
130	cnh	1.1
131			RETURN
132			END