/[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.9 - (hide annotations) (download)
Mon Sep 10 13:09:04 2001 UTC (22 years, 8 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint40pre9
Changes since 1.8: +5 -7 lines
Added third dimension for DST method.

1 adcroft 1.9 C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_calc_rhs.F,v 1.8 2001/09/10 01:22:48 adcroft 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 adcroft 1.9 CALL GAD_DST3_ADV_R(
277     & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
278 adcroft 1.4 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
279 adcroft 1.9 CALL GAD_DST3FL_ADV_R(
280     & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
281 jmc 1.2 ELSE
282 adcroft 1.3 STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
283 jmc 1.2 ENDIF
284     C- Surface "correction" term at k>1 :
285     DO j=1-Oly,sNy+Oly
286     DO i=1-Olx,sNx+Olx
287     af(i,j) = af(i,j)
288     & + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*
289     & rTrans(i,j)*Tracer(i,j,k,bi,bj)
290     ENDDO
291     ENDDO
292 adcroft 1.1 ELSE
293 jmc 1.2 C- Surface "correction" term at k=1 :
294     DO j=1-Oly,sNy+Oly
295     DO i=1-Olx,sNx+Olx
296     af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)
297     ENDDO
298     ENDDO
299 adcroft 1.1 ENDIF
300 jmc 1.2 C- add the advective flux to fVerT
301 adcroft 1.5 DO j=1-Oly,sNy+Oly
302     DO i=1-Olx,sNx+Olx
303 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)
304     ENDDO
305     ENDDO
306 adcroft 1.8 ENDIF
307 adcroft 1.1
308     C- Diffusive flux in R
309     C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
310     C boundary condition.
311     IF (implicitDiffusion) THEN
312 adcroft 1.5 DO j=1-Oly,sNy+Oly
313     DO i=1-Olx,sNx+Olx
314 adcroft 1.1 df(i,j) = 0.
315     ENDDO
316     ENDDO
317     ELSE
318     CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
319     ENDIF
320 adcroft 1.5 c DO j=1-Oly,sNy+Oly
321     c DO i=1-Olx,sNx+Olx
322 adcroft 1.1 c fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
323     c ENDDO
324     c ENDDO
325    
326     #ifdef ALLOW_GMREDI
327     C- GM/Redi flux in R
328     IF (useGMRedi) THEN
329     C *note* should update GMREDI_RTRANSPORT to set df *aja*
330     CALL GMREDI_RTRANSPORT(
331     I iMin,iMax,jMin,jMax,bi,bj,K,
332 adcroft 1.6 I Tracer,
333 adcroft 1.1 U df,
334     I myThid)
335 adcroft 1.5 c DO j=1-Oly,sNy+Oly
336     c DO i=1-Olx,sNx+Olx
337 adcroft 1.1 c fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
338     c ENDDO
339     c ENDDO
340     ENDIF
341     #endif
342    
343 adcroft 1.5 DO j=1-Oly,sNy+Oly
344     DO i=1-Olx,sNx+Olx
345 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
346     ENDDO
347     ENDDO
348    
349     #ifdef ALLOW_KPP
350     C- Add non local KPP transport term (ghat) to diffusive T flux.
351     IF (useKPP) THEN
352 adcroft 1.5 DO j=1-Oly,sNy+Oly
353     DO i=1-Olx,sNx+Olx
354 adcroft 1.1 df(i,j) = 0.
355     ENDDO
356     ENDDO
357     IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
358     C *note* should update KPP_TRANSPORT_T to set df *aja*
359     CALL KPP_TRANSPORT_T(
360     I iMin,iMax,jMin,jMax,bi,bj,k,km1,
361     I KappaRT,
362     U df )
363     ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
364     CALL KPP_TRANSPORT_S(
365     I iMin,iMax,jMin,jMax,bi,bj,k,km1,
366     I KappaRT,
367     U df )
368     ELSE
369     STOP 'GAD_CALC_RHS: Ooops'
370     ENDIF
371 adcroft 1.5 DO j=1-Oly,sNy+Oly
372     DO i=1-Olx,sNx+Olx
373 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
374     ENDDO
375     ENDDO
376     ENDIF
377     #endif
378    
379     C-- Divergence of fluxes
380 adcroft 1.5 DO j=1-Oly,sNy+Oly
381     DO i=1-Olx,sNx+Olx
382 adcroft 1.8 gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
383 adcroft 1.1 & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
384     & *recip_rA(i,j,bi,bj)
385     & *(
386     & +( fZon(i+1,j)-fZon(i,j) )
387     & +( fMer(i,j+1)-fMer(i,j) )
388     & +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
389     & )
390     ENDDO
391     ENDDO
392    
393     RETURN
394     END

  ViewVC Help
Powered by ViewVC 1.1.22