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

Contents of /MITgcm_contrib/gael/pkg/smooth/smooth_rhs.F

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


Revision 1.2 - (show annotations) (download)
Fri Oct 16 03:36:34 2009 UTC (15 years, 9 months ago) by gforget
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +5 -5 lines
bring pkg/smooth up to date

1 #include "CPP_OPTIONS.h"
2 #include "GAD_OPTIONS.h"
3
4 C !INTERFACE: ==========================================================
5 SUBROUTINE smooth_rhs(fld_in,gt_in,myThid)
6
7 C !DESCRIPTION:
8 C Calculates the tendency of a tracer due to advection and diffusion.
9 C It calculates the fluxes in each direction indepentently and then
10 C sets the tendency to the divergence of these fluxes. The advective
11 C fluxes are only calculated here when using the linear advection schemes
12 C otherwise only the diffusive and parameterized fluxes are calculated.
13 C
14 C Contributions to the flux are calculated and added:
15 C \begin{equation*}
16 C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
17 C \end{equation*}
18 C
19 C The tendency is the divergence of the fluxes:
20 C \begin{equation*}
21 C G_\theta = G_\theta + \nabla \cdot {\bf F}
22 C \end{equation*}
23 C
24 C The tendency is assumed to contain data on entry.
25
26 C !USES: ===============================================================
27 IMPLICIT NONE
28 #include "SIZE.h"
29 #include "EEPARAMS.h"
30 #include "PARAMS.h"
31 #include "GRID.h"
32 #include "SURFACE.h"
33 #include "GAD.h"
34 #include "GMREDI.h"
35
36 #ifdef ALLOW_AUTODIFF_TAMC
37 #include "tamc.h"
38 #include "tamc_keys.h"
39 #endif /* ALLOW_AUTODIFF_TAMC */
40
41 #include "smooth.h"
42
43 C !INPUT PARAMETERS: ===================================================
44
45 c ... xA, yA, maskUp => in commons ?determined at init?
46 c dTdy...: peuvent devenir (bi,bj)
47
48 INTEGER bi,bj,iMin,iMax,jMin,jMax
49 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
50 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51 _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
52 _RL dTdz (nSx,nSy)
53 _RL dTdx (nSx,nSy)
54 _RL dTdy (nSx,nSy)
55 INTEGER myThid
56 INTEGER i,j,k
57
58 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
59 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
60 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
61 _RL df (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
62 _RL fld_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
63 _RL gt_in(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
64 integer wc01_di, wc01_dj, wc01_dk, ii, jj, kk
65 _RL wc01_tmp(2)
66
67 DO bj=myByLo(myThid),myByHi(myThid)
68 DO bi=myBxLo(myThid),myBxHi(myThid)
69
70
71 c 1rst k loop: initialization
72 DO k=1,Nr
73 DO j=1-OLy,sNy+OLy
74 DO i=1-OLx,sNx+OLx
75 fZon(i,j,k,bi,bj) = 0. _d 0
76 fMer(i,j,k,bi,bj) = 0. _d 0
77 fVerT(i,j,k,bi,bj) = 0. _d 0
78 gt_in(i,j,k,bi,bj) = 0. _d 0
79 ENDDO
80 ENDDO
81 ENDDO
82
83 cnot needed et !!!!!!!!bibj _EXCH_XYZ_RL( gt_in , myThid )
84
85 iMin = 1-OLx+1
86 iMax = sNx+OLx-1
87 jMin = 1-OLy+1
88 jMax = sNy+OLy-1
89
90 c 2nd k loop: flux computation
91 DO k=1,Nr
92
93 DO j=1-OLy,sNy+OLy
94 DO i=1-OLx,sNx+OLx
95 df(i,j,bi,bj) = 0. _d 0
96 xA(i,j,bi,bj) = _dyG(i,j,bi,bj)
97 & *drF(k)*_hFacW(i,j,k,bi,bj)
98 yA(i,j,bi,bj) = _dxG(i,j,bi,bj)
99 & *drF(k)*_hFacS(i,j,k,bi,bj)
100 IF (K .EQ. 1) THEN
101 maskUp(i,j,bi,bj) = 0.
102 ELSE
103 maskUp(i,j,bi,bj) =
104 & maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
105 ENDIF
106 ENDDO
107 ENDDO
108
109
110 c ///gmredi_xtr///
111
112 DO j=jMin,jMax
113 DO i=iMin,iMax
114 df(i,j,bi,bj) = df(i,j,bi,bj)
115 & -xA(i,j,bi,bj)
116 & *wc01_Kux(i,j,k,bi,bj)
117 & *recip_dxC(i,j,bi,bj)
118 & *(fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj))
119 ENDDO
120 ENDDO
121
122 DO j=jMin,jMax
123 DO i=iMin,iMax
124 dTdz(bi,bj) = op5*(
125 & +op5*recip_drC(k)*
126 & ( maskC(i-1,j,k,bi,bj)*
127 & (fld_in(i-1,j, MAX(k-1,1) ,bi,bj)-fld_in(i-1,j,k,bi,bj))
128 & +maskC( i ,j,k,bi,bj)*
129 & (fld_in( i ,j, MAX(k-1,1) ,bi,bj)-fld_in( i ,j,k,bi,bj))
130 & )
131 & +op5*recip_drC(MIN(k+1,Nr))*
132 & ( maskC(i-1,j,MIN(k+1,Nr),bi,bj)*
133 & (fld_in(i-1,j,k,bi,bj)-fld_in(i-1,j,MIN(k+1,Nr),bi,bj))
134 & +maskC( i ,j,MIN(k+1,Nr),bi,bj)*
135 & (fld_in( i ,j,k,bi,bj)-fld_in( i ,j,MIN(k+1,Nr),bi,bj))
136 & ) )
137 df(i,j,bi,bj) = df(i,j,bi,bj)
138 & - xA(i,j,bi,bj)*wc01_Kuz(i,j,k,bi,bj)*dTdz(bi,bj)
139 ENDDO
140 ENDDO
141
142 DO j=jMin,jMax
143 DO i=iMin,iMax
144 dTdy(bi,bj) = op5*(
145 & +op5*(maskS(i,j,k,bi,bj)
146 & *recip_dyC(i,j,bi,bj)*
147 & (fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj))
148 & +maskS(i,j+1,k,bi,bj)
149 & *recip_dyC(i,j+1,bi,bj)*
150 & (fld_in(i,j+1,k,bi,bj)-fld_in(i,j,k,bi,bj)))
151 & +op5*(maskS(i-1,j,k,bi,bj)
152 & *recip_dyC(i,j,bi,bj)*
153 & (fld_in(i-1,j,k,bi,bj)-fld_in(i-1,j-1,k,bi,bj))
154 & +maskS(i-1,j+1,k,bi,bj)
155 & *recip_dyC(i,j+1,bi,bj)*
156 & (fld_in(i-1,j+1,k,bi,bj)-fld_in(i-1,j,k,bi,bj)))
157 & )
158 df(i,j,bi,bj) = df(i,j,bi,bj)
159 & - xA(i,j,bi,bj)*wc01_Kuy(i,j,k,bi,bj)*dTdy(bi,bj)
160 ENDDO
161 ENDDO
162
163
164 c /// end for x ///
165
166 DO j=jMin,jMax
167 DO i=iMin,iMax
168 fZon(i,j,k,bi,bj) = fZon(i,j,k,bi,bj) + df(i,j,bi,bj)
169 ENDDO
170 ENDDO
171
172 DO j=jMin,jMax
173 DO i=iMin,iMax
174 df(i,j,bi,bj) = 0.
175 ENDDO
176 ENDDO
177
178 c ///gmredi_ytr///
179
180 DO j=jMin,jMax
181 DO i=iMin,iMax
182 df(i,j,bi,bj) = df(i,j,bi,bj)
183 & -yA(i,j,bi,bj)
184 & *wc01_Kvy(i,j,k,bi,bj)
185 & *recip_dyC(i,j,bi,bj)
186 & *(fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj))
187 ENDDO
188 ENDDO
189
190 DO j=jMin,jMax
191 DO i=iMin,iMax
192 dTdz(bi,bj) = op5*(
193 & +op5*recip_drC(k)*
194 & ( maskC(i,j-1,k,bi,bj)*
195 & (fld_in(i,j-1,MAX(k-1,1),bi,bj)-fld_in(i,j-1,k,bi,bj))
196 & +maskC(i, j ,k,bi,bj)*
197 & (fld_in(i, j ,MAX(k-1,1),bi,bj)-fld_in(i, j ,k,bi,bj))
198 & )
199 & +op5*recip_drC(MIN(k+1,Nr))*
200 & ( maskC(i,j-1,MIN(k+1,Nr),bi,bj)*
201 & (fld_in(i,j-1,k,bi,bj)-fld_in(i,j-1,MIN(k+1,Nr),bi,bj))
202 & +maskC(i, j ,MIN(k+1,Nr),bi,bj)*
203 & (fld_in(i, j ,k,bi,bj)-fld_in(i, j ,MIN(k+1,Nr),bi,bj))
204 & ) )
205 df(i,j,bi,bj) = df(i,j,bi,bj)
206 & - yA(i,j,bi,bj)*wc01_Kvz(i,j,k,bi,bj)*dTdz(bi,bj)
207 ENDDO
208 ENDDO
209
210 DO j=jMin,jMax
211 DO i=iMin,iMax
212 dTdx(bi,bj) = op5*(
213 & +op5*(maskW(i+1,j,k,bi,bj)
214 & *recip_dxC(i+1,j,bi,bj)*
215 & (fld_in(i+1,j,k,bi,bj)-fld_in(i,j,k,bi,bj))
216 & +maskW(i,j,k,bi,bj)
217 & *recip_dxC(i,j,bi,bj)*
218 & (fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj)))
219 & +op5*(maskW(i+1,j-1,k,bi,bj)
220 & *recip_dxC(i+1,j,bi,bj)*
221 & (fld_in(i+1,j-1,k,bi,bj)-fld_in(i,j-1,k,bi,bj))
222 & +maskW(i,j-1,k,bi,bj)
223 & *recip_dxC(i,j,bi,bj)*
224 & (fld_in(i,j-1,k,bi,bj)-fld_in(i-1,j-1,k,bi,bj)))
225 & )
226 df(i,j,bi,bj) = df(i,j,bi,bj)
227 & - yA(i,j,bi,bj)*wc01_Kvx(i,j,k,bi,bj)*dTdx(bi,bj)
228 ENDDO
229 ENDDO
230
231 c /// end for y ///
232
233 DO j=jMin,jMax
234 DO i=iMin,iMax
235 fMer(i,j,k,bi,bj) = fMer(i,j,k,bi,bj) + df(i,j,bi,bj)
236 ENDDO
237 ENDDO
238
239 DO j=jMin,jMax
240 DO i=iMin,iMax
241 df(i,j,bi,bj) = 0.
242 ENDDO
243 ENDDO
244
245 c /// GAD_DIFF_R ///
246
247 if (.NOT. wc01_implicitDiffusion ) then
248
249 IF (k.gt.1) then
250 DO j=jMin,jMax
251 DO i=iMin,iMax
252 df(i,j,bi,bj) =
253 & -_rA(i,j,bi,bj)
254 & *KappaRwc01(i,j,k,bi,bj)*recip_drC(k)
255 & *(fld_in(i,j,k,bi,bj)
256 & -fld_in(i,j,k-1,bi,bj))*rkSign
257 ENDDO
258 ENDDO
259 ENDIF
260
261 endif
262
263 c ///gmredi rtrans///
264
265 IF (K.GT.1) THEN
266 DO j=jMin,jMax
267 DO i=iMin,iMax
268 dTdx(bi,bj) = op5*(
269 & +op5*(maskW(i+1,j,k,bi,bj)
270 & *recip_dxC(i+1,j,bi,bj)*
271 & (fld_in(i+1,j,k,bi,bj)-fld_in(i,j,k,bi,bj))
272 & +maskW(i,j,k,bi,bj)
273 & *recip_dxC(i,j,bi,bj)*
274 & (fld_in(i,j,k,bi,bj)-fld_in(i-1,j,k,bi,bj)))
275 & +op5*(maskW(i+1,j,k-1,bi,bj)
276 & *recip_dxC(i+1,j,bi,bj)*
277 & (fld_in(i+1,j,k-1,bi,bj)-fld_in(i,j,k-1,bi,bj))
278 & +maskW(i,j,k-1,bi,bj)
279 & *recip_dxC(i,j,bi,bj)*
280 & (fld_in(i,j,k-1,bi,bj)-fld_in(i-1,j,k-1,bi,bj)))
281 & )
282
283 dTdy(bi,bj) = op5*(
284 & +op5*(maskS(i,j,k,bi,bj)
285 & *recip_dyC(i,j,bi,bj)*
286 & (fld_in(i,j,k,bi,bj)-fld_in(i,j-1,k,bi,bj))
287 & +maskS(i,j+1,k,bi,bj)
288 & *recip_dyC(i,j+1,bi,bj)*
289 & (fld_in(i,j+1,k,bi,bj)-fld_in(i,j,k,bi,bj)))
290 & +op5*(maskS(i,j,k-1,bi,bj)
291 & *recip_dyC(i,j,bi,bj)*
292 & (fld_in(i,j,k-1,bi,bj)-fld_in(i,j-1,k-1,bi,bj))
293 & +maskS(i,j+1,k-1,bi,bj)
294 & *recip_dyC(i,j+1,bi,bj)*
295 & (fld_in(i,j+1,k-1,bi,bj)-fld_in(i,j,k-1,bi,bj)))
296 & )
297
298 df(i,j,bi,bj) = df(i,j,bi,bj)
299 & - rA(i,j,bi,bj)
300 & *( wc01_Kwx(i,j,k,bi,bj)*dTdx(bi,bj)
301 & +wc01_Kwy(i,j,k,bi,bj)*dTdy(bi,bj) )
302
303 ENDDO
304 ENDDO
305
306 ENDIF
307
308
309 c /// end for r ///
310
311 IF (K.GT.1) THEN
312 DO j=jMin,jMax
313 DO i=iMin,iMax
314 fVerT(i,j,k-1,bi,bj) = fVerT(i,j,k-1,bi,bj) +
315 & df(i,j,bi,bj)*maskUp(i,j,bi,bj)
316 ENDDO
317 ENDDO
318 ENDIF
319
320 DO j=jMin,jMax
321 DO i=iMin,iMax
322 df(i,j,bi,bj) = 0.
323 ENDDO
324 ENDDO
325
326 ENDDO
327
328 ENDDO
329 ENDDO
330
331 c these exchanges are crucial... two bibj loops
332 _EXCH_XYZ_RL ( fZon , myThid )
333 _EXCH_XYZ_RL ( fMer , myThid )
334 _EXCH_XYZ_RL ( fVerT, myThid )
335
336 DO bj=myByLo(myThid),myByHi(myThid)
337 DO bi=myBxLo(myThid),myBxHi(myThid)
338 c 3rd k loop: Divergence of fluxes
339 DO k=1,Nr
340 IF (K.GT.1) THEN
341 DO j=jMin,jMax
342 DO i=iMin,iMax
343 gt_in(i,j,k,bi,bj)=gt_in(i,j,k,bi,bj)
344 & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
345 & *( (fZon(i+1,j,k,bi,bj)-fZon(i,j,k,bi,bj))
346 & +(fMer(i,j+1,k,bi,bj)-fMer(i,j,k,bi,bj))
347 & +(fVerT(i,j,k,bi,bj)-fVerT(i,j,k-1,bi,bj))*rkSign
348 & )
349 ENDDO
350 ENDDO
351 ELSE
352 DO j=jMin,jMax
353 DO i=iMin,iMax
354 gt_in(i,j,k,bi,bj)=gt_in(i,j,k,bi,bj)
355 & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
356 & *( (fZon(i+1,j,k,bi,bj)-fZon(i,j,k,bi,bj))
357 & +(fMer(i,j+1,k,bi,bj)-fMer(i,j,k,bi,bj))
358 & +(fVerT(i,j,k,bi,bj))*rkSign
359 & )
360 ENDDO
361 ENDDO
362 ENDIF
363 ENDDO
364 ENDDO
365 ENDDO
366
367 _EXCH_XYZ_RL ( gt_in , myThid )
368
369 END
370

  ViewVC Help
Powered by ViewVC 1.1.22