model/src/calc_div_ghat.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_div_ghat.F,v 1.11 2001/02/04 14:38:45 cnh 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,
     I        K,
     I        xA,yA,
     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"
#include "CG2D.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     myThid - Instance number for this innvocation of CALC_MOM_RHS
      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)
      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 ( k .EQ. Nr ) THEN
C      Initialise source term on first pass
       DO j=1,sNy
        DO i=1,sNx
C Note: The source term containing cg2d_x and cg3d_x (at k=1)
C       has been moved to solve_for_pressure.F for convenience.
         cg2d_b(i,j,bi,bj) = 0.
#ifdef USE_NATURAL_BCS
     &     + freeSurfFac*_rA(i,j,bi,bj)*horiVertRatio*
     &       EmPmR(I,J,bi,bj)/deltaTMom
#endif
        ENDDO
       ENDDO
      ENDIF

      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	C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_div_ghat.F,v 1.11 2001/02/04 14:38:45 cnh Exp $
2	C $Name: $
3
4	#include "CPP_OPTIONS.h"
5
6	C /==========================================================\
7	C \| S/R CALC_DIV_GHAT \|
8	C \| o Form the right hand-side of the surface pressure eqn. \|
9	C \|==========================================================\|
10	C \==========================================================/
11	SUBROUTINE CALC_DIV_GHAT(
12	I bi,bj,iMin,iMax,jMin,jMax,
13	I K,
14	I xA,yA,
15	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	#include "CG2D.h"
27	#ifdef ALLOW_NONHYDROSTATIC
28	#include "CG3D.h"
29	#endif
30
31	C == Routine arguments ==
32	C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
33	C results will be set.
34	C k - Index of layer.
35	C xA, yA - Cell face areas
36	C myThid - Instance number for this innvocation of CALC_MOM_RHS
37	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	INTEGER myThid
42
43	C == Local variables ==
44	INTEGER i,j
45	_RL pf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46
47	C-- Pressure equation source term
48	C Term is the vertical integral of the divergence of the
49	C time tendency terms along with a relaxation term that
50	C pulls div(U) + dh/dt back toward zero.
51
52	IF ( k .EQ. Nr ) THEN
53	C Initialise source term on first pass
54	DO j=1,sNy
55	DO i=1,sNx
56	C Note: The source term containing cg2d_x and cg3d_x (at k=1)
57	C has been moved to solve_for_pressure.F for convenience.
58	cg2d_b(i,j,bi,bj) = 0.
59	#ifdef USE_NATURAL_BCS
60	& + freeSurfFac_rA(i,j,bi,bj)horiVertRatio*
61	& EmPmR(I,J,bi,bj)/deltaTMom
62	#endif
63	ENDDO
64	ENDDO
65	ENDIF
66
67	IF (implicDiv2Dflow.EQ.1.) then
68	C Fully Implicit treatment of the Barotropic Flow Divergence
69	DO j=1,sNy
70	DO i=1,sNx+1
71	pf(i,j) = xA(i,j)*gUNm1(i,j,k,bi,bj) / deltaTmom
72	ENDDO
73	ENDDO
74	ELSE
75	C Explicit+Implicit part of the Barotropic Flow Divergence
76	C => Filtering of uVel,vVel is necessary
77	#ifdef ALLOW_ZONAL_FILT
78	IF (zonal_filt_lat.LT.90.) THEN
79	CALL ZONAL_FILTER(
80	& uVel, hFacW, 1-1, sNy+1, k, k, bi, bj, 1, myThid)
81	CALL ZONAL_FILTER(
82	& vVel, hFacS, 1-1, sNy+1, k, k, bi, bj, 2, myThid)
83	ENDIF
84	#endif
85	DO j=1,sNy
86	DO i=1,sNx+1
87	pf(i,j) = ( implicDiv2Dflow * gUNm1(i,j,k,bi,bj)
88	& + (1.-implicDiv2Dflow) * uVel(i,j,k,bi,bj)
89	& ) * xA(i,j) / deltaTmom
90	ENDDO
91	ENDDO
92	ENDIF
93	DO j=1,sNy
94	DO i=1,sNx
95	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
96	& pf(i+1,j)-pf(i,j)
97	ENDDO
98	ENDDO
99
100	#ifdef ALLOW_NONHYDROSTATIC
101	IF (nonHydrostatic) THEN
102	DO j=1,sNy
103	DO i=1,sNx
104	cg3d_b(i,j,k,bi,bj) = pf(i+1,j)-pf(i,j)
105	ENDDO
106	ENDDO
107	ENDIF
108	#endif
109
110	IF (implicDiv2Dflow.EQ.1.) then
111	C Fully Implicit treatment of the Barotropic Flow Divergence
112	DO j=1,sNy+1
113	DO i=1,sNx
114	pf(i,j) = yA(i,j)*gVNm1(i,j,k,bi,bj) / deltatmom
115	ENDDO
116	ENDDO
117	ELSE
118	C Explicit+Implicit part of the Barotropic Flow Divergence
119	DO j=1,sNy+1
120	DO i=1,sNx
121	pf(i,j) = ( implicDiv2Dflow * gVNm1(i,j,k,bi,bj)
122	& + (1.-implicDiv2Dflow) * vVel(i,j,k,bi,bj)
123	& ) * yA(i,j) / deltaTmom
124	ENDDO
125	ENDDO
126	ENDIF
127	DO j=1,sNy
128	DO i=1,sNx
129	cg2d_b(i,j,bi,bj) = cg2d_b(i,j,bi,bj) +
130	& pf(i,j+1)-pf(i,j)
131	ENDDO
132	ENDDO
133
134	#ifdef ALLOW_NONHYDROSTATIC
135	IF (nonHydrostatic) THEN
136	DO j=1,sNy
137	DO i=1,sNx
138	cg3d_b(i,j,k,bi,bj) = cg3d_b(i,j,k,bi,bj) +
139	& pf(i,j+1)-pf(i,j)
140	ENDDO
141	ENDDO
142	ENDIF
143	#endif
144
145	RETURN
146	END