/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F
ViewVC logotype

Contents of /MITgcm/pkg/generic_advdiff/gad_calc_rhs.F

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


Revision 1.14 - (show annotations) (download)
Sat Jun 15 03:31:17 2002 UTC (21 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint46l_post, checkpoint46g_pre, checkpoint46f_post, checkpoint46b_post, checkpoint46l_pre, checkpoint46d_pre, checkpoint45d_post, checkpoint46j_pre, checkpoint46a_post, checkpoint46j_post, checkpoint46k_post, checkpoint46e_pre, checkpoint46b_pre, checkpoint46c_pre, checkpoint46, checkpoint46h_pre, checkpoint46m_post, checkpoint46a_pre, checkpoint46g_post, checkpoint46i_post, checkpoint46c_post, checkpoint46e_post, checkpoint46h_post, checkpoint46d_post
Changes since 1.13: +8 -28 lines
for each tracer, define an internal flag for multiDimAdvection & A-B

1 C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F,v 1.13 2002/03/24 02:12:50 heimbach Exp $
2 C $Name: $
3
4 #include "GAD_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: GAD_CALC_RHS
8
9 C !INTERFACE: ==========================================================
10 SUBROUTINE GAD_CALC_RHS(
11 I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12 I xA,yA,uTrans,vTrans,rTrans,maskUp,
13 I diffKh, diffK4, KappaRT, Tracer,
14 I tracerIdentity, advectionScheme, calcAdvection,
15 U fVerT, gTracer,
16 I myThid )
17
18 C !DESCRIPTION:
19 C Calculates the tendancy of a tracer due to advection and diffusion.
20 C It calculates the fluxes in each direction indepentently and then
21 C sets the tendancy to the divergence of these fluxes. The advective
22 C fluxes are only calculated here when using the linear advection schemes
23 C otherwise only the diffusive and parameterized fluxes are calculated.
24 C
25 C Contributions to the flux are calculated and added:
26 C \begin{equation*}
27 C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
28 C \end{equation*}
29 C
30 C The tendancy is the divergence of the fluxes:
31 C \begin{equation*}
32 C G_\theta = G_\theta + \nabla \cdot {\bf F}
33 C \end{equation*}
34 C
35 C The tendancy is assumed to contain data on entry.
36
37 C !USES: ===============================================================
38 IMPLICIT NONE
39 #include "SIZE.h"
40 #include "EEPARAMS.h"
41 #include "PARAMS.h"
42 #include "GRID.h"
43 #include "DYNVARS.h"
44 #include "GAD.h"
45
46 #ifdef ALLOW_AUTODIFF_TAMC
47 #include "tamc.h"
48 #include "tamc_keys.h"
49 #endif /* ALLOW_AUTODIFF_TAMC */
50
51 C !INPUT PARAMETERS: ===================================================
52 C bi,bj :: tile indices
53 C iMin,iMax,jMin,jMax :: loop range for called routines
54 C kup :: index into 2 1/2D array, toggles between 1 and 2
55 C kdown :: index into 2 1/2D array, toggles between 2 and 1
56 C kp1 :: =k+1 for k<Nr, =Nr for k=Nr
57 C xA,yA :: areas of X and Y face of tracer cells
58 C uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points
59 C maskUp :: 2-D array for mask at W points
60 C diffKh :: horizontal diffusion coefficient
61 C diffK4 :: bi-harmonic diffusion coefficient
62 C KappaRT :: 3-D array for vertical diffusion coefficient
63 C Tracer :: tracer field
64 C tracerIdentity :: identifier for the tracer (required only for KPP)
65 C advectionScheme :: advection scheme to use
66 C calcAdvection :: =False if Advec terms computed with multiDim scheme
67 C myThid :: thread number
68 INTEGER bi,bj,iMin,iMax,jMin,jMax
69 INTEGER k,kUp,kDown,kM1
70 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 _RL diffKh, diffK4
77 _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
78 _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
79 INTEGER tracerIdentity
80 INTEGER advectionScheme
81 LOGICAL calcAdvection
82 INTEGER myThid
83
84 C !OUTPUT PARAMETERS: ==================================================
85 C gTracer :: tendancy array
86 C fVerT :: 2 1/2D arrays for vertical advective flux
87 _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
88 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
89
90 C !LOCAL VARIABLES: ====================================================
91 C i,j :: loop indices
92 C df4 :: used for storing del^2 T for bi-harmonic term
93 C fZon :: zonal flux
94 C fmer :: meridional flux
95 C af :: advective flux
96 C df :: diffusive flux
97 C localT :: local copy of tracer field
98 INTEGER i,j
99 _RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102 _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103 _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105 CEOP
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 fVerT(1,1,kDown) = fVerT(1,1,kDown)
111 #endif
112
113 DO j=1-OLy,sNy+OLy
114 DO i=1-OLx,sNx+OLx
115 fZon(i,j) = 0. _d 0
116 fMer(i,j) = 0. _d 0
117 fVerT(i,j,kUp) = 0. _d 0
118 df(i,j) = 0. _d 0
119 df4(i,j) = 0. _d 0
120 localT(i,j) = 0. _d 0
121 ENDDO
122 ENDDO
123
124 C-- Make local copy of tracer array
125 DO j=1-OLy,sNy+OLy
126 DO i=1-OLx,sNx+OLx
127 localT(i,j)=tracer(i,j,k,bi,bj)
128 ENDDO
129 ENDDO
130
131 C-- Unless we have already calculated the advection terms we initialize
132 C the tendency to zero.
133 IF (calcAdvection) THEN
134 DO j=1-Oly,sNy+Oly
135 DO i=1-Olx,sNx+Olx
136 gTracer(i,j,k,bi,bj)=0. _d 0
137 ENDDO
138 ENDDO
139 ENDIF
140
141 C-- Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
142 IF (diffK4 .NE. 0.) THEN
143 CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
144 CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
145 CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
146 ENDIF
147
148 C-- Initialize net flux in X direction
149 DO j=1-Oly,sNy+Oly
150 DO i=1-Olx,sNx+Olx
151 fZon(i,j) = 0. _d 0
152 ENDDO
153 ENDDO
154
155 C- Advective flux in X
156 IF (calcAdvection) THEN
157 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
158 CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
159 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
160 CALL GAD_FLUXLIMIT_ADV_X(
161 & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
162 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
163 CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
164 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
165 CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
166 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
167 CALL GAD_DST3_ADV_X(
168 & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
169 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
170 CALL GAD_DST3FL_ADV_X(
171 & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
172 ELSE
173 STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
174 ENDIF
175 DO j=1-Oly,sNy+Oly
176 DO i=1-Olx,sNx+Olx
177 fZon(i,j) = fZon(i,j) + af(i,j)
178 ENDDO
179 ENDDO
180 ENDIF
181
182 C- Diffusive flux in X
183 IF (diffKh.NE.0.) THEN
184 CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
185 ELSE
186 DO j=1-Oly,sNy+Oly
187 DO i=1-Olx,sNx+Olx
188 df(i,j) = 0. _d 0
189 ENDDO
190 ENDDO
191 ENDIF
192
193 #ifdef ALLOW_GMREDI
194 C- GM/Redi flux in X
195 IF (useGMRedi) THEN
196 C *note* should update GMREDI_XTRANSPORT to use localT and set df *aja*
197 CALL GMREDI_XTRANSPORT(
198 I iMin,iMax,jMin,jMax,bi,bj,K,
199 I xA,Tracer,
200 U df,
201 I myThid)
202 ENDIF
203 #endif
204 DO j=1-Oly,sNy+Oly
205 DO i=1-Olx,sNx+Olx
206 fZon(i,j) = fZon(i,j) + df(i,j)
207 ENDDO
208 ENDDO
209
210 C- Bi-harmonic duffusive flux in X
211 IF (diffK4 .NE. 0.) THEN
212 CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
213 DO j=1-Oly,sNy+Oly
214 DO i=1-Olx,sNx+Olx
215 fZon(i,j) = fZon(i,j) + df(i,j)
216 ENDDO
217 ENDDO
218 ENDIF
219
220 C-- Initialize net flux in Y direction
221 DO j=1-Oly,sNy+Oly
222 DO i=1-Olx,sNx+Olx
223 fMer(i,j) = 0. _d 0
224 ENDDO
225 ENDDO
226
227 C- Advective flux in Y
228 IF (calcAdvection) THEN
229 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
230 CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
231 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
232 CALL GAD_FLUXLIMIT_ADV_Y(
233 & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
234 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
235 CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
236 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
237 CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
238 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
239 CALL GAD_DST3_ADV_Y(
240 & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
241 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
242 CALL GAD_DST3FL_ADV_Y(
243 & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
244 ELSE
245 STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
246 ENDIF
247 DO j=1-Oly,sNy+Oly
248 DO i=1-Olx,sNx+Olx
249 fMer(i,j) = fMer(i,j) + af(i,j)
250 ENDDO
251 ENDDO
252 ENDIF
253
254 C- Diffusive flux in Y
255 IF (diffKh.NE.0.) THEN
256 CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
257 ELSE
258 DO j=1-Oly,sNy+Oly
259 DO i=1-Olx,sNx+Olx
260 df(i,j) = 0. _d 0
261 ENDDO
262 ENDDO
263 ENDIF
264
265 #ifdef ALLOW_GMREDI
266 C- GM/Redi flux in Y
267 IF (useGMRedi) THEN
268 C *note* should update GMREDI_YTRANSPORT to use localT and set df *aja*
269 CALL GMREDI_YTRANSPORT(
270 I iMin,iMax,jMin,jMax,bi,bj,K,
271 I yA,Tracer,
272 U df,
273 I myThid)
274 ENDIF
275 #endif
276 DO j=1-Oly,sNy+Oly
277 DO i=1-Olx,sNx+Olx
278 fMer(i,j) = fMer(i,j) + df(i,j)
279 ENDDO
280 ENDDO
281
282 C- Bi-harmonic flux in Y
283 IF (diffK4 .NE. 0.) THEN
284 CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
285 DO j=1-Oly,sNy+Oly
286 DO i=1-Olx,sNx+Olx
287 fMer(i,j) = fMer(i,j) + df(i,j)
288 ENDDO
289 ENDDO
290 ENDIF
291
292 C- Advective flux in R
293 IF (calcAdvection) THEN
294 C Note: wVel needs to be masked
295 IF (K.GE.2) THEN
296 C- Compute vertical advective flux in the interior:
297 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
298 CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
299 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
300 CALL GAD_FLUXLIMIT_ADV_R(
301 & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
302 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
303 CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
304 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
305 CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
306 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
307 CALL GAD_DST3_ADV_R(
308 & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
309 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
310 CALL GAD_DST3FL_ADV_R(
311 & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
312 ELSE
313 STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
314 ENDIF
315 C- Surface "correction" term at k>1 :
316 DO j=1-Oly,sNy+Oly
317 DO i=1-Olx,sNx+Olx
318 af(i,j) = af(i,j)
319 & + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*
320 & rTrans(i,j)*Tracer(i,j,k,bi,bj)
321 ENDDO
322 ENDDO
323 ELSE
324 C- Surface "correction" term at k=1 :
325 DO j=1-Oly,sNy+Oly
326 DO i=1-Olx,sNx+Olx
327 af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)
328 ENDDO
329 ENDDO
330 ENDIF
331 C- add the advective flux to fVerT
332 DO j=1-Oly,sNy+Oly
333 DO i=1-Olx,sNx+Olx
334 fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
335 ENDDO
336 ENDDO
337 ENDIF
338
339 C- Diffusive flux in R
340 C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
341 C boundary condition.
342 IF (implicitDiffusion) THEN
343 DO j=1-Oly,sNy+Oly
344 DO i=1-Olx,sNx+Olx
345 df(i,j) = 0. _d 0
346 ENDDO
347 ENDDO
348 ELSE
349 CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
350 ENDIF
351
352 #ifdef ALLOW_GMREDI
353 C- GM/Redi flux in R
354 IF (useGMRedi) THEN
355 C *note* should update GMREDI_RTRANSPORT to set df *aja*
356 CALL GMREDI_RTRANSPORT(
357 I iMin,iMax,jMin,jMax,bi,bj,K,
358 I Tracer,
359 U df,
360 I myThid)
361 ENDIF
362 #endif
363
364 DO j=1-Oly,sNy+Oly
365 DO i=1-Olx,sNx+Olx
366 fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
367 ENDDO
368 ENDDO
369
370 #ifdef ALLOW_KPP
371 C- Add non local KPP transport term (ghat) to diffusive T flux.
372 IF (useKPP) THEN
373 DO j=1-Oly,sNy+Oly
374 DO i=1-Olx,sNx+Olx
375 df(i,j) = 0. _d 0
376 ENDDO
377 ENDDO
378 IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
379 C *note* should update KPP_TRANSPORT_T to set df *aja*
380 CALL KPP_TRANSPORT_T(
381 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
382 I KappaRT,
383 U df )
384 ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
385 CALL KPP_TRANSPORT_S(
386 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
387 I KappaRT,
388 U df )
389 ELSE
390 STOP 'GAD_CALC_RHS: Ooops'
391 ENDIF
392 DO j=1-Oly,sNy+Oly
393 DO i=1-Olx,sNx+Olx
394 fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
395 ENDDO
396 ENDDO
397 ENDIF
398 #endif
399
400 C-- Divergence of fluxes
401 DO j=1-Oly,sNy+Oly-1
402 DO i=1-Olx,sNx+Olx-1
403 gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
404 & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
405 & *recip_rA(i,j,bi,bj)
406 & *(
407 & +( fZon(i+1,j)-fZon(i,j) )
408 & +( fMer(i,j+1)-fMer(i,j) )
409 & +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
410 & )
411 ENDDO
412 ENDDO
413
414 RETURN
415 END

  ViewVC Help
Powered by ViewVC 1.1.22