/[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.8 - (hide annotations) (download)
Mon Sep 10 01:22:48 2001 UTC (22 years, 8 months ago) by adcroft
Branch: MAIN
Changes since 1.7: +29 -2 lines
Added multi-dimensional form of advection
 o available only for single step schemes (ie. can't be used with ABII)
 o stable for max(cfl_u,cfl_v,cfl_w)<=1  (without cfl_u+cfl_v+cfl_w <=1)
 o selected using multiDimAdvection=.T.  (default)
 o had to hack some existing routines to work on local arrays

1 adcroft 1.8 C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_calc_rhs.F,v 1.7 2001/09/05 17:46:03 heimbach Exp $
2 jmc 1.2 C $Name: $
3 adcroft 1.1
4     #include "GAD_OPTIONS.h"
5    
6     SUBROUTINE GAD_CALC_RHS(
7     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
8     I xA,yA,uTrans,vTrans,rTrans,maskUp,
9     I diffKh, diffK4, KappaRT, Tracer,
10 adcroft 1.3 I tracerIdentity, advectionScheme,
11 adcroft 1.1 U fVerT, gTracer,
12     I myThid )
13     C /==========================================================\
14     C | SUBROUTINE GAD_CALC_RHS |
15     C |==========================================================|
16     C \==========================================================/
17     IMPLICIT NONE
18    
19     C == GLobal variables ==
20     #include "SIZE.h"
21     #include "EEPARAMS.h"
22     #include "PARAMS.h"
23     #include "GRID.h"
24     #include "DYNVARS.h"
25     #include "GAD.h"
26    
27     C == Routine arguments ==
28     INTEGER k,kUp,kDown,kM1
29     INTEGER bi,bj,iMin,iMax,jMin,jMax
30     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32     _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33     _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
34     _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
35     _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
36     _RL diffKh, diffK4
37     _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
38     _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
39     INTEGER tracerIdentity
40 adcroft 1.3 INTEGER advectionScheme
41 adcroft 1.1 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
42     _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
43     INTEGER myThid
44    
45     C == Local variables ==
46     C I, J, K - Loop counters
47     INTEGER i,j
48     LOGICAL TOP_LAYER
49     _RL afFacT, dfFacT
50     _RL df4 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51     _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53     _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54     _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55     _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56    
57     #ifdef ALLOW_AUTODIFF_TAMC
58     C-- only the kUp part of fverT is set in this subroutine
59     C-- the kDown is still required
60     fVerT(1,1,kDown) = fVerT(1,1,kDown)
61     #endif
62     DO j=1-OLy,sNy+OLy
63     DO i=1-OLx,sNx+OLx
64     fZon(i,j) = 0.0
65     fMer(i,j) = 0.0
66     fVerT(i,j,kUp) = 0.0
67     ENDDO
68     ENDDO
69    
70     afFacT = 1. _d 0
71     dfFacT = 1. _d 0
72     TOP_LAYER = K .EQ. 1
73    
74     C-- Make local copy of tracer array
75     DO j=1-OLy,sNy+OLy
76     DO i=1-OLx,sNx+OLx
77     localT(i,j)=tracer(i,j,k,bi,bj)
78     ENDDO
79     ENDDO
80    
81 adcroft 1.8 C-- Unless we have already calculated the advection terms we initialize
82     C the tendency to zero.
83     IF (.NOT. multiDimAdvection .OR.
84     & advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
85     & advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
86     & advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
87     DO j=1-Oly,sNy+Oly
88     DO i=1-Olx,sNx+Olx
89     gTracer(i,j,k,bi,bj)=0.
90     ENDDO
91     ENDDO
92     ENDIF
93 adcroft 1.1
94     C-- Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
95     IF (diffK4 .NE. 0.) THEN
96     CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
97     CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
98     CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
99     ENDIF
100    
101     C-- Initialize net flux in X direction
102     DO j=1-Oly,sNy+Oly
103     DO i=1-Olx,sNx+Olx
104     fZon(i,j) = 0.
105     ENDDO
106     ENDDO
107    
108     C- Advective flux in X
109 adcroft 1.8 IF (.NOT. multiDimAdvection .OR.
110     & advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
111     & advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
112     & advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
113 adcroft 1.3 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
114 adcroft 1.1 CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
115 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
116 adcroft 1.1 CALL GAD_FLUXLIMIT_ADV_X(
117     & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
118 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
119 jmc 1.2 CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
120 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
121 adcroft 1.1 CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
122 adcroft 1.4 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
123     CALL GAD_DST3_ADV_X(
124     & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
125     ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
126     CALL GAD_DST3FL_ADV_X(
127     & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
128 adcroft 1.1 ELSE
129 adcroft 1.3 STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
130 adcroft 1.1 ENDIF
131 adcroft 1.5 DO j=1-Oly,sNy+Oly
132     DO i=1-Olx,sNx+Olx
133 adcroft 1.1 fZon(i,j) = fZon(i,j) + af(i,j)
134     ENDDO
135     ENDDO
136 adcroft 1.8 ENDIF
137 adcroft 1.1
138     C- Diffusive flux in X
139     IF (diffKh.NE.0.) THEN
140     CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
141     ELSE
142 adcroft 1.5 DO j=1-Oly,sNy+Oly
143     DO i=1-Olx,sNx+Olx
144 adcroft 1.1 df(i,j) = 0.
145     ENDDO
146     ENDDO
147     ENDIF
148    
149     #ifdef ALLOW_GMREDI
150     C- GM/Redi flux in X
151     IF (useGMRedi) THEN
152     C *note* should update GMREDI_XTRANSPORT to use localT and set df *aja*
153     CALL GMREDI_XTRANSPORT(
154     I iMin,iMax,jMin,jMax,bi,bj,K,
155     I xA,Tracer,
156     U df,
157     I myThid)
158     ENDIF
159     #endif
160 adcroft 1.5 DO j=1-Oly,sNy+Oly
161     DO i=1-Olx,sNx+Olx
162 adcroft 1.1 fZon(i,j) = fZon(i,j) + df(i,j)
163     ENDDO
164     ENDDO
165    
166     C- Bi-harmonic duffusive flux in X
167     IF (diffK4 .NE. 0.) THEN
168     CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
169 adcroft 1.5 DO j=1-Oly,sNy+Oly
170     DO i=1-Olx,sNx+Olx
171 adcroft 1.1 fZon(i,j) = fZon(i,j) + df(i,j)
172     ENDDO
173     ENDDO
174     ENDIF
175    
176     C-- Initialize net flux in Y direction
177     DO j=1-Oly,sNy+Oly
178     DO i=1-Olx,sNx+Olx
179     fMer(i,j) = 0.
180     ENDDO
181     ENDDO
182    
183     C- Advective flux in Y
184 adcroft 1.8 IF (.NOT. multiDimAdvection .OR.
185     & advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
186     & advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
187     & advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
188 adcroft 1.3 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
189 adcroft 1.1 CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
190 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
191 adcroft 1.1 CALL GAD_FLUXLIMIT_ADV_Y(
192     & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
193 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
194 jmc 1.2 CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
195 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
196 adcroft 1.1 CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
197 adcroft 1.4 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
198     CALL GAD_DST3_ADV_Y(
199     & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
200     ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
201     CALL GAD_DST3FL_ADV_Y(
202     & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
203 adcroft 1.1 ELSE
204 adcroft 1.3 STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
205 adcroft 1.1 ENDIF
206     DO j=1-Oly,sNy+Oly
207     DO i=1-Olx,sNx+Olx
208     fMer(i,j) = fMer(i,j) + af(i,j)
209     ENDDO
210     ENDDO
211 adcroft 1.8 ENDIF
212 adcroft 1.1
213     C- Diffusive flux in Y
214     IF (diffKh.NE.0.) THEN
215     CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
216     ELSE
217     DO j=1-Oly,sNy+Oly
218     DO i=1-Olx,sNx+Olx
219     df(i,j) = 0.
220     ENDDO
221     ENDDO
222     ENDIF
223    
224     #ifdef ALLOW_GMREDI
225     C- GM/Redi flux in Y
226     IF (useGMRedi) THEN
227 heimbach 1.7 C *note* should update GMREDI_YTRANSPORT to use localT and set df *aja*
228 adcroft 1.1 CALL GMREDI_YTRANSPORT(
229     I iMin,iMax,jMin,jMax,bi,bj,K,
230     I yA,Tracer,
231     U df,
232     I myThid)
233     ENDIF
234     #endif
235     DO j=1-Oly,sNy+Oly
236     DO i=1-Olx,sNx+Olx
237     fMer(i,j) = fMer(i,j) + df(i,j)
238     ENDDO
239     ENDDO
240    
241     C- Bi-harmonic flux in Y
242     IF (diffK4 .NE. 0.) THEN
243     CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
244     DO j=1-Oly,sNy+Oly
245     DO i=1-Olx,sNx+Olx
246     fMer(i,j) = fMer(i,j) + df(i,j)
247     ENDDO
248     ENDDO
249     ENDIF
250    
251     C-- Initialize net flux in R
252 adcroft 1.5 DO j=1-Oly,sNy+Oly
253     DO i=1-Olx,sNx+Olx
254 adcroft 1.1 fVerT(i,j,kUp) = 0.
255     ENDDO
256     ENDDO
257    
258     C- Advective flux in R
259 adcroft 1.8 IF (.NOT. multiDimAdvection .OR.
260     & advectionScheme.EQ.ENUM_CENTERED_2ND .OR.
261     & advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
262     & advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN
263 jmc 1.2 C Note: wVel needs to be masked
264     IF (K.GE.2) THEN
265     C- Compute vertical advective flux in the interior:
266 adcroft 1.3 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
267 jmc 1.2 CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
268 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
269 jmc 1.2 CALL GAD_FLUXLIMIT_ADV_R(
270 adcroft 1.1 & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
271 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
272 jmc 1.2 CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
273 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
274 jmc 1.2 CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
275 adcroft 1.4 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
276     c CALL GAD_DST3_ADV_R(
277     c & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
278     STOP 'GAD_CALC_RHS: GAD_DST3_ADV_R not coded yet'
279     ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
280     c CALL GAD_DST3FL_ADV_R(
281     c & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
282     STOP 'GAD_CALC_RHS: GAD_DST3FL_ADV_R not coded yet'
283 jmc 1.2 ELSE
284 adcroft 1.3 STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
285 jmc 1.2 ENDIF
286     C- Surface "correction" term at k>1 :
287     DO j=1-Oly,sNy+Oly
288     DO i=1-Olx,sNx+Olx
289     af(i,j) = af(i,j)
290     & + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*
291     & rTrans(i,j)*Tracer(i,j,k,bi,bj)
292     ENDDO
293     ENDDO
294 adcroft 1.1 ELSE
295 jmc 1.2 C- Surface "correction" term at k=1 :
296     DO j=1-Oly,sNy+Oly
297     DO i=1-Olx,sNx+Olx
298     af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)
299     ENDDO
300     ENDDO
301 adcroft 1.1 ENDIF
302 jmc 1.2 C- add the advective flux to fVerT
303 adcroft 1.5 DO j=1-Oly,sNy+Oly
304     DO i=1-Olx,sNx+Olx
305 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)
306     ENDDO
307     ENDDO
308 adcroft 1.8 ENDIF
309 adcroft 1.1
310     C- Diffusive flux in R
311     C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
312     C boundary condition.
313     IF (implicitDiffusion) THEN
314 adcroft 1.5 DO j=1-Oly,sNy+Oly
315     DO i=1-Olx,sNx+Olx
316 adcroft 1.1 df(i,j) = 0.
317     ENDDO
318     ENDDO
319     ELSE
320     CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
321     ENDIF
322 adcroft 1.5 c DO j=1-Oly,sNy+Oly
323     c DO i=1-Olx,sNx+Olx
324 adcroft 1.1 c fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
325     c ENDDO
326     c ENDDO
327    
328     #ifdef ALLOW_GMREDI
329     C- GM/Redi flux in R
330     IF (useGMRedi) THEN
331     C *note* should update GMREDI_RTRANSPORT to set df *aja*
332     CALL GMREDI_RTRANSPORT(
333     I iMin,iMax,jMin,jMax,bi,bj,K,
334 adcroft 1.6 I Tracer,
335 adcroft 1.1 U df,
336     I myThid)
337 adcroft 1.5 c DO j=1-Oly,sNy+Oly
338     c DO i=1-Olx,sNx+Olx
339 adcroft 1.1 c fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
340     c ENDDO
341     c ENDDO
342     ENDIF
343     #endif
344    
345 adcroft 1.5 DO j=1-Oly,sNy+Oly
346     DO i=1-Olx,sNx+Olx
347 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
348     ENDDO
349     ENDDO
350    
351     #ifdef ALLOW_KPP
352     C- Add non local KPP transport term (ghat) to diffusive T flux.
353     IF (useKPP) THEN
354 adcroft 1.5 DO j=1-Oly,sNy+Oly
355     DO i=1-Olx,sNx+Olx
356 adcroft 1.1 df(i,j) = 0.
357     ENDDO
358     ENDDO
359     IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
360     C *note* should update KPP_TRANSPORT_T to set df *aja*
361     CALL KPP_TRANSPORT_T(
362     I iMin,iMax,jMin,jMax,bi,bj,k,km1,
363     I KappaRT,
364     U df )
365     ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
366     CALL KPP_TRANSPORT_S(
367     I iMin,iMax,jMin,jMax,bi,bj,k,km1,
368     I KappaRT,
369     U df )
370     ELSE
371     STOP 'GAD_CALC_RHS: Ooops'
372     ENDIF
373 adcroft 1.5 DO j=1-Oly,sNy+Oly
374     DO i=1-Olx,sNx+Olx
375 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
376     ENDDO
377     ENDDO
378     ENDIF
379     #endif
380    
381     C-- Divergence of fluxes
382 adcroft 1.5 DO j=1-Oly,sNy+Oly
383     DO i=1-Olx,sNx+Olx
384 adcroft 1.8 gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
385 adcroft 1.1 & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
386     & *recip_rA(i,j,bi,bj)
387     & *(
388     & +( fZon(i+1,j)-fZon(i,j) )
389     & +( fMer(i,j+1)-fMer(i,j) )
390     & +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
391     & )
392     ENDDO
393     ENDDO
394    
395     RETURN
396     END

  ViewVC Help
Powered by ViewVC 1.1.22