/[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.27 - (hide annotations) (download)
Mon Nov 13 16:32:57 2000 UTC (23 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: branch-atmos-merge-start, checkpoint33, checkpoint32, checkpoint34, branch-atmos-merge-phase1
Branch point for: branch-atmos-merge
Changes since 1.26: +1 -2 lines
Rescaling of forcing fields done immediately after reading fields.

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

  ViewVC Help
Powered by ViewVC 1.1.22