/[MITgcm]/MITgcm/model/src/calc_gt.F
ViewVC logotype

Annotation of /MITgcm/model/src/calc_gt.F

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


Revision 1.29 - (hide annotations) (download)
Sun Feb 4 14:38:46 2001 UTC (23 years, 5 months ago) by cnh
Branch: MAIN
Changes since 1.28: +2 -1 lines
Made sure each .F and .h file had
the CVS keywords Header and Name at its start.
Most had header but very few currently have Name, so
lots of changes!

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

  ViewVC Help
Powered by ViewVC 1.1.22