model/src/calc_gw.F

C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.13 2004/04/06 00:31:54 jmc Exp $
C     !DESCRIPTION: \bv
C $Name:  $

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

CBOP
C     !ROUTINE: CALC_GW
C     !INTERFACE:
      SUBROUTINE CALC_GW(     
     I        myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R CALC_GW                                               
C     | o Calculate vert. velocity tendency terms ( NH, QH only ) 
C     *==========================================================*
C     | In NH and QH, the vertical momentum tendency must be
C     | calculated explicitly and included as a source term 
C     | for a 3d pressure eqn. Calculate that term here.
C     | This routine is not used in HYD calculations.
C     *==========================================================*
C     \ev 

C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "GW.h"
#include "CG3D.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid - Instance number for this innvocation of CALC_GW
      INTEGER myThid

#ifdef ALLOW_NONHYDROSTATIC

C     !LOCAL VARIABLES:
C     == Local variables ==
C     bi, bj,      :: Loop counters
C     iMin, iMax,
C     jMin, jMax
C     flx_NS       :: Temp. used for fVol meridional terms.
C     flx_EW       :: Temp. used for fVol zonal terms.
C     flx_Up       :: Temp. used for fVol vertical terms.
C     flx_Dn       :: Temp. used for fVol vertical terms.
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      _RL    flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL    flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL    flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL    flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     I,J,K - Loop counters
      INTEGER i,j,k, kP1, kUp
      _RL  wOverride
      _RS  hFacROpen
      _RS  hFacRClosed
      _RL ab15,ab05
      _RL slipSideFac
      _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ

      _RL  Half
      _RL  One
      PARAMETER(Half=0.5D0)
      PARAMETER(One=0.5D0)
CEOP

ceh3 needs an IF ( useNONHYDROSTATIC ) THEN

      iMin = 1
      iMax = sNx
      jMin = 1
      jMax = sNy

C     Adams-Bashforth timestepping weights
      ab15 =  1.5 _d 0 + abeps
      ab05 = -0.5 _d 0 - abeps

C     Lateral friction (no-slip, free slip, or half slip):
      IF ( no_slip_sides ) THEN
        slipSideFac = -One
      ELSE
        slipSideFac =  One 
      ENDIF
C-    half slip was used before ; keep it for now.
C        slipSideFac = 0. _d 0

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO K=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj) 
           gW(i,j,k,bi,bj) = 0.
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C Catch barotropic mode
      IF ( Nr .LT. 2 ) RETURN

C For each tile
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

C Boundaries condition at top
        DO J=jMin,jMax
         DO I=iMin,iMax
          Flx_Dn(I,J,bi,bj)=0.
         ENDDO
        ENDDO

C Sweep down column
        DO K=2,Nr
         Kp1=K+1
         wOverRide=1.
         if (K.EQ.Nr) then
          Kp1=Nr
          wOverRide=0.
         endif
C Flux on Southern face
         DO J=jMin,jMax+1
          DO I=iMin,iMax
           tmp_VbarZ=Half*(
     &          _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)
     &         +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
           Flx_NS(I,J,bi,bj)=
     &     tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))
     &    -viscAh*_recip_dyC(I,J,bi,bj)  
     &      *(1. _d 0 + slipSideFac*
     &         (maskS(I,J,K-1,bi,bj)*maskS(I,J,K,bi,bj)-2. _d 0))
     &        *(max(hFacS(I,J,K-1,bi,bj)-Half,0) 
     &         +min(hFacS(I,J,K,bi,bj),Half))
     &                   *(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))
          ENDDO
         ENDDO
