model/src/calc_div_ghat.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_div_ghat.F,v 1.13 2001/03/06 16:51:02 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"
#undef ALLOW_ZONAL_FILT
#undef ALLOW_SHAP_FILT

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