Ticket #649: mpdata_adiff.F

File mpdata_adiff.F, 52.4 KB (added by jcwarner, 11 years ago)

mpdata_adiff corrected ranges

Line 
1#include "cppdefs.h"
2#define MPDATA_HOT
3
4 MODULE mpdata_adiff_mod
5#if defined NONLINEAR && defined TS_MPDATA && defined SOLVE3D
6!
7!svn $Id: mpdata_adiff.F 732 2008-09-07 01:55:51Z jcwarner $
8!================================================== Hernan G. Arango ===
9! Copyright (c) 2002-2014 The ROMS/TOMS Group John C. Warner !
10! Licensed under a MIT/X style license !
11! See License_ROMS.txt !
12!=======================================================================
13! !
14! This routine computes anti-diffusive velocities to correct tracer !
15! advection using MPDATA Recursive method. !
16! !
17! On Output: !
18! !
19! Ua Andi-diffusive velocity in the XI-direction (m/s). !
20! Va Anti-diffusive velocity in the ETA-direction (m/s). !
21! Wa Anti-diffusive velocity in the S-direction (m3/s). !
22! !
23! Reference: !
24! !
25! Margolin, L. and P.K. Smolarkiewicz, 1998: Antidiffusive !
26! velocities for multipass donor cell advection, SIAM J. !
27! Sci. Comput., 907-929. !
28! !
29!=======================================================================
30!
31 implicit none
32
33 PUBLIC :: mpdata_adiff_tile
34
35 CONTAINS
36!
37!***********************************************************************
38 SUBROUTINE mpdata_adiff_tile (ng, tile, &
39 & LBi, UBi, LBj, UBj, &
40 & IminS, ImaxS, JminS, JmaxS, &
41# ifdef MASKING
42 & rmask, umask, vmask, &
43# endif
44# ifdef WET_DRY
45 & rmask_wet, umask_wet, vmask_wet, &
46# endif
47 & pm, pn, omn, om_u, on_v, &
48 & z_r, oHz, &
49 & Huon, Hvom, w, &
50# ifdef WEC_VF
51 & W_stokes, &
52# endif
53 & t, Ta, Ua, Va, Wa)
54!***********************************************************************
55!
56 USE mod_param
57 USE mod_ncparam
58 USE mod_scalars
59!
60! Imported variable declarations.
61!
62 integer, intent(in) :: ng, tile
63 integer, intent(in) :: LBi, UBi, LBj, UBj
64 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
65!
66# ifdef ASSUMED_SHAPE
67# ifdef MASKING
68 real(r8), intent(in) :: rmask(LBi:,LBj:)
69 real(r8), intent(in) :: umask(LBi:,LBj:)
70 real(r8), intent(in) :: vmask(LBi:,LBj:)
71# endif
72# ifdef WET_DRY
73 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
74 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
75 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
76# endif
77 real(r8), intent(in) :: pm(LBi:,LBj:)
78 real(r8), intent(in) :: pn(LBi:,LBj:)
79 real(r8), intent(in) :: omn(LBi:,LBj:)
80 real(r8), intent(in) :: om_u(LBi:,LBj:)
81 real(r8), intent(in) :: on_v(LBi:,LBj:)
82 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
83 real(r8), intent(in) :: oHz(IminS:,JminS:,:)
84 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
85 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
86 real(r8), intent(in) :: t(LBi:,LBj:,:)
87 real(r8), intent(in) :: w(LBi:,LBj:,0:)
88# ifdef WEC_VF
89 real(r8), intent(in) :: W_stokes(LBi:,LBj:,0:)
90# endif
91 real(r8), intent(inout) :: Ta(IminS:,JminS:,:)
92
93 real(r8), intent(out) :: Ua(IminS:,JminS:,:)
94 real(r8), intent(out) :: Va(IminS:,JminS:,:)
95 real(r8), intent(out) :: Wa(IminS:,JminS:,0:)
96# else
97# ifdef MASKING
98 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
99 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
100 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
101# endif
102# ifdef WET_DRY
103 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
104 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
105 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
106# endif
107 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
108 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
109 real(r8), intent(in) :: omn(LBi:UBi,LBj:UBj)
110 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
111 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
112 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
113 real(r8), intent(in) :: oHz(IminS:ImaxS,JminS:JmaxS,N(ng))
114 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
115 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
116 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng))
117 real(r8), intent(in) :: w(LBi:UBi,LBj:UBj,0:N(ng))
118# ifdef WEC_VF
119 real(r8), intent(in) :: W_stokes(LBi:UBi,LBj:UBj,0:N(ng))
120# endif
121 real(r8), intent(inout) :: Ta(IminS:ImaxS,JminS:JmaxS,N(ng))
122
123 real(r8), intent(out) :: Ua(IminS:ImaxS,JminS:JmaxS,N(ng))
124 real(r8), intent(out) :: Va(IminS:ImaxS,JminS:JmaxS,N(ng))
125 real(r8), intent(out) :: Wa(IminS:ImaxS,JminS:JmaxS,0:N(ng))
126# endif
127!
128! Local variable declarations.
129!
130 integer :: i, is, j, k
131
132 real(r8), parameter :: Large = 1.0E+20_r8
133 real(r8), parameter :: eps = 1.0E-18_r8
134 real(r8), parameter :: eps2= 1.0E-10_r8
135
136 real(r8) :: A, B, Tmax, Tmin, Um, Vm, X, Y, Z
137 real(r8) :: cff, cff1, cff2, sig_alfa, fac
138# ifdef MPDATA_HOT
139 real(r8) :: AA, BB, CC, AB, AC, BC
140 real(r8) :: XX, YY, ZZ, XY, XZ, YZ
141 real(r8) :: sig_beta, sig_gama
142 real(r8) :: sig_a, sig_b, sig_c, sig_d, sig_e, sig_f
143# endif
144
145 real(r8), dimension(IminS:ImaxS,N(ng)) :: C
146 real(r8), dimension(IminS:ImaxS,N(ng)) :: Wm
147
148 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: mask_dn
149 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: mask_up
150
151 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: beta_dn
152 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: beta_up
153 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: odz
154
155# include "set_bounds.h"
156# ifdef TS_MPDATA_LIMIT
157 fac=0.25_r8
158# else
159 fac=1.0_r8
160# endif
161!
162!-----------------------------------------------------------------------
163! Compute anti-diffusive MPDATA velocities (Ua,Va,Wa).
164!-----------------------------------------------------------------------
165!
166! Set boundary conditions to advected tracer, Ta.
167!
168 IF (.not.EWperiodic(ng)) THEN
169 IF (DOMAIN(ng)%Western_Edge(tile)) THEN
170 DO k=1,N(ng)
171 DO j=JstrVm2,Jendp2i
172 Ta(Istr-1,j,k)=Ta(Istr,j,k)
173 END DO
174 END DO
175 END IF
176 IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
177 DO k=1,N(ng)
178 DO j=JstrVm2,Jendp2i
179 Ta(Iend+1,j,k)=Ta(Iend,j,k)
180 END DO
181 END DO
182 END IF
183 END IF
184
185 IF (.not.NSperiodic(ng)) THEN
186 IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
187 DO k=1,N(ng)
188 DO i=IstrUm2,Iendp2i
189 Ta(i,Jstr-1,k)=Ta(i,Jstr,k)
190 END DO
191 END DO
192 END IF
193 IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
194 DO k=1,N(ng)
195 DO i=IstrUm2,Iendp2i
196 Ta(i,Jend+1,k)=Ta(i,Jend,k)
197 END DO
198 END DO
199 END IF
200 END IF
201
202 IF (.not.(EWperiodic(ng).or.NSperiodic(ng))) THEN
203 IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
204 DO k=1,N(ng)
205 Ta(Istr-1,Jstr-1,k)=0.5_r8*(Ta(Istr ,Jstr-1,k)+ &
206 & Ta(Istr-1,Jstr ,k))
207 END DO
208 END IF
209 IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
210 DO k=1,N(ng)
211 Ta(Iend+1,Jstr-1,k)=0.5_r8*(Ta(Iend+1,Jstr ,k)+ &
212 & Ta(Iend ,Jstr-1,k))
213 END DO
214 END IF
215 IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
216 DO k=1,N(ng)
217 Ta(Istr-1,Jend+1,k)=0.5_r8*(Ta(Istr-1,Jend ,k)+ &
218 & Ta(Istr ,Jend+1,k))
219 END DO
220 END IF
221 IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
222 DO k=1,N(ng)
223 Ta(Iend+1,Jend+1,k)=0.5_r8*(Ta(Iend+1,Jend ,k)+ &
224 & Ta(Iend ,Jend+1,k))
225 END DO
226 END IF
227 END IF
228!
229! Compute inverse vertical grid spacing at W-points.
230!
231 DO k=1,N(ng)-1
232 DO j=Jstrm2,Jendp2
233 DO i=Istrm2,Iendp2
234 odz(i,j,k)=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
235 END DO
236 END DO
237 END DO
238 cff=1.0_r8/dt(ng)
239!
240! Compute nondimensional U-antidiffusive velocities, Ua. If applicable,
241! retain up to third-order terms of the power series.
242!
243 DO j=JstrV-1,Jendp1
244 k=1
245! DO i=IstrU-1,Iendp2
246 DO i=IstrUm1,Iendp2
247 C(i,k)=0.25_r8* &
248 & ((Ta(i ,j,k+1)-Ta(i ,j,k ))*odz(i ,j,k )+ &
249 & (Ta(i-1,j,k+1)-Ta(i-1,j,k ))*odz(i-1,j,k ))* &
250 & (z_r(i ,j,k+1)-z_r(i ,j,k)+ &
251 & z_r(i-1,j,k+1)-z_r(i-1,j,k))/ &
252 & (Ta(i-1,j,k)+Ta(i,j,k)+eps)
253 Wm(i,k)=0.25_r8*dt(ng)* &
254 & (w(i-1,j,k )*odz(i-1,j,k)*pm(i-1,j)*pn(i-1,j)+ &
255 & w(i ,j,k )*odz(i ,j,k)*pm(i ,j)*pn(i ,j))
256# ifdef WEC_VF
257 Wm(i,k)=Wm(i,k)+ &
258 & 0.25_r8*dt(ng)* &
259 & (W_stokes(i-1,j,k )*odz(i-1,j,k)*pm(i-1,j)*pn(i-1,j)+&
260 & W_stokes(i ,j,k )*odz(i ,j,k)*pm(i ,j)*pn(i ,j))
261# endif
262 END DO
263 DO k=2,N(ng)-1
264 DO i=IstrU-1,Iendp2
265 C(i,k)=0.0625_r8* &
266 & ((Ta(i ,j,k+1)-Ta(i ,j,k ))*odz(i ,j,k )+ &
267 & (Ta(i ,j,k )-Ta(i ,j,k-1))*odz(i ,j,k-1)+ &
268 & (Ta(i-1,j,k+1)-Ta(i-1,j,k ))*odz(i-1,j,k )+ &
269 & (Ta(i-1,j,k )-Ta(i-1,j,k-1))*odz(i-1,j,k-1))* &
270 & (z_r(i ,j,k+1)-z_r(i ,j,k-1)+ &
271 & z_r(i-1,j,k+1)-z_r(i-1,j,k-1))/ &
272 & (Ta(i-1,j,k)+Ta(i,j,k)+eps)
273 Wm(i,k)=0.25_r8*dt(ng)* &
274 & ((w(i-1,j,k-1)*odz(i-1,j,k-1)+ &
275 & w(i-1,j,k )*odz(i-1,j,k ))*pm(i-1,j)*pn(i-1,j)+ &
276 & (w(i ,j,k )*odz(i ,j,k )+ &
277 & w(i ,j,k-1)*odz(i ,j,k-1))*pm(i ,j)*pn(i ,j))
278# ifdef WEC_VF
279 Wm(i,k)=Wm(i,k)+ &
280 & 0.25_r8*dt(ng)* &
281 & ((W_stokes(i-1,j,k-1)*odz(i-1,j,k-1)+ &
282 & W_stokes(i-1,j,k )*odz(i-1,j,k ))* &
283 & pm(i-1,j)*pn(i-1,j)+ &
284 & (W_stokes(i ,j,k )*odz(i ,j,k )+ &
285 & W_stokes(i ,j,k-1)*odz(i ,j,k-1))* &
286 & pm(i ,j)*pn(i ,j))
287# endif
288 END DO
289 END DO
290 k=N(ng)
291 DO i=IstrU-1,Iendp2
292 C(i,k)=0.25_r8* &
293 & ((Ta(i ,j,k )-Ta(i ,j,k-1))*odz(i ,j,k-1)+ &
294 & (Ta(i-1,j,k )-Ta(i-1,j,k-1))*odz(i-1,j,k-1))* &
295 & (z_r(i ,j,k )-z_r(i ,j,k-1)+ &
296 & z_r(i-1,j,k )-z_r(i-1,j,k-1))/ &
297 & (Ta(i-1,j,k)+Ta(i,j,k)+eps)
298 Wm(i,k)=0.25_r8*dt(ng)* &
299 & (w(i-1,j,k-1)*odz(i-1,j,k-1)*pm(i-1,j)*pn(i-1,j)+ &
300 & w(i ,j,k-1)*odz(i ,j,k-1)*pm(i ,j)*pn(i ,j))
301# ifdef WEC_VF
302 Wm(i,k)=Wm(i,k)+ &
303 & 0.25_r8*dt(ng)* &
304 & (W_stokes(i-1,j,k-1)*odz(i-1,j,k-1)* &
305 & pm(i-1,j)*pn(i-1,j)+ &
306 & W_stokes(i ,j,k-1)*odz(i ,j,k-1)* &
307 & pm(i ,j)*pn(i ,j))
308# endif
309 END DO
310 DO k=1,N(ng)
311 DO i=IstrU-1,Iendp2
312 IF ((Ta(i-1,j,k).le.0.0_r8).or. &
313 & (Ta(i ,j,k).le.0.0_r8).or. &
314 & (ABS(Ta(i-1,j,k)-Ta(i,j,k)).le.eps2)) THEN
315 Ua(i,j,k)=0.0_r8
316 ELSE
317 A=(Ta(i,j,k)-Ta(i-1,j,k))/ &
318 & (Ta(i,j,k)+Ta(i-1,j,k)+eps)
319# ifdef MASKING
320 B=0.03125_r8* &
321 & ((Ta(i ,j+1,k)-Ta(i ,j ,k))* &
322 & (pn(i ,j )+pn(i ,j+1))*vmask(i ,j+1)+ &
323 & (Ta(i ,j ,k)-Ta(i ,j-1,k))* &
324 & (pn(i ,j-1)+pn(i ,j ))*vmask(i ,j )+ &
325 & (Ta(i-1,j+1,k)-Ta(i-1,j ,k))* &
326 & (pn(i-1,j )+pn(i-1,j+1))*vmask(i-1,j+1)+ &
327 & (Ta(i-1,j ,k)-Ta(i-1,j-1,k))* &
328 & (pn(i-1,j-1)+pn(i-1,j ))*vmask(i-1,j ))
329# else
330 B=0.03125_r8* &
331 & ((Ta(i ,j+1,k)-Ta(i ,j ,k))* &
332 & (pn(i ,j )+pn(i ,j+1))+ &
333 & (Ta(i ,j ,k)-Ta(i ,j-1,k))* &
334 & (pn(i ,j-1)+pn(i ,j ))+ &
335 & (Ta(i-1,j+1,k)-Ta(i-1,j ,k))* &
336 & (pn(i-1,j )+pn(i-1,j+1))+ &
337 & (Ta(i-1,j ,k)-Ta(i-1,j-1,k))* &
338 & (pn(i-1,j-1)+pn(i-1,j )))
339# endif
340 B=B*(on_v(i ,j )+on_v(i ,j+1)+ &
341 & on_v(i-1,j )+on_v(i-1,j+1))/ &
342 & (Ta(i-1,j,k)+Ta(i,j,k)+eps)
343!
344 Um=0.125_r8*Huon(i,j,k)* &
345 & dt(ng)*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))* &
346 & (oHz(i-1,j,k)+oHz(i,j,k))
347 Vm=0.03125_r8*dt(ng)* &
348 & (Hvom(i-1,j ,k)*(pm(i-1,j)+pm(i-1,j-1))* &
349 & (pn(i-1,j)+pn(i-1,j-1))* &
350 & (oHz(i-1,j, k)+oHz(i-1,j-1,k))+ &
351 & Hvom(i-1,j+1,k)*(pm(i-1,j+1)+pm(i-1,j))* &
352 & (pn(i-1,j+1)+pn(i-1,j))* &
353 & (oHz(i-1,j+1,k)+oHz(i-1,j ,k))+ &
354 & Hvom(i ,j ,k)*(pm(i ,j)+pm(i ,j-1))* &
355 & (pn(i ,j)+pn(i ,j-1))* &
356 & (oHz(i ,j ,k)+oHz(i ,j-1,k))+ &
357 & Hvom(i ,j+1,k)*(pm(i ,j+1)+pm(i ,j))* &
358 & (pn(i ,j+1)+pn(i ,j))* &
359 & (oHz(i ,j+1,k)+oHz(i ,j ,k)))
360!
361 X=(ABS(Um)-Um*Um)*A-B*Um*Vm-C(i,k)*Um*Wm(i,k)
362 Y=(ABS(Vm)-Vm*Vm)*B-A*Um*Vm-C(i,k)*Vm*Wm(i,k)
363 Z=(ABS(Wm(i,k))-Wm(i,k)*Wm(i,k))*C(i,k)- &
364 & A*Um*Wm(i,k)-B*Vm*Wm(i,k)
365
366# ifdef MPDATA_HOT
367!
368 AA=A*A
369 BB=B*B
370 CC=C(i,k)*C(i,k)
371 AB=A*B
372 AC=A*C(i,k)
373 BC=B*C(i,k)
374
375 XX=X*X
376 YY=Y*Y
377 ZZ=Z*Z
378 XY=X*Y
379 XZ=X*Z
380 YZ=Y*Z
381# endif
382!
383 sig_alfa=1.0_r8/(1.0_r8-ABS(A)+eps)
384# ifdef MPDATA_HOT
385 sig_beta=-A/((1.0_r8-ABS(A))* &
386 & (1.0_r8-AA)+eps)
387 sig_gama=2.0_r8*ABS(AA*A)/((1.0_r8-ABS(A))* &
388 & (1.0_r8-AA)* &
389 & (1.0_r8-ABS(AA*A))+eps)
390 sig_a=-B/((1.0_r8-ABS(A))* &
391 & (1.0_r8-ABS(AB))+eps)
392 sig_b=AB/((1.0_r8-ABS(A))* &
393 & (1.0_r8-AA*ABS(B))+eps)* &
394 & (ABS(B)/(1.0_r8-ABS(AB)+eps)+ &
395 & 2.0_r8*A/(1.0_r8-AA+eps))
396 sig_c=ABS(A)*BB/((1.0_r8-ABS(A))* &
397 & (1.0_r8-BB*ABS(A))* &
398 & (1.0_r8-ABS(AB))+eps)
399 sig_d=-C(i,k)/((1.0_r8-ABS(A))* &
400 & (1.0_r8-ABS(AC))+eps)
401 sig_e=AC/((1.0_r8-ABS(A))* &
402 & (1.0_r8-AA*ABS(C(i,k)))+eps)* &
403 & (ABS(C(i,k))/(1.0_r8-ABS(AC)+eps)+ &
404 & 2.0_r8*A/(1.0_r8-AA+eps))
405 sig_f=ABS(A)*CC/((1.0_r8-ABS(A))* &
406 & (1.0_r8-CC*ABS(A))* &
407 & (1.0_r8-ABS(AC))+eps)
408
409 Ua(i,j,k)=sig_alfa*X+ &
410 & sig_beta*XX+ &
411 & sig_gama*XX*X+ &
412 & sig_a*XY+ &
413 & sig_b*XX*Y+ &
414 & sig_c*X*YY+ &
415 & sig_d*XZ+ &
416 & sig_e*XX*Z+ &
417 & sig_f*X*ZZ
418# else
419 Ua(i,j,k)=sig_alfa*X
420# endif
421!
422! Limit by physical velocity.
423!
424 Ua(i,j,k)=MIN(ABS(Ua(i,j,k)),fac*ABS(Um))* &
425 & SIGN(1.0_r8,Ua(i,j,k))
426# ifdef MASKING
427 Ua(i,j,k)=Ua(i,j,k)*umask(i,j)
428# endif
429# ifdef WET_DRY
430 Ua(i,j,k)=Ua(i,j,k)*umask_wet(i,j)
431# endif
432 END IF
433 END DO
434 END DO
435 END DO
436!
437! Compute nondimensional V-antidiffusive velocities, Va. If applicable,
438! retain up to third-order terms of the power series.
439!
440! DO j=JstrV-1,Jendp2
441 DO j=JstrVm1,Jendp2
442 k=1
443 DO i=IstrU-1,Iendp1
444 C(i,k)=0.25_r8* &
445 & ((Ta(i,j ,k+1)-Ta(i,j ,k ))*odz(i,j ,k )+ &
446 & (Ta(i,j-1,k+1)-Ta(i,j-1,k ))*odz(i,j-1,k ))* &
447 & (z_r(i,j ,k+1)-z_r(i,j ,k)+ &
448 & z_r(i,j-1,k+1)-z_r(i,j-1,k))/ &
449 & (Ta(i,j-1,k)+Ta(i,j,k)+eps)
450 Wm(i,k)=0.25_r8*dt(ng)* &
451 & (w(i,j-1,k )*odz(i,j-1,k )*pm(i,j-1)*pn(i,j-1)+ &
452 & w(i,j ,k )*odz(i,j ,k )*pm(i,j )*pn(i,j ))
453# ifdef WEC_VF
454 Wm(i,k)=Wm(i,k)+ &
455 & 0.25_r8*dt(ng)* &
456 & (W_stokes(i,j-1,k )*odz(i,j-1,k )*pm(i,j-1)*pn(i,j-1)+&
457 & W_stokes(i,j ,k )*odz(i,j ,k )*pm(i,j )*pn(i,j ))
458# endif
459 END DO
460 DO k=2,N(ng)-1
461 DO i=IstrU-1,Iendp1
462 C(i,k)=0.0625_r8* &
463 & ((Ta(i,j ,k+1)-Ta(i,j ,k ))*odz(i,j ,k )+ &
464 & (Ta(i,j ,k )-Ta(i,j ,k-1))*odz(i,j ,k-1)+ &
465 & (Ta(i,j-1,k+1)-Ta(i,j-1,k ))*odz(i,j-1,k )+ &
466 & (Ta(i,j-1,k )-Ta(i,j-1,k-1))*odz(i,j-1,k-1))* &
467 & (z_r(i,j ,k+1)-z_r(i,j ,k-1)+ &
468 & z_r(i,j-1,k+1)-z_r(i,j-1,k-1))/ &
469 & (Ta(i,j-1,k)+Ta(i,j,k)+eps)
470 Wm(i,k)=0.25_r8*dt(ng)* &
471 & ((w(i,j-1,k-1)*odz(i,j-1,k-1)+ &
472 & w(i,j-1,k )*odz(i,j-1,k ))*pm(i,j-1)*pn(i,j-1)+ &
473 & (w(i,j ,k )*odz(i,j ,k )+ &
474 & w(i,j ,k-1)*odz(i,j ,k-1))*pm(i,j )*pn(i,j ))
475# ifdef WEC_VF
476 Wm(i,k)=Wm(i,k)+ &
477 & 0.25_r8*dt(ng)* &
478 & ((W_stokes(i,j-1,k-1)*odz(i,j-1,k-1)+ &
479 & W_stokes(i,j-1,k )*odz(i,j-1,k ))* &
480 & pm(i,j-1)*pn(i,j-1)+ &
481 & (W_stokes(i,j ,k )*odz(i,j ,k )+ &
482 & W_stokes(i,j ,k-1)*odz(i,j ,k-1))* &
483 & pm(i,j )*pn(i,j ))
484# endif
485 END DO
486 END DO
487 k=N(ng)
488 DO i=IstrU-1,Iendp1
489 C(i,k)=0.25_r8* &
490 & ((Ta(i,j ,k )-Ta(i,j ,k-1))*odz(i,j ,k-1)+ &
491 & (Ta(i,j-1,k )-Ta(i,j-1,k-1))*odz(i,j-1,k-1))* &
492 & (z_r(i,j ,k )-z_r(i,j ,k-1)+ &
493 & z_r(i,j-1,k )-z_r(i,j-1,k-1))/ &
494 & (Ta(i,j-1,k)+Ta(i,j,k)+eps)
495 Wm(i,k)=0.25_r8*dt(ng)* &
496 & (w(i,j-1,k-1)*odz(i,j-1,k-1)*pm(i,j-1)*pn(i,j-1)+ &
497 & w(i,j ,k-1)*odz(i,j ,k-1)*pm(i,j )*pn(i,j ))
498# ifdef WEC_VF
499 Wm(i,k)=Wm(i,k)+ &
500 & 0.25_r8*dt(ng)* &
501 & (W_stokes(i,j-1,k-1)*odz(i,j-1,k-1)* &
502 & pm(i,j-1)*pn(i,j-1)+ &
503 & W_stokes(i,j ,k-1)*odz(i,j ,k-1)* &
504 & pm(i,j )*pn(i,j ))
505# endif
506 END DO
507 DO k=1,N(ng)
508 DO i=IstrU-1,Iendp1
509 IF ((Ta(i,j-1,k).le.0.0_r8).or. &
510 & (Ta(i,j ,k).le.0.0_r8).or. &
511 & (ABS(Ta(i,j-1,k)-Ta(i,j,k)).le.eps2)) THEN
512 Va(i,j,k)=0.0_r8
513 ELSE
514# ifdef MASKING
515 A=0.03125_r8* &
516 & ((Ta(i+1,j ,k)-Ta(i ,j ,k))* &
517 & (pm(i+1,j )+pm(i ,j ))*umask(i+1,j )+ &
518 & (Ta(i ,j ,k)-Ta(i-1,j ,k))* &
519 & (pm(i-1,j )+pm(i ,j ))*umask(i ,j )+ &
520 & (Ta(i+1,j-1,k)-Ta(i ,j-1,k))* &
521 & (pm(i+1,j-1)+pm(i ,j-1))*umask(i+1,j-1)+ &
522 & (Ta(i ,j-1,k)-Ta(i-1,j-1,k))* &
523 & (pm(i-1,j-1)+pm(i ,j-1))*umask(i ,j-1))
524# else
525 A=0.03125_r8* &
526 & ((Ta(i+1,j ,k)-Ta(i ,j ,k))* &
527 & (pm(i+1,j )+pm(i ,j ))+ &
528 & (Ta(i ,j ,k)-Ta(i-1,j ,k))* &
529 & (pm(i-1,j )+pm(i ,j ))+ &
530 & (Ta(i+1,j-1,k)-Ta(i ,j-1,k))* &
531 & (pm(i+1,j-1)+pm(i ,j-1))+ &
532 & (Ta(i ,j-1,k)-Ta(i-1,j-1,k))* &
533 & (pm(i-1,j-1)+pm(i ,j-1)))
534# endif
535 A=A*(om_u(i ,j )+om_u(i+1,j )+ &
536 & om_u(i ,j-1)+om_u(i+1,j-1))/ &
537 & (Ta(i ,j-1,k)+Ta(i,j,k)+eps)
538
539 B=(Ta(i,j,k)-Ta(i,j-1,k))/ &
540 & (Ta(i,j,k)+Ta(i,j-1,k)+eps)
541!
542 Um=0.03125_r8*dt(ng)* &
543 & (Huon(i+1,j ,k)*(pm(i+1,j)+pm(i,j))* &
544 & (pn(i+1,j)+pn(i,j))* &
545 & (oHz(i+1,j ,k)+oHz(i,j ,k))+ &
546 & Huon(i+1,j-1,k)*(pm(i+1,j-1)+pm(i,j-1))* &
547 & (pn(i+1,j-1)+pn(i,j-1))* &
548 & (oHz(i+1,j-1,k)+oHz(i,j-1,k))+ &
549 & Huon(i ,j ,k)*(pm(i-1,j)+pm(i,j))* &
550 & (pn(i-1,j)+pn(i,j))* &
551 & (oHz(i-1,j ,k)+oHz(i,j ,k))+ &
552 & Huon(i ,j-1,k)*(pm(i-1,j-1)+pm(i,j-1))* &
553 & (pn(i-1,j-1)+pn(i,j-1))* &
554 & (oHz(i-1,j-1,k)+oHz(i,j-1,k)))
555 Vm=0.125_r8*Hvom(i,j,k)* &
556 & dt(ng)*(pn(i,j-1)+pn(i,j))*(pm(i,j-1)+pm(i,j))* &
557 & (oHz(i,j-1,k)+oHz(i,j,k))
558!
559 X=(ABS(Um)-Um*Um)*A-B*Um*Vm-C(i,k)*Um*Wm(i,k)
560 Y=(ABS(Vm)-Vm*Vm)*B-A*Um*Vm-C(i,k)*Vm*Wm(i,k)
561 Z=(ABS(Wm(i,k))-Wm(i,k)*Wm(i,k))*C(i,k)- &
562 & A*Um*Wm(i,k)-B*Vm*Wm(i,k)
563
564# ifdef MPDATA_HOT
565!
566 AA=A*A
567 BB=B*B
568 CC=C(i,k)*C(i,k)
569 AB=A*B
570 AC=A*C(i,k)
571 BC=B*C(i,k)
572
573 XX=X*X
574 YY=Y*Y
575 ZZ=Z*Z
576 XY=X*Y
577 XZ=X*Z
578 YZ=Y*Z
579# endif
580!
581 sig_alfa=1.0_r8/(1.0_r8-ABS(B)+eps)
582# ifdef MPDATA_HOT
583 sig_beta=-B/((1.0_r8-ABS(B))* &
584 & (1.0_r8-BB)+eps)
585 sig_gama=2.0_r8*ABS(BB*B)/((1.0_r8-ABS(B))* &
586 & (1.0_r8-BB)* &
587 & (1.0_r8-ABS(BB*B))+eps)
588 sig_a=-A/((1.0_r8-ABS(B))* &
589 & (1.0_r8-ABS(AB))+eps)
590 sig_b=AB/((1.0_r8-ABS(B))* &
591 & (1.0_r8-BB*ABS(A))+eps)* &
592 & (ABS(A)/(1.0_r8-ABS(AB)+eps)+ &
593 & 2.0_r8*B/(1.0_r8-BB+eps))
594 sig_c=ABS(B)*AA/((1.0_r8-ABS(B))* &
595 & (1.0_r8-AA*ABS(B))* &
596 & (1.0_r8-ABS(AB))+eps)
597 sig_d=-C(i,k)/((1.0_r8-ABS(B))* &
598 & (1.0_r8-ABS(BC))+eps)
599 sig_e=BC/((1.0_r8-ABS(B))* &
600 & (1.0_r8-BB*ABS(C(i,k)))+eps)* &
601 & (ABS(C(i,k))/(1.0_r8-ABS(BC)+eps)+ &
602 & 2.0_r8*B/(1.0_r8-BB+eps))
603 sig_f=ABS(B)*CC/((1.0_r8-ABS(B))* &
604 & (1.0_r8-CC*ABS(B))* &
605 & (1.0_r8-ABS(BC))+eps)
606
607 Va(i,j,k)=sig_alfa*Y+ &
608 & sig_beta*YY+ &
609 & sig_gama*YY*Y+ &
610 & sig_a*XY+ &
611 & sig_b*Y*XX+ &
612 & sig_c*YY*X+ &
613 & sig_d*YZ+ &
614 & sig_e*YY*Z+ &
615 & sig_f*Y*ZZ
616# else
617 Va(i,j,k)=sig_alfa*Y
618# endif
619!
620! Limit by physical velocity.
621!
622 Va(i,j,k)=MIN(ABS(Va(i,j,k)),fac*ABS(Vm))* &
623 & SIGN(1.0_r8,Va(i,j,k))
624# ifdef MASKING
625 Va(i,j,k)=Va(i,j,k)*vmask(i,j)
626# endif
627# ifdef WET_DRY
628 Va(i,j,k)=Va(i,j,k)*vmask_wet(i,j)
629# endif
630 END IF
631 END DO
632 END DO
633 END DO
634!
635! Apply boundary conditions to anti-diffusive velocities.
636!
637 IF (.not.EWperiodic(ng)) THEN
638 IF (DOMAIN(ng)%Western_Edge(tile)) THEN
639 IF (LBC(iwest,isBu3d,ng)%closed) THEN
640 DO k=1,N(ng)
641! DO j=Jstr,Jend
642 DO j=Jstrm1,Jendp1
643 Ua(Istr,j,k)=0.0_r8
644 END DO
645 END DO
646 ELSE
647 DO k=1,N(ng)
648! DO j=Jstr,Jend
649 DO j=Jstrm1,Jendp1
650 Ua(Istr,j,k)=Ua(Istr+1,j,k)
651 END DO
652 END DO
653 END IF
654 END IF
655
656 IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
657 IF (LBC(ieast,isBu3d,ng)%closed) THEN
658 DO k=1,N(ng)
659! DO j=Jstr,Jend
660 DO j=Jstrm1,Jendp1
661 Ua(Iend+1,j,k)=0.0_r8
662 END DO
663 END DO
664 ELSE
665 DO k=1,N(ng)
666! DO j=Jstr,Jend
667 DO j=Jstrm1,Jendp1
668 Ua(Iend+1,j,k)=Ua(Iend,j,k)
669 END DO
670 END DO
671 END IF
672 END IF
673 END IF
674
675 IF (.not.NSperiodic(ng)) THEN
676 IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
677 IF (LBC(isouth,isBv3d,ng)%closed) THEN
678 DO k=1,N(ng)
679! DO i=Istr,Iend
680 DO i=Istrm1,Iendp1
681 Va(i,Jstr,k)=0.0_r8
682 END DO
683 END DO
684 ELSE
685 DO k=1,N(ng)
686! DO i=Istr,Iend
687 DO i=Istrm1,Iendp1
688 Va(i,Jstr,k)=Va(i,Jstr+1,k)
689 END DO
690 END DO
691 END IF
692 END IF
693
694 IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
695 IF (LBC(inorth,isBv3d,ng)%closed) THEN
696 DO k=1,N(ng)
697! DO i=Istr,Iend
698 DO i=Istrm1,Iendp1
699 Va(i,Jend+1,k)=0.0_r8
700 END DO
701 END DO
702 ELSE
703 DO k=1,N(ng)
704! DO i=Istr,Iend
705 DO i=Istrm1,Iendp1
706 Va(i,Jend+1,k)=Va(i,Jend,k)
707 END DO
708 END DO
709 END IF
710 END IF
711 END IF
712
713!
714! Compute nondimensional W-antidiffusive velocities, Wa. If applicable,
715! retain up to third-order terms of the power series.
716!
717 DO j=JstrV-1,Jendp1
718 DO k=1,N(ng)-1
719 DO i=IstrU-1,Iendp1
720 IF ((Ta(i,j,k ).le.0.0_r8).or. &
721 & (Ta(i,j,k+1).le.0.0_r8).or. &
722 & (ABS(Ta(i,j,k)-Ta(i,j,k+1)).le.eps2)) THEN
723 Wa(i,j,k)=0.0_r8
724 ELSE
725 C(i,k)=(Ta(i,j,k+1)-Ta(i,j,k))/ &
726 & (Ta(i,j,k+1)+Ta(i,j,k)+eps)
727# ifdef MASKING
728 A=0.0625_r8* &
729 & ((Ta(i+1,j,k+1)-Ta(i ,j,k+1))* &
730 & (pm(i+1,j )+pm(i ,j ))*umask(i+1,j)+ &
731 & (Ta(i ,j,k+1)-Ta(i-1,j,k+1))* &
732 & (pm(i ,j )+pm(i-1,j ))*umask(i ,j)+ &
733 & (Ta(i+1,j,k )-Ta(i ,j,k ))* &
734 & (pm(i+1,j )+pm(i ,j ))*umask(i+1,j)+ &
735 & (Ta(i ,j,k )-Ta(i-1,j,k ))* &
736 & (pm(i ,j )+pm(i-1,j ))*umask(i ,j))
737 B=0.0625_r8* &
738 & ((Ta(i,j+1,k+1)-Ta(i,j ,k+1))* &
739 & (pn(i ,j+1)+pn(i ,j ))*vmask(i,j+1)+ &
740 & (Ta(i,j ,k+1)-Ta(i,j-1,k+1))* &
741 & (pn(i ,j )+pn(i ,j-1))*vmask(i,j )+ &
742 & (Ta(i,j+1,k )-Ta(i,j ,k ))* &
743 & (pn(i ,j+1)+pn(i ,j ))*vmask(i,j+1)+ &
744 & (Ta(i,j ,k )-Ta(i,j-1,k ))* &
745 & (pn(i ,j )+pn(i ,j-1))*vmask(i,j ))
746# else
747 A=0.0625_r8* &
748 & ((Ta(i+1,j,k+1)-Ta(i ,j,k+1))* &
749 & (pm(i+1,j )+pm(i ,j ))+ &
750 & (Ta(i ,j,k+1)-Ta(i-1,j,k+1))* &
751 & (pm(i ,j )+pm(i-1,j ))+ &
752 & (Ta(i+1,j,k )-Ta(i ,j,k ))* &
753 & (pm(i+1,j )+pm(i ,j ))+ &
754 & (Ta(i ,j,k )-Ta(i-1,j,k ))* &
755 & (pm(i ,j )+pm(i-1,j )))
756 B=0.0625_r8* &
757 & ((Ta(i,j+1,k+1)-Ta(i,j ,k+1))* &
758 & (pn(i ,j+1)+pn(i ,j ))+ &
759 & (Ta(i,j ,k+1)-Ta(i,j-1,k+1))* &
760 & (pn(i ,j )+pn(i ,j-1))+ &
761 & (Ta(i,j+1,k )-Ta(i,j ,k ))* &
762 & (pn(i ,j+1)+pn(i ,j ))+ &
763 & (Ta(i,j ,k )-Ta(i,j-1,k ))* &
764 & (pn(i ,j )+pn(i ,j-1)))
765# endif
766 A=A*(om_u(i+1,j)+om_u(i ,j))/ &
767 & (Ta(i,j,k+1)+Ta(i,j,k)+eps)
768 B=B*(on_v(i,j+1)+on_v(i,j ))/ &
769 & (Ta(i,j,k+1)+Ta(i,j,k)+eps)
770!
771 Um=0.03125_r8*dt(ng)* &
772 & (Huon(i ,j,k )*(pm(i,j)+pm(i-1,j))* &
773 & (pn(i,j)+pn(i-1,j))* &
774 & (oHz(i,j,k )+oHz(i-1,j,k ))+ &
775 & Huon(i ,j,k+1)*(pm(i,j)+pm(i-1,j))* &
776 & (pn(i,j)+pn(i-1,j))* &
777 & (oHz(i,j,k+1)+oHz(i-1,j,k+1))+ &
778 & Huon(i+1,j,k )*(pm(i,j)+pm(i+1,j))* &
779 & (pn(i,j)+pn(i+1,j))* &
780 & (oHz(i,j,k )+oHz(i+1,j,k ))+ &
781 & Huon(i+1,j,k+1)*(pm(i,j)+pm(i+1,j))* &
782 & (pn(i,j)+pn(i+1,j))* &
783 & (oHz(i,j,k+1)+oHz(i+1,j,k+1)))
784 Vm=0.03125_r8*dt(ng)* &
785 & (Hvom(i,j ,k )*(pm(i,j)+pm(i,j-1))* &
786 & (pn(i,j)+pn(i,j-1))* &
787 & (oHz(i,j,k )+oHz(i,j-1,k ))+ &
788 & Hvom(i,j ,k+1)*(pm(i,j)+pm(i,j-1))* &
789 & (pn(i,j)+pn(i,j-1))* &
790 & (oHz(i,j,k+1)+oHz(i,j-1,k+1))+ &
791 & Hvom(i,j+1,k )*(pm(i,j)+pm(i,j+1))* &
792 & (pn(i,j)+pn(i,j+1))* &
793 & (oHz(i,j,k )+oHz(i,j+1,k ))+ &
794 & Hvom(i,j+1,k+1)*(pm(i,j)+pm(i,j+1))* &
795 & (pn(i,j)+pn(i,j+1))* &
796 & (oHz(i,j,k+1)+oHz(i,j+1,k+1)))
797
798 Wm(i,k)=w(i,j,k)*odz(i,j,k)*pm(i,j)*pn(i,j)*dt(ng)
799# ifdef WEC_VF
800 Wm(i,k)=Wm(i,k)+W_stokes(i,j,k)*odz(i,j,k)* &
801 & pm(i,j)*pn(i,j)*dt(ng)
802# endif
803!
804 X=(ABS(Um)-Um*Um)*A-B*Um*Vm-C(i,k)*Um*Wm(i,k)
805 Y=(ABS(Vm)-Vm*Vm)*B-A*Um*Vm-C(i,k)*Vm*Wm(i,k)
806 Z=(ABS(Wm(i,k))-Wm(i,k)*Wm(i,k))*C(i,k)- &
807 & A*Um*Wm(i,k)-B*Vm*Wm(i,k)
808
809# ifdef MPDATA_HOT
810!
811 AA=A*A
812 BB=B*B
813 CC=C(i,k)*C(i,k)
814 AB=A*B
815 AC=A*C(i,k)
816 BC=B*C(i,k)
817
818 XX=X*X
819 YY=Y*Y
820 ZZ=Z*Z
821 XY=X*Y
822 XZ=X*Z
823 YZ=Y*Z
824# endif
825!
826 sig_alfa=1.0_r8/(1.0_r8-ABS(C(i,k))+eps)
827# ifdef MPDATA_HOT
828 sig_beta=-C(i,k)/((1.0_r8-ABS(C(i,k)))* &
829 & (1.0_r8-CC)+eps)
830 sig_gama=2.0_r8*ABS(CC*C(i,k))/ &
831 & ((1.0_r8-ABS(C(i,k)))* &
832 & (1.0_r8-CC)* &
833 & (1.0_r8-ABS(CC*C(i,k)))+eps)
834 sig_a=-B/((1.0_r8-ABS(C(i,k)))* &
835 & (1.0_r8-ABS(BC))+eps)
836 sig_b=BC/((1.0_r8-ABS(C(i,k)))* &
837 & (1.0_r8-CC*ABS(B))+eps)* &
838 & (ABS(B)/(1.0_r8-ABS(BC)+eps)+ &
839 & 2.0_r8*C(i,k)/(1.0_r8-CC+eps))
840 sig_c=ABS(C(i,k))*BB/((1.0_r8-ABS(C(i,k)))* &
841 & (1.0_r8-B*B*ABS(C(i,k)))* &
842 & (1.0_r8-ABS(BC))+eps)
843 sig_d=-A/((1.0_r8-ABS(C(i,k)))* &
844 & (1.0_r8-ABS(AC))+eps)
845 sig_e=AC/((1.0_r8-ABS(C(i,k)))* &
846 & (1.0_r8-CC*ABS(A))+eps)* &
847 & (ABS(A)/(1.0_r8-ABS(AC)+eps)+ &
848 & 2.0_r8*C(i,k)/(1.0_r8-CC+eps))
849 sig_f=ABS(C(i,k))*AA/((1.0_r8-ABS(C(i,k)))* &
850 & (1.0_r8-AA*ABS(C(i,k)))* &
851 & (1.0_r8-ABS(AC))+eps)
852
853 Wa(i,j,k)=sig_alfa*Z+ &
854 & sig_beta*ZZ+ &
855 & sig_gama*ZZ*Z+ &
856 & sig_a*YZ+ &
857 & sig_b*ZZ*Y+ &
858 & sig_c*Z*YY+ &
859 & sig_d*XZ+ &
860 & sig_e*ZZ*X+ &
861 & sig_f*Z*XX
862# else
863 Wa(i,j,k)=sig_alfa*Z
864# endif
865!
866! Limit by physical velocity.
867!
868 Wa(i,j,k)=MIN(ABS(Wa(i,j,k)),fac*ABS(Wm(i,k)))* &
869 & SIGN(1.0_r8,Wa(i,j,k))
870# ifdef MASKING
871 Wa(i,j,k)=Wa(i,j,k)*rmask(i,j)
872# endif
873# ifdef WET_DRY
874 Wa(i,j,k)=Wa(i,j,k)*rmask_wet(i,j)
875# endif
876 END IF
877 END DO
878 END DO
879 DO i=IstrU-1,Iendp1
880 Wa(i,j,0)=0.0_r8
881 Wa(i,j,N(ng))=0.0_r8
882 END DO
883 END DO
884!
885!-----------------------------------------------------------------------
886! Supress false oscillations in the solution by imposing appropriate
887! limits on the transport fluxes. Compute the UP and DOWN beta-ratios
888! described in Smolarkiewicz and Grabowski (1990).
889!-----------------------------------------------------------------------
890!
891! Build special mask array used to limit the UP and DOWN beta-ratios
892! to avoid including land/sea masking when computing limiting Tmin
893! and Tmax values. Notice that a zero Tmin due to land mask is not
894! considered.
895!
896 DO j=Jstrm2,Jendp2
897 DO i=Istrm2,Iendp2
898# ifdef MASKING
899 mask_up(i,j)=rmask(i,j)
900 mask_dn(i,j)=MAX(1.0_r8,MIN(Large,(1.0_r8-rmask(i,j))*Large))
901# else
902 mask_up(i,j)=1.0_r8
903 mask_dn(i,j)=1.0_r8
904# endif
905 END DO
906 END DO
907!
908! Compute UP and DOWN beta-ratios.
909!
910 DO j=JstrV-1,Jendp1
911 k=1
912 DO i=IstrU-1,Iendp1
913 Tmax=MAX(Ta(i-1,j ,k )*mask_up(i-1,j ), &
914 & t (i-1,j ,k )*mask_up(i-1,j ), &
915 & Ta(i ,j ,k )*mask_up(i ,j ), &
916 & t (i ,j ,k )*mask_up(i ,j ), &
917 & Ta(i+1,j ,k )*mask_up(i+1,j ), &
918 & t (i+1,j ,k )*mask_up(i+1,j ), &
919 & Ta(i ,j-1,k )*mask_up(i ,j-1), &
920 & t (i ,j-1,k )*mask_up(i ,j-1), &
921 & Ta(i ,j+1,k )*mask_up(i ,j+1), &
922 & t (i ,j+1,k )*mask_up(i ,j+1), &
923 & Ta(i ,j ,k+1)*mask_up(i ,j ), &
924 & t (i ,j ,k+1)*mask_up(i ,j ))
925 cff1=Ta(i-1,j ,k )*MAX(0.0_r8,Ua(i ,j ,k ))- &
926 & Ta(i+1,j ,k )*MIN(0.0_r8,Ua(i+1,j ,k ))+ &
927 & Ta(i ,j-1,k )*MAX(0.0_r8,Va(i ,j ,k ))- &
928 & Ta(i ,j+1,k )*MIN(0.0_r8,Va(i ,j+1,k ))- &
929 & Ta(i ,j ,k+1)*MIN(0.0_r8,Wa(i ,j ,k ))
930 beta_up(i,j,k)=(Tmax-Ta(i,j,k))/(cff1+eps)
931!
932 Tmin=MIN(Ta(i-1,j ,k )*mask_dn(i-1,j ), &
933 & t (i-1,j ,k )*mask_dn(i-1,j ), &
934 & Ta(i ,j ,k )*mask_dn(i ,j ), &
935 & t (i ,j ,k )*mask_dn(i ,j ), &
936 & Ta(i+1,j ,k )*mask_dn(i+1,j ), &
937 & t (i+1,j ,k )*mask_dn(i+1,j ), &
938 & Ta(i ,j-1,k )*mask_dn(i ,j-1), &
939 & t (i ,j-1,k )*mask_dn(i ,j-1), &
940 & Ta(i ,j+1,k )*mask_dn(i ,j+1), &
941 & t (i ,j+1,k )*mask_dn(i ,j+1), &
942 & Ta(i ,j ,k+1)*mask_dn(i ,j ), &
943 & t (i ,j ,k+1)*mask_dn(i ,j ))
944 cff2=Ta(i ,j ,k )*MAX(0.0_r8,Ua(i+1,j ,k ))- &
945 & Ta(i ,j ,k )*MIN(0.0_r8,Ua(i ,j ,k ))+ &
946 & Ta(i ,j ,k )*MAX(0.0_r8,Va(i ,j+1,k ))- &
947 & Ta(i ,j ,k )*MIN(0.0_r8,Va(i ,j ,k ))+ &
948 & Ta(i ,j ,k )*MAX(0.0_r8,Wa(i ,j ,k ))
949 beta_dn(i,j,k)=(Ta(i,j,k)-Tmin)/(cff2+eps)
950 END DO
951 DO k=2,N(ng)-1
952 DO i=IstrU-1,Iendp1
953 Tmax=MAX(Ta(i-1,j ,k )*mask_up(i-1,j ), &
954 & t (i-1,j ,k )*mask_up(i-1,j ), &
955 & Ta(i ,j ,k )*mask_up(i ,j ), &
956 & t (i ,j ,k )*mask_up(i ,j ), &
957 & Ta(i+1,j ,k )*mask_up(i+1,j ), &
958 & t (i+1,j ,k )*mask_up(i+1,j ), &
959 & Ta(i ,j-1,k )*mask_up(i ,j-1), &
960 & t (i ,j-1,k )*mask_up(i ,j-1), &
961 & Ta(i ,j+1,k )*mask_up(i ,j+1), &
962 & t (i ,j+1,k )*mask_up(i ,j+1), &
963 & Ta(i ,j ,k-1)*mask_up(i ,j ), &
964 & t (i ,j ,k-1)*mask_up(i ,j ), &
965 & Ta(i ,j ,k+1)*mask_up(i ,j ), &
966 & t (i ,j ,k+1)*mask_up(i ,j ))
967 cff1=Ta(i-1,j ,k )*MAX(0.0_r8,Ua(i ,j ,k ))- &
968 & Ta(i+1,j ,k )*MIN(0.0_r8,Ua(i+1,j ,k ))+ &
969 & Ta(i ,j-1,k )*MAX(0.0_r8,Va(i ,j ,k ))- &
970 & Ta(i ,j+1,k )*MIN(0.0_r8,Va(i ,j+1,k ))+ &
971 & Ta(i ,j ,k-1)*MAX(0.0_r8,Wa(i ,j ,k-1))- &
972 & Ta(i ,j ,k+1)*MIN(0.0_r8,Wa(i ,j ,k ))
973 beta_up(i,j,k)=(Tmax-Ta(i,j,k))/(cff1+eps)
974!
975 Tmin=MIN(Ta(i-1,j ,k )*mask_dn(i-1,j ), &
976 & t (i-1,j ,k )*mask_dn(i-1,j ), &
977 & Ta(i ,j ,k )*mask_dn(i ,j ), &
978 & t (i ,j ,k )*mask_dn(i ,j ), &
979 & Ta(i+1,j ,k )*mask_dn(i+1,j ), &
980 & t (i+1,j ,k )*mask_dn(i+1,j ), &
981 & Ta(i ,j-1,k )*mask_dn(i ,j-1), &
982 & t (i ,j-1,k )*mask_dn(i ,j-1), &
983 & Ta(i ,j+1,k )*mask_dn(i ,j+1), &
984 & t (i ,j+1,k )*mask_dn(i ,j+1), &
985 & Ta(i ,j ,k-1)*mask_dn(i ,j ), &
986 & t (i ,j ,k-1)*mask_dn(i ,j ), &
987 & Ta(i ,j ,k+1)*mask_dn(i ,j ), &
988 & t (i ,j ,k+1)*mask_dn(i ,j ))
989 cff2=Ta(i ,j ,k )*MAX(0.0_r8,Ua(i+1,j ,k ))- &
990 & Ta(i ,j ,k )*MIN(0.0_r8,Ua(i ,j ,k ))+ &
991 & Ta(i ,j ,k )*MAX(0.0_r8,Va(i ,j+1,k ))- &
992 & Ta(i ,j ,k )*MIN(0.0_r8,Va(i ,j ,k ))+ &
993 & Ta(i ,j ,k )*MAX(0.0_r8,Wa(i ,j ,k ))- &
994 & Ta(i ,j ,k )*MIN(0.0_r8,Wa(i ,j ,k-1))
995 beta_dn(i,j,k)=(Ta(i,j,k)-Tmin)/(cff2+eps)
996 END DO
997 END DO
998 k=N(ng)
999 DO i=IstrU-1,Iendp1
1000 Tmax=MAX(Ta(i-1,j ,k )*mask_up(i-1,j ), &
1001 & t (i-1,j ,k )*mask_up(i-1,j ), &
1002 & Ta(i ,j ,k )*mask_up(i ,j ), &
1003 & t (i ,j ,k )*mask_up(i ,j ), &
1004 & Ta(i+1,j ,k )*mask_up(i+1,j ), &
1005 & t (i+1,j ,k )*mask_up(i+1,j ), &
1006 & Ta(i ,j-1,k )*mask_up(i ,j-1), &
1007 & t (i ,j-1,k )*mask_up(i ,j-1), &
1008 & Ta(i ,j+1,k )*mask_up(i ,j+1), &
1009 & t (i ,j+1,k )*mask_up(i ,j+1), &
1010 & Ta(i ,j ,k-1)*mask_up(i ,j ), &
1011 & t (i ,j ,k-1)*mask_up(i ,j ))
1012 cff1=Ta(i-1,j ,k )*MAX(0.0_r8,Ua(i ,j ,k ))- &
1013 & Ta(i+1,j ,k )*MIN(0.0_r8,Ua(i+1,j ,k ))+ &
1014 & Ta(i ,j-1,k )*MAX(0.0_r8,Va(i ,j ,k ))- &
1015 & Ta(i ,j+1,k )*MIN(0.0_r8,Va(i ,j+1,k ))+ &
1016 & Ta(i ,j ,k-1)*MAX(0.0_r8,Wa(i ,j ,k-1))
1017 beta_up(i,j,k)=(Tmax-Ta(i,j,k))/(cff1+eps)
1018!
1019 Tmin=MIN(Ta(i-1,j ,k )*mask_dn(i-1,j ), &
1020 & t (i-1,j ,k )*mask_dn(i-1,j ), &
1021 & Ta(i ,j ,k )*mask_dn(i ,j ), &
1022 & t (i ,j ,k )*mask_dn(i ,j ), &
1023 & Ta(i+1,j ,k )*mask_dn(i+1,j ), &
1024 & t (i+1,j ,k )*mask_dn(i+1,j ), &
1025 & Ta(i ,j-1,k )*mask_dn(i ,j-1), &
1026 & t (i ,j-1,k )*mask_dn(i ,j-1), &
1027 & Ta(i ,j+1,k )*mask_dn(i ,j+1), &
1028 & t (i ,j+1,k )*mask_dn(i ,j+1), &
1029 & Ta(i ,j ,k-1)*mask_dn(i ,j ), &
1030 & t (i ,j ,k-1)*mask_dn(i ,j ))
1031 cff2=Ta(i ,j ,k )*MAX(0.0_r8,Ua(i+1,j ,k ))- &
1032 & Ta(i ,j ,k )*MIN(0.0_r8,Ua(i ,j ,k ))+ &
1033 & Ta(i ,j ,k )*MAX(0.0_r8,Va(i ,j+1,k ))- &
1034 & Ta(i ,j ,k )*MIN(0.0_r8,Va(i ,j ,k ))- &
1035 & Ta(i ,j ,k )*MIN(0.0_r8,Wa(i ,j ,k-1))
1036 beta_dn(i,j,k)=(Ta(i,j,k)-Tmin)/(cff2+eps)
1037 END DO
1038 END DO
1039 DO k=1,N(ng)
1040 DO j=JstrV-1,Jendp1
1041 DO i=IstrU-1,Iendp1
1042 IF (mask_up(i,j).eq.0.0_r8) THEN
1043 beta_up(i,j,k)=2.0_r8
1044 beta_dn(i,j,k)=2.0_r8
1045 END IF
1046 END DO
1047 END DO
1048 END DO
1049!
1050! Calculate monotonic velocities. Scale back to dimensional units.
1051!
1052 DO k=1,N(ng)
1053 DO j=Jstr,Jend
1054 DO i=IstrU,Iendp1
1055 cff1=MIN(beta_dn(i-1,j,k),beta_up(i,j,k),1.0_r8)
1056 cff2=MIN(beta_up(i-1,j,k),beta_dn(i,j,k),1.0_r8)
1057 Ua(i,j,k)=(cff1*MAX(0.0_r8,Ua(i,j,k))+ &
1058 & cff2*MIN(0.0_r8,Ua(i,j,k)))* &
1059 & cff*om_u(i,j)
1060# ifdef MASKING
1061 Ua(i,j,k)=Ua(i,j,k)*umask(i,j)
1062# endif
1063# ifdef WET_DRY
1064 Ua(i,j,k)=Ua(i,j,k)*umask_wet(i,j)
1065# endif
1066 END DO
1067 END DO
1068 DO j=JstrV,Jendp1
1069 DO i=Istr,Iend
1070 cff1=MIN(beta_dn(i,j-1,k),beta_up(i,j,k),1.0_r8)
1071 cff2=MIN(beta_up(i,j-1,k),beta_dn(i,j,k),1.0_r8)
1072 Va(i,j,k)=(cff1*MAX(0.0_r8,Va(i,j,k))+ &
1073 & cff2*MIN(0.0_r8,Va(i,j,k)))* &
1074 & cff*on_v(i,j)
1075# ifdef MASKING
1076 Va(i,j,k)=Va(i,j,k)*vmask(i,j)
1077# endif
1078# ifdef WET_DRY
1079 Va(i,j,k)=Va(i,j,k)*vmask_wet(i,j)
1080# endif
1081 END DO
1082 END DO
1083 IF (k.lt.N(ng)) THEN
1084 DO j=Jstr,Jend
1085 DO i=Istr,Iend
1086 cff1=MIN(beta_dn(i,j,k),beta_up(i,j,k+1),1.0_r8)
1087 cff2=MIN(beta_up(i,j,k),beta_dn(i,j,k+1),1.0_r8)
1088 Wa(i,j,k)=(cff1*MAX(0.0_r8,Wa(i,j,k))+ &
1089 & cff2*MIN(0.0_r8,Wa(i,j,k)))* &
1090 & cff*omn(i,j)*(z_r(i,j,k+1)-z_r(i,j,k))
1091# ifdef MASKING
1092 Wa(i,j,k)=Wa(i,j,k)*rmask(i,j)
1093# endif
1094# ifdef WET_DRY
1095 Wa(i,j,k)=Wa(i,j,k)*rmask_wet(i,j)
1096# endif
1097 END DO
1098 END DO
1099 END IF
1100 END DO
1101!
1102! Apply boundary conditions to anti-diffusive velocities.
1103!
1104 IF (.not.EWperiodic(ng)) THEN
1105 IF (DOMAIN(ng)%Western_Edge(tile)) THEN
1106 IF (LBC(iwest,isBu3d,ng)%closed) THEN
1107 DO k=1,N(ng)
1108 DO j=Jstr,Jend
1109 Ua(Istr,j,k)=0.0_r8
1110 END DO
1111 END DO
1112 ELSE
1113 DO k=1,N(ng)
1114 DO j=Jstr,Jend
1115 Ua(Istr,j,k)=Ua(Istr+1,j,k)
1116 END DO
1117 END DO
1118 END IF
1119 END IF
1120
1121 IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
1122 IF (LBC(ieast,isBu3d,ng)%closed) THEN
1123 DO k=1,N(ng)
1124 DO j=Jstr,Jend
1125 Ua(Iend+1,j,k)=0.0_r8
1126 END DO
1127 END DO
1128 ELSE
1129 DO k=1,N(ng)
1130 DO j=Jstr,Jend
1131 Ua(Iend+1,j,k)=Ua(Iend,j,k)
1132 END DO
1133 END DO
1134 END IF
1135 END IF
1136 END IF
1137
1138 IF (.not.NSperiodic(ng)) THEN
1139 IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
1140 IF (LBC(isouth,isBv3d,ng)%closed) THEN
1141 DO k=1,N(ng)
1142 DO i=Istr,Iend
1143 Va(i,Jstr,k)=0.0_r8
1144 END DO
1145 END DO
1146 ELSE
1147 DO k=1,N(ng)
1148 DO i=Istr,Iend
1149 Va(i,Jstr,k)=Va(i,Jstr+1,k)
1150 END DO
1151 END DO
1152 END IF
1153 END IF
1154
1155 IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
1156 IF (LBC(inorth,isBv3d,ng)%closed) THEN
1157 DO k=1,N(ng)
1158 DO i=Istr,Iend
1159 Va(i,Jend+1,k)=0.0_r8
1160 END DO
1161 END DO
1162 ELSE
1163 DO k=1,N(ng)
1164 DO i=Istr,Iend
1165 Va(i,Jend+1,k)=Va(i,Jend,k)
1166 END DO
1167 END DO
1168 END IF
1169 END IF
1170 END IF
1171
1172 RETURN
1173 END SUBROUTINE mpdata_adiff_tile
1174#endif
1175 END MODULE mpdata_adiff_mod