/[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.21 - (hide annotations) (download)
Wed Jun 21 19:15:26 2000 UTC (23 years, 11 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint29, checkpoint30
Changes since 1.20: +43 -86 lines
Moved all terms associated with KPP and GM/Redi to separate routines
in respective packages. Terms are included through calls to these
routines.

1 adcroft 1.21 C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gs.F,v 1.20 2000/06/09 14:26:30 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 adcroft 1.21 I KappaRS,
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 cnh 1.1
47     C == Routine arguments ==
48     C fZon - Work array for flux of temperature in the east-west
49     C direction at the west face of a cell.
50     C fMer - Work array for flux of temperature in the north-south
51     C direction at the south face of a cell.
52 adcroft 1.7 C fVerS - Flux of salt (S) in the vertical
53 cnh 1.1 C direction at the upper(U) and lower(D) faces of a cell.
54     C maskUp - Land mask used to denote base of the domain.
55 adcroft 1.10 C maskC - Land mask for salt cells (used in TOP_LAYER only)
56 cnh 1.1 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 adcroft 1.18 C rTrans - Vertical volume transport through cell face
61 cnh 1.1 C af - Advective flux component work array
62     C df - Diffusive flux component work array
63     C bi, bj, iMin, iMax, jMin, jMax - Range of points for which calculation
64     C results will be set.
65 adcroft 1.7 C myThid - Instance number for this innvocation of CALC_GT
66 cnh 1.1 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
69     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72     _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 cnh 1.12 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 cnh 1.1 _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 adcroft 1.10 _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 cnh 1.14 _RL KappaRS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
77 cnh 1.1 _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78     _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 adcroft 1.7 INTEGER k,kUp,kDown,kM1
80 cnh 1.1 INTEGER bi,bj,iMin,iMax,jMin,jMax
81 adcroft 1.18 _RL myCurrentTime
82 cnh 1.1 INTEGER myThid
83     CEndOfInterface
84    
85     C == Local variables ==
86     C I, J, K - Loop counters
87 adcroft 1.7 INTEGER i,j
88     LOGICAL TOP_LAYER
89     _RL afFacS, dfFacS
90     _RL dSdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91     _RL dSdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 adcroft 1.19 _RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 heimbach 1.20
94     #ifdef ALLOW_AUTODIFF_TAMC
95     C-- only the kUp part of fverS is set in this subroutine
96     C-- the kDown is still required
97    
98     fVerS(1,1,kDown) = fVerS(1,1,kDown)
99     DO j=1-OLy,sNy+OLy
100     DO i=1-OLx,sNx+OLx
101     fZon(i,j) = 0.0
102     fMer(i,j) = 0.0
103     fVerS(i,j,kUp) = 0.0
104     ENDDO
105     ENDDO
106     #endif
107 cnh 1.1
108     afFacS = 1. _d 0
109     dfFacS = 1. _d 0
110 adcroft 1.7 TOP_LAYER = K .EQ. 1
111 cnh 1.1
112     C--- Calculate advective and diffusive fluxes between cells.
113    
114 adcroft 1.19 #ifdef INCLUDE_T_DIFFUSION_CODE
115     C o Zonal tracer gradient
116     DO j=1-Oly,sNy+Oly
117     DO i=1-Olx+1,sNx+Olx
118     dSdx(i,j) = _recip_dxC(i,j,bi,bj)*
119     & (salt(i,j,k,bi,bj)-salt(i-1,j,k,bi,bj))
120     ENDDO
121     ENDDO
122     C o Meridional tracer gradient
123     DO j=1-Oly+1,sNy+Oly
124     DO i=1-Olx,sNx+Olx
125     dSdy(i,j) = _recip_dyC(i,j,bi,bj)*
126     & (salt(i,j,k,bi,bj)-salt(i,j-1,k,bi,bj))
127     ENDDO
128     ENDDO
129    
130     C-- del^2 of S, needed for bi-harmonic (del^4) term
131     IF (diffK4S .NE. 0.) THEN
132     DO j=1-Oly+1,sNy+Oly-1
133     DO i=1-Olx+1,sNx+Olx-1
134     df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
135     & *recip_drF(k)/_rA(i,j,bi,bj)
136     & *(
137     & +( xA(i+1,j)*dSdx(i+1,j)-xA(i,j)*dSdx(i,j) )
138     & +( yA(i,j+1)*dSdy(i,j+1)-yA(i,j)*dSdy(i,j) )
139     & )
140     ENDDO
141     ENDDO
142     ENDIF
143     #endif
144    
145 cnh 1.1 C-- Zonal flux (fZon is at west face of "salt" cell)
146     C Advective component of zonal flux
147     DO j=jMin,jMax
148     DO i=iMin,iMax
149     af(i,j) =
150     & uTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj))*0.5 _d 0
151     ENDDO
152     ENDDO
153 adcroft 1.19 C o Diffusive component of zonal flux
154 cnh 1.1 DO j=jMin,jMax
155     DO i=iMin,iMax
156 adcroft 1.21 df(i,j) = -diffKhS*xA(i,j)*dSdx(i,j)
157 cnh 1.1 ENDDO
158     ENDDO
159 adcroft 1.21 #ifdef ALLOW_GMREDI
160     IF (use_GMRedi) CALL GMREDI_XTRANSPORT(
161     I iMin,iMax,jMin,jMax,bi,bj,K,
162     I xA,salt,
163     U df,
164     I myThid)
165     #endif
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.21 df(i,j) = -diffKhS*yA(i,j)*dSdy(i,j)
195 cnh 1.1 ENDDO
196     ENDDO
197 adcroft 1.21 #ifdef ALLOW_GMREDI
198     IF (use_GMRedi) CALL GMREDI_YTRANSPORT(
199     I iMin,iMax,jMin,jMax,bi,bj,K,
200     I yA,salt,
201     U df,
202     I myThid)
203     #endif
204 adcroft 1.19 C o Add the bi-harmonic contribution
205     IF (diffK4S .NE. 0.) THEN
206     DO j=jMin,jMax
207     DO i=iMin,iMax
208     df(i,j) = df(i,j) + yA(i,j)*
209     & diffK4S*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
210     ENDDO
211     ENDDO
212     ENDIF
213    
214 cnh 1.1 C Net meridional flux
215     DO j=jMin,jMax
216     DO i=iMin,iMax
217     fMer(i,j) = afFacS*af(i,j) + dfFacS*df(i,j)
218     ENDDO
219     ENDDO
220    
221     C-- Vertical flux (fVerS) above
222     C Advective component of vertical flux
223 adcroft 1.7 C Note: For K=1 then KM1=1 this gives a barZ(T) = T
224     C (this plays the role of the free-surface correction)
225 cnh 1.1 DO j=jMin,jMax
226     DO i=iMin,iMax
227     af(i,j) =
228 cnh 1.12 & rTrans(i,j)*(salt(i,j,k,bi,bj)+salt(i,j,kM1,bi,bj))*0.5 _d 0
229 cnh 1.1 ENDDO
230     ENDDO
231 adcroft 1.21 C o Diffusive component of vertical flux
232     C Note: For K=1 then KM1=1 and this gives a dS/dr = 0 upper
233 adcroft 1.7 C boundary condition.
234 adcroft 1.21 IF (implicitDiffusion) THEN
235     DO j=jMin,jMax
236     DO i=iMin,iMax
237     df(i,j) = 0.
238     ENDDO
239 cnh 1.1 ENDDO
240 adcroft 1.21 ELSE
241 adcroft 1.7 DO j=jMin,jMax
242     DO i=iMin,iMax
243 adcroft 1.21 df(i,j) = - _rA(i,j,bi,bj)*(
244     & KappaRS(i,j,k)*recip_drC(k)
245 cnh 1.13 & *(salt(i,j,kM1,bi,bj)-salt(i,j,k,bi,bj))*rkFac
246 adcroft 1.7 & )
247     ENDDO
248     ENDDO
249     ENDIF
250 adcroft 1.21
251     #ifdef ALLOW_GMREDI
252     IF (use_GMRedi) CALL GMREDI_RTRANSPORT(
253     I iMin,iMax,jMin,jMax,bi,bj,K,
254     I maskUp,salt,
255     U df,
256     I myThid)
257     #endif
258    
259 adcroft 1.18 #ifdef ALLOW_KPP
260 adcroft 1.21 C-- Add non-local KPP transport term (ghat) to diffusive salt flux.
261     IF (use_KPPmixing) CALL KPP_TRANSPORT_S(
262     I iMin,iMax,jMin,jMax,bi,bj,k,km1,
263     I maskC,KappaRS,
264     U df )
265     #endif
266 adcroft 1.18
267 cnh 1.1 C Net vertical flux
268     DO j=jMin,jMax
269     DO i=iMin,iMax
270 adcroft 1.7 fVerS(i,j,kUp) = ( afFacS*af(i,j)+ dfFacS*df(i,j) )*maskUp(i,j)
271 cnh 1.1 ENDDO
272     ENDDO
273 adcroft 1.7 IF ( TOP_LAYER ) THEN
274     DO j=jMin,jMax
275     DO i=iMin,iMax
276     fVerS(i,j,kUp) = afFacS*af(i,j)*freeSurfFac
277     ENDDO
278     ENDDO
279     ENDIF
280 cnh 1.1
281     C-- Tendency is minus divergence of the fluxes.
282     C Note. Tendency terms will only be correct for range
283     C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
284     C will contain valid floating point numbers but
285     C they are not algorithmically correct. These points
286     C are not used.
287     DO j=jMin,jMax
288     DO i=iMin,iMax
289 cnh 1.15 #define _recip_VolS1(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
290     #define _recip_VolS2(i,j,k,bi,bj) /_rA(i,j,bi,bj)
291 cnh 1.1 gS(i,j,k,bi,bj)=
292 cnh 1.15 & -_recip_VolS1(i,j,k,bi,bj)
293     & _recip_VolS2(i,j,k,bi,bj)
294 cnh 1.1 & *(
295     & +( fZon(i+1,j)-fZon(i,j) )
296     & +( fMer(i,j+1)-fMer(i,j) )
297 cnh 1.12 & +( fVerS(i,j,kUp)-fVerS(i,j,kDown) )*rkFac
298 cnh 1.1 & )
299     ENDDO
300     ENDDO
301    
302 cnh 1.16 C-- External forcing term(s)
303     CALL EXTERNAL_FORCING_S(
304     I iMin,iMax,jMin,jMax,bi,bj,k,
305 cnh 1.17 I maskC,
306 cnh 1.16 I myCurrentTime,myThid)
307    
308 cnh 1.17 #ifdef INCLUDE_LAT_CIRC_FFT_FILTER_CODE
309 cnh 1.16 C--
310     CALL FILTER_LATCIRCS_FFT_APPLY( gS, 1, sNy, k, k, bi, bj, 1, myThid)
311     #endif
312 cnh 1.1
313     RETURN
314     END

  ViewVC Help
Powered by ViewVC 1.1.22