/[MITgcm]/MITgcm/verification/aim.5l_LatLon/code/calc_gs.F
ViewVC logotype

Contents of /MITgcm/verification/aim.5l_LatLon/code/calc_gs.F

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


Revision 1.2 - (show annotations) (download)
Tue May 29 14:01:48 2001 UTC (23 years ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint40pre3, checkpoint40pre1, checkpoint40pre6, checkpoint40pre9, checkpoint40pre8, checkpoint40pre5, checkpoint40pre7, checkpoint40pre4, checkpoint40pre2
Changes since 1.1: +473 -0 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 C $Header: /u/gcmpack/models/MITgcmUV/verification/aim.5l_LatLon/code/Attic/calc_gs.F,v 1.1.2.2 2001/04/17 01:38:31 jmc Exp $
2 C $Name: pre38-close $
3
4 #include "CPP_OPTIONS.h"
5
6 #define COSINEMETH_III
7 #undef ISOTROPIC_COS_SCALING
8 #define USE_3RD_O_ADVEC
9
10 CStartOfInterFace
11 SUBROUTINE CALC_GS(
12 I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
13 I xA,yA,uTrans,vTrans,rTrans,maskUp,
14 I KappaRS,
15 U fVerS,
16 I myCurrentTime, myThid )
17 C /==========================================================\
18 C | SUBROUTINE CALC_GS |
19 C | o Calculate the salt tendency terms. |
20 C |==========================================================|
21 C | A procedure called EXTERNAL_FORCING_S is called from |
22 C | here. These procedures can be used to add per problem |
23 C | E-P flux source terms. |
24 C | Note: Although it is slightly counter-intuitive the |
25 C | EXTERNAL_FORCING routine is not the place to put |
26 C | file I/O. Instead files that are required to |
27 C | calculate the external source terms are generally |
28 C | read during the model main loop. This makes the |
29 C | logisitics of multi-processing simpler and also |
30 C | makes the adjoint generation simpler. It also |
31 C | allows for I/O to overlap computation where that |
32 C | is supported by hardware. |
33 C | Aside from the problem specific term the code here |
34 C | forms the tendency terms due to advection and mixing |
35 C | The baseline implementation here uses a centered |
36 C | difference form for the advection term and a tensorial |
37 C | divergence of a flux form for the diffusive term. The |
38 C | diffusive term is formulated so that isopycnal mixing and|
39 C | GM-style subgrid-scale terms can be incorporated b simply|
40 C | setting the diffusion tensor terms appropriately. |
41 C \==========================================================/
42 IMPLICIT NONE
43
44 C == GLobal variables ==
45 #include "SIZE.h"
46 #include "DYNVARS.h"
47 #include "EEPARAMS.h"
48 #include "PARAMS.h"
49 #include "GRID.h"
50 #include "FFIELDS.h"
51
52 C == Routine arguments ==
53 C fVerS - Flux of salt (S) in the vertical
54 C direction at the upper(U) and lower(D) faces of a cell.
55 C maskUp - Land mask used to denote base of the domain.
56 C xA - Tracer cell face area normal to X
57 C yA - Tracer cell face area normal to X
58 C uTrans - Zonal volume transport through cell face
59 C vTrans - Meridional volume transport through cell face
60 C rTrans - Vertical volume transport through cell face
61 C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
62 C results will be set.
63 C myThid - Instance number for this innvocation of CALC_GT
64 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
65 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
72 INTEGER k,kUp,kDown,kM1
73 INTEGER bi,bj,iMin,iMax,jMin,jMax
74 _RL myCurrentTime
75 INTEGER myThid
76 CEndOfInterface
77
78 C == Local variables ==
79 C I, J, K - Loop counters
80 C tauUpwH - Horizontal upwind weight : 1=Upwind ; 0=Centered
81 C tauUpwV - Vertical upwind weight : 1=Upwind ; 0=Centered
82 INTEGER i,j
83 LOGICAL TOP_LAYER
84 _RL afFacS, dfFacS
85 _RL tauUpwH, tauUpwV
86 _RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 c_jmc:
92 _RL ddx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 _RL d2dx2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 _RL ddy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95 _RL d2dy2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96 _RL phiLo(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97 _RL phiHi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL dist
99 c_jmc.
100
101 #ifdef ALLOW_AUTODIFF_TAMC
102 C-- only the kUp part of fverS is set in this subroutine
103 C-- the kDown is still required
104 fVerS(1,1,kDown) = fVerS(1,1,kDown)
105 #endif
106 DO j=1-OLy,sNy+OLy
107 DO i=1-OLx,sNx+OLx
108 fZon(i,j) = 0.0
109 fMer(i,j) = 0.0
110 fVerS(i,j,kUp) = 0.0
111 ENDDO
112 ENDDO
113
114 afFacS = 1. _d 0
115 dfFacS = 1. _d 0
116 tauUpwH = 1. _d 0
117 tauUpwV = 1. _d 0
118 TOP_LAYER = K .EQ. 1
119
120 C--- Calculate advective and diffusive fluxes between cells.
121
122 C o Zonal tracer gradient
123 DO j=1-Oly,sNy+Oly
124 DO i=1-Olx+1,sNx+Olx
125 fZon(i,j) = _recip_dxC(i,j,bi,bj)*xA(i,j)
126 & *(salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
127 #ifdef COSINEMETH_III
128 & *sqCosFacU(j,bi,bj)
129 #endif
130 ENDDO
131 ENDDO
132 C o Meridional tracer gradient
133 DO j=1-Oly+1,sNy+Oly
134 DO i=1-Olx,sNx+Olx
135 fMer(i,j) = _recip_dyC(i,j,bi,bj)*yA(i,j)
136 & *(salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
137 #ifdef ISOTROPIC_COS_SCALING
138 #ifdef COSINEMETH_III
139 & *sqCosFacV(j,bi,bj)
140 #endif
141 #endif
142 ENDDO
143 ENDDO
144
145 C-- del^2 of S, needed for bi-harmonic (del^4) term
146 IF (diffK4S .NE. 0.) THEN
147 DO j=1-Oly+1,sNy+Oly-1
148 DO i=1-Olx+1,sNx+Olx-1
149 df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
150 & *recip_drF(k)/_rA(i,j,bi,bj)
151 & *(
152 & +( fZon(i+1,j)-fZon(i,j) )
153 & +( fMer(i,j+1)-fMer(i,j) )
154 & )
155 ENDDO
156 ENDDO
157 ENDIF
158
159 C-- Zonal flux (fZon is at west face of "salt" cell)
160 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
161 #ifdef USE_3RD_O_ADVEC
162 C o Advective component of zonal flux, 3rd order Advec Scheme
163 DO j=jMin,jMax
164 DO i=1-OLx+1,sNx+OLx
165 ddx(i,j) = (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
166 & *_recip_dxC(i,j,bi,bj)*_maskW(i,j,k,bi,bj)
167 ENDDO
168 ENDDO
169 DO j=jMin,jMax
170 DO i=1-OLx,sNx+OLx-1
171 d2dx2(i,j) = ( ddx(i+1,j)-ddx(i,j) )
172 & *_recip_dxF(i,j,bi,bj)*maskC(i,j,k,bi,bj)
173 ENDDO
174 ENDDO
175 DO j=jMin,jMax
176 DO i=1-OLx+1,sNx+OLx
177 dist = _dxF(i-1,j,bi,bj)*0.5 _d 0
178 phiLo(i,j) = salt(i-1,j,k,bi,bj)
179 & +dist
180 & *( ddx(i ,j)+ddx(i-1,j) )*0.5 _d 0
181 & +0.5 _d 0*dist*dist
182 & *d2dx2(i-1,j)
183 dist = -_dxF(i,j,bi,bj)*0.5 _d 0
184 phiHi(i,j) = salt(i,j,k,bi,bj)
185 & +dist
186 & *( ddx(i+1,j)+ddx(i ,j) )*0.5 _d 0
187 & +0.5 _d 0*dist*dist
188 & *d2dx2(i,j)
189 ENDDO
190 ENDDO
191 DO j=jMin,jMax
192 DO i=1-OLx,sNx+OLx
193 IF ( uTrans(i,j) .GT. 0. ) THEN
194 af(i,j) = uTrans(i,j)*phiLo(i,j)
195 ELSE
196 af(i,j) = uTrans(i,j)*phiHi(i,j)
197 ENDIF
198 ENDDO
199 ENDDO
200 #else
201 C o Advective component of zonal flux, 1rst & 2nd order Advec Scheme
202 IF (tauUpwH.EQ.0. _d 0) THEN
203 C Centered scheme :
204 DO j=jMin,jMax
205 DO i=iMin,iMax
206 af(i,j) = uTrans(i,j)*
207 & (salt(i-1,j,k,bi,bj)+salt(i,j,k,bi,bj))*0.5 _d 0
208 ENDDO
209 ENDDO
210 ELSE
211 C Upwind weighted scheme :
212 DO j=jMin,jMax
213 DO i=iMin,iMax
214 af(i,j) = uTrans(i,j)*
215 & (salt(i-1,j,k,bi,bj)+salt(i,j,k,bi,bj))*0.5 _d 0
216 & +tauUpwH*abs(uTrans(i,j))*
217 & (salt(i-1,j,k,bi,bj)-salt(i,j,k,bi,bj))*0.5 _d 0
218 ENDDO
219 ENDDO
220 ENDIF
221 #endif
222 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
223 C o Diffusive component of zonal flux
224 DO j=jMin,jMax
225 DO i=iMin,iMax
226 df(i,j) = -diffKhS*xA(i,j)*_recip_dxC(i,j,bi,bj)*
227 & (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
228 & *CosFacU(j,bi,bj)
229 ENDDO
230 ENDDO
231 #ifdef ALLOW_GMREDI
232 IF (useGMRedi) CALL GMREDI_XTRANSPORT(
233 I iMin,iMax,jMin,jMax,bi,bj,K,
234 I xA,salt,
235 U df,
236 I myThid)
237 #endif
238 C o Add the bi-harmonic contribution
239 IF (diffK4S .NE. 0.) THEN
240 DO j=jMin,jMax
241 DO i=iMin,iMax
242 df(i,j) = df(i,j) + xA(i,j)*
243 & diffK4S*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)
244 #ifdef COSINEMETH_III
245 & *sqCosFacU(j,bi,bj)
246 #else
247 & *CosFacU(j,bi,bj)
248 #endif
249 ENDDO
250 ENDDO
251 ENDIF
252 C Net zonal flux
253 DO j=jMin,jMax
254 DO i=iMin,iMax
255 fZon(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
256 ENDDO
257 ENDDO
258
259 C-- Meridional flux (fMer is at south face of "salt" cell)
260 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
261 #ifdef USE_3RD_O_ADVEC
262 C o Advective component of meridional flux, 3rd order Advec Scheme
263 DO j=jMin,jMax
264 DO i=iMin,iMax
265 ddy(i,j) = (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
266 & *_recip_dyC(i,j,bi,bj)*_maskS(i,j,k,bi,bj)
267 ENDDO
268 ENDDO
269 DO j=jMin,jMax
270 DO i=iMin,iMax
271 d2dy2(i,j) = ( ddy(i,j+1)-ddy(i,j) )
272 & *_recip_dyF(i,j,bi,bj)*maskC(i,j,k,bi,bj)
273 ENDDO
274 ENDDO
275 DO j=jMin,jMax
276 DO i=iMin,iMax
277 dist = _dyF(i,j-1,bi,bj)*0.5 _d 0
278 phiLo(i,j) = salt(i,j-1,k,bi,bj)
279 & +dist
280 & *( ddy(i ,j)+ddy(i,j-1) )*0.5 _d 0
281 & +0.5 _d 0*dist*dist
282 & *d2dy2(i,j-1)
283 dist = -_dyF(i,j,bi,bj)*0.5 _d 0
284 phiHi(i,j) = salt(i,j,k,bi,bj)
285 & +dist
286 & *( ddy(i,j+1)+ddy(i ,j) )*0.5 _d 0
287 & +0.5 _d 0*dist*dist
288 & *d2dy2(i,j)
289 ENDDO
290 ENDDO
291 DO j=jMin,jMax
292 DO i=iMin,iMax
293 IF ( vTrans(i,j) .GT. 0. ) THEN
294 af(i,j) = vTrans(i,j)*phiLo(i,j)
295 ELSE
296 af(i,j) = vTrans(i,j)*phiHi(i,j)
297 ENDIF
298 ENDDO
299 ENDDO
300 #else
301 C o Advective component of meridional flux, 1rst & 2nd order Advec Scheme
302 IF (tauUpwH.EQ.0. _d 0) THEN
303 C Centered scheme :
304 DO j=jMin,jMax
305 DO i=iMin,iMax
306 af(i,j) = vTrans(i,j)*
307 & (salt(i,j-1,k,bi,bj)+salt(i,j,k,bi,bj))*0.5 _d 0
308 ENDDO
309 ENDDO
310 ELSE
311 C Upwind weighted scheme :
312 DO j=jMin,jMax
313 DO i=iMin,iMax
314 af(i,j) = vTrans(i,j)*
315 & (salt(i,j-1,k,bi,bj)+salt(i,j,k,bi,bj))*0.5 _d 0
316 & +tauUpwH*abs(vTrans(i,j))*
317 & (salt(i,j-1,k,bi,bj)-salt(i,j,k,bi,bj))*0.5 _d 0
318 ENDDO
319 ENDDO
320 ENDIF
321 #endif
322 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
323 C o Diffusive component of meridional flux
324 DO j=jMin,jMax
325 DO i=iMin,iMax
326 df(i,j) = -diffKhS*yA(i,j)*_recip_dyC(i,j,bi,bj)*
327 & (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
328 & *CosFacV(j,bi,bj)
329 ENDDO
330 ENDDO
331 #ifdef ALLOW_GMREDI
332 IF (useGMRedi) CALL GMREDI_YTRANSPORT(
333 I iMin,iMax,jMin,jMax,bi,bj,K,
334 I yA,salt,
335 U df,
336 I myThid)
337 #endif
338 C o Add the bi-harmonic contribution
339 IF (diffK4S .NE. 0.) THEN
340 DO j=jMin,jMax
341 DO i=iMin,iMax
342 df(i,j) = df(i,j) + yA(i,j)*
343 & diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
344 #ifdef ISOTROPIC_COS_SCALING
345 #ifdef COSINEMETH_III
346 & *sqCosFacV(j,bi,bj)
347 #else
348 & *CosFacV(j,bi,bj)
349 #endif
350 #endif
351 ENDDO
352 ENDDO
353 ENDIF
354
355 C Net meridional flux
356 DO j=jMin,jMax
357 DO i=iMin,iMax
358 fMer(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
359 ENDDO
360 ENDDO
361
362 C-- Vertical flux ( fVerS(,,kUp) is at upper face of "Tracer" cell )
363 C o Advective component of vertical flux : assume W_bottom=0 (mask)
364 C Note: For K=1 then KM1=1 this gives a barZ(S) = S
365 C (this plays the role of the free-surface correction for k=1)
366 IF ( rigidLid .AND. TOP_LAYER) THEN
367 DO j=jMin,jMax
368 DO i=iMin,iMax
369 af(i,j) = 0.
370 ENDDO
371 ENDDO
372 ELSE
373 IF (tauUpwV.EQ.0. _d 0) THEN
374 C Centered scheme :
375 DO j=jMin,jMax
376 DO i=iMin,iMax
377 af(i,j) = rTrans(i,j)*
378 & (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
379 ENDDO
380 ENDDO
381 ELSE
382 C Upwind weighted scheme :
383 DO j=jMin,jMax
384 DO i=iMin,iMax
385 af(i,j) = rTrans(i,j)*
386 & (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
387 & +tauUpwV*abs(rTrans(i,j))*
388 & (salt(i,j,k,bi,bj)-salt(i,j,kM1,bi,bj))*0.5 _d 0
389 ENDDO
390 ENDDO
391 ENDIF
392 IF (.NOT.rigidLid ) THEN
393 C free-surface correction for k > 1
394 DO j=jMin,jMax
395 DO i=iMin,iMax
396 af(i,j) = af(i,j)*maskC(i,j,kM1,bi,bj)
397 & +rTrans(i,j)*(maskC(i,j,k,bi,bj)-maskC(i,j,kM1,bi,bj))*
398 & salt(i,j,k,bi,bj)
399 ENDDO
400 ENDDO
401 ENDIF
402 ENDIF
403 C o Diffusive component of vertical flux
404 C Note: For K=1 then KM1=1 and this gives a dS/dr = 0 upper
405 C boundary condition.
406 IF (implicitDiffusion) THEN
407 DO j=jMin,jMax
408 DO i=iMin,iMax
409 df(i,j) = 0.
410 ENDDO
411 ENDDO
412 ELSE
413 DO j=jMin,jMax
414 DO i=iMin,iMax
415 df(i,j) = - _rA(i,j,bi,bj)*(
416 & KappaRS(i,j,k)*recip_drC(k)
417 & *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))*rkFac
418 & )
419 ENDDO
420 ENDDO
421 ENDIF
422
423 #ifdef ALLOW_GMREDI
424 IF (useGMRedi) CALL GMREDI_RTRANSPORT(
425 I iMin,iMax,jMin,jMax,bi,bj,K,
426 I maskUp,salt,
427 U df,
428 I myThid)
429 #endif
430
431 #ifdef ALLOW_KPP
432 C-- Add non-local KPP transport term (ghat) to diffusive salt flux.
433 IF (useKPP) CALL KPP_TRANSPORT_S(
434 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
435 I KappaRS,
436 U df )
437 #endif
438
439 C Net vertical flux
440 DO j=jMin,jMax
441 DO i=iMin,iMax
442 fVerS(i,j,kUp) = afFacS*af(i,j) + dfFacS*df(i,j)*maskUp(i,j)
443 ENDDO
444 ENDDO
445
446 C-- Tendency is minus divergence of the fluxes.
447 C Note. Tendency terms will only be correct for range
448 C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
449 C will contain valid floating point numbers but
450 C they are not algorithmically correct. These points
451 C are not used.
452 DO j=jMin,jMax
453 DO i=iMin,iMax
454 #define _recip_VolS1(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
455 #define _recip_VolS2(i,j,k,bi,bj) /_rA(i,j,bi,bj)
456 gS(i,j,k,bi,bj)=
457 & -_recip_VolS1(i,j,k,bi,bj)
458 & _recip_VolS2(i,j,k,bi,bj)
459 & *(
460 & +( fZon(i+1,j)-fZon(i,j) )
461 & +( fMer(i,j+1)-fMer(i,j) )
462 & +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )*rkFac
463 & )
464 ENDDO
465 ENDDO
466
467 C-- External forcing term(s)
468 CALL EXTERNAL_FORCING_S(
469 I iMin,iMax,jMin,jMax,bi,bj,k,
470 I myCurrentTime,myThid)
471
472 RETURN
473 END

  ViewVC Help
Powered by ViewVC 1.1.22