ROMS
Loading...
Searching...
No Matches
Macros | Functions/Subroutines
mpdata_adiff.F File Reference
#include "cppdefs.h"
#include "set_bounds.h"
Include dependency graph for mpdata_adiff.F:

Go to the source code of this file.

Macros

#define MPDATA_HOT
 

Functions/Subroutines

program __mpdata_adiff_f__
 
subroutine mpdata_adiff_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, rmask, umask, vmask, rmask_wet, umask_wet, vmask_wet, pm, pn, omn, om_u, on_v, z_r, ohz, huon, hvom, w, t, ta, ua, va, wa)
 

Macro Definition Documentation

◆ MPDATA_HOT

#define MPDATA_HOT

Function/Subroutine Documentation

◆ __mpdata_adiff_f__()

program __mpdata_adiff_f__

Definition at line 2 of file mpdata_adiff.F.

References mpdata_adiff_tile().

Here is the call graph for this function:

◆ mpdata_adiff_tile()

subroutine __mpdata_adiff_f__::mpdata_adiff_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
real(r8), dimension(lbi:,lbj:), intent(in) rmask,
real(r8), dimension(lbi:,lbj:), intent(in) umask,
real(r8), dimension(lbi:,lbj:), intent(in) vmask,
real(r8), dimension(lbi:,lbj:), intent(in) rmask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) umask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) vmask_wet,
real(r8), dimension(lbi:,lbj:), intent(in) pm,
real(r8), dimension(lbi:,lbj:), intent(in) pn,
real(r8), dimension(lbi:,lbj:), intent(in) omn,
real(r8), dimension(lbi:,lbj:), intent(in) om_u,
real(r8), dimension(lbi:,lbj:), intent(in) on_v,
real(r8), dimension(lbi:,lbj:,:), intent(in) z_r,
real(r8), dimension(imins:,jmins:,:), intent(in) ohz,
real(r8), dimension(lbi:,lbj:,:), intent(in) huon,
real(r8), dimension(lbi:,lbj:,:), intent(in) hvom,
real(r8), dimension(lbi:,lbj:,0:), intent(in) w,
real(r8), dimension(lbi:,lbj:,:), intent(in) t,
real(r8), dimension(imins:,jmins:,:), intent(inout) ta,
real(r8), dimension(imins:,jmins:,:), intent(out) ua,
real(r8), dimension(imins:,jmins:,:), intent(out) va,
real(r8), dimension(imins:,jmins:,0:), intent(out) wa )
private

Definition at line 38 of file mpdata_adiff.F.

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

References mod_param::domain, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::ieast, mod_scalars::inorth, mod_ncparam::isbu3d, mod_ncparam::isbv3d, mod_scalars::isouth, mod_scalars::iwest, mod_scalars::large, mod_param::lbc, and mod_scalars::nsperiodic.

Referenced by __mpdata_adiff_f__(), and step3d_t_mod::step3d_t_tile().

Here is the caller graph for this function: