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

Annotation of /MITgcm/verification/aim.5l_cs/code/calc_gs.F

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


Revision 1.1 - (hide annotations) (download)
Mon Jun 18 17:40:06 2001 UTC (22 years, 11 months ago) by cnh
Branch: MAIN
CVS Tags: checkpoint40pre3, checkpoint40pre1, checkpoint40pre6, checkpoint40pre9, checkpoint40pre8, checkpoint40pre5, checkpoint40pre7, checkpoint40pre4, checkpoint40pre2
Add to main branch of
  o CS atmos with AIM physics
  o Multi-threaded AIM physics for LatLon and CS tests
  o Tidied up monitor() output

1 cnh 1.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: $
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