/[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.2 - (hide annotations) (download)
Tue Mar 14 17:47:25 2000 UTC (24 years, 2 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint25
Changes since 1.1: +0 -41 lines 
Various updates for OBCs and Non-hydrostatic routines.
 o OBCs now fits into time-stepping properly
 o div.G has been moved to solve_for_pressure()
1 adcroft 1.1 #include "CPP_OPTIONS.h"
2    
3     SUBROUTINE CALC_GW(
4     I myThid)
5     C /==========================================================\
6     C | S/R CALC_GW |
7     C \==========================================================/
8     IMPLICIT NONE
9    
10     C == Global variables ==
11     #include "SIZE.h"
12     #include "DYNVARS.h"
13     #include "FFIELDS.h"
14     #include "EEPARAMS.h"
15     #include "PARAMS.h"
16     #include "GRID.h"
17     #include "CG2D.h"
18     #include "GW.h"
19     #include "CG3D.h"
20    
21     C == Routine arguments ==
22     C myThid - Instance number for this innvocation of CALC_GW
23     INTEGER myThid
24    
25     C == Local variables ==
26     INTEGER bi,bj,iMin,iMax,jMin,jMax
27     _RL aF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
28     _RL vF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
29     _RL v4F(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30     _RL cF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31     _RL mT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32     _RL pF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33     _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
34     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
35    
36     _RL flx_NS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
37     _RL flx_EW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
38     _RL flx_Dn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
39     _RL flx_Up(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
40     C I,J,K - Loop counters
41     INTEGER i,j,k, kP1
42     _RL wOverride
43     _RS hFacROpen
44     _RS hFacRClosed
45     _RL ab15,ab05
46     _RL tmp_VbarZ, tmp_UbarZ, tmp_WbarZ
47    
48     _RL Half
49     PARAMETER(Half=0.5D0)
50    
51     #ifdef ALLOW_NONHYDROSTATIC
52    
53     #define I0 1
54     #define In sNx
55     #define J0 1
56     #define Jn sNy
57    
58     C Adams-Bashforth timestepping weights
59     ab15=1.5+abeps
60     ab05=-0.5-abeps
61    
62     DO bj=myByLo(myThid),myByHi(myThid)
63     DO bi=myBxLo(myThid),myBxHi(myThid)
64     DO K=1,Nr
65     DO j=1-OLy,sNy+OLy
66     DO i=1-OLx,sNx+OLx
67     gW(i,j,k,bi,bj) = 0.
68     ENDDO
69     ENDDO
70     ENDDO
71     ENDDO
72     ENDDO
73    
74     C Catch barotropic mode
75     IF ( Nr .LT. 2 ) RETURN
76    
77     C For each tile
78     DO bj=myByLo(myThid),myByHi(myThid)
79     DO bi=myBxLo(myThid),myBxHi(myThid)
80    
81     C Boundaries condition at top
82     DO J=J0,Jn
83     DO I=I0,In
84     Flx_Dn(I,J,bi,bj)=0.
85     ENDDO
86     ENDDO
87    
88     C Sweep down column
89     DO K=2,Nr
90     Kp1=K+1
91     wOverRide=1.
92     if (K.EQ.Nr) then
93     Kp1=Nr
94     wOverRide=0.
95     endif
96     C Flux on Southern face
97     DO J=J0,Jn+1
98     DO I=I0,In
99     tmp_VbarZ=Half*(
100     & _hFacS(I,J,K-1,bi,bj)*vVel( I ,J,K-1,bi,bj)
101     & +_hFacS(I,J, K ,bi,bj)*vVel( I ,J, K ,bi,bj))
102     Flx_NS(I,J,bi,bj)=
103     & tmp_VbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I,J-1,K,bi,bj))
104     & -viscAh*_recip_dyC(I,J,bi,bj)*(
105     & wVel(I,J,K,bi,bj)-wVel(I,J-1,K,bi,bj) )
106     ENDDO
107     ENDDO
108     C Flux on Western face
109     DO J=J0,Jn
110     DO I=I0,In+1
111     tmp_UbarZ=Half*(
112     & _hFacW(I,J,K-1,bi,bj)*uVel( I ,J,K-1,bi,bj)
113     & +_hFacW(I,J, K ,bi,bj)*uVel( I ,J, K ,bi,bj))
114     Flx_EW(I,J,bi,bj)=
115     & tmp_UbarZ*Half*(wVel(I,J,K,bi,bj)+wVel(I-1,J,K,bi,bj))
116     & -viscAh*_recip_dxC(I,J,bi,bj)*(
117     & wVel(I,J,K,bi,bj)-wVel(I-1,J,K,bi,bj) )
118     ENDDO
119     ENDDO
120     C Flux on Lower face
121     DO J=J0,Jn
122     DO I=I0,In
123     Flx_Up(I,J,bi,bj)=Flx_Dn(I,J,bi,bj)
124     tmp_WbarZ=Half*(wVel(I,J,K,bi,bj)+wVel(I,J,Kp1,bi,bj))
125     Flx_Dn(I,J,bi,bj)=
126     & tmp_WbarZ*tmp_WbarZ
127     & -viscAr*recip_drF(K)*( wVel(I,J,K,bi,bj)
128     & -wOverRide*wVel(I,J,Kp1,bi,bj) )
129     ENDDO
130     ENDDO
131     C Divergence of fluxes
132     DO J=J0,Jn
133     DO I=I0,In
134     gW(I,J,K,bi,bj) = 0.
135     & -(
136     & +_recip_dxF(I,J,bi,bj)*(
137     & Flx_EW(I+1,J,bi,bj)-Flx_EW(I,J,bi,bj) )
138     & +_recip_dyF(I,J,bi,bj)*(
139     & Flx_NS(I,J+1,bi,bj)-Flx_NS(I,J,bi,bj) )
140     & +recip_drC(K) *(
141     & Flx_Up(I,J,bi,bj) -Flx_Dn(I,J,bi,bj) )
142     & )
143     caja * recip_hFacU(I,J,K,bi,bj)
144     caja NOTE: This should be included
145     caja but we need an hFacUW (above U points)
146     caja and an hFacUS (above V points) too...
147     ENDDO
148     ENDDO
149     ENDDO
150     ENDDO
151     ENDDO
152    
153    
154     DO bj=myByLo(myThid),myByHi(myThid)
155     DO bi=myBxLo(myThid),myBxHi(myThid)
156     DO K=2,Nr
157     DO j=J0,Jn
158     DO i=I0,In
159     wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
160     & +deltatMom*( ab15*gW(i,j,k,bi,bj)
161     & +ab05*gWNM1(i,j,k,bi,bj) )
162     IF (hFacC(I,J,K,bi,bj).EQ.0.) wVel(i,j,k,bi,bj)=0.
163     ENDDO
164     ENDDO
165     ENDDO
166     ENDDO
167     ENDDO
168     DO bj=myByLo(myThid),myByHi(myThid)
169     DO bi=myBxLo(myThid),myBxHi(myThid)
170     DO K=1,Nr
171     DO j=J0,Jn
172     DO i=I0,In
173     gWNM1(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
174     ENDDO
175     ENDDO
176     ENDDO
177     ENDDO
178     ENDDO
179    
180     #endif /* ALLOW_NONHYDROSTATIC */
181    
182     RETURN
183     END
184    
  ViewVC Help
Powered by ViewVC 1.1.22