/[MITgcm]/MITgcm/model/src/calc_gw.F
ViewVC logotype

Annotation of /MITgcm/model/src/calc_gw.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.13 - (hide annotations) (download)
Tue Apr 6 00:31:54 2004 UTC (20 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint53, checkpoint52m_post, checkpoint53a_post, checkpoint52n_post, checkpoint53b_pre, checkpoint53b_post
Changes since 1.12: +1 -2 lines
include FFIELDS.h not needed ;

1 jmc 1.13 C $Header: /u/gcmpack/MITgcm/model/src/calc_gw.F,v 1.12 2004/04/05 21:46:18 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     PARAMETER(Half=0.5D0)
68    
69     #define I0 1
70     #define In sNx
71     #define J0 1
72     #define Jn sNy
73 cnh 1.9 CEOP
74 edhill 1.10
75     ceh3 needs an IF ( useNONHYDROSTATIC ) THEN
76 adcroft 1.1
77     C Adams-Bashforth timestepping weights
78 jmc 1.11 ab15 = 1.5 _d 0 + abeps
79     ab05 = -0.5 _d 0 - abeps
80    
81     C Lateral friction (no-slip, free slip, or half slip):
82     IF ( no_slip_sides ) THEN
83     slipSideFac = -Half
84     ELSE
85     slipSideFac = Half
86     ENDIF
87     C- half slip was used before ; keep it for now.
88     slipSideFac = 0. _d 0
89    
90 adcroft 1.1 DO bj=myByLo(myThid),myByHi(myThid)
91     DO bi=myBxLo(myThid),myBxHi(myThid)
92     DO K=1,Nr
93     DO j=1-OLy,sNy+OLy
94     DO i=1-OLx,sNx+OLx
95 jmc 1.8 gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
96 adcroft 1.1 gW(i,j,k,bi,bj) = 0.
97     ENDDO
98     ENDDO
99     ENDDO
100     ENDDO
101     ENDDO
102    
103     C Catch barotropic mode
104     IF ( Nr .LT. 2 ) RETURN
105    
106     C For each tile
107     DO bj=myByLo(myThid),myByHi(myThid)
108     DO bi=myBxLo(myThid),myBxHi(myThid)
109    
110     C Boundaries condition at top
111     DO J=J0,Jn
112     DO I=I0,In
113     Flx_Dn(I,J,bi,bj)=0.
114     ENDDO
115     ENDDO
116    
117     C Sweep down column
118     DO K=2,Nr
119     Kp1=K+1
120     wOverRide=1.
121     if (K.EQ.Nr) then
122     Kp1=Nr
123     wOverRide=0.
124     endif
125     C Flux on Southern face
126     DO J=J0,Jn+1
127     DO I=I0,In
128     tmp_VbarZ=Half*(
129     & _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)
130     & +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
131     Flx_NS(I,J,bi,bj)=
132     & tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))
133 jmc 1.11 & -viscAh*_recip_dyC(I,J,bi,bj)
134     & *(1. _d 0 + slipSideFac*
135     & (maskS(I,J,K-1,bi,bj)+maskS(I,J,K,bi,bj)-2. _d 0))
136     & *(wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj))
137 adcroft 1.1 ENDDO
138     ENDDO
139     C Flux on Western face
140     DO J=J0,Jn
141     DO I=I0,In+1
142     tmp_UbarZ=Half*(
143     & _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)
144     & +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
145     Flx_EW(I,J,bi,bj)=
146     & tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))
147 jmc 1.11 & -viscAh*_recip_dxC(I,J,bi,bj)
148     & *(1. _d 0 + slipSideFac*
149     & (maskW(I,J,K-1,bi,bj)+maskW(I,J,K,bi,bj)-2. _d 0))
150     & *(wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj))
151 adcroft 1.1 ENDDO
152     ENDDO
153     C Flux on Lower face
154     DO J=J0,Jn
155     DO I=I0,In
156     Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
157     tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+wVel(I,J,Kp1,bi,bj))
158     Flx_Dn(I,J,bi,bj)=
159     & tmp_WbarZ*tmp_WbarZ
160 jmc 1.11 & -viscAr*recip_drF(K)
161     & *( wVel(I,J,K,bi,bj)-wOverRide*wVel(I,J,Kp1,bi,bj) )
162 adcroft 1.1 ENDDO
163     ENDDO
164     C Divergence of fluxes
165     DO J=J0,Jn
166     DO I=I0,In
167     gW(I,J,K,bi,bj) = 0.
168     & -(
169     & +_recip_dxF(I,J,bi,bj)*(
170     & Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )
171     & +_recip_dyF(I,J,bi,bj)*(
172     & Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )
173     & +recip_drC(K) *(
174     & Flx_Up(I,J,bi,bj) -Flx_Dn(I,J,bi,bj) )
175     & )
176     caja * recip_hFacU(I,J,K,bi,bj)
177     caja NOTE: This should be included
178     caja but we need an hFacUW (above U points)
179     caja and an hFacUS (above V points) too...
180     ENDDO
181     ENDDO
182     ENDDO
183     ENDDO
184     ENDDO
185    
186    
187     DO bj=myByLo(myThid),myByHi(myThid)
188     DO bi=myBxLo(myThid),myBxHi(myThid)
189     DO K=2,Nr
190     DO j=J0,Jn
191     DO i=I0,In
192     wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
193     & +deltatMom*( ab15*gW(i,j,k,bi,bj)
194     & +ab05*gWNM1(i,j,k,bi,bj) )
195     IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.
196     ENDDO
197     ENDDO
198 adcroft 1.4 ENDDO
199     ENDDO
200     ENDDO
201    
202 adcroft 1.5 #ifdef ALLOW_OBCS
203     IF (useOBCS) THEN
204 adcroft 1.4 C-- This call is aesthetic: it makes the W field
205     C consistent with the OBs but this has no algorithmic
206     C impact. This is purely for diagnostic purposes.
207 adcroft 1.5 DO bj=myByLo(myThid),myByHi(myThid)
208     DO bi=myBxLo(myThid),myBxHi(myThid)
209     DO K=1,Nr
210     CALL OBCS_APPLY_W( bi, bj, K, wVel, myThid )
211     ENDDO
212 adcroft 1.1 ENDDO
213     ENDDO
214 adcroft 1.5 ENDIF
215     #endif /* ALLOW_OBCS */
216 adcroft 1.1
217     #endif /* ALLOW_NONHYDROSTATIC */
218    
219     RETURN
220     END

  ViewVC Help
Powered by ViewVC 1.1.22