C Flux on Western face
         DO J=jMin,jMax
          DO I=iMin,iMax+1
           tmp_UbarZ=Half*(
     &         _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)
     &        +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
           Flx_EW(I,J,bi,bj)=
     &     tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))
     &    -viscAh*_recip_dxC(I,J,bi,bj)
     &      *(1. _d 0 + slipSideFac*
     &         (maskW(I,J,K-1,bi,bj)*maskW(I,J,K,bi,bj)-1. _d 0))
     &        *(max(hFacW(I,J,K-1,bi,bj)-Half,0) 
     &         +min(hFacW(I,J,K,bi,bj),Half))
     &                   *(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))
          ENDDO
         ENDDO
C Flux on Lower face
         DO J=jMin,jMax
          DO I=iMin,iMax
           Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
           tmp_WbarZ=Half*(wVel(I,J,K,bi,bj) 
     &         +wOverRide*wVel(I,J,Kp1,bi,bj))
           Flx_Dn(I,J,bi,bj)=
     &     tmp_WbarZ*tmp_WbarZ
     &    -viscAr*recip_drF(K) 
     &       *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )
          ENDDO
         ENDDO
C        Divergence of fluxes
         DO J=jMin,jMax
          DO I=iMin,iMax
           gW(I,J,K,bi,bj) = 0.
     &      -(
     &        +_recip_dxF(I,J,bi,bj)*(
     &              Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )
     &        +_recip_dyF(I,J,bi,bj)*(
     &              Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )
     &        +recip_drC(K)         *(
     &              Flx_Up(I,J,bi,bj)  -Flx_Dn(I,J,bi,bj) )
     &       )
caja    * recip_hFacU(I,J,K,bi,bj)
caja   NOTE:  This should be included   
caja   but we need an hFacUW (above U points)
caja           and an hFacUS (above V points) too...
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

     
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO K=2,Nr
         DO j=jMin,jMax
          DO i=iMin,iMax
           wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj) 
     &     +deltatMom*( ab15*gW(i,j,k,bi,bj)
     &                 +ab05*gWNM1(i,j,k,bi,bj) )
           IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_OBCS
      IF (useOBCS) THEN
C--   This call is aesthetic: it makes the W field
C     consistent with the OBs but this has no algorithmic
C     impact. This is purely for diagnostic purposes.
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO K=1,Nr
          CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )
         ENDDO
        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_OBCS */

#endif /* ALLOW_NONHYDROSTATIC */

      RETURN
      END
