model/src/calc_div_ghat.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_div_ghat.F,v 1.26 2009/11/30 16:26:48 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CALC_DIV_GHAT
C     !INTERFACE:
      SUBROUTINE CALC_DIV_GHAT(
     I                bi,bj,k,
     U                cg2d_b, cg3d_b,
     I                myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R CALC_DIV_GHAT
C     | o Form the right hand-side of the surface pressure eqn.
C     *==========================================================*
C     | Right hand side of pressure equation is divergence
C     | of veclocity tendency (GHAT) term along with a relaxation
C     | term equal to the barotropic flow field divergence.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_ADDFLUID
# include "FFIELDS.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     bi, bj  :: tile indices
C     k       :: Index of layer.
C     cg2d_b  :: Conjugate Gradient 2-D solver : Right-hand side vector
C     cg3d_b  :: Conjugate Gradient 3-D solver : Right-hand side vector
C     myThid  :: Instance number for this call of CALC_DIV_GHAT
      INTEGER bi,bj
      INTEGER k
      _RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#ifdef ALLOW_NONHYDROSTATIC
      _RL cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
#else
      _RL cg3d_b(1)
#endif
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     i,j    :: Loop counters
C     xA, yA :: Cell vertical face areas
C     pf     :: Intermediate array for building RHS source term.
      INTEGER i,j
      _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

C     Calculate vertical face areas
      DO j=1,sNy+1
        DO i=1,sNx+1
          xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
     &              *drF(k)*_hFacW(i,j,k,bi,bj)*rhoFacC(k)
          yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
     &              *drF(k)*_hFacS(i,j,k,bi,bj)*rhoFacC(k)
        ENDDO
      ENDDO

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)*gU(i,j,k,bi,bj) / deltaTMom
         ENDDO
        ENDDO
      ELSEIF (exactConserv) THEN
