/[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.3 - (show annotations) (download)
Wed Sep 19 02:43:27 2001 UTC (22 years, 8 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint40, checkpoint41
Changes since 1.2: +23 -3 lines
Re-arranged sequence of operations for Adams-Bashforth
 o this does not change numbers
 o this makes it very easy to extract forcing/diffusion out of ABII
   by changing calling sequence in calc_gt, calc_gs,...

Key modifications:
 o new s/r: ADAMS_BASHFORTH2  gT=3/2*gT-1/2*gTnm1
 o changed TIMESTEP_TRACER from gTnm1=t+dt*(3/2*gT-1/2*gTnm1)
   to  gT=T+dt*gT
 o changed CALC_GT,CALC_GS & CALC_GTR1 to calcuate "gT" defined
   by new timestep_tracer (ie. including forcing, ABII, N-L F-S, etc...)
   now calls ADAMS_BASHFORTH2 and FREESURF_RESCALE_G
 o changed CYCLE_TRACER appropriately  T=gT only

Other associated mods:
 o new s/r: FREESURF_RESCALE_G applies non-linear free-surface term
   this used to be in TIMESTEP_TRACER
 o added myIter as argument to CALC_GS,CALC_GT,CALC_GTR1

1 C $Header: /u/gcmpack/models/MITgcmUV/verification/aim.5l_LatLon/code/calc_gs.F,v 1.2 2001/05/29 14:01:48 adcroft Exp $
2 C $Name: $
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,myIter,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 #include "GAD.h"
52
53 C == Routine arguments ==
54 C fVerS - Flux of salt (S) in the vertical
55 C direction at the upper(U) and lower(D) faces of a cell.
56 C maskUp - Land mask used to denote base of the domain.
57 C xA - Tracer cell face area normal to X
58 C yA - Tracer cell face area normal to X
59 C uTrans - Zonal volume transport through cell face
60 C vTrans - Meridional volume transport through cell face
61 C rTrans - Vertical volume transport through cell face
62 C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
63 C results will be set.
64 C myThid - Instance number for this innvocation of CALC_GT
65 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
66 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72 _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
73 INTEGER k,kUp,kDown,kM1
74 INTEGER bi,bj,iMin,iMax,jMin,jMax
75 _RL myCurrentTime
76 INTEGER myIter
77 INTEGER myThid
78 CEndOfInterface
79
80 C == Local variables ==
81 C I, J, K - Loop counters
82 C tauUpwH - Horizontal upwind weight : 1=Upwind ; 0=Centered
83 C tauUpwV - Vertical upwind weight : 1=Upwind ; 0=Centered
84 INTEGER i,j
85 LOGICAL TOP_LAYER
86 _RL afFacS, dfFacS
87 _RL tauUpwH, tauUpwV
88 _RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 c_jmc:
94 _RL ddx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95 _RL d2dx2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96 _RL ddy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97 _RL d2dy2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL phiLo(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99 _RL phiHi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 _RL dist
101 c_jmc.
102
103 #ifdef ALLOW_AUTODIFF_TAMC
104 C-- only the kUp part of fverS is set in this subroutine
105 C-- the kDown is still required
106 fVerS(1,1,kDown) = fVerS(1,1,kDown)
107 #endif
108 DO j=1-OLy,sNy+OLy
109 DO i=1-OLx,sNx+OLx
110 fZon(i,j) = 0.0
111 fMer(i,j) = 0.0
112 fVerS(i,j,kUp) = 0.0
113 ENDDO
114 ENDDO
115
116 afFacS = 1. _d 0
117 dfFacS = 1. _d 0
118 tauUpwH = 1. _d 0
119 tauUpwV = 1. _d 0
120 TOP_LAYER = K .EQ. 1
121
122 C--- Calculate advective and diffusive fluxes between cells.
123
124 C o Zonal tracer gradient
125 DO j=1-Oly,sNy+Oly
126 DO i=1-Olx+1,sNx+Olx
127 fZon(i,j) = _recip_dxC(i,j,bi,bj)*xA(i,j)
128 & *(salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
129 #ifdef COSINEMETH_III
130 & *sqCosFacU(j,bi,bj)
131 #endif
132 ENDDO
133 ENDDO
134 C o Meridional tracer gradient
135 DO j=1-Oly+1,sNy+Oly
136 DO i=1-Olx,sNx+Olx
137 fMer(i,j) = _recip_dyC(i,j,bi,bj)*yA(i,j)
138 & *(salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
139 #ifdef ISOTROPIC_COS_SCALING
140 #ifdef COSINEMETH_III
141 & *sqCosFacV(j,bi,bj)
142 #endif
143 #endif
144 ENDDO
145 ENDDO
146
147 C-- del^2 of S, needed for bi-harmonic (del^4) term
148 IF (diffK4S .NE. 0.) THEN
149 DO j=1-Oly+1,sNy+Oly-1
150 DO i=1-Olx+1,sNx+Olx-1
151 df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
152 & *recip_drF(k)/_rA(i,j,bi,bj)
153 & *(
154 & +( fZon(i+1,j)-fZon(i,j) )
155 & +( fMer(i,j+1)-fMer(i,j) )
156 & )
157 ENDDO
158 ENDDO
159 ENDIF
160
161 C-- Zonal flux (fZon is at west face of "salt" cell)
162 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
163 #ifdef USE_3RD_O_ADVEC
164 C o Advective component of zonal flux, 3rd order Advec Scheme
165 DO j=jMin,jMax
166 DO i=1-OLx+1,sNx+OLx
167 ddx(i,j) = (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
168 & *_recip_dxC(i,j,bi,bj)*_maskW(i,j,k,bi,bj)
169 ENDDO
170 ENDDO
171 DO j=jMin,jMax
172 DO i=1-OLx,sNx+OLx-1
173 d2dx2(i,j) = ( ddx(i+1,j)-ddx(i,j) )
174 & *_recip_dxF(i,j,bi,bj)*maskC(i,j,k,bi,bj)
175 ENDDO
176 ENDDO
177 DO j=jMin,jMax
178 DO i=1-OLx+1,sNx+OLx
179 dist = _dxF(i-1,j,bi,bj)*0.5 _d 0
180 phiLo(i,j) = salt(i-1,j,k,bi,bj)
181 & +dist
182 & *( ddx(i ,j)+ddx(i-1,j) )*0.5 _d 0
183 & +0.5 _d 0*dist*dist
184 & *d2dx2(i-1,j)
185 dist = -_dxF(i,j,bi,bj)*0.5 _d 0
186 phiHi(i,j) = salt(i,j,k,bi,bj)
187 & +dist
188 & *( ddx(i+1,j)+ddx(i ,j) )*0.5 _d 0
189 & +0.5 _d 0*dist*dist
190 & *d2dx2(i,j)
191 ENDDO
192 ENDDO
193 DO j=jMin,jMax
194 DO i=1-OLx,sNx+OLx
195 IF ( uTrans(i,j) .GT. 0. ) THEN
196 af(i,j) = uTrans(i,j)*phiLo(i,j)
197 ELSE
198 af(i,j) = uTrans(i,j)*phiHi(i,j)
199 ENDIF
200 ENDDO
201 ENDDO
202 #else
203 C o Advective component of zonal flux, 1rst & 2nd order Advec Scheme
204 IF (tauUpwH.EQ.0. _d 0) THEN
205 C Centered scheme :
206 DO j=jMin,jMax
207 DO i=iMin,iMax
208 af(i,j) = uTrans(i,j)*
209 & (salt(i-1,j,k,bi,bj)+salt(i,j,k,bi,bj))*0.5 _d 0
210 ENDDO
211 ENDDO
212 ELSE
213 C Upwind weighted scheme :
214 DO j=jMin,jMax
215 DO i=iMin,iMax
216 af(i,j) = uTrans(i,j)*
217 & (salt(i-1,j,k,bi,bj)+salt(i,j,k,bi,bj))*0.5 _d 0
218 & +tauUpwH*abs(uTrans(i,j))*
219 & (salt(i-1,j,k,bi,bj)-salt(i,j,k,bi,bj))*0.5 _d 0
220 ENDDO
221 ENDDO
222 ENDIF
223 #endif
224 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
225 C o Diffusive component of zonal flux
226 DO j=jMin,jMax
227 DO i=iMin,iMax
228 df(i,j) = -diffKhS*xA(i,j)*_recip_dxC(i,j,bi,bj)*
229 & (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
230 & *CosFacU(j,bi,bj)
231 ENDDO
232 ENDDO
233 #ifdef ALLOW_GMREDI
234 IF (useGMRedi) CALL GMREDI_XTRANSPORT(
235 I iMin,iMax,jMin,jMax,bi,bj,K,
236 I xA,salt,
237 U df,
238 I myThid)
239 #endif
240 C o Add the bi-harmonic contribution
241 IF (diffK4S .NE. 0.) THEN
242 DO j=jMin,jMax
243 DO i=iMin,iMax
244 df(i,j) = df(i,j) + xA(i,j)*
245 & diffK4S*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)
246 #ifdef COSINEMETH_III
247 & *sqCosFacU(j,bi,bj)
248 #else
249 & *CosFacU(j,bi,bj)
250 #endif
251 ENDDO
252 ENDDO
253 ENDIF
254 C Net zonal flux
255 DO j=jMin,jMax
256 DO i=iMin,iMax
257 fZon(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
258 ENDDO
259 ENDDO
260
261 C-- Meridional flux (fMer is at south face of "salt" cell)
262 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
263 #ifdef USE_3RD_O_ADVEC
264 C o Advective component of meridional flux, 3rd order Advec Scheme
265 DO j=jMin,jMax
266 DO i=iMin,iMax
267 ddy(i,j) = (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
268 & *_recip_dyC(i,j,bi,bj)*_maskS(i,j,k,bi,bj)
269 ENDDO
270 ENDDO
271 DO j=jMin,jMax
272 DO i=iMin,iMax
273 d2dy2(i,j) = ( ddy(i,j+1)-ddy(i,j) )
274 & *_recip_dyF(i,j,bi,bj)*maskC(i,j,k,bi,bj)
275 ENDDO
276 ENDDO
277 DO j=jMin,jMax
278 DO i=iMin,iMax
279 dist = _dyF(i,j-1,bi,bj)*0.5 _d 0
280 phiLo(i,j) = salt(i,j-1,k,bi,bj)
281 & +dist
282 & *( ddy(i ,j)+ddy(i,j-1) )*0.5 _d 0
283 & +0.5 _d 0*dist*dist
284 & *d2dy2(i,j-1)
285 dist = -_dyF(i,j,bi,bj)*0.5 _d 0
286 phiHi(i,j) = salt(i,j,k,bi,bj)
287 & +dist
288 & *( ddy(i,j+1)+ddy(i ,j) )*0.5 _d 0
289 & +0.5 _d 0*dist*dist
290 & *d2dy2(i,j)
291 ENDDO
292 ENDDO
293 DO j=jMin,jMax
294 DO i=iMin,iMax
295 IF ( vTrans(i,j) .GT. 0. ) THEN
296 af(i,j) = vTrans(i,j)*phiLo(i,j)
297 ELSE
298 af(i,j) = vTrans(i,j)*phiHi(i,j)
299 ENDIF
300 ENDDO
301 ENDDO
302 #else
303 C o Advective component of meridional flux, 1rst & 2nd order Advec Scheme
304 IF (tauUpwH.EQ.0. _d 0) THEN
305 C Centered scheme :
306 DO j=jMin,jMax
307 DO i=iMin,iMax
308 af(i,j) = vTrans(i,j)*
309 & (salt(i,j-1,k,bi,bj)+salt(i,j,k,bi,bj))*0.5 _d 0
310 ENDDO
311 ENDDO
312 ELSE
313 C Upwind weighted scheme :
314 DO j=jMin,jMax
315 DO i=iMin,iMax
316 af(i,j) = vTrans(i,j)*
317 & (salt(i,j-1,k,bi,bj)+salt(i,j,k,bi,bj))*0.5 _d 0
318 & +tauUpwH*abs(vTrans(i,j))*
319 & (salt(i,j-1,k,bi,bj)-salt(i,j,k,bi,bj))*0.5 _d 0
320 ENDDO
321 ENDDO
322 ENDIF
323 #endif
324 c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
325 C o Diffusive component of meridional flux
326 DO j=jMin,jMax
327 DO i=iMin,iMax
328 df(i,j) = -diffKhS*yA(i,j)*_recip_dyC(i,j,bi,bj)*
329 & (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
330 & *CosFacV(j,bi,bj)
331 ENDDO
332 ENDDO
333 #ifdef ALLOW_GMREDI
334 IF (useGMRedi) CALL GMREDI_YTRANSPORT(
335 I iMin,iMax,jMin,jMax,bi,bj,K,
336 I yA,salt,
337 U df,
338 I myThid)
339 #endif
340 C o Add the bi-harmonic contribution
341 IF (diffK4S .NE. 0.) THEN
342 DO j=jMin,jMax
343 DO i=iMin,iMax
344 df(i,j) = df(i,j) + yA(i,j)*
345 & diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
346 #ifdef ISOTROPIC_COS_SCALING
347 #ifdef COSINEMETH_III
348 & *sqCosFacV(j,bi,bj)
349 #else
350 & *CosFacV(j,bi,bj)
351 #endif
352 #endif
353 ENDDO
354 ENDDO
355 ENDIF
356
357 C Net meridional flux
358 DO j=jMin,jMax
359 DO i=iMin,iMax
360 fMer(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
361 ENDDO
362 ENDDO
363
364 C-- Vertical flux ( fVerS(,,kUp) is at upper face of "Tracer" cell )
365 C o Advective component of vertical flux : assume W_bottom=0 (mask)
366 C Note: For K=1 then KM1=1 this gives a barZ(S) = S
367 C (this plays the role of the free-surface correction for k=1)
368 IF ( rigidLid .AND. TOP_LAYER) THEN
369 DO j=jMin,jMax
370 DO i=iMin,iMax
371 af(i,j) = 0.
372 ENDDO
373 ENDDO
374 ELSE
375 IF (tauUpwV.EQ.0. _d 0) THEN
376 C Centered scheme :
377 DO j=jMin,jMax
378 DO i=iMin,iMax
379 af(i,j) = rTrans(i,j)*
380 & (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
381 ENDDO
382 ENDDO
383 ELSE
384 C Upwind weighted scheme :
385 DO j=jMin,jMax
386 DO i=iMin,iMax
387 af(i,j) = rTrans(i,j)*
388 & (salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
389 & +tauUpwV*abs(rTrans(i,j))*
390 & (salt(i,j,k,bi,bj)-salt(i,j,kM1,bi,bj))*0.5 _d 0
391 ENDDO
392 ENDDO
393 ENDIF
394 IF (.NOT.rigidLid ) THEN
395 C free-surface correction for k > 1
396 DO j=jMin,jMax
397 DO i=iMin,iMax
398 af(i,j) = af(i,j)*maskC(i,j,kM1,bi,bj)
399 & +rTrans(i,j)*(maskC(i,j,k,bi,bj)-maskC(i,j,kM1,bi,bj))*
400 & salt(i,j,k,bi,bj)
401 ENDDO
402 ENDDO
403 ENDIF
404 ENDIF
405 C o Diffusive component of vertical flux
406 C Note: For K=1 then KM1=1 and this gives a dS/dr = 0 upper
407 C boundary condition.
408 IF (implicitDiffusion) THEN
409 DO j=jMin,jMax
410 DO i=iMin,iMax
411 df(i,j) = 0.
412 ENDDO
413 ENDDO
414 ELSE
415 DO j=jMin,jMax
416 DO i=iMin,iMax
417 df(i,j) = - _rA(i,j,bi,bj)*(
418 & KappaRS(i,j,k)*recip_drC(k)
419 & *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))*rkFac
420 & )
421 ENDDO
422 ENDDO
423 ENDIF
424
425 #ifdef ALLOW_GMREDI
426 IF (useGMRedi) CALL GMREDI_RTRANSPORT(
427 I iMin,iMax,jMin,jMax,bi,bj,K,
428 I maskUp,salt,
429 U df,
430 I myThid)
431 #endif
432
433 #ifdef ALLOW_KPP
434 C-- Add non-local KPP transport term (ghat) to diffusive salt flux.
435 IF (useKPP) CALL KPP_TRANSPORT_S(
436 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
437 I KappaRS,
438 U df )
439 #endif
440
441 C Net vertical flux
442 DO j=jMin,jMax
443 DO i=iMin,iMax
444 fVerS(i,j,kUp) = afFacS*af(i,j) + dfFacS*df(i,j)*maskUp(i,j)
445 ENDDO
446 ENDDO
447
448 C-- Tendency is minus divergence of the fluxes.
449 C Note. Tendency terms will only be correct for range
450 C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
451 C will contain valid floating point numbers but
452 C they are not algorithmically correct. These points
453 C are not used.
454 DO j=jMin,jMax
455 DO i=iMin,iMax
456 #define _recip_VolS1(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
457 #define _recip_VolS2(i,j,k,bi,bj) /_rA(i,j,bi,bj)
458 gS(i,j,k,bi,bj)=
459 & -_recip_VolS1(i,j,k,bi,bj)
460 & _recip_VolS2(i,j,k,bi,bj)
461 & *(
462 & +( fZon(i+1,j)-fZon(i,j) )
463 & +( fMer(i,j+1)-fMer(i,j) )
464 & +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )*rkFac
465 & )
466 ENDDO
467 ENDDO
468
469 C-- External forcing term(s)
470 CALL EXTERNAL_FORCING_S(
471 I iMin,iMax,jMin,jMax,bi,bj,k,
472 I myCurrentTime,myThid)
473
474 IF ( saltAdvScheme.EQ.ENUM_CENTERED_2ND
475 & .OR.saltAdvScheme.EQ.ENUM_UPWIND_3RD
476 & .OR.saltAdvScheme.EQ.ENUM_CENTERED_4TH ) THEN
477 CALL ADAMS_BASHFORTH2(
478 I bi, bj, K,
479 U gS, gSnm1,
480 I myIter, myThid )
481 ENDIF
482
483 #ifdef NONLIN_FRSURF
484 IF (nonlinFreeSurf.GT.0) THEN
485 CALL FREESURF_RESCALE_G(
486 I bi, bj, K,
487 U gS,
488 I myThid )
489 ENDIF
490 #endif /* NONLIN_FRSURF */
491
492 RETURN
493 END

  ViewVC Help
Powered by ViewVC 1.1.22