/[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.20 - (hide annotations) (download)
Fri Jun 9 14:26:30 2000 UTC (23 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint28
Changes since 1.19: +15 -1 lines
Included initialisations required for TAMC compatibility.

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

  ViewVC Help
Powered by ViewVC 1.1.22