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

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

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


Revision 1.11 - (hide annotations) (download)
Wed Sep 19 20:45:09 2001 UTC (22 years, 8 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint41
Changes since 1.10: +63 -24 lines
Added comments in form compatible with "protex".

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

  ViewVC Help
Powered by ViewVC 1.1.22