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

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

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

revision 1.1 by jmc, Mon Apr 9 18:10:55 2001 UTC revision 1.2 by adcroft, Tue May 29 14:01:48 2001 UTC
# Line 0  Line 1 
1    C $Header$
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, 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
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22