model/src/calc_div_ghat.F

C $Header: /u/u3/gcmpack/MITgcm/model/src/calc_div_ghat.F,v 1.18 2003/04/17 13:40:06 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CALC_DIV_GHAT
C     !INTERFACE:
      SUBROUTINE CALC_DIV_GHAT(
     I        bi,bj,iMin,iMax,jMin,jMax,K,
     I        xA,yA,
     U        cg2d_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 "DYNVARS.h"
#include "FFIELDS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SOLVE_FOR_PRESSURE3D.h"

C     !INPUT/OUTPUT PARAMETERS:
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:
C     == Local variables ==
C     i,j :: Loop counters
C     pf  :: Intermediate array for building RHS source term.
      INTEGER i,j
      _RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

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