1	mlosch	1.14	C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.13 2004/04/06 00:31:54 jmc Exp $
2	cnh	1.9	C !DESCRIPTION: \bv
3	jmc	1.7	C $Name: $
4	edhill	1.10
5			#include "PACKAGES_CONFIG.h"
6	adcroft	1.1	#include "CPP_OPTIONS.h"
7
8	cnh	1.9	CBOP
9			C !ROUTINE: CALC_GW
10			C !INTERFACE:
11	adcroft	1.1	SUBROUTINE CALC_GW(
12			I myThid)
13	cnh	1.9	C !DESCRIPTION: \bv
14			C ==========================================================
15			C \| S/R CALC_GW
16			C \| o Calculate vert. velocity tendency terms ( NH, QH only )
17			C ==========================================================
18			C \| In NH and QH, the vertical momentum tendency must be
19			C \| calculated explicitly and included as a source term
20			C \| for a 3d pressure eqn. Calculate that term here.
21			C \| This routine is not used in HYD calculations.
22			C ==========================================================
23			C \ev
24
25			C !USES:
26	adcroft	1.1	IMPLICIT NONE
27			C == Global variables ==
28			#include "SIZE.h"
29			#include "DYNVARS.h"
30			#include "EEPARAMS.h"
31			#include "PARAMS.h"
32			#include "GRID.h"
33			#include "GW.h"
34			#include "CG3D.h"
35
36	cnh	1.9	C !INPUT/OUTPUT PARAMETERS:
37	adcroft	1.1	C == Routine arguments ==
38			C myThid - Instance number for this innvocation of CALC_GW
39			INTEGER myThid
40
41	adcroft	1.3	#ifdef ALLOW_NONHYDROSTATIC
42
43	cnh	1.9	C !LOCAL VARIABLES:
44	adcroft	1.1	C == Local variables ==
45	cnh	1.9	C bi, bj, :: Loop counters
46			C iMin, iMax,
47			C jMin, jMax
48			C flx_NS :: Temp. used for fVol meridional terms.
49			C flx_EW :: Temp. used for fVol zonal terms.
50			C flx_Up :: Temp. used for fVol vertical terms.
51			C flx_Dn :: Temp. used for fVol vertical terms.
52	adcroft	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
53			_RL flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54			_RL flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55			_RL flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56			_RL flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57			C I,J,K - Loop counters
58	adcroft	1.4	INTEGER i,j,k, kP1, kUp
59	adcroft	1.1	_RL wOverride
60			_RS hFacROpen
61			_RS hFacRClosed
62			_RL ab15,ab05
63	jmc	1.12	_RL slipSideFac
64	adcroft	1.1	_RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
65
66			_RL Half
67	mlosch	1.14	_RL One
68	adcroft	1.1	PARAMETER(Half=0.5D0)
69	mlosch	1.14	PARAMETER(One=0.5D0)
70	cnh	1.9	CEOP
71	edhill	1.10
72			ceh3 needs an IF ( useNONHYDROSTATIC ) THEN
73	adcroft	1.1
74	mlosch	1.14	iMin = 1
75			iMax = sNx
76			jMin = 1
77			jMax = sNy
78
79	adcroft	1.1	C Adams-Bashforth timestepping weights
80	jmc	1.11	ab15 = 1.5 _d 0 + abeps
81			ab05 = -0.5 _d 0 - abeps
82
83			C Lateral friction (no-slip, free slip, or half slip):
84			IF ( no_slip_sides ) THEN
85	mlosch	1.14	slipSideFac = -One
86	jmc	1.11	ELSE
87	mlosch	1.14	slipSideFac = One
88	jmc	1.11	ENDIF
89			C- half slip was used before ; keep it for now.
90	mlosch	1.14	C slipSideFac = 0. _d 0
91	jmc	1.11
92	adcroft	1.1	DO bj=myByLo(myThid),myByHi(myThid)
93			DO bi=myBxLo(myThid),myBxHi(myThid)
94			DO K=1,Nr
95			DO j=1-OLy,sNy+OLy
96			DO i=1-OLx,sNx+OLx
97	jmc	1.8	gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
98	adcroft	1.1	gW(i,j,k,bi,bj) = 0.
99			ENDDO
100			ENDDO
101			ENDDO
102			ENDDO
103			ENDDO
104
105			C Catch barotropic mode
106			IF ( Nr .LT. 2 ) RETURN
107
108			C For each tile
109			DO bj=myByLo(myThid),myByHi(myThid)
110			DO bi=myBxLo(myThid),myBxHi(myThid)
111
112			C Boundaries condition at top
113	mlosch	1.14	DO J=jMin,jMax
114			DO I=iMin,iMax
115	adcroft	1.1	Flx_Dn(I,J,bi,bj)=0.
116			ENDDO
117			ENDDO
118
119			C Sweep down column
120			DO K=2,Nr
121			Kp1=K+1
122			wOverRide=1.
123			if (K.EQ.Nr) then
124			Kp1=Nr
125			wOverRide=0.
126			endif
127			C Flux on Southern face
128	mlosch	1.14	DO J=jMin,jMax+1
129			DO I=iMin,iMax
130	adcroft	1.1	tmp_VbarZ=Half*(
131			& _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)
132			& +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
133			Flx_NS(I,J,bi,bj)=
134			& tmp_VbarZHalf(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))
135	jmc	1.11	& -viscAh*_recip_dyC(I,J,bi,bj)
136			& (1. _d 0 + slipSideFac
137	mlosch	1.14	& (maskS(I,J,K-1,bi,bj)*maskS(I,J,K,bi,bj)-2. _d 0))
138			& *(max(hFacS(I,J,K-1,bi,bj)-Half,0)
139			& +min(hFacS(I,J,K,bi,bj),Half))
140	jmc	1.11	& *(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))
141	adcroft	1.1	ENDDO
142			ENDDO
143			C Flux on Western face
144	mlosch	1.14	DO J=jMin,jMax
145			DO I=iMin,iMax+1
146	adcroft	1.1	tmp_UbarZ=Half*(
147			& _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)
148			& +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
149			Flx_EW(I,J,bi,bj)=
150			& tmp_UbarZHalf(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))
151	jmc	1.11	& -viscAh*_recip_dxC(I,J,bi,bj)
152			& (1. _d 0 + slipSideFac
153	mlosch	1.14	& (maskW(I,J,K-1,bi,bj)*maskW(I,J,K,bi,bj)-1. _d 0))
154			& *(max(hFacW(I,J,K-1,bi,bj)-Half,0)
155			& +min(hFacW(I,J,K,bi,bj),Half))
156	jmc	1.11	& *(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))
157	adcroft	1.1	ENDDO
158			ENDDO
159			C Flux on Lower face
160	mlosch	1.14	DO J=jMin,jMax
161			DO I=iMin,iMax
162	adcroft	1.1	Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
163	mlosch	1.14	tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)
164			& +wOverRide*wVel(I,J,Kp1,bi,bj))
165	adcroft	1.1	Flx_Dn(I,J,bi,bj)=
166			& tmp_WbarZ*tmp_WbarZ
167	jmc	1.11	& -viscAr*recip_drF(K)
168			& ( wVel(I,J,K,bi,bj)-wOverRidewVel(I,J,Kp1,bi,bj) )
169	adcroft	1.1	ENDDO
170			ENDDO
171			C Divergence of fluxes
172	mlosch	1.14	DO J=jMin,jMax
173			DO I=iMin,iMax
174	adcroft	1.1	gW(I,J,K,bi,bj) = 0.
175			& -(
176			& +_recip_dxF(I,J,bi,bj)*(
177			& Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )
178			& +_recip_dyF(I,J,bi,bj)*(
179			& Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )
180			& +recip_drC(K) *(
181			& Flx_Up(I,J,bi,bj) -Flx_Dn(I,J,bi,bj) )
182			& )
183			caja * recip_hFacU(I,J,K,bi,bj)
184			caja NOTE: This should be included
185			caja but we need an hFacUW (above U points)
186			caja and an hFacUS (above V points) too...
187			ENDDO
188			ENDDO
189			ENDDO
190			ENDDO
191			ENDDO
192
193
194			DO bj=myByLo(myThid),myByHi(myThid)
195			DO bi=myBxLo(myThid),myBxHi(myThid)
196			DO K=2,Nr
197	mlosch	1.14	DO j=jMin,jMax
198			DO i=iMin,iMax
199	adcroft	1.1	wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
200			& +deltatMom( ab15gW(i,j,k,bi,bj)
201			& +ab05*gWNM1(i,j,k,bi,bj) )
202			IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.
203			ENDDO
204			ENDDO
205	adcroft	1.4	ENDDO
206			ENDDO
207			ENDDO
208
209	adcroft	1.5	#ifdef ALLOW_OBCS
210			IF (useOBCS) THEN
211	adcroft	1.4	C-- This call is aesthetic: it makes the W field
212			C consistent with the OBs but this has no algorithmic
213			C impact. This is purely for diagnostic purposes.
214	adcroft	1.5	DO bj=myByLo(myThid),myByHi(myThid)
215			DO bi=myBxLo(myThid),myBxHi(myThid)
216			DO K=1,Nr
217			CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )
218			ENDDO
219	adcroft	1.1	ENDDO
220			ENDDO
221	adcroft	1.5	ENDIF
222			#endif /* ALLOW_OBCS */
223	adcroft	1.1
224			#endif /* ALLOW_NONHYDROSTATIC */
225
226			RETURN
227			END