/[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.31 - (hide annotations) (download)
Wed Feb 7 16:27:29 2001 UTC (23 years, 4 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint38, c37_adj, checkpoint39, checkpoint37, checkpoint36, checkpoint35
Branch point for: pre38
Changes since 1.30: +1 -12 lines
Applied same fix for DEC in both calc_gs and calc_gt. Keeping these
to routines close together will make it easier to replace them
with a generic routine later.

1 adcroft 1.31 C $Header: /u/gcmpack/models/MITgcmUV/model/src/calc_gt.F,v 1.30 2001/02/06 03:08:59 cnh Exp $
2 cnh 1.30 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     fVerT(1,1,kDown) = fVerT(1,1,kDown)
95 adcroft 1.20 #endif
96 cnh 1.30 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 cnh 1.1
104     afFacT = 1. _d 0
105     dfFacT = 1. _d 0
106 cnh 1.10 TOP_LAYER = K .EQ. 1
107 cnh 1.1
108     C--- Calculate advective and diffusive fluxes between cells.
109    
110 adcroft 1.22 #ifdef INCLUDE_T_DIFFUSION_CODE
111     C o Zonal tracer gradient
112     DO j=1-Oly,sNy+Oly
113     DO i=1-Olx+1,sNx+Olx
114     dTdx(i,j) = _recip_dxC(i,j,bi,bj)*
115     & (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj))
116     ENDDO
117     ENDDO
118     C o Meridional tracer gradient
119     DO j=1-Oly+1,sNy+Oly
120     DO i=1-Olx,sNx+Olx
121     dTdy(i,j) = _recip_dyC(i,j,bi,bj)*
122     & (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
123     ENDDO
124     ENDDO
125    
126     C-- del^2 of T, needed for bi-harmonic (del^4) term
127     IF (diffK4T .NE. 0.) THEN
128     DO j=1-Oly+1,sNy+Oly-1
129     DO i=1-Olx+1,sNx+Olx-1
130     df4(i,j)= _recip_hFacC(i,j,k,bi,bj)
131     & *recip_drF(k)/_rA(i,j,bi,bj)
132     & *(
133     & +( xA(i+1,j)*dTdx(i+1,j)-xA(i,j)*dTdx(i,j) )
134     & +( yA(i,j+1)*dTdy(i,j+1)-yA(i,j)*dTdy(i,j) )
135     & )
136     ENDDO
137     ENDDO
138     ENDIF
139     #endif
140    
141 cnh 1.1 C-- Zonal flux (fZon is at west face of "theta" cell)
142 cnh 1.19 #ifdef INCLUDE_T_ADVECTION_CODE
143     C o Advective component of zonal flux
144 cnh 1.1 DO j=jMin,jMax
145     DO i=iMin,iMax
146     af(i,j) =
147     & uTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj))*0.5 _d 0
148     ENDDO
149     ENDDO
150 cnh 1.19 #endif /* INCLUDE_T_ADVECTION_CODE */
151     #ifdef INCLUDE_T_DIFFUSION_CODE
152     C o Diffusive component of zonal flux
153 cnh 1.1 DO j=jMin,jMax
154     DO i=iMin,iMax
155 adcroft 1.25 df(i,j) = -diffKhT*xA(i,j)*dTdx(i,j)
156 cnh 1.1 ENDDO
157     ENDDO
158 adcroft 1.25 #ifdef ALLOW_GMREDI
159 heimbach 1.26 IF (useGMRedi) CALL GMREDI_XTRANSPORT(
160 adcroft 1.25 I iMin,iMax,jMin,jMax,bi,bj,K,
161     I xA,theta,
162     U df,
163     I myThid)
164     #endif
165 adcroft 1.22 C o Add the bi-harmonic contribution
166     IF (diffK4T .NE. 0.) THEN
167     DO j=jMin,jMax
168     DO i=iMin,iMax
169     df(i,j) = df(i,j) + xA(i,j)*
170     & diffK4T*(df4(i,j)-df4(i-1,j))*_recip_dxC(i,j,bi,bj)
171     ENDDO
172     ENDDO
173     ENDIF
174 cnh 1.19 #endif /* INCLUDE_T_DIFFUSION_CODE */
175     C o Net zonal flux
176 cnh 1.1 DO j=jMin,jMax
177     DO i=iMin,iMax
178 cnh 1.19 fZon(i,j) = 0.
179 adcroft 1.23 & _ADT( + afFacT*af(i,j) )
180     & _LPT( + dfFacT*df(i,j) )
181 cnh 1.1 ENDDO
182     ENDDO
183    
184     C-- Meridional flux (fMer is at south face of "theta" cell)
185 cnh 1.19 #ifdef INCLUDE_T_ADVECTION_CODE
186     C o Advective component of meridional flux
187 cnh 1.1 DO j=jMin,jMax
188     DO i=iMin,iMax
189     af(i,j) =
190     & vTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))*0.5 _d 0
191     ENDDO
192     ENDDO
193 cnh 1.19 #endif /* INCLUDE_T_ADVECTION_CODE */
194     #ifdef INCLUDE_T_DIFFUSION_CODE
195     C o Diffusive component of meridional flux
196 cnh 1.1 DO j=jMin,jMax
197     DO i=iMin,iMax
198 adcroft 1.25 df(i,j) = -diffKhT*yA(i,j)*dTdy(i,j)
199 cnh 1.1 ENDDO
200     ENDDO
201 adcroft 1.25 #ifdef ALLOW_GMREDI
202 heimbach 1.26 IF (useGMRedi) CALL GMREDI_YTRANSPORT(
203 adcroft 1.25 I iMin,iMax,jMin,jMax,bi,bj,K,
204     I yA,theta,
205     U df,
206     I myThid)
207     #endif
208 adcroft 1.22 C o Add the bi-harmonic contribution
209     IF (diffK4T .NE. 0.) THEN
210     DO j=jMin,jMax
211     DO i=iMin,iMax
212     df(i,j) = df(i,j) + yA(i,j)*
213     & diffK4T*(df4(i,j)-df4(i,j-1))*_recip_dyC(i,j,bi,bj)
214     ENDDO
215     ENDDO
216     ENDIF
217 cnh 1.19 #endif /* INCLUDE_T_DIFFUSION_CODE */
218     C o Net meridional flux
219 cnh 1.1 DO j=jMin,jMax
220     DO i=iMin,iMax
221 cnh 1.19 fMer(i,j) = 0.
222 adcroft 1.23 & _ADT( + afFacT*af(i,j) )
223     & _LPT( + dfFacT*df(i,j) )
224 cnh 1.1 ENDDO
225     ENDDO
226    
227 cnh 1.19 #ifdef INCLUDE_T_DIFFUSION_CODE
228     C-- Terms that diffusion tensor projects onto z
229 adcroft 1.3 DO j=jMin,jMax
230     DO i=iMin,iMax
231     dTdx(i,j) = 0.5*(
232 cnh 1.17 & +0.5*(_maskW(i+1,j,k,bi,bj)
233     & *_recip_dxC(i+1,j,bi,bj)*
234 adcroft 1.3 & (theta(i+1,j,k,bi,bj)-theta(i,j,k,bi,bj))
235 cnh 1.17 & +_maskW(i,j,k,bi,bj)
236     & *_recip_dxC(i,j,bi,bj)*
237 adcroft 1.3 & (theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj)))
238 cnh 1.17 & +0.5*(_maskW(i+1,j,km1,bi,bj)
239     & *_recip_dxC(i+1,j,bi,bj)*
240 adcroft 1.3 & (theta(i+1,j,km1,bi,bj)-theta(i,j,km1,bi,bj))
241 cnh 1.17 & +_maskW(i,j,km1,bi,bj)
242     & *_recip_dxC(i,j,bi,bj)*
243 adcroft 1.3 & (theta(i,j,km1,bi,bj)-theta(i-1,j,km1,bi,bj)))
244     & )
245     ENDDO
246     ENDDO
247     DO j=jMin,jMax
248     DO i=iMin,iMax
249     dTdy(i,j) = 0.5*(
250 cnh 1.17 & +0.5*(_maskS(i,j,k,bi,bj)
251     & *_recip_dyC(i,j,bi,bj)*
252 adcroft 1.3 & (theta(i,j,k,bi,bj)-theta(i,j-1,k,bi,bj))
253 cnh 1.17 & +_maskS(i,j+1,k,bi,bj)
254     & *_recip_dyC(i,j+1,bi,bj)*
255 adcroft 1.3 & (theta(i,j+1,k,bi,bj)-theta(i,j,k,bi,bj)))
256 cnh 1.17 & +0.5*(_maskS(i,j,km1,bi,bj)
257     & *_recip_dyC(i,j,bi,bj)*
258 adcroft 1.3 & (theta(i,j,km1,bi,bj)-theta(i,j-1,km1,bi,bj))
259 cnh 1.17 & +_maskS(i,j+1,km1,bi,bj)
260     & *_recip_dyC(i,j+1,bi,bj)*
261 adcroft 1.3 & (theta(i,j+1,km1,bi,bj)-theta(i,j,km1,bi,bj)))
262     & )
263     ENDDO
264     ENDDO
265 cnh 1.19 #endif /* INCLUDE_T_DIFFUSION_CODE */
266 adcroft 1.3
267 cnh 1.19 C-- Vertical flux ( fVerT(,,kUp) is at upper face of "theta" cell )
268     #ifdef INCLUDE_T_ADVECTION_CODE
269     C o Advective component of vertical flux
270 adcroft 1.3 C Note: For K=1 then KM1=1 this gives a barZ(T) = T
271     C (this plays the role of the free-surface correction)
272 cnh 1.1 DO j=jMin,jMax
273     DO i=iMin,iMax
274     af(i,j) =
275 cnh 1.14 & rTrans(i,j)*(theta(i,j,k,bi,bj)+theta(i,j,kM1,bi,bj))*0.5 _d 0
276 cnh 1.1 ENDDO
277     ENDDO
278 cnh 1.19 #endif /* INCLUDE_T_ADVECTION_CODE */
279     #ifdef INCLUDE_T_DIFFUSION_CODE
280     C o Diffusive component of vertical flux
281     C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
282 adcroft 1.3 C boundary condition.
283 adcroft 1.25 IF (implicitDiffusion) THEN
284     DO j=jMin,jMax
285     DO i=iMin,iMax
286     df(i,j) = 0.
287     ENDDO
288 cnh 1.1 ENDDO
289 adcroft 1.25 ELSE
290 adcroft 1.9 DO j=jMin,jMax
291     DO i=iMin,iMax
292 adcroft 1.25 df(i,j) = - _rA(i,j,bi,bj)*(
293     & KappaRT(i,j,k)*recip_drC(k)
294 cnh 1.15 & *(theta(i,j,kM1,bi,bj)-theta(i,j,k,bi,bj))*rkFac
295 adcroft 1.9 & )
296     ENDDO
297     ENDDO
298     ENDIF
299 cnh 1.19 #endif /* INCLUDE_T_DIFFUSION_CODE */
300 adcroft 1.20
301 adcroft 1.25 #ifdef ALLOW_GMREDI
302 heimbach 1.26 IF (useGMRedi) CALL GMREDI_RTRANSPORT(
303 adcroft 1.25 I iMin,iMax,jMin,jMax,bi,bj,K,
304     I maskUp,theta,
305     U df,
306     I myThid)
307     #endif
308    
309 adcroft 1.20 #ifdef ALLOW_KPP
310 adcroft 1.25 C-- Add non local KPP transport term (ghat) to diffusive T flux.
311 heimbach 1.26 IF (useKPP) CALL KPP_TRANSPORT_T(
312 adcroft 1.25 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
313     I maskC,KappaRT,
314     U df )
315     #endif
316 adcroft 1.20
317 cnh 1.19 C o Net vertical flux
318 cnh 1.1 DO j=jMin,jMax
319     DO i=iMin,iMax
320 cnh 1.19 fVerT(i,j,kUp) = 0.
321 adcroft 1.23 & _ADT( +afFacT*af(i,j)*maskUp(i,j) )
322     & _LPT( +dfFacT*df(i,j)*maskUp(i,j) )
323 cnh 1.1 ENDDO
324     ENDDO
325 cnh 1.19 #ifdef INCLUDE_T_ADVECTION_CODE
326 cnh 1.10 IF ( TOP_LAYER ) THEN
327     DO j=jMin,jMax
328     DO i=iMin,iMax
329     fVerT(i,j,kUp) = afFacT*af(i,j)*freeSurfFac
330     ENDDO
331     ENDDO
332     ENDIF
333 cnh 1.19 #endif /* INCLUDE_T_ADVECTION_CODE */
334 cnh 1.1
335     C-- Tendency is minus divergence of the fluxes.
336     C Note. Tendency terms will only be correct for range
337     C i=iMin+1:iMax-1, j=jMin+1:jMax-1. Edge points
338     C will contain valid floating point numbers but
339     C they are not algorithmically correct. These points
340     C are not used.
341     DO j=jMin,jMax
342     DO i=iMin,iMax
343 cnh 1.17 #define _recip_VolT1(i,j,k,bi,bj) _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
344     #define _recip_VolT2(i,j,k,bi,bj) /_rA(i,j,bi,bj)
345 cnh 1.1 gT(i,j,k,bi,bj)=
346 cnh 1.17 & -_recip_VolT1(i,j,k,bi,bj)
347     & _recip_VolT2(i,j,k,bi,bj)
348 cnh 1.1 & *(
349     & +( fZon(i+1,j)-fZon(i,j) )
350     & +( fMer(i,j+1)-fMer(i,j) )
351 cnh 1.14 & +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
352 cnh 1.1 & )
353     ENDDO
354     ENDDO
355    
356 cnh 1.19 #ifdef INCLUDE_T_FORCING_CODE
357 cnh 1.1 C-- External thermal forcing term(s)
358 cnh 1.19 CALL EXTERNAL_FORCING_T(
359     I iMin,iMax,jMin,jMax,bi,bj,k,
360     I maskC,
361     I myCurrentTime,myThid)
362     #endif /* INCLUDE_T_FORCING_CODE */
363 cnh 1.1
364     RETURN
365     END

  ViewVC Help
Powered by ViewVC 1.1.22