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

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

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


Revision 1.1.2.1 - (hide annotations) (download)
Thu Aug 25 16:22:17 2005 UTC (18 years, 7 months ago) by dimitri
Branch: release1_50yr
Changes since 1.1: +440 -0 lines
adding ecco1x1 verification/global_with_CFC11 experiment

1 dimitri 1.1.2.1 C $Header: /u/gcmpack/MITgcm/verification/global_with_CFC11/code50yr/Attic/gad_calc_rhs.F,v 1.1.2.1 2003/05/04 23:19:18 dimitri Exp $
2     C $Name: release1_50yr $
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

  ViewVC Help
Powered by ViewVC 1.1.22