/[MITgcm]/MITgcm/verification/global_with_CFC11/code1x1/gad_calc_rhs.F
ViewVC logotype

Contents of /MITgcm/verification/global_with_CFC11/code1x1/gad_calc_rhs.F

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


Revision 1.1.2.1 - (show annotations) (download)
Thu Aug 25 16:22:17 2005 UTC (18 years, 7 months ago) by dimitri
Branch: release1_50yr
Changes since 1.1: +440 -0 lines
adding ecco1x1 verification/global_with_CFC11 experiment

1 C $Header: /u/gcmpack/MITgcm/verification/global_with_CFC11/code50yr/Attic/gad_calc_rhs.F,v 1.1.2.1 2003/05/04 23:19:18 dimitri Exp $
2 C $Name: release1_50yr $
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,
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 myThid :: thread number
67 INTEGER bi,bj,iMin,iMax,jMin,jMax
68 INTEGER k,kUp,kDown,kM1
69 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74 _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 _RL diffKh, diffK4
76 _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
77 _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
78 INTEGER tracerIdentity
79 INTEGER advectionScheme
80 INTEGER myThid
81
82 C !OUTPUT PARAMETERS: ==================================================
83 C gTracer :: tendancy array
84 C fVerT :: 2 1/2D arrays for vertical advective flux
85 _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
86 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
87
88 C !LOCAL VARIABLES: ====================================================
89 C i,j :: loop indices
90 C df4 :: used for storing del^2 T for bi-harmonic term
91 C fZon :: zonal flux
92 C fmer :: meridional flux
93 C af :: advective flux
94 C df :: diffusive flux
95 C localT :: local copy of tracer field
96 INTEGER i,j
97 _RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101 _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102 _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103 CEOP
104
105 #ifdef ALLOW_AUTODIFF_TAMC
106 C-- only the kUp part of fverT is set in this subroutine
107 C-- the kDown is still required
108 fVerT(1,1,kDown) = fVerT(1,1,kDown)
109 #endif
110
111 DO j=1-OLy,sNy+OLy
112 DO i=1-OLx,sNx+OLx
113 fZon(i,j) = 0. _d 0
114 fMer(i,j) = 0. _d 0
115 fVerT(i,j,kUp) = 0. _d 0
116 df(i,j) = 0. _d 0
117 df4(i,j) = 0. _d 0
118 localT(i,j) = 0. _d 0
119 ENDDO
120 ENDDO
121
122 C-- Make local copy of tracer array
123 DO j=1-OLy,sNy+OLy
124 DO i=1-OLx,sNx+OLx
125 localT(i,j)=tracer(i,j,k,bi,bj)
126 ENDDO
127 ENDDO
128
129 C-- Unless we have already calculated the advection terms we initialize
130 C the tendency to zero.
131 IF (.NOT. multiDimAdvection .OR.
132 & advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
133 & advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
134 & advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
135 DO j=1-Oly,sNy+Oly
136 DO i=1-Olx,sNx+Olx
137 gTracer(i,j,k,bi,bj)=0. _d 0
138 ENDDO
139 ENDDO
140 ENDIF
141
142 C-- Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
143 IF (diffK4 .NE. 0.) THEN
144 CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
145 CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
146 CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
147 ENDIF
148
149 C-- Initialize net flux in X direction
150 DO j=1-Oly,sNy+Oly
151 DO i=1-Olx,sNx+Olx
152 fZon(i,j) = 0. _d 0
153 ENDDO
154 ENDDO
155
156 C- Advective flux in X
157 IF (.NOT. multiDimAdvection .OR.
158 & advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
159 & advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
160 & advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
161 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
162 CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
163 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
164 CALL GAD_FLUXLIMIT_ADV_X(
165 & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
166 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
167 CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
168 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
169 CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
170 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
171 CALL GAD_DST3_ADV_X(
172 & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
173 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
174 CALL GAD_DST3FL_ADV_X(
175 & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
176 ELSE
177 STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
178 ENDIF
179 DO j=1-Oly,sNy+Oly
180 DO i=1-Olx,sNx+Olx
181 fZon(i,j) = fZon(i,j) + af(i,j)
182 ENDDO
183 ENDDO
184 ENDIF
185
186 C- Diffusive flux in X
187 IF (diffKh.NE.0.) THEN
188 CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
189 ELSE
190 DO j=1-Oly,sNy+Oly
191 DO i=1-Olx,sNx+Olx
192 df(i,j) = 0. _d 0
193 ENDDO
194 ENDDO
195 ENDIF
196
197 #ifdef ALLOW_GMREDI
198 C- GM/Redi flux in X
199 IF (useGMRedi) THEN
200 C *note* should update GMREDI_XTRANSPORT to use localT and set df *aja*
201 CALL GMREDI_XTRANSPORT(
202 I iMin,iMax,jMin,jMax,bi,bj,K,
203 I xA,Tracer,tracerIdentity,
204 U df,
205 I myThid)
206 ENDIF
207 #endif
208 DO j=1-Oly,sNy+Oly
209 DO i=1-Olx,sNx+Olx
210 fZon(i,j) = fZon(i,j) + df(i,j)
211 ENDDO
212 ENDDO
213
214 C- Bi-harmonic duffusive flux in X
215 IF (diffK4 .NE. 0.) THEN
216 CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
217 DO j=1-Oly,sNy+Oly
218 DO i=1-Olx,sNx+Olx
219 fZon(i,j) = fZon(i,j) + df(i,j)
220 ENDDO
221 ENDDO
222 ENDIF
223
224 C-- Initialize net flux in Y direction
225 DO j=1-Oly,sNy+Oly
226 DO i=1-Olx,sNx+Olx
227 fMer(i,j) = 0. _d 0
228 ENDDO
229 ENDDO
230
231 C- Advective flux in Y
232 IF (.NOT. multiDimAdvection .OR.
233 & advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
234 & advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
235 & advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
236 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
237 CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
238 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
239 CALL GAD_FLUXLIMIT_ADV_Y(
240 & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
241 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
242 CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
243 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
244 CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
245 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
246 CALL GAD_DST3_ADV_Y(
247 & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
248 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
249 CALL GAD_DST3FL_ADV_Y(
250 & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
251 ELSE
252 STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
253 ENDIF
254 DO j=1-Oly,sNy+Oly
255 DO i=1-Olx,sNx+Olx
256 fMer(i,j) = fMer(i,j) + af(i,j)
257 ENDDO
258 ENDDO
259 ENDIF
260
261 C- Diffusive flux in Y
262 IF (diffKh.NE.0.) THEN
263 CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
264 ELSE
265 DO j=1-Oly,sNy+Oly
266 DO i=1-Olx,sNx+Olx
267 df(i,j) = 0. _d 0
268 ENDDO
269 ENDDO
270 ENDIF
271
272 #ifdef ALLOW_GMREDI
273 C- GM/Redi flux in Y
274 IF (useGMRedi) THEN
275 C *note* should update GMREDI_YTRANSPORT to use localT and set df *aja*
276 CALL GMREDI_YTRANSPORT(
277 I iMin,iMax,jMin,jMax,bi,bj,K,
278 I yA,Tracer,tracerIdentity,
279 U df,
280 I myThid)
281 ENDIF
282 #endif
283 DO j=1-Oly,sNy+Oly
284 DO i=1-Olx,sNx+Olx
285 fMer(i,j) = fMer(i,j) + df(i,j)
286 ENDDO
287 ENDDO
288
289 C- Bi-harmonic flux in Y
290 IF (diffK4 .NE. 0.) THEN
291 CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
292 DO j=1-Oly,sNy+Oly
293 DO i=1-Olx,sNx+Olx
294 fMer(i,j) = fMer(i,j) + df(i,j)
295 ENDDO
296 ENDDO
297 ENDIF
298
299 C- Advective flux in R
300 IF (.NOT. multiDimAdvection .OR.
301 & advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
302 & advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
303 & advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
304 C Note: wVel needs to be masked
305 IF (K.GE.2) THEN
306 C- Compute vertical advective flux in the interior:
307 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
308 CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
309 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
310 CALL GAD_FLUXLIMIT_ADV_R(
311 & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
312 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
313 CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
314 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
315 CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
316 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
317 CALL GAD_DST3_ADV_R(
318 & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
319 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
320 CALL GAD_DST3FL_ADV_R(
321 & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
322 ELSE
323 STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
324 ENDIF
325 C- Surface "correction" term at k>1 :
326 DO j=1-Oly,sNy+Oly
327 DO i=1-Olx,sNx+Olx
328 af(i,j) = af(i,j)
329 & + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*
330 & rTrans(i,j)*Tracer(i,j,k,bi,bj)
331 ENDDO
332 ENDDO
333 ELSE
334 C- Surface "correction" term at k=1 :
335 DO j=1-Oly,sNy+Oly
336 DO i=1-Olx,sNx+Olx
337 af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)
338 ENDDO
339 ENDDO
340 ENDIF
341 C- add the advective flux to fVerT
342 DO j=1-Oly,sNy+Oly
343 DO i=1-Olx,sNx+Olx
344 fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
345 ENDDO
346 ENDDO
347 ENDIF
348
349 C- Diffusive flux in R
350 C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
351 C boundary condition.
352 IF (implicitDiffusion) THEN
353 DO j=1-Oly,sNy+Oly
354 DO i=1-Olx,sNx+Olx
355 df(i,j) = 0. _d 0
356 ENDDO
357 ENDDO
358 ELSE
359 CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
360 ENDIF
361 c DO j=1-Oly,sNy+Oly
362 c DO i=1-Olx,sNx+Olx
363 c fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
364 c ENDDO
365 c ENDDO
366
367 #ifdef ALLOW_GMREDI
368 C- GM/Redi flux in R
369 IF (useGMRedi) THEN
370 C *note* should update GMREDI_RTRANSPORT to set df *aja*
371 CALL GMREDI_RTRANSPORT(
372 I iMin,iMax,jMin,jMax,bi,bj,K,
373 I Tracer,tracerIdentity,
374 U df,
375 I myThid)
376 c DO j=1-Oly,sNy+Oly
377 c DO i=1-Olx,sNx+Olx
378 c fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
379 c ENDDO
380 c ENDDO
381 ENDIF
382 #endif
383
384 DO j=1-Oly,sNy+Oly
385 DO i=1-Olx,sNx+Olx
386 fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
387 ENDDO
388 ENDDO
389
390 #ifdef ALLOW_KPP
391 C- Add non local KPP transport term (ghat) to diffusive T flux.
392 IF (useKPP) THEN
393 DO j=1-Oly,sNy+Oly
394 DO i=1-Olx,sNx+Olx
395 df(i,j) = 0. _d 0
396 ENDDO
397 ENDDO
398 IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
399 C *note* should update KPP_TRANSPORT_T to set df *aja*
400 CALL KPP_TRANSPORT_T(
401 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
402 I KappaRT,
403 U df )
404 ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
405 CALL KPP_TRANSPORT_S(
406 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
407 I KappaRT,
408 U df )
409 ELSEIF (tracerIdentity.EQ.GAD_TR1) THEN
410 CALL KPP_TRANSPORT_TR1(
411 I iMin,iMax,jMin,jMax,bi,bj,k,km1,
412 I KappaRT,
413 U df )
414 ELSE
415 STOP 'GAD_CALC_RHS: Ooops'
416 ENDIF
417 DO j=1-Oly,sNy+Oly
418 DO i=1-Olx,sNx+Olx
419 fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
420 ENDDO
421 ENDDO
422 ENDIF
423 #endif
424
425 C-- Divergence of fluxes
426 DO j=1-Oly,sNy+Oly-1
427 DO i=1-Olx,sNx+Olx-1
428 gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
429 & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
430 & *recip_rA(i,j,bi,bj)
431 & *(
432 & +( fZon(i+1,j)-fZon(i,j) )
433 & +( fMer(i,j+1)-fMer(i,j) )
434 & +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
435 & )
436 ENDDO
437 ENDDO
438
439 RETURN
440 END

  ViewVC Help
Powered by ViewVC 1.1.22