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

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

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


Revision 1.24 - (show annotations) (download)
Fri Jun 9 02:45:04 2000 UTC (23 years, 11 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint28
Changes since 1.23: +15 -1 lines
Modifications to include TAMC directives, tape key computations
and initialisations to make code TAMC compatible.
Routines the_model_main.F and initialise_fixed.F
are left unchanged for the moment. (P.H.)

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

  ViewVC Help
Powered by ViewVC 1.1.22