model/src/calc_div_ghat.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_div_ghat.F,v 1.23 2008/08/24 21:38:18 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"

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