/[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.6 - (hide annotations) (download)
Tue Sep 4 17:16:11 2001 UTC (22 years, 9 months ago) by adcroft
Branch: MAIN
Changes since 1.5: +2 -2 lines
maskUp was unused by gmrei_rtransport. Deleted from argument list.

1 adcroft 1.6 C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_calc_rhs.F,v 1.5 2001/09/04 17:00: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    
82     C-- Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
83     IF (diffK4 .NE. 0.) THEN
84     CALL GAD_GRAD_X(bi,bj,k,xA,localT,fZon,myThid)
85     CALL GAD_GRAD_Y(bi,bj,k,yA,localT,fMer,myThid)
86     CALL GAD_DEL2(bi,bj,k,fZon,fMer,df4,myThid)
87     ENDIF
88    
89     C-- Initialize net flux in X direction
90     DO j=1-Oly,sNy+Oly
91     DO i=1-Olx,sNx+Olx
92     fZon(i,j) = 0.
93     ENDDO
94     ENDDO
95    
96     C- Advective flux in X
97 adcroft 1.3 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
98 adcroft 1.1 CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
99 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
100 adcroft 1.1 CALL GAD_FLUXLIMIT_ADV_X(
101     & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
102 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
103 jmc 1.2 CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
104 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
105 adcroft 1.1 CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
106 adcroft 1.4 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
107     CALL GAD_DST3_ADV_X(
108     & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
109     ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
110     CALL GAD_DST3FL_ADV_X(
111     & bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
112 adcroft 1.1 ELSE
113 adcroft 1.3 STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
114 adcroft 1.1 ENDIF
115 adcroft 1.5 DO j=1-Oly,sNy+Oly
116     DO i=1-Olx,sNx+Olx
117 adcroft 1.1 fZon(i,j) = fZon(i,j) + af(i,j)
118     ENDDO
119     ENDDO
120    
121     C- Diffusive flux in X
122     IF (diffKh.NE.0.) THEN
123     CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
124     ELSE
125 adcroft 1.5 DO j=1-Oly,sNy+Oly
126     DO i=1-Olx,sNx+Olx
127 adcroft 1.1 df(i,j) = 0.
128     ENDDO
129     ENDDO
130     ENDIF
131    
132     #ifdef ALLOW_GMREDI
133     C- GM/Redi flux in X
134     IF (useGMRedi) THEN
135     C *note* should update GMREDI_XTRANSPORT to use localT and set df *aja*
136     CALL GMREDI_XTRANSPORT(
137     I iMin,iMax,jMin,jMax,bi,bj,K,
138     I xA,Tracer,
139     U df,
140     I myThid)
141     ENDIF
142     #endif
143 adcroft 1.5 DO j=1-Oly,sNy+Oly
144     DO i=1-Olx,sNx+Olx
145 adcroft 1.1 fZon(i,j) = fZon(i,j) + df(i,j)
146     ENDDO
147     ENDDO
148    
149     C- Bi-harmonic duffusive flux in X
150     IF (diffK4 .NE. 0.) THEN
151     CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
152 adcroft 1.5 DO j=1-Oly,sNy+Oly
153     DO i=1-Olx,sNx+Olx
154 adcroft 1.1 fZon(i,j) = fZon(i,j) + df(i,j)
155     ENDDO
156     ENDDO
157     ENDIF
158    
159     C-- Initialize net flux in Y direction
160     DO j=1-Oly,sNy+Oly
161     DO i=1-Olx,sNx+Olx
162     fMer(i,j) = 0.
163     ENDDO
164     ENDDO
165    
166     C- Advective flux in Y
167 adcroft 1.3 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
168 adcroft 1.1 CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
169 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
170 adcroft 1.1 CALL GAD_FLUXLIMIT_ADV_Y(
171     & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
172 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
173 jmc 1.2 CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
174 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
175 adcroft 1.1 CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
176 adcroft 1.4 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
177     CALL GAD_DST3_ADV_Y(
178     & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
179     ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
180     CALL GAD_DST3FL_ADV_Y(
181     & bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
182 adcroft 1.1 ELSE
183 adcroft 1.3 STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
184 adcroft 1.1 ENDIF
185     DO j=1-Oly,sNy+Oly
186     DO i=1-Olx,sNx+Olx
187     fMer(i,j) = fMer(i,j) + af(i,j)
188     ENDDO
189     ENDDO
190    
191     C- Diffusive flux in Y
192     IF (diffKh.NE.0.) THEN
193     CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
194     ELSE
195     DO j=1-Oly,sNy+Oly
196     DO i=1-Olx,sNx+Olx
197     df(i,j) = 0.
198     ENDDO
199     ENDDO
200     ENDIF
201    
202     #ifdef ALLOW_GMREDI
203     C- GM/Redi flux in Y
204     IF (useGMRedi) THEN
205     CALL GMREDI_YTRANSPORT(
206     C *note* should update GMREDI_YTRANSPORT to use localT and set df *aja*
207     I iMin,iMax,jMin,jMax,bi,bj,K,
208     I yA,Tracer,
209     U df,
210     I myThid)
211     ENDIF
212     #endif
213     DO j=1-Oly,sNy+Oly
214     DO i=1-Olx,sNx+Olx
215     fMer(i,j) = fMer(i,j) + df(i,j)
216     ENDDO
217     ENDDO
218    
219     C- Bi-harmonic flux in Y
220     IF (diffK4 .NE. 0.) THEN
221     CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
222     DO j=1-Oly,sNy+Oly
223     DO i=1-Olx,sNx+Olx
224     fMer(i,j) = fMer(i,j) + df(i,j)
225     ENDDO
226     ENDDO
227     ENDIF
228    
229     C-- Initialize net flux in R
230 adcroft 1.5 DO j=1-Oly,sNy+Oly
231     DO i=1-Olx,sNx+Olx
232 adcroft 1.1 fVerT(i,j,kUp) = 0.
233     ENDDO
234     ENDDO
235    
236     C- Advective flux in R
237 jmc 1.2 C Note: wVel needs to be masked
238     IF (K.GE.2) THEN
239     C- Compute vertical advective flux in the interior:
240 adcroft 1.3 IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
241 jmc 1.2 CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
242 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
243 jmc 1.2 CALL GAD_FLUXLIMIT_ADV_R(
244 adcroft 1.1 & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
245 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
246 jmc 1.2 CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
247 adcroft 1.3 ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
248 jmc 1.2 CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
249 adcroft 1.4 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
250     c CALL GAD_DST3_ADV_R(
251     c & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
252     STOP 'GAD_CALC_RHS: GAD_DST3_ADV_R not coded yet'
253     ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
254     c CALL GAD_DST3FL_ADV_R(
255     c & bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
256     STOP 'GAD_CALC_RHS: GAD_DST3FL_ADV_R not coded yet'
257 jmc 1.2 ELSE
258 adcroft 1.3 STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
259 jmc 1.2 ENDIF
260     C- Surface "correction" term at k>1 :
261     DO j=1-Oly,sNy+Oly
262     DO i=1-Olx,sNx+Olx
263     af(i,j) = af(i,j)
264     & + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*
265     & rTrans(i,j)*Tracer(i,j,k,bi,bj)
266     ENDDO
267     ENDDO
268 adcroft 1.1 ELSE
269 jmc 1.2 C- Surface "correction" term at k=1 :
270     DO j=1-Oly,sNy+Oly
271     DO i=1-Olx,sNx+Olx
272     af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)
273     ENDDO
274     ENDDO
275 adcroft 1.1 ENDIF
276 jmc 1.2 C- add the advective flux to fVerT
277 adcroft 1.5 DO j=1-Oly,sNy+Oly
278     DO i=1-Olx,sNx+Olx
279 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)
280     ENDDO
281     ENDDO
282    
283     C- Diffusive flux in R
284     C Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
285     C boundary condition.
286     IF (implicitDiffusion) THEN
287 adcroft 1.5 DO j=1-Oly,sNy+Oly
288     DO i=1-Olx,sNx+Olx
289 adcroft 1.1 df(i,j) = 0.
290     ENDDO
291     ENDDO
292     ELSE
293     CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
294     ENDIF
295 adcroft 1.5 c DO j=1-Oly,sNy+Oly
296     c DO i=1-Olx,sNx+Olx
297 adcroft 1.1 c fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
298     c ENDDO
299     c ENDDO
300    
301     #ifdef ALLOW_GMREDI
302     C- GM/Redi flux in R
303     IF (useGMRedi) THEN
304     C *note* should update GMREDI_RTRANSPORT to set df *aja*
305     CALL GMREDI_RTRANSPORT(
306     I iMin,iMax,jMin,jMax,bi,bj,K,
307 adcroft 1.6 I Tracer,
308 adcroft 1.1 U df,
309     I myThid)
310 adcroft 1.5 c DO j=1-Oly,sNy+Oly
311     c DO i=1-Olx,sNx+Olx
312 adcroft 1.1 c fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
313     c ENDDO
314     c ENDDO
315     ENDIF
316     #endif
317    
318 adcroft 1.5 DO j=1-Oly,sNy+Oly
319     DO i=1-Olx,sNx+Olx
320 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
321     ENDDO
322     ENDDO
323    
324     #ifdef ALLOW_KPP
325     C- Add non local KPP transport term (ghat) to diffusive T flux.
326     IF (useKPP) THEN
327 adcroft 1.5 DO j=1-Oly,sNy+Oly
328     DO i=1-Olx,sNx+Olx
329 adcroft 1.1 df(i,j) = 0.
330     ENDDO
331     ENDDO
332     IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
333     C *note* should update KPP_TRANSPORT_T to set df *aja*
334     CALL KPP_TRANSPORT_T(
335     I iMin,iMax,jMin,jMax,bi,bj,k,km1,
336     I KappaRT,
337     U df )
338     ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
339     CALL KPP_TRANSPORT_S(
340     I iMin,iMax,jMin,jMax,bi,bj,k,km1,
341     I KappaRT,
342     U df )
343     ELSE
344     STOP 'GAD_CALC_RHS: Ooops'
345     ENDIF
346 adcroft 1.5 DO j=1-Oly,sNy+Oly
347     DO i=1-Olx,sNx+Olx
348 adcroft 1.1 fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)
349     ENDDO
350     ENDDO
351     ENDIF
352     #endif
353    
354     C-- Divergence of fluxes
355 adcroft 1.5 DO j=1-Oly,sNy+Oly
356     DO i=1-Olx,sNx+Olx
357 adcroft 1.1 gTracer(i,j,k,bi,bj)=
358     & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
359     & *recip_rA(i,j,bi,bj)
360     & *(
361     & +( fZon(i+1,j)-fZon(i,j) )
362     & +( fMer(i,j+1)-fMer(i,j) )
363     & +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac
364     & )
365     ENDDO
366     ENDDO
367    
368     RETURN
369     END

  ViewVC Help
Powered by ViewVC 1.1.22