/[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.19 - (hide annotations) (download)
Wed May 26 20:26:41 1999 UTC (24 years, 11 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint23, checkpoint24, checkpoint25, checkpoint27, checkpoint26
Changes since 1.18: +53 -16 lines
Added bi-harmonic diffusion to calc_gs (Salt) and calc_gt (Temperature).

1 adcroft 1.19 C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gs.F,v 1.18 1999/05/18 18:01:12 adcroft 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 cnh 1.1
100     afFacS = 1. _d 0
101     dfFacS = 1. _d 0
102 adcroft 1.7 TOP_LAYER = K .EQ. 1
103 cnh 1.1
104     C--- Calculate advective and diffusive fluxes between cells.
105    
106 adcroft 1.19 #ifdef INCLUDE_T_DIFFUSION_CODE
107     C o Zonal tracer gradient
108     DO j=1-Oly,sNy+Oly
109     DO i=1-Olx+1,sNx+Olx
110     dSdx(i,j) = _recip_dxC(i,j,bi,bj)*
111     & (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
112     ENDDO
113     ENDDO
114     C o Meridional tracer gradient
115     DO j=1-Oly+1,sNy+Oly
116     DO i=1-Olx,sNx+Olx
117     dSdy(i,j) = _recip_dyC(i,j,bi,bj)*
118     & (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
119     ENDDO
120     ENDDO
121    
122     C-- del^2 of S, needed for bi-harmonic (del^4) term
123     IF (diffK4S .NE. 0.) THEN
124     DO j=1-Oly+1,sNy+Oly-1
125     DO i=1-Olx+1,sNx+Olx-1
126     df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
127     & *recip_drF(k)/_rA(i,j,bi,bj)
128     & *(
129     & +( xA(i+1,j)*dSdx(i+1,j)-xA(i,j)*dSdx(i,j) )
130     & +( yA(i,j+1)*dSdy(i,j+1)-yA(i,j)*dSdy(i,j) )
131     & )
132     ENDDO
133     ENDDO
134     ENDIF
135     #endif
136    
137 cnh 1.1 C-- Zonal flux (fZon is at west face of "salt" cell)
138     C Advective component of zonal flux
139     DO j=jMin,jMax
140     DO i=iMin,iMax
141     af(i,j) =
142     & uTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj))*0.5 _d 0
143     ENDDO
144     ENDDO
145 adcroft 1.19 C o Diffusive component of zonal flux
146 cnh 1.1 DO j=jMin,jMax
147     DO i=iMin,iMax
148 adcroft 1.7 df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i-1,j)))*
149     & xA(i,j)*dSdx(i,j)
150 cnh 1.1 ENDDO
151     ENDDO
152 adcroft 1.19 C o Add the bi-harmonic contribution
153     IF (diffK4S .NE. 0.) THEN
154     DO j=jMin,jMax
155     DO i=iMin,iMax
156     df(i,j) = df(i,j) + xA(i,j)*
157     & diffK4S*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)
158     ENDDO
159     ENDDO
160     ENDIF
161 cnh 1.1 C Net zonal flux
162     DO j=jMin,jMax
163     DO i=iMin,iMax
164     fZon(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
165     ENDDO
166     ENDDO
167    
168     C-- Meridional flux (fMer is at south face of "salt" cell)
169     C Advective component of meridional flux
170     DO j=jMin,jMax
171     DO i=iMin,iMax
172     C Advective component of meridional flux
173     af(i,j) =
174     & vTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj))*0.5 _d 0
175     ENDDO
176     ENDDO
177     C Diffusive component of meridional flux
178     DO j=jMin,jMax
179     DO i=iMin,iMax
180 adcroft 1.7 df(i,j) = -(diffKhS+0.5*(KapGM(i,j)+KapGM(i,j-1)))*
181     & yA(i,j)*dSdy(i,j)
182 cnh 1.1 ENDDO
183     ENDDO
184 adcroft 1.19 C o Add the bi-harmonic contribution
185     IF (diffK4S .NE. 0.) THEN
186     DO j=jMin,jMax
187     DO i=iMin,iMax
188     df(i,j) = df(i,j) + yA(i,j)*
189     & diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
190     ENDDO
191     ENDDO
192     ENDIF
193    
194 cnh 1.1 C Net meridional flux
195     DO j=jMin,jMax
196     DO i=iMin,iMax
197     fMer(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
198     ENDDO
199     ENDDO
200    
201 adcroft 1.7 C-- Interpolate terms for Redi/GM scheme
202     DO j=jMin,jMax
203     DO i=iMin,iMax
204     dSdx(i,j) = 0.5*(
205 cnh 1.15 & +0.5*(_maskW(i+1,j,k,bi,bj)
206     & *_recip_dxC(i+1,j,bi,bj)*
207 adcroft 1.7 & (salt(i+1,j,k,bi,bj)-salt(i,j,k,bi,bj))
208 cnh 1.15 & +_maskW(i,j,k,bi,bj)
209     & *_recip_dxC(i,j,bi,bj)*
210 adcroft 1.7 & (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj)))
211 cnh 1.15 & +0.5*(_maskW(i+1,j,km1,bi,bj)
212     & *_recip_dxC(i+1,j,bi,bj)*
213 adcroft 1.7 & (salt(i+1,j,km1,bi,bj)-salt(i,j,km1,bi,bj))
214 cnh 1.15 & +_maskW(i,j,km1,bi,bj)
215     & *_recip_dxC(i,j,bi,bj)*
216 adcroft 1.7 & (salt(i,j,km1,bi,bj)-salt(i-1,j,km1,bi,bj)))
217     & )
218     ENDDO
219     ENDDO
220     DO j=jMin,jMax
221     DO i=iMin,iMax
222     dSdy(i,j) = 0.5*(
223 cnh 1.15 & +0.5*(_maskS(i,j,k,bi,bj)
224     & *_recip_dyC(i,j,bi,bj)*
225 adcroft 1.7 & (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
226 cnh 1.15 & +_maskS(i,j+1,k,bi,bj)
227     & *_recip_dyC(i,j+1,bi,bj)*
228 adcroft 1.7 & (salt(i,j+1,k,bi,bj)-salt(i,j,k,bi,bj)))
229 cnh 1.15 & +0.5*(_maskS(i,j,km1,bi,bj)
230     & *_recip_dyC(i,j,bi,bj)*
231 adcroft 1.7 & (salt(i,j,km1,bi,bj)-salt(i,j-1,km1,bi,bj))
232 cnh 1.15 & +_maskS(i,j+1,km1,bi,bj)
233     & *_recip_dyC(i,j+1,bi,bj)*
234 adcroft 1.7 & (salt(i,j+1,km1,bi,bj)-salt(i,j,km1,bi,bj)))
235     & )
236     ENDDO
237     ENDDO
238    
239 cnh 1.1 C-- Vertical flux (fVerS) above
240     C Advective component of vertical flux
241 adcroft 1.7 C Note: For K=1 then KM1=1 this gives a barZ(T) = T
242     C (this plays the role of the free-surface correction)
243 cnh 1.1 DO j=jMin,jMax
244     DO i=iMin,iMax
245     af(i,j) =
246 cnh 1.12 & rTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
247 cnh 1.1 ENDDO
248     ENDDO
249     C Diffusive component of vertical flux
250 adcroft 1.7 C Note: For K=1 then KM1=1 this gives a dS/dz = 0 upper
251     C boundary condition.
252 cnh 1.1 DO j=jMin,jMax
253     DO i=iMin,iMax
254 cnh 1.13 df(i,j) = _rA(i,j,bi,bj)*(
255 adcroft 1.7 & -KapGM(i,j)*K13(i,j,k)*dSdx(i,j)
256     & -KapGM(i,j)*K23(i,j,k)*dSdy(i,j)
257     & )
258 cnh 1.1 ENDDO
259     ENDDO
260 adcroft 1.7 IF (.NOT.implicitDiffusion) THEN
261     DO j=jMin,jMax
262     DO i=iMin,iMax
263 cnh 1.12 df(i,j) = df(i,j) + _rA(i,j,bi,bj)*(
264     & -KappaRS(i,j,k)*recip_drC(k)
265 cnh 1.13 & *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))*rkFac
266 adcroft 1.7 & )
267     ENDDO
268     ENDDO
269     ENDIF
270 adcroft 1.18 #ifdef ALLOW_KPP
271     IF (usingKPPmixing) THEN
272     C-- Add non local transport coefficient (ghat term) to right-hand-side
273     C The nonlocal transport term is noNrero only for scalars in unstable
274     C (convective) forcing conditions.
275     IF ( TOP_LAYER ) THEN
276     DO j=jMin,jMax
277     DO i=iMin,iMax
278     df(i,j) = df(i,j) - _rA(i,j,bi,bj) *
279     & EmPmR(i,j,bi,bj) * delZ(1) *
280     & ( KappaRS(i,j,k) * KPPghat(i,j,k,bi,bj) )
281     ENDDO
282     ENDDO
283     ELSE
284     DO j=jMin,jMax
285     DO i=iMin,iMax
286     df(i,j) = df(i,j) - _rA(i,j,bi,bj) *
287     & EmPmR(i,j,bi,bj) * delZ(1) *
288     & ( KappaRS(i,j,k) * KPPghat(i,j,k,bi,bj)
289     & - KappaRS(i,j,k-1) * KPPghat(i,j,k-1,bi,bj) )
290     ENDDO
291     ENDDO
292     ENDIF
293     ENDIF
294     #endif /* ALLOW_KPP */
295    
296 cnh 1.1 C Net vertical flux
297     DO j=jMin,jMax
298     DO i=iMin,iMax
299 adcroft 1.7 fVerS(i,j,kUp) = ( afFacS*af(i,j)+ dfFacS*df(i,j) )*maskUp(i,j)
300 cnh 1.1 ENDDO
301     ENDDO
302 adcroft 1.7 IF ( TOP_LAYER ) THEN
303     DO j=jMin,jMax
304     DO i=iMin,iMax
305     fVerS(i,j,kUp) = afFacS*af(i,j)*freeSurfFac
306     ENDDO
307     ENDDO
308     ENDIF
309 cnh 1.1
310     C-- Tendency is minus divergence of the fluxes.
311     C Note. Tendency terms will only be correct for range
312     C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
313     C will contain valid floating point numbers but
314     C they are not algorithmically correct. These points
315     C are not used.
316     DO j=jMin,jMax
317     DO i=iMin,iMax
318 cnh 1.15 #define _recip_VolS1(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
319     #define _recip_VolS2(i,j,k,bi,bj) /_rA(i,j,bi,bj)
320 cnh 1.1 gS(i,j,k,bi,bj)=
321 cnh 1.15 & -_recip_VolS1(i,j,k,bi,bj)
322     & _recip_VolS2(i,j,k,bi,bj)
323 cnh 1.1 & *(
324     & +( fZon(i+1,j)-fZon(i,j) )
325     & +( fMer(i,j+1)-fMer(i,j) )
326 cnh 1.12 & +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )*rkFac
327 cnh 1.1 & )
328     ENDDO
329     ENDDO
330    
331 cnh 1.16 C-- External forcing term(s)
332     CALL EXTERNAL_FORCING_S(
333     I iMin,iMax,jMin,jMax,bi,bj,k,
334 cnh 1.17 I maskC,
335 cnh 1.16 I myCurrentTime,myThid)
336    
337 cnh 1.17 #ifdef INCLUDE_LAT_CIRC_FFT_FILTER_CODE
338 cnh 1.16 C--
339     CALL FILTER_LATCIRCS_FFT_APPLY( gS, 1, sNy, k, k, bi, bj, 1, myThid)
340     #endif
341 cnh 1.1
342     RETURN
343     END

  ViewVC Help
Powered by ViewVC 1.1.22