/[MITgcm]/MITgcm/verification/global_with_CFC11/code1x1/gad_calc_rhs.F
ViewVC logotype

Diff of /MITgcm/verification/global_with_CFC11/code1x1/gad_calc_rhs.F

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

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