c     ELSEIF (nonlinFreeSurf.GT.0) THEN
C     Implicit treatment of the Barotropic Flow Divergence
        DO j=1,sNy
         DO i=1,sNx+1
          pf(i,j) = implicDiv2Dflow
     &             *xA(i,j)*gU(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
C-- Now the filter are applied in the_correction_step().
C   We have left this code here to indicate where the filters used to be
C   in the algorithm before JMC moved them to after the pressure solver.
c#ifdef ALLOW_ZONAL_FILT
c        IF (zonal_filt_lat.LT.90.) THEN
c          CALL ZONAL_FILTER(
c    U                     uVel( 1-OLx,1-OLy,k,bi,bj),
c    I                     hFacW(1-OLx,1-OLy,k,bi,bj),
c    I                     0, sNy+1, 1, bi, bj, 1, myThid )
c          CALL ZONAL_FILTER(
c    U                     vVel( 1-OLx,1-OLy,k,bi,bj),
c    I                     hFacS(1-OLx,1-OLy,k,bi,bj),
c    I                     0, sNy+1, 1, bi, bj, 2, myThid )
c        ENDIF
c#endif
        DO j=1,sNy
         DO i=1,sNx+1
          pf(i,j) = ( implicDiv2Dflow * gU(i,j,k,bi,bj)
     &     + (1. _d 0-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 (use3Dsolver) 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)*gV(i,j,k,bi,bj) / deltatmom
         ENDDO
        ENDDO
      ELSEIF (exactConserv) THEN
c     ELSEIF (nonlinFreeSurf.GT.0) THEN
C     Implicit treatment of the Barotropic Flow Divergence
        DO j=1,sNy+1
         DO i=1,sNx
          pf(i,j) = implicDiv2Dflow
     &             *yA(i,j)*gV(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 * gV(i,j,k,bi,bj)
     &     + (1. _d 0-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 (use3Dsolver) 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

#ifdef ALLOW_ADDFLUID
      IF ( selectAddFluid.GE.1 ) THEN
        DO j=1,sNy
         DO i=1,sNx
          cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
     &         - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTMom
         ENDDO
        ENDDO
#ifdef ALLOW_NONHYDROSTATIC
       IF (use3Dsolver) THEN
        DO j=1,sNy
         DO i=1,sNx
          cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
     &         - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTMom
         ENDDO
        ENDDO
       ENDIF
#endif
      ENDIF
#endif /* ALLOW_ADDFLUID */

      RETURN
      END
1	jmc	1.27	C $Header: /u/gcmpack/MITgcm/model/src/calc_div_ghat.F,v 1.26 2009/11/30 16:26:48 jmc Exp $
2	jmc	1.12	C $Name: $
3	cnh	1.1
4	cnh	1.6	#include "CPP_OPTIONS.h"
5	edhill	1.19
6	cnh	1.17	CBOP
7			C !ROUTINE: CALC_DIV_GHAT
8			C !INTERFACE:
9			SUBROUTINE CALC_DIV_GHAT(
10	jmc	1.22	I bi,bj,k,
11			U cg2d_b, cg3d_b,
12			I myThid)
13	cnh	1.17	C !DESCRIPTION: \bv
14			C ==========================================================
15	jmc	1.21	C \| S/R CALC_DIV_GHAT
16			C \| o Form the right hand-side of the surface pressure eqn.
17	cnh	1.17	C ==========================================================
18			C \| Right hand side of pressure equation is divergence
19			C \| of veclocity tendency (GHAT) term along with a relaxation
20			C \| term equal to the barotropic flow field divergence.
21			C ==========================================================
22			C \ev
23	cnh	1.1
24	cnh	1.17	C !USES:
25	cnh	1.1	IMPLICIT NONE
26			C == Global variables ==
27			#include "SIZE.h"
28			#include "EEPARAMS.h"
29			#include "PARAMS.h"
30			#include "GRID.h"
31	jmc	1.22	#include "DYNVARS.h"
32	jmc	1.27	#ifdef ALLOW_ADDFLUID
33			# include "FFIELDS.h"
34			#endif
35	cnh	1.1
36	cnh	1.17	C !INPUT/OUTPUT PARAMETERS:
37	cnh	1.1	C == Routine arguments ==
38	jmc	1.21	C bi, bj :: tile indices
39			C k :: Index of layer.
40			C cg2d_b :: Conjugate Gradient 2-D solver : Right-hand side vector
41	jmc	1.22	C cg3d_b :: Conjugate Gradient 3-D solver : Right-hand side vector
42	jmc	1.21	C myThid :: Instance number for this call of CALC_DIV_GHAT
43	jmc	1.22	INTEGER bi,bj
44	jmc	1.21	INTEGER k
45	jmc	1.13	_RL cg2d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
46	jmc	1.24	#ifdef ALLOW_NONHYDROSTATIC
47	jmc	1.22	_RL cg3d_b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
48	jmc	1.24	#else
49			_RL cg3d_b(1)
50			#endif
51	cnh	1.1	INTEGER myThid
52
53	cnh	1.17	C !LOCAL VARIABLES:
54	cnh	1.1	C == Local variables ==
55	jmc	1.22	C i,j :: Loop counters
56			C xA, yA :: Cell vertical face areas
57			C pf :: Intermediate array for building RHS source term.
58	cnh	1.1	INTEGER i,j
59	jmc	1.22	_RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60			_RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61	cnh	1.1	_RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62	cnh	1.17	CEOP
63	cnh	1.1
64	jmc	1.22	C Calculate vertical face areas
65			DO j=1,sNy+1
66			DO i=1,sNx+1
67			xA(i,j) = _dyG(i,j,bi,bj)*deepFacC(k)
68			& drF(k)_hFacW(i,j,k,bi,bj)*rhoFacC(k)
69			yA(i,j) = _dxG(i,j,bi,bj)*deepFacC(k)
70			& drF(k)_hFacS(i,j,k,bi,bj)*rhoFacC(k)
71			ENDDO
72			ENDDO
73
74	cnh	1.1	C-- Pressure equation source term
75	jmc	1.21	C Term is the vertical integral of the divergence of the
76	cnh	1.1	C time tendency terms along with a relaxation term that
77			C pulls div(U) + dh/dt back toward zero.
78
79	jmc	1.16	IF (implicDiv2Dflow.EQ.1.) THEN
80	jmc	1.12	C Fully Implicit treatment of the Barotropic Flow Divergence
81			DO j=1,sNy
82			DO i=1,sNx+1
83	jmc	1.27	pf(i,j) = xA(i,j)*gU(i,j,k,bi,bj) / deltaTMom
84	jmc	1.12	ENDDO
85			ENDDO
86	jmc	1.16	ELSEIF (exactConserv) THEN
87			c ELSEIF (nonlinFreeSurf.GT.0) THEN
88			C Implicit treatment of the Barotropic Flow Divergence
89			DO j=1,sNy
90			DO i=1,sNx+1
91	jmc	1.21	pf(i,j) = implicDiv2Dflow
92	jmc	1.27	& xA(i,j)gU(i,j,k,bi,bj) / deltaTMom
93	jmc	1.16	ENDDO
94			ENDDO
95	jmc	1.12	ELSE
96			C Explicit+Implicit part of the Barotropic Flow Divergence
97			C => Filtering of uVel,vVel is necessary
98	jmc	1.16	C-- Now the filter are applied in the_correction_step().
99			C We have left this code here to indicate where the filters used to be
100	jmc	1.21	C in the algorithm before JMC moved them to after the pressure solver.
101	jmc	1.25	c#ifdef ALLOW_ZONAL_FILT
102			c IF (zonal_filt_lat.LT.90.) THEN
103			c CALL ZONAL_FILTER(
104			c U uVel( 1-OLx,1-OLy,k,bi,bj),
105			c I hFacW(1-OLx,1-OLy,k,bi,bj),
106			c I 0, sNy+1, 1, bi, bj, 1, myThid )
107			c CALL ZONAL_FILTER(
108			c U vVel( 1-OLx,1-OLy,k,bi,bj),
109			c I hFacS(1-OLx,1-OLy,k,bi,bj),
110			c I 0, sNy+1, 1, bi, bj, 2, myThid )
111			c ENDIF
112			c#endif
113	jmc	1.12	DO j=1,sNy
114			DO i=1,sNx+1
115	jmc	1.18	pf(i,j) = ( implicDiv2Dflow * gU(i,j,k,bi,bj)
116	jmc	1.22	& + (1. _d 0-implicDiv2Dflow)* uVel(i,j,k,bi,bj)
117	jmc	1.27	& ) * xA(i,j) / deltaTMom
118	jmc	1.12	ENDDO
119			ENDDO
120			ENDIF
121	cnh	1.1	DO j=1,sNy
122			DO i=1,sNx
123	adcroft	1.8	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
124	cnh	1.1	& pf(i+1,j)-pf(i,j)
125			ENDDO
126			ENDDO
127
128	adcroft	1.8	#ifdef ALLOW_NONHYDROSTATIC
129	jmc	1.21	IF (use3Dsolver) THEN
130	adcroft	1.8	DO j=1,sNy
131			DO i=1,sNx
132	jmc	1.26	cg3d_b(i,j,k,bi,bj) = ( pf(i+1,j)-pf(i,j) )
133	adcroft	1.8	ENDDO
134			ENDDO
135			ENDIF
136			#endif
137
138	jmc	1.16	IF (implicDiv2Dflow.EQ.1.) THEN
139	jmc	1.12	C Fully Implicit treatment of the Barotropic Flow Divergence
140			DO j=1,sNy+1
141			DO i=1,sNx
142	jmc	1.18	pf(i,j) = yA(i,j)*gV(i,j,k,bi,bj) / deltatmom
143	jmc	1.16	ENDDO
144			ENDDO
145			ELSEIF (exactConserv) THEN
146			c ELSEIF (nonlinFreeSurf.GT.0) THEN
147			C Implicit treatment of the Barotropic Flow Divergence
148			DO j=1,sNy+1
149			DO i=1,sNx
150			pf(i,j) = implicDiv2Dflow
151	jmc	1.18	& yA(i,j)gV(i,j,k,bi,bj) / deltatmom
152	jmc	1.12	ENDDO
153			ENDDO
154			ELSE
155			C Explicit+Implicit part of the Barotropic Flow Divergence
156			DO j=1,sNy+1
157			DO i=1,sNx
158	jmc	1.18	pf(i,j) = ( implicDiv2Dflow * gV(i,j,k,bi,bj)
159	jmc	1.22	& + (1. _d 0-implicDiv2Dflow)* vVel(i,j,k,bi,bj)
160	jmc	1.27	& ) * yA(i,j) / deltaTMom
161	jmc	1.12	ENDDO
162			ENDDO
163			ENDIF
164	cnh	1.1	DO j=1,sNy
165			DO i=1,sNx
166	adcroft	1.8	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
167	cnh	1.1	& pf(i,j+1)-pf(i,j)
168			ENDDO
169			ENDDO
170	cnh	1.4
171	adcroft	1.8	#ifdef ALLOW_NONHYDROSTATIC
172	jmc	1.21	IF (use3Dsolver) THEN
173	adcroft	1.8	DO j=1,sNy
174			DO i=1,sNx
175	jmc	1.26	cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
176			& + ( pf(i,j+1)-pf(i,j) )
177	adcroft	1.8	ENDDO
178			ENDDO
179			ENDIF
180			#endif
181	cnh	1.1
182	jmc	1.23	#ifdef ALLOW_ADDFLUID
183			IF ( selectAddFluid.GE.1 ) THEN
184			DO j=1,sNy
185			DO i=1,sNx
186			cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj)
187	jmc	1.27	& - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTMom
188	jmc	1.23	ENDDO
189			ENDDO
190			#ifdef ALLOW_NONHYDROSTATIC
191			IF (use3Dsolver) THEN
192			DO j=1,sNy
193			DO i=1,sNx
194			cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj)
195	jmc	1.27	& - addMass(i,j,k,bi,bj)*mass2rUnit/deltaTMom
196	jmc	1.23	ENDDO
197			ENDDO
198			ENDIF
199			#endif
200			ENDIF
201			#endif /* ALLOW_ADDFLUID */
202
203	cnh	1.1	RETURN
204			END