/[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.16 - (hide annotations) (download)
Tue Dec 10 03:02:00 2002 UTC (21 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint47e_post, checkpoint47c_post, checkpoint48c_pre, checkpoint47d_pre, checkpoint47i_post, checkpoint47d_post, checkpoint47g_post, checkpoint48a_post, checkpoint47j_post, branch-exfmods-tag, checkpoint48c_post, checkpoint47f_post, checkpoint48, checkpoint47h_post
Branch point for: branch-exfmods-curt
Changes since 1.15: +17 -1 lines
 * OCEANICP & realFreshWater: include P-E direct effect on wVel ;
   NOTES: requires option NONLIN_FRSURF to be "#define".

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

  ViewVC Help
Powered by ViewVC 1.1.22