/[MITgcm]/MITgcm_contrib/gael/pkg/smooth2/smooth_rhs.F
ViewVC logotype

Annotation of /MITgcm_contrib/gael/pkg/smooth2/smooth_rhs.F

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


Revision 1.2 - (hide annotations) (download)
Mon Feb 15 19:18:55 2010 UTC (15 years, 5 months ago) by gforget
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -2 lines
bug fix

1 gforget 1.1 C $Header: /u/gcmpack/MITgcm_contrib/gael/pkg/smooth2/smooth_rhs.F,v 1.1 2009/10/24 23:27:24 gforget Exp $
2     C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     C !INTERFACE: ==========================================================
7     SUBROUTINE smooth_rhs(fld_in,gt_in,myThid)
8    
9     C *==========================================================*
10     C | SUBROUTINE smooth_rhs
11     C | o As part of smooth_diff3D, this routine computes the
12     C | right hand side of the tendency equation (see below).
13     C | It is made of bits from model/src and pkg/generic_advdiff
14     C | pieced togheter.
15     C *==========================================================*
16    
17    
18     C !DESCRIPTION:
19     C Calculates the tendency of a tracer due to advection and diffusion.
20     C It calculates the fluxes in each direction indepentently and then
21     C sets the tendency 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 tendency 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 tendency 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 "SURFACE.h"
44    
45     #ifdef ALLOW_AUTODIFF_TAMC
46     #include "tamc.h"
47     #include "tamc_keys.h"
48     #endif /* ALLOW_AUTODIFF_TAMC */
49    
50     #include "smooth.h"
51    
52     C !INPUT PARAMETERS: ===================================================
53    
54    
55     INTEGER bi,bj,iMin,iMax,jMin,jMax
56     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
58     _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
59     _RL dTdz (nSx,nSy)
60     _RL dTdx (nSx,nSy)
61     _RL dTdy (nSx,nSy)
62     INTEGER myThid
63     INTEGER i,j,k
64    
65     _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
66     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
67     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
68     _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
69     _RL fld_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
70     _RL gt_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
71    
72     DO bj=myByLo(myThid),myByHi(myThid)
73     DO bi=myBxLo(myThid),myBxHi(myThid)
74    
75    
76     c 1rst k loop: initialization
77     DO k=1,Nr
78     DO j=1-OLy,sNy+OLy
79     DO i=1-OLx,sNx+OLx
80     fZon(i,j,k,bi,bj) = 0. _d 0
81     fMer(i,j,k,bi,bj) = 0. _d 0
82     fVerT(i,j,k,bi,bj) = 0. _d 0
83     gt_in(i,j,k,bi,bj) = 0. _d 0
84     ENDDO
85     ENDDO
86     ENDDO
87    
88    
89     iMin = 1-OLx+1
90     iMax = sNx+OLx-1
91     jMin = 1-OLy+1
92     jMax = sNy+OLy-1
93    
94     c 2nd k loop: flux computation
95     DO k=1,Nr
96    
97     DO j=1-OLy,sNy+OLy
98     DO i=1-OLx,sNx+OLx
99     df(i,j,bi,bj) = 0. _d 0
100     xA(i,j,bi,bj) = _dyG(i,j,bi,bj)
101     & *drF(k)*_hFacW(i,j,k,bi,bj)
102     yA(i,j,bi,bj) = _dxG(i,j,bi,bj)
103     & *drF(k)*_hFacS(i,j,k,bi,bj)
104     IF (K .EQ. 1) THEN
105     maskUp(i,j,bi,bj) = 0.
106     ELSE
107     maskUp(i,j,bi,bj) =
108     & maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
109     ENDIF
110     ENDDO
111     ENDDO
112    
113    
114     c ///gmredi_xtr///
115    
116     DO j=jMin,jMax
117     DO i=iMin,iMax
118     df(i,j,bi,bj) = df(i,j,bi,bj)
119     & -xA(i,j,bi,bj)
120     & *smooth3D_Kux(i,j,k,bi,bj)
121     & *recip_dxC(i,j,bi,bj)
122     & *(fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj))
123     ENDDO
124     ENDDO
125    
126     DO j=jMin,jMax
127     DO i=iMin,iMax
128     dTdz(bi,bj) = 0.5*(
129     & +0.5*recip_drC(k)*
130     & ( maskC(i-1,j,k,bi,bj)*
131     & (fld_in(i-1,j, MAX(k-1,1) ,bi,bj)-fld_in(i-1,j,k,bi,bj))
132     & +maskC( i ,j,k,bi,bj)*
133     & (fld_in( i ,j, MAX(k-1,1) ,bi,bj)-fld_in( i ,j,k,bi,bj))
134     & )
135     & +0.5*recip_drC(MIN(k+1,Nr))*
136     & ( maskC(i-1,j,MIN(k+1,Nr),bi,bj)*
137     & (fld_in(i-1,j,k,bi,bj)-fld_in(i-1,j,MIN(k+1,Nr),bi,bj))
138     & +maskC( i ,j,MIN(k+1,Nr),bi,bj)*
139     & (fld_in( i ,j,k,bi,bj)-fld_in( i ,j,MIN(k+1,Nr),bi,bj))
140     & ) )
141     df(i,j,bi,bj) = df(i,j,bi,bj)
142     & - xA(i,j,bi,bj)*smooth3D_Kuz(i,j,k,bi,bj)*dTdz(bi,bj)
143     ENDDO
144     ENDDO
145    
146     DO j=jMin,jMax
147     DO i=iMin,iMax
148     dTdy(bi,bj) = 0.5*(
149     & +0.5*(maskS(i,j,k,bi,bj)
150     & *recip_dyC(i,j,bi,bj)*
151     & (fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj))
152     & +maskS(i,j+1,k,bi,bj)
153     & *recip_dyC(i,j+1,bi,bj)*
154     & (fld_in(i,j+1,k,bi,bj)-fld_in(i,j,k,bi,bj)))
155     & +0.5*(maskS(i-1,j,k,bi,bj)
156     & *recip_dyC(i,j,bi,bj)*
157     & (fld_in(i-1,j,k,bi,bj)-fld_in(i-1,j-1,k,bi,bj))
158     & +maskS(i-1,j+1,k,bi,bj)
159     & *recip_dyC(i,j+1,bi,bj)*
160     & (fld_in(i-1,j+1,k,bi,bj)-fld_in(i-1,j,k,bi,bj)))
161     & )
162     df(i,j,bi,bj) = df(i,j,bi,bj)
163     & - xA(i,j,bi,bj)*smooth3D_Kuy(i,j,k,bi,bj)*dTdy(bi,bj)
164     ENDDO
165     ENDDO
166    
167    
168     c /// end for x ///
169    
170     DO j=jMin,jMax
171     DO i=iMin,iMax
172     fZon(i,j,k,bi,bj) = fZon(i,j,k,bi,bj) + df(i,j,bi,bj)
173     ENDDO
174     ENDDO
175    
176     DO j=jMin,jMax
177     DO i=iMin,iMax
178     df(i,j,bi,bj) = 0.
179     ENDDO
180     ENDDO
181    
182     c ///gmredi_ytr///
183    
184     DO j=jMin,jMax
185     DO i=iMin,iMax
186     df(i,j,bi,bj) = df(i,j,bi,bj)
187     & -yA(i,j,bi,bj)
188     & *smooth3D_Kvy(i,j,k,bi,bj)
189     & *recip_dyC(i,j,bi,bj)
190     & *(fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj))
191     ENDDO
192     ENDDO
193    
194     DO j=jMin,jMax
195     DO i=iMin,iMax
196     dTdz(bi,bj) = 0.5*(
197     & +0.5*recip_drC(k)*
198     & ( maskC(i,j-1,k,bi,bj)*
199     & (fld_in(i,j-1,MAX(k-1,1),bi,bj)-fld_in(i,j-1,k,bi,bj))
200     & +maskC(i, j ,k,bi,bj)*
201     & (fld_in(i, j ,MAX(k-1,1),bi,bj)-fld_in(i, j ,k,bi,bj))
202     & )
203     & +0.5*recip_drC(MIN(k+1,Nr))*
204     & ( maskC(i,j-1,MIN(k+1,Nr),bi,bj)*
205     & (fld_in(i,j-1,k,bi,bj)-fld_in(i,j-1,MIN(k+1,Nr),bi,bj))
206     & +maskC(i, j ,MIN(k+1,Nr),bi,bj)*
207     & (fld_in(i, j ,k,bi,bj)-fld_in(i, j ,MIN(k+1,Nr),bi,bj))
208     & ) )
209     df(i,j,bi,bj) = df(i,j,bi,bj)
210     & - yA(i,j,bi,bj)*smooth3D_Kvz(i,j,k,bi,bj)*dTdz(bi,bj)
211     ENDDO
212     ENDDO
213    
214     DO j=jMin,jMax
215     DO i=iMin,iMax
216     dTdx(bi,bj) = 0.5*(
217     & +0.5*(maskW(i+1,j,k,bi,bj)
218     & *recip_dxC(i+1,j,bi,bj)*
219     & (fld_in(i+1,j,k,bi,bj)-fld_in(i,j,k,bi,bj))
220     & +maskW(i,j,k,bi,bj)
221     & *recip_dxC(i,j,bi,bj)*
222     & (fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj)))
223     & +0.5*(maskW(i+1,j-1,k,bi,bj)
224     & *recip_dxC(i+1,j,bi,bj)*
225     & (fld_in(i+1,j-1,k,bi,bj)-fld_in(i,j-1,k,bi,bj))
226     & +maskW(i,j-1,k,bi,bj)
227     & *recip_dxC(i,j,bi,bj)*
228     & (fld_in(i,j-1,k,bi,bj)-fld_in(i-1,j-1,k,bi,bj)))
229     & )
230     df(i,j,bi,bj) = df(i,j,bi,bj)
231     & - yA(i,j,bi,bj)*smooth3D_Kvx(i,j,k,bi,bj)*dTdx(bi,bj)
232     ENDDO
233     ENDDO
234    
235     c /// end for y ///
236    
237     DO j=jMin,jMax
238     DO i=iMin,iMax
239     fMer(i,j,k,bi,bj) = fMer(i,j,k,bi,bj) + df(i,j,bi,bj)
240     ENDDO
241     ENDDO
242    
243     DO j=jMin,jMax
244     DO i=iMin,iMax
245     df(i,j,bi,bj) = 0.
246     ENDDO
247     ENDDO
248    
249     c /// GAD_DIFF_R ///
250    
251     if (.NOT. smooth3DdoImpldiff ) then
252    
253     IF (k.gt.1) then
254     DO j=jMin,jMax
255     DO i=iMin,iMax
256     df(i,j,bi,bj) =
257     & -_rA(i,j,bi,bj)
258     & *smooth3D_kappaR(i,j,k,bi,bj)*recip_drC(k)
259     & *(fld_in(i,j,k,bi,bj)
260     & -fld_in(i,j,k-1,bi,bj))*rkSign
261     ENDDO
262     ENDDO
263     ENDIF
264    
265     endif
266    
267     c ///gmredi rtrans///
268    
269     IF (K.GT.1) THEN
270     DO j=jMin,jMax
271     DO i=iMin,iMax
272     dTdx(bi,bj) = 0.5*(
273     & +0.5*(maskW(i+1,j,k,bi,bj)
274     & *recip_dxC(i+1,j,bi,bj)*
275     & (fld_in(i+1,j,k,bi,bj)-fld_in(i,j,k,bi,bj))
276     & +maskW(i,j,k,bi,bj)
277     & *recip_dxC(i,j,bi,bj)*
278     & (fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj)))
279     & +0.5*(maskW(i+1,j,k-1,bi,bj)
280     & *recip_dxC(i+1,j,bi,bj)*
281     & (fld_in(i+1,j,k-1,bi,bj)-fld_in(i,j,k-1,bi,bj))
282     & +maskW(i,j,k-1,bi,bj)
283     & *recip_dxC(i,j,bi,bj)*
284     & (fld_in(i,j,k-1,bi,bj)-fld_in(i-1,j,k-1,bi,bj)))
285     & )
286    
287     dTdy(bi,bj) = 0.5*(
288     & +0.5*(maskS(i,j,k,bi,bj)
289     & *recip_dyC(i,j,bi,bj)*
290     & (fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj))
291     & +maskS(i,j+1,k,bi,bj)
292     & *recip_dyC(i,j+1,bi,bj)*
293     & (fld_in(i,j+1,k,bi,bj)-fld_in(i,j,k,bi,bj)))
294     & +0.5*(maskS(i,j,k-1,bi,bj)
295     & *recip_dyC(i,j,bi,bj)*
296     & (fld_in(i,j,k-1,bi,bj)-fld_in(i,j-1,k-1,bi,bj))
297     & +maskS(i,j+1,k-1,bi,bj)
298     & *recip_dyC(i,j+1,bi,bj)*
299     & (fld_in(i,j+1,k-1,bi,bj)-fld_in(i,j,k-1,bi,bj)))
300     & )
301    
302     df(i,j,bi,bj) = df(i,j,bi,bj)
303     & - rA(i,j,bi,bj)
304     & *( smooth3D_Kwx(i,j,k,bi,bj)*dTdx(bi,bj)
305     & +smooth3D_Kwy(i,j,k,bi,bj)*dTdy(bi,bj) )
306    
307     ENDDO
308     ENDDO
309    
310     ENDIF
311    
312    
313     c /// end for r ///
314    
315     IF (K.GT.1) THEN
316     DO j=jMin,jMax
317     DO i=iMin,iMax
318     fVerT(i,j,k-1,bi,bj) = fVerT(i,j,k-1,bi,bj) +
319     & df(i,j,bi,bj)*maskUp(i,j,bi,bj)
320     ENDDO
321     ENDDO
322     ENDIF
323    
324     DO j=jMin,jMax
325     DO i=iMin,iMax
326     df(i,j,bi,bj) = 0.
327     ENDDO
328     ENDDO
329    
330     ENDDO
331    
332     ENDDO
333     ENDDO
334    
335     c these exchanges are crucial:
336 gforget 1.2 CALL EXCH_UV_XYZ_RL(fZon,fMer,.TRUE.,myThid)
337 gforget 1.1 _EXCH_XYZ_RL ( fVerT, myThid )
338    
339     DO bj=myByLo(myThid),myByHi(myThid)
340     DO bi=myBxLo(myThid),myBxHi(myThid)
341     c 3rd k loop: Divergence of fluxes
342     DO k=1,Nr
343     IF (K.GT.1) THEN
344     DO j=jMin,jMax
345     DO i=iMin,iMax
346     gt_in(i,j,k,bi,bj)=gt_in(i,j,k,bi,bj)
347     & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
348     & *( (fZon(i+1,j,k,bi,bj)-fZon(i,j,k,bi,bj))
349     & +(fMer(i,j+1,k,bi,bj)-fMer(i,j,k,bi,bj))
350     & +(fVerT(i,j,k,bi,bj)-fVerT(i,j,k-1,bi,bj))*rkSign
351     & )
352     ENDDO
353     ENDDO
354     ELSE
355     DO j=jMin,jMax
356     DO i=iMin,iMax
357     gt_in(i,j,k,bi,bj)=gt_in(i,j,k,bi,bj)
358     & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
359     & *( (fZon(i+1,j,k,bi,bj)-fZon(i,j,k,bi,bj))
360     & +(fMer(i,j+1,k,bi,bj)-fMer(i,j,k,bi,bj))
361     & +(fVerT(i,j,k,bi,bj))*rkSign
362     & )
363     ENDDO
364     ENDDO
365     ENDIF
366     ENDDO
367     ENDDO
368     ENDDO
369    
370     _EXCH_XYZ_RL ( gt_in , myThid )
371    
372     END
373    

  ViewVC Help
Powered by ViewVC 1.1.22