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

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

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


Revision 1.28 - (hide annotations) (download)
Tue May 29 14:01:36 2001 UTC (23 years ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint40pre3, checkpoint40pre1, checkpoint40pre7, checkpoint40pre6, checkpoint40pre2, checkpoint40pre4, checkpoint40pre5
Changes since 1.27: +69 -42 lines
Merge from branch pre38:
 o essential mods for cubed sphere
 o debugged atmosphere, dynamcis + physics (aim)
 o new packages (mom_vecinv, mom_fluxform, ...)

1 adcroft 1.28 C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gs.F,v 1.27.2.4 2001/04/09 18:52:45 adcroft Exp $
2 cnh 1.25 C $Name: $
3 cnh 1.1
4 cnh 1.17 #include "CPP_OPTIONS.h"
5 cnh 1.1
6 adcroft 1.28 #define COSINEMETH_III
7     #undef ISOTROPIC_COS_SCALING
8    
9 cnh 1.1 CStartOfInterFace
10     SUBROUTINE CALC_GS(
11     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12 adcroft 1.28 I xA,yA,uTrans,vTrans,rTrans,maskUp,
13 adcroft 1.21 I KappaRS,
14 adcroft 1.23 U fVerS,
15 cnh 1.16 I myCurrentTime, myThid )
16 cnh 1.1 C /==========================================================\
17     C | SUBROUTINE CALC_GS |
18 adcroft 1.7 C | o Calculate the salt tendency terms. |
19 cnh 1.1 C |==========================================================|
20     C | A procedure called EXTERNAL_FORCING_S is called from |
21     C | here. These procedures can be used to add per problem |
22 adcroft 1.7 C | E-P flux source terms. |
23 cnh 1.1 C | Note: Although it is slightly counter-intuitive the |
24     C | EXTERNAL_FORCING routine is not the place to put |
25     C | file I/O. Instead files that are required to |
26     C | calculate the external source terms are generally |
27     C | read during the model main loop. This makes the |
28     C | logisitics of multi-processing simpler and also |
29     C | makes the adjoint generation simpler. It also |
30     C | allows for I/O to overlap computation where that |
31     C | is supported by hardware. |
32     C | Aside from the problem specific term the code here |
33     C | forms the tendency terms due to advection and mixing |
34     C | The baseline implementation here uses a centered |
35     C | difference form for the advection term and a tensorial |
36     C | divergence of a flux form for the diffusive term. The |
37     C | diffusive term is formulated so that isopycnal mixing and|
38     C | GM-style subgrid-scale terms can be incorporated b simply|
39     C | setting the diffusion tensor terms appropriately. |
40     C \==========================================================/
41     IMPLICIT NONE
42    
43     C == GLobal variables ==
44     #include "SIZE.h"
45     #include "DYNVARS.h"
46     #include "EEPARAMS.h"
47     #include "PARAMS.h"
48     #include "GRID.h"
49 cnh 1.8 #include "FFIELDS.h"
50 cnh 1.1
51     C == Routine arguments ==
52 adcroft 1.7 C fVerS - Flux of salt (S) in the vertical
53 cnh 1.1 C direction at the upper(U) and lower(D) faces of a cell.
54     C maskUp - Land mask used to denote base of the domain.
55     C xA - Tracer cell face area normal to X
56     C yA - Tracer cell face area normal to X
57     C uTrans - Zonal volume transport through cell face
58     C vTrans - Meridional volume transport through cell face
59 adcroft 1.18 C rTrans - Vertical volume transport through cell face
60 cnh 1.1 C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
61     C results will be set.
62 adcroft 1.7 C myThid - Instance number for this innvocation of CALC_GT
63 cnh 1.1 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
64     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66     _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 cnh 1.12 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69 cnh 1.1 _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 cnh 1.14 _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
71 adcroft 1.7 INTEGER k,kUp,kDown,kM1
72 cnh 1.1 INTEGER bi,bj,iMin,iMax,jMin,jMax
73 adcroft 1.18 _RL myCurrentTime
74 cnh 1.1 INTEGER myThid
75     CEndOfInterface
76    
77     C == Local variables ==
78     C I, J, K - Loop counters
79 adcroft 1.7 INTEGER i,j
80     LOGICAL TOP_LAYER
81     _RL afFacS, dfFacS
82 adcroft 1.19 _RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83 adcroft 1.23 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85     _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86     _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87 heimbach 1.20
88     #ifdef ALLOW_AUTODIFF_TAMC
89     C-- only the kUp part of fverS is set in this subroutine
90     C-- the kDown is still required
91     fVerS(1,1,kDown) = fVerS(1,1,kDown)
92 adcroft 1.27 #endif
93 heimbach 1.20 DO j=1-OLy,sNy+OLy
94     DO i=1-OLx,sNx+OLx
95     fZon(i,j) = 0.0
96     fMer(i,j) = 0.0
97     fVerS(i,j,kUp) = 0.0
98     ENDDO
99     ENDDO
100 cnh 1.1
101     afFacS = 1. _d 0
102     dfFacS = 1. _d 0
103 adcroft 1.7 TOP_LAYER = K .EQ. 1
104 cnh 1.1
105     C--- Calculate advective and diffusive fluxes between cells.
106    
107 adcroft 1.19 C o Zonal tracer gradient
108     DO j=1-Oly,sNy+Oly
109     DO i=1-Olx+1,sNx+Olx
110 adcroft 1.28 fZon(i,j) = _recip_dxC(i,j,bi,bj)*xA(i,j)
111     & *(salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
112     #ifdef COSINEMETH_III
113     & *sqCosFacU(j,bi,bj)
114     #endif
115 adcroft 1.19 ENDDO
116     ENDDO
117     C o Meridional tracer gradient
118     DO j=1-Oly+1,sNy+Oly
119     DO i=1-Olx,sNx+Olx
120 adcroft 1.28 fMer(i,j) = _recip_dyC(i,j,bi,bj)*yA(i,j)
121     & *(salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
122     #ifdef ISOTROPIC_COS_SCALING
123     #ifdef COSINEMETH_III
124     & *sqCosFacV(j,bi,bj)
125     #endif
126     #endif
127 adcroft 1.19 ENDDO
128     ENDDO
129    
130     C-- del^2 of S, needed for bi-harmonic (del^4) term
131     IF (diffK4S .NE. 0.) THEN
132     DO j=1-Oly+1,sNy+Oly-1
133     DO i=1-Olx+1,sNx+Olx-1
134     df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
135     & *recip_drF(k)/_rA(i,j,bi,bj)
136     & *(
137 adcroft 1.28 & +( fZon(i+1,j)-fZon(i,j) )
138     & +( fMer(i,j+1)-fMer(i,j) )
139 adcroft 1.19 & )
140     ENDDO
141     ENDDO
142     ENDIF
143    
144 cnh 1.1 C-- Zonal flux (fZon is at west face of "salt" cell)
145     C Advective component of zonal flux
146     DO j=jMin,jMax
147     DO i=iMin,iMax
148     af(i,j) =
149     & uTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj))*0.5 _d 0
150     ENDDO
151     ENDDO
152 adcroft 1.19 C o Diffusive component of zonal flux
153 cnh 1.1 DO j=jMin,jMax
154     DO i=iMin,iMax
155 adcroft 1.28 df(i,j) = -diffKhS*xA(i,j)*_recip_dxC(i,j,bi,bj)*
156     & (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
157     & *CosFacU(j,bi,bj)
158 cnh 1.1 ENDDO
159     ENDDO
160 adcroft 1.21 #ifdef ALLOW_GMREDI
161 heimbach 1.22 IF (useGMRedi) CALL GMREDI_XTRANSPORT(
162 adcroft 1.21 I iMin,iMax,jMin,jMax,bi,bj,K,
163     I xA,salt,
164     U df,
165     I myThid)
166     #endif
167 adcroft 1.19 C o Add the bi-harmonic contribution
168     IF (diffK4S .NE. 0.) THEN
169     DO j=jMin,jMax
170     DO i=iMin,iMax
171     df(i,j) = df(i,j) + xA(i,j)*
172     & diffK4S*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)
173 adcroft 1.28 #ifdef COSINEMETH_III
174     & *sqCosFacU(j,bi,bj)
175     #else
176     & *CosFacU(j,bi,bj)
177     #endif
178 adcroft 1.19 ENDDO
179     ENDDO
180     ENDIF
181 cnh 1.1 C Net zonal flux
182     DO j=jMin,jMax
183     DO i=iMin,iMax
184     fZon(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
185     ENDDO
186     ENDDO
187    
188     C-- Meridional flux (fMer is at south face of "salt" cell)
189     C Advective component of meridional flux
190     DO j=jMin,jMax
191     DO i=iMin,iMax
192     C Advective component of meridional flux
193     af(i,j) =
194     & vTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj))*0.5 _d 0
195     ENDDO
196     ENDDO
197     C Diffusive component of meridional flux
198     DO j=jMin,jMax
199     DO i=iMin,iMax
200 adcroft 1.28 df(i,j) = -diffKhS*yA(i,j)*_recip_dyC(i,j,bi,bj)*
201     & (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
202     & *CosFacV(j,bi,bj)
203 cnh 1.1 ENDDO
204     ENDDO
205 adcroft 1.21 #ifdef ALLOW_GMREDI
206 heimbach 1.22 IF (useGMRedi) CALL GMREDI_YTRANSPORT(
207 adcroft 1.21 I iMin,iMax,jMin,jMax,bi,bj,K,
208     I yA,salt,
209     U df,
210     I myThid)
211     #endif
212 adcroft 1.19 C o Add the bi-harmonic contribution
213     IF (diffK4S .NE. 0.) THEN
214     DO j=jMin,jMax
215     DO i=iMin,iMax
216     df(i,j) = df(i,j) + yA(i,j)*
217     & diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
218 adcroft 1.28 #ifdef ISOTROPIC_COS_SCALING
219     #ifdef COSINEMETH_III
220     & *sqCosFacV(j,bi,bj)
221     #else
222     & *CosFacV(j,bi,bj)
223     #endif
224     #endif
225 adcroft 1.19 ENDDO
226     ENDDO
227     ENDIF
228    
229 cnh 1.1 C Net meridional flux
230     DO j=jMin,jMax
231     DO i=iMin,iMax
232     fMer(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
233     ENDDO
234     ENDDO
235    
236 adcroft 1.28 C-- Vertical flux ( fVerS(,,kUp) is at upper face of "Tracer" cell )
237     C o Advective component of vertical flux : assume W_bottom=0 (mask)
238     C Note: For K=1 then KM1=1 this gives a barZ(S) = S
239 adcroft 1.7 C (this plays the role of the free-surface correction)
240 adcroft 1.28 IF ( rigidLid .AND. TOP_LAYER) THEN
241     DO j=jMin,jMax
242     DO i=iMin,iMax
243     af(i,j) = 0.
244     ENDDO
245     ENDDO
246     ELSEIF ( rigidLid ) THEN
247     DO j=jMin,jMax
248     DO i=iMin,iMax
249     af(i,j) = rTrans(i,j)*
250     & (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
251     ENDDO
252     ENDDO
253     ELSE
254     DO j=jMin,jMax
255     DO i=iMin,iMax
256     af(i,j) = rTrans(i,j)*(
257     & maskC(i,j,kM1,bi,bj)*
258     & (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
259     & +(maskC(i,j,k,bi,bj)-maskC(i,j,kM1,bi,bj))*
260     & salt(i,j,k,bi,bj) )
261     ENDDO
262 cnh 1.1 ENDDO
263 adcroft 1.28 ENDIF
264 adcroft 1.21 C o Diffusive component of vertical flux
265     C Note: For K=1 then KM1=1 and this gives a dS/dr = 0 upper
266 adcroft 1.7 C boundary condition.
267 adcroft 1.21 IF (implicitDiffusion) THEN
268     DO j=jMin,jMax
269     DO i=iMin,iMax
270     df(i,j) = 0.
271     ENDDO
272 cnh 1.1 ENDDO
273 adcroft 1.21 ELSE
274 adcroft 1.7 DO j=jMin,jMax
275     DO i=iMin,iMax
276 adcroft 1.21 df(i,j) = - _rA(i,j,bi,bj)*(
277     & KappaRS(i,j,k)*recip_drC(k)
278 cnh 1.13 & *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))*rkFac
279 adcroft 1.7 & )
280     ENDDO
281     ENDDO
282     ENDIF
283 adcroft 1.21
284     #ifdef ALLOW_GMREDI
285 heimbach 1.22 IF (useGMRedi) CALL GMREDI_RTRANSPORT(
286 adcroft 1.21 I iMin,iMax,jMin,jMax,bi,bj,K,
287     I maskUp,salt,
288     U df,
289     I myThid)
290     #endif
291    
292 adcroft 1.18 #ifdef ALLOW_KPP
293 adcroft 1.21 C-- Add non-local KPP transport term (ghat) to diffusive salt flux.
294 heimbach 1.22 IF (useKPP) CALL KPP_TRANSPORT_S(
295 adcroft 1.21 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
296 adcroft 1.28 I KappaRS,
297 adcroft 1.21 U df )
298     #endif
299 adcroft 1.18
300 cnh 1.1 C Net vertical flux
301     DO j=jMin,jMax
302     DO i=iMin,iMax
303 adcroft 1.28 fVerS(i,j,kUp) = afFacS*af(i,j) + dfFacS*df(i,j)*maskUp(i,j)
304 cnh 1.1 ENDDO
305     ENDDO
306    
307     C-- Tendency is minus divergence of the fluxes.
308     C Note. Tendency terms will only be correct for range
309     C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
310     C will contain valid floating point numbers but
311     C they are not algorithmically correct. These points
312     C are not used.
313 adcroft 1.28 DO j=jMin,jMax
314     DO i=iMin,iMax
315 cnh 1.1 gS(i,j,k,bi,bj)=
316 adcroft 1.28 & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
317     & *recip_rA(i,j,bi,bj)
318 cnh 1.1 & *(
319     & +( fZon(i+1,j)-fZon(i,j) )
320     & +( fMer(i,j+1)-fMer(i,j) )
321 cnh 1.12 & +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )*rkFac
322 cnh 1.1 & )
323     ENDDO
324     ENDDO
325    
326 cnh 1.16 C-- External forcing term(s)
327     CALL EXTERNAL_FORCING_S(
328     I iMin,iMax,jMin,jMax,bi,bj,k,
329     I myCurrentTime,myThid)
330 cnh 1.1
331     RETURN
332     END

  ViewVC Help
Powered by ViewVC 1.1.22