ROMS
Loading...
Searching...
No Matches
uv3dmix4_s.h
Go to the documentation of this file.
1 MODULE uv3dmix4_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! This subroutine computes biharmonic mixing of momentum, along !
11! constant S-surfaces, from the horizontal divergence of the !
12! stress tensor. A transverse isotropy is assumed so the stress !
13! tensor is split into vertical and horizontal subtensors. !
14! !
15! Reference: !
16! !
17! Wajsowicz, R.C, 1993: A consistent formulation of the !
18! anisotropic stress tensor for use in models of the !
19! large-scale ocean circulation, JCP, 105, 333-338. !
20! !
21! Sadourny, R. and K. Maynard, 1997: Formulations of !
22! lateral diffusion in geophysical fluid dynamics !
23! models, In Numerical Methods of Atmospheric and !
24! Oceanic Modelling. Lin, Laprise, and Ritchie, !
25! Eds., NRC Research Press, 547-556. !
26! !
27! Griffies, S.M. and R.W. Hallberg, 2000: Biharmonic !
28! friction with a Smagorinsky-like viscosity for !
29! use in large-scale eddy-permitting ocean models, !
30! Monthly Weather Rev., 128, 8, 2935-2946. !
31! !
32!=======================================================================
33!
34 implicit none
35!
36 PRIVATE
37 PUBLIC uv3dmix4
38!
39 CONTAINS
40!
41!***********************************************************************
42 SUBROUTINE uv3dmix4 (ng, tile)
43!***********************************************************************
44!
45 USE mod_param
46 USE mod_coupling
47#ifdef DIAGNOSTICS_UV
48 USE mod_diags
49#endif
50 USE mod_grid
51 USE mod_mixing
52 USE mod_ocean
53 USE mod_stepping
54!
55! Imported variable declarations.
56!
57 integer, intent(in) :: ng, tile
58!
59! Local variable declarations.
60!
61 character (len=*), parameter :: MyFile = &
62 & __FILE__
63!
64#include "tile.h"
65!
66#ifdef PROFILE
67 CALL wclock_on (ng, inlm, 32, __line__, myfile)
68#endif
69 CALL uv3dmix4_s_tile (ng, tile, &
70 & lbi, ubi, lbj, ubj, &
71 & imins, imaxs, jmins, jmaxs, &
72 & nrhs(ng), nnew(ng), &
73#ifdef MASKING
74 & grid(ng) % pmask, &
75#endif
76#ifdef WET_DRY
77 & grid(ng) % pmask_wet, &
78#endif
79 & grid(ng) % Hz, &
80 & grid(ng) % om_p, &
81 & grid(ng) % om_r, &
82 & grid(ng) % on_p, &
83 & grid(ng) % on_r, &
84 & grid(ng) % pm, &
85 & grid(ng) % pmon_p, &
86 & grid(ng) % pmon_r, &
87 & grid(ng) % pn, &
88 & grid(ng) % pnom_p, &
89 & grid(ng) % pnom_r, &
90#ifdef VISC_3DCOEF
91# ifdef UV_U3ADV_SPLIT
92 & mixing(ng) % Uvis3d_r, &
93 & mixing(ng) % Vvis3d_r, &
94# else
95 & mixing(ng) % visc3d_r, &
96# endif
97#else
98 & mixing(ng) % visc4_p, &
99 & mixing(ng) % visc4_r, &
100#endif
101#ifdef DIAGNOSTICS_UV
102 & diags(ng) % DiaRUfrc, &
103 & diags(ng) % DiaRVfrc, &
104 & diags(ng) % DiaU3wrk, &
105 & diags(ng) % DiaV3wrk, &
106#endif
107 & coupling(ng) % rufrc, &
108 & coupling(ng) % rvfrc, &
109 & ocean(ng) % u, &
110 & ocean(ng) % v)
111#ifdef PROFILE
112 CALL wclock_off (ng, inlm, 32, __line__, myfile)
113#endif
114!
115 RETURN
116 END SUBROUTINE uv3dmix4
117!
118!***********************************************************************
119 SUBROUTINE uv3dmix4_s_tile (ng, tile, &
120 & LBi, UBi, LBj, UBj, &
121 & IminS, ImaxS, JminS, JmaxS, &
122 & nrhs, nnew, &
123#ifdef MASKING
124 & pmask, &
125#endif
126#ifdef WET_DRY
127 & pmask_wet, &
128#endif
129 & Hz, &
130 & om_p, om_r, on_p, on_r, &
131 & pm, pmon_p, pmon_r, &
132 & pn, pnom_p, pnom_r, &
133#ifdef VISC_3DCOEF
134# ifdef UV_U3ADV_SPLIT
135 & Uvis3d_r, Vvis3d_r, &
136# else
137 & visc3d_r, &
138# endif
139#else
140 & visc4_p, visc4_r, &
141#endif
142#ifdef DIAGNOSTICS_UV
143 & DiaRUfrc, DiaRVfrc, &
144 & DiaU3wrk, DiaV3wrk, &
145#endif
146 & rufrc, rvfrc, u, v)
147!***********************************************************************
148!
149 USE mod_param
150 USE mod_ncparam
151 USE mod_scalars
152!
153! Imported variable declarations.
154!
155 integer, intent(in) :: ng, tile
156 integer, intent(in) :: LBi, UBi, LBj, UBj
157 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
158 integer, intent(in) :: nrhs, nnew
159
160#ifdef ASSUMED_SHAPE
161# ifdef MASKING
162 real(r8), intent(in) :: pmask(LBi:,LBj:)
163# endif
164# ifdef WET_DRY
165 real(r8), intent(in) :: pmask_wet(LBi:,LBj:)
166# endif
167 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
168 real(r8), intent(in) :: om_p(LBi:,LBj:)
169 real(r8), intent(in) :: om_r(LBi:,LBj:)
170 real(r8), intent(in) :: on_p(LBi:,LBj:)
171 real(r8), intent(in) :: on_r(LBi:,LBj:)
172 real(r8), intent(in) :: pm(LBi:,LBj:)
173 real(r8), intent(in) :: pmon_p(LBi:,LBj:)
174 real(r8), intent(in) :: pmon_r(LBi:,LBj:)
175 real(r8), intent(in) :: pn(LBi:,LBj:)
176 real(r8), intent(in) :: pnom_p(LBi:,LBj:)
177 real(r8), intent(in) :: pnom_r(LBi:,LBj:)
178# ifdef VISC_3DCOEF
179# ifdef UV_U3ADV_SPLIT
180 real(r8), intent(in) :: Uvis3d_r(LBi:,LBj:,:)
181 real(r8), intent(in) :: Vvis3d_r(LBi:,LBj:,:)
182# else
183 real(r8), intent(in) :: visc3d_r(LBi:,LBj:,:)
184# endif
185# else
186 real(r8), intent(in) :: visc4_p(LBi:,LBj:)
187 real(r8), intent(in) :: visc4_r(LBi:,LBj:)
188# endif
189# ifdef DIAGNOSTICS_UV
190 real(r8), intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
191 real(r8), intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
192 real(r8), intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
193 real(r8), intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
194# endif
195 real(r8), intent(inout) :: rufrc(LBi:,LBj:)
196 real(r8), intent(inout) :: rvfrc(LBi:,LBj:)
197 real(r8), intent(inout) :: u(LBi:,LBj:,:,:)
198 real(r8), intent(inout) :: v(LBi:,LBj:,:,:)
199#else
200
201# ifdef MASKING
202 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
203# endif
204# ifdef WET_DRY
205 real(r8), intent(in) :: pmask_wet(LBi:UBi,LBj:UBj)
206# endif
207 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
208 real(r8), intent(in) :: om_p(LBi:UBi,LBj:UBj)
209 real(r8), intent(in) :: om_r(LBi:UBi,LBj:UBj)
210 real(r8), intent(in) :: on_p(LBi:UBi,LBj:UBj)
211 real(r8), intent(in) :: on_r(LBi:UBi,LBj:UBj)
212 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
213 real(r8), intent(in) :: pmon_p(LBi:UBi,LBj:UBj)
214 real(r8), intent(in) :: pmon_r(LBi:UBi,LBj:UBj)
215 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
216 real(r8), intent(in) :: pnom_p(LBi:UBi,LBj:UBj)
217 real(r8), intent(in) :: pnom_r(LBi:UBi,LBj:UBj)
218# ifdef VISC_3DCOEF
219# ifdef UV_U3ADV_SPLIT
220 real(r8), intent(in) :: Uvis3d_r(LBi:UBi,LBj:UBj,N(ng))
221 real(r8), intent(in) :: Vvis3d_r(LBi:UBi,LBj:UBj,N(ng))
222# else
223 real(r8), intent(in) :: visc3d_r(LBi:UBi,LBj:UBj,N(ng))
224# endif
225# else
226 real(r8), intent(in) :: visc4_p(LBi:UBi,LBj:UBj)
227 real(r8), intent(in) :: visc4_r(LBi:UBi,LBj:UBj)
228# endif
229# ifdef DIAGNOSTICS_UV
230 real(r8), intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
231 real(r8), intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
232 real(r8), intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
233 real(r8), intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
234# endif
235 real(r8), intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
236 real(r8), intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
237 real(r8), intent(inout) :: u(LBi:UBi,LBj:UBj,N(ng),2)
238 real(r8), intent(inout) :: v(LBi:UBi,LBj:UBj,N(ng),2)
239#endif
240!
241! Local variable declarations.
242!
243 integer :: IminU, IminV, ImaxU, ImaxV
244 integer :: JminU, JminV, JmaxU, JmaxV
245 integer :: i, j, k
246
247 real(r8) :: cff, cff1, cff2, cff3
248#ifdef VISC_3DCOEF
249 real(r8) :: Uvis_p, Vvis_p, visc_p
250#endif
251 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LapU
252 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LapV
253 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
254 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
255 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
256 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
257
258#include "set_bounds.h"
259!
260!-----------------------------------------------------------------------
261! Compute horizontal biharmonic viscosity along constant S-surfaces.
262! The biharmonic operator is computed by applying the harmonic
263! operator twice.
264!-----------------------------------------------------------------------
265!
266! Set local I- and J-ranges.
267!
268 IF (ewperiodic(ng)) THEN
269 iminu=istr-1
270 imaxu=iend+1
271 iminv=istr-1
272 imaxv=iend+1
273 ELSE
274 iminu=max(2,istru-1)
275 imaxu=min(iend+1,lm(ng))
276 iminv=max(1,istr-1)
277 imaxv=min(iend+1,lm(ng))
278 END IF
279 IF (nsperiodic(ng)) THEN
280 jminu=jstr-1
281 jmaxu=jend+1
282 jminv=jstr-1
283 jmaxv=jend+1
284 ELSE
285 jminu=max(1,jstr-1)
286 jmaxu=min(jend+1,mm(ng))
287 jminv=max(2,jstrv-1)
288 jmaxv=min(jend+1,mm(ng))
289 END IF
290!
291! Compute flux-components of the horizontal divergence of the stress
292! tensor (m4 s^-3/2) in XI- and ETA-directions. It is assumed here
293! that mixing coefficients are the squared root of the biharmonic
294! viscosity coefficient. For momentum balance purposes, the
295! thickness "Hz" appears only when computing the second harmonic
296! operator.
297!
298 k_loop : DO k=1,n(ng)
299 DO j=jminv-1,jmaxv
300 DO i=iminu-1,imaxu
301 cff=0.5_r8* &
302 & (pmon_r(i,j)* &
303 & ((pn(i ,j)+pn(i+1,j))*u(i+1,j,k,nrhs)- &
304 & (pn(i-1,j)+pn(i ,j))*u(i ,j,k,nrhs))- &
305 & pnom_r(i,j)* &
306 & ((pm(i,j )+pm(i,j+1))*v(i,j+1,k,nrhs)- &
307 & (pm(i,j-1)+pm(i,j ))*v(i,j ,k,nrhs)))
308#ifdef VISC_3DCOEF
309# ifdef UV_U3ADV_SPLIT
310 ufx(i,j)=on_r(i,j)*on_r(i,j)*uvis3d_r(i,j,k)*cff
311 vfe(i,j)=om_r(i,j)*om_r(i,j)*vvis3d_r(i,j,k)*cff
312# else
313 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc3d_r(i,j,k)*cff
314 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc3d_r(i,j,k)*cff
315# endif
316#else
317 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc4_r(i,j)*cff
318 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc4_r(i,j)*cff
319#endif
320 END DO
321 END DO
322 DO j=jminu,jmaxu+1
323 DO i=iminv,imaxv+1
324 cff=0.5_r8* &
325 & (pmon_p(i,j)* &
326 & ((pn(i ,j-1)+pn(i ,j))*v(i ,j,k,nrhs)- &
327 & (pn(i-1,j-1)+pn(i-1,j))*v(i-1,j,k,nrhs))+ &
328 & pnom_p(i,j)* &
329 & ((pm(i-1,j )+pm(i,j ))*u(i,j ,k,nrhs)- &
330 & (pm(i-1,j-1)+pm(i,j-1))*u(i,j-1,k,nrhs)))
331#ifdef MASKING
332 cff=cff*pmask(i,j)
333#endif
334#ifdef WET_DRY
335 cff=cff*pmask_wet(i,j)
336#endif
337#ifdef VISC_3DCOEF
338# ifdef UV_U3ADV_SPLIT
339 uvis_p=0.25_r8*(uvis3d_r(i-1,j-1,k)+uvis3d_r(i-1,j,k)+ &
340 & uvis3d_r(i ,j-1,k)+uvis3d_r(i ,j,k))
341 vvis_p=0.25_r8*(vvis3d_r(i-1,j-1,k)+vvis3d_r(i-1,j,k)+ &
342 & vvis3d_r(i ,j-1,k)+vvis3d_r(i ,j,k))
343 ufe(i,j)=om_p(i,j)*om_p(i,j)*uvis_p*cff
344 vfx(i,j)=on_p(i,j)*on_p(i,j)*vvis_p*cff
345# else
346 visc_p=0.25_r8*(visc3d_r(i-1,j-1,k)+visc3d_r(i-1,j,k)+ &
347 & visc3d_r(i ,j-1,k)+visc3d_r(i ,j,k))
348 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc_p*cff
349 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc_p*cff
350# endif
351#else
352 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc4_p(i,j)*cff
353 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc4_p(i,j)*cff
354#endif
355 END DO
356 END DO
357!
358! Compute first harmonic operator (m s^-3/2).
359!
360 DO j=jminu,jmaxu
361 DO i=iminu,imaxu
362 lapu(i,j)=0.125_r8* &
363 & (pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))* &
364 & ((pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))+ &
365 & (pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j)))
366 END DO
367 END DO
368 DO j=jminv,jmaxv
369 DO i=iminv,imaxv
370 lapv(i,j)=0.125_r8* &
371 & (pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))* &
372 & ((pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))- &
373 & (pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1)))
374 END DO
375 END DO
376!
377! Apply boundary conditions (other than periodic) to the first
378! harmonic operator. These are gradient or closed (free slip or
379! no slip) boundary conditions.
380!
381 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
382 IF (domain(ng)%Western_Edge(tile)) THEN
383 IF (lbc(iwest,isuvel,ng)%closed) THEN
384 DO j=jminu,jmaxu
385 lapu(istr,j)=0.0_r8
386 END DO
387 ELSE
388 DO j=jminu,jmaxu
389 lapu(istr,j)=lapu(istr+1,j)
390 END DO
391 END IF
392 IF (lbc(iwest,isvvel,ng)%closed) THEN
393 DO j=jminv,jmaxv
394 lapv(istr-1,j)=gamma2(ng)*lapv(istr,j)
395 END DO
396 ELSE
397 DO j=jminv,jmaxv
398 lapv(istr-1,j)=0.0_r8
399 END DO
400 END IF
401 END IF
402 END IF
403!
404 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
405 IF (domain(ng)%Eastern_Edge(tile)) THEN
406 IF (lbc(ieast,isuvel,ng)%closed) THEN
407 DO j=jminu,jmaxu
408 lapu(iend+1,j)=0.0_r8
409 END DO
410 ELSE
411 DO j=jminu,jmaxu
412 lapu(iend+1,j)=lapu(iend,j)
413 END DO
414 END IF
415 IF (lbc(ieast,isvvel,ng)%closed) THEN
416 DO j=jminv,jmaxv
417 lapv(iend+1,j)=gamma2(ng)*lapv(iend,j)
418 END DO
419 ELSE
420 DO j=jminv,jmaxv
421 lapv(iend+1,j)=0.0_r8
422 END DO
423 END IF
424 END IF
425 END IF
426!
427 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
428 IF (domain(ng)%Southern_Edge(tile)) THEN
429 IF (lbc(isouth,isuvel,ng)%closed) THEN
430 DO i=iminu,imaxu
431 lapu(i,jstr-1)=gamma2(ng)*lapu(i,jstr)
432 END DO
433 ELSE
434 DO i=iminu,imaxu
435 lapu(i,jstr-1)=0.0_r8
436 END DO
437 END IF
438 IF (lbc(isouth,isvvel,ng)%closed) THEN
439 DO i=iminv,imaxv
440 lapv(i,jstr)=0.0_r8
441 END DO
442 ELSE
443 DO i=iminv,imaxv
444 lapv(i,jstr)=lapv(i,jstr+1)
445 END DO
446 END IF
447 END IF
448 END IF
449!
450 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
451 IF (domain(ng)%Northern_Edge(tile)) THEN
452 IF (lbc(inorth,isuvel,ng)%closed) THEN
453 DO i=iminu,imaxu
454 lapu(i,jend+1)=gamma2(ng)*lapu(i,jend)
455 END DO
456 ELSE
457 DO i=iminu,imaxu
458 lapu(i,jend+1)=0.0_r8
459 END DO
460 END IF
461 IF (lbc(inorth,isvvel,ng)%closed) THEN
462 DO i=iminv,imaxv
463 lapv(i,jend+1)=0.0_r8
464 END DO
465 ELSE
466 DO i=iminv,imaxv
467 lapv(i,jend+1)=lapv(i,jend)
468 END DO
469 END IF
470 END IF
471 END IF
472!
473 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
474 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
475 IF (domain(ng)%SouthWest_Corner(tile)) THEN
476 lapu(istr ,jstr-1)=0.5_r8* &
477 & (lapu(istr+1,jstr-1)+ &
478 & lapu(istr ,jstr ))
479 lapv(istr-1,jstr )=0.5_r8* &
480 & (lapv(istr-1,jstr+1)+ &
481 & lapv(istr ,jstr ))
482 END IF
483 END IF
484
485 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
486 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
487 IF (domain(ng)%SouthEast_Corner(tile)) THEN
488 lapu(iend+1,jstr-1)=0.5_r8* &
489 & (lapu(iend ,jstr-1)+ &
490 & lapu(iend+1,jstr ))
491 lapv(iend+1,jstr )=0.5_r8* &
492 & (lapv(iend ,jstr )+ &
493 & lapv(iend+1,jstr+1))
494 END IF
495 END IF
496
497 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
498 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
499 IF (domain(ng)%NorthWest_Corner(tile)) THEN
500 lapu(istr ,jend+1)=0.5_r8* &
501 & (lapu(istr+1,jend+1)+ &
502 & lapu(istr ,jend ))
503 lapv(istr-1,jend+1)=0.5_r8* &
504 & (lapv(istr ,jend+1)+ &
505 & lapv(istr-1,jend ))
506 END IF
507 END IF
508
509 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
510 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
511 IF (domain(ng)%NorthEast_Corner(tile)) THEN
512 lapu(iend+1,jend+1)=0.5_r8* &
513 & (lapu(iend ,jend+1)+ &
514 & lapu(iend+1,jend ))
515 lapv(iend+1,jend+1)=0.5_r8* &
516 & (lapv(iend ,jend+1)+ &
517 & lapv(iend+1,jend ))
518 END IF
519 END IF
520!
521! Compute flux-components of the horizontal divergence of the
522! harmonic stress tensor (m4/s2) in XI- and ETA-directions.
523!
524 DO j=jstrv-1,jend
525 DO i=istru-1,iend
526 cff=hz(i,j,k)*0.5_r8* &
527 & (pmon_r(i,j)* &
528 & ((pn(i ,j)+pn(i+1,j))*lapu(i+1,j)- &
529 & (pn(i-1,j)+pn(i ,j))*lapu(i ,j))- &
530 & pnom_r(i,j)* &
531 & ((pm(i,j )+pm(i,j+1))*lapv(i,j+1)- &
532 & (pm(i,j-1)+pm(i,j ))*lapv(i,j )))
533#ifdef VISC_3DCOEF
534# ifdef UV_U3ADV_SPLIT
535 ufx(i,j)=on_r(i,j)*on_r(i,j)*uvis3d_r(i,j,k)*cff
536 vfe(i,j)=om_r(i,j)*om_r(i,j)*vvis3d_r(i,j,k)*cff
537# else
538 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc3d_r(i,j,k)*cff
539 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc3d_r(i,j,k)*cff
540# endif
541#else
542 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc4_r(i,j)*cff
543 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc4_r(i,j)*cff
544#endif
545 END DO
546 END DO
547 DO j=jstr,jend+1
548 DO i=istr,iend+1
549 cff=0.125_r8*(hz(i-1,j ,k)+hz(i,j ,k)+ &
550 & hz(i-1,j-1,k)+hz(i,j-1,k))* &
551 & (pmon_p(i,j)* &
552 & ((pn(i ,j-1)+pn(i ,j))*lapv(i ,j)- &
553 & (pn(i-1,j-1)+pn(i-1,j))*lapv(i-1,j))+ &
554 & pnom_p(i,j)* &
555 & ((pm(i-1,j )+pm(i,j ))*lapu(i,j )- &
556 & (pm(i-1,j-1)+pm(i,j-1))*lapu(i,j-1)))
557#ifdef MASKING
558 cff=cff*pmask(i,j)
559#endif
560#ifdef WET_DRY
561 cff=cff*pmask_wet(i,j)
562#endif
563#ifdef VISC_3DCOEF
564# ifdef UV_U3ADV_SPLIT
565 uvis_p=0.25_r8*(uvis3d_r(i-1,j-1,k)+uvis3d_r(i-1,j,k)+ &
566 & uvis3d_r(i ,j-1,k)+uvis3d_r(i ,j,k))
567 vvis_p=0.25_r8*(vvis3d_r(i-1,j-1,k)+vvis3d_r(i-1,j,k)+ &
568 & vvis3d_r(i ,j-1,k)+vvis3d_r(i ,j,k))
569 ufe(i,j)=om_p(i,j)*om_p(i,j)*uvis_p*cff
570 vfx(i,j)=on_p(i,j)*on_p(i,j)*vvis_p*cff
571# else
572 visc_p=0.25_r8*(visc3d_r(i-1,j-1,k)+visc3d_r(i-1,j,k)+ &
573 & visc3d_r(i ,j-1,k)+visc3d_r(i ,j,k))
574 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc_p*cff
575 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc_p*cff
576# endif
577#else
578 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc4_p(i,j)*cff
579 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc4_p(i,j)*cff
580#endif
581 END DO
582 END DO
583!
584! Time-step biharmonic, S-surfaces viscosity term. Notice that
585! momentum at this stage is HzU and HzV and has units m2/s. Add
586! contribution for barotropic forcing terms.
587!
588 DO j=jstr,jend
589 DO i=istru,iend
590 cff=dt(ng)*0.25_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
591 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
592 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
593 cff3=cff*(cff1+cff2)
594 rufrc(i,j)=rufrc(i,j)-cff1-cff2
595 u(i,j,k,nnew)=u(i,j,k,nnew)-cff3
596#ifdef DIAGNOSTICS_UV
597 diarufrc(i,j,3,m2hvis)=diarufrc(i,j,3,m2hvis)-cff1-cff2
598 diarufrc(i,j,3,m2xvis)=diarufrc(i,j,3,m2xvis)-cff1
599 diarufrc(i,j,3,m2yvis)=diarufrc(i,j,3,m2yvis)-cff2
600 diau3wrk(i,j,k,m3hvis)=-cff3
601 diau3wrk(i,j,k,m3xvis)=-cff*cff1
602 diau3wrk(i,j,k,m3yvis)=-cff*cff2
603#endif
604 END DO
605 END DO
606 DO j=jstrv,jend
607 DO i=istr,iend
608 cff=dt(ng)*0.25_r8*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
609 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
610 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
611 cff3=cff*(cff1-cff2)
612 rvfrc(i,j)=rvfrc(i,j)-cff1+cff2
613 v(i,j,k,nnew)=v(i,j,k,nnew)-cff3
614#ifdef DIAGNOSTICS_UV
615 diarvfrc(i,j,3,m2hvis)=diarvfrc(i,j,3,m2hvis)-cff1+cff2
616 diarvfrc(i,j,3,m2xvis)=diarvfrc(i,j,3,m2xvis)-cff1
617 diarvfrc(i,j,3,m2yvis)=diarvfrc(i,j,3,m2yvis)+cff2
618 diav3wrk(i,j,k,m3hvis)=-cff3
619 diav3wrk(i,j,k,m3xvis)=-cff*cff1
620 diav3wrk(i,j,k,m3yvis)= cff*cff2
621#endif
622 END DO
623 END DO
624 END DO k_loop
625!
626 RETURN
627 END SUBROUTINE uv3dmix4_s_tile
628
629 END MODULE uv3dmix4_mod
mod_coupling
Definition mod_coupling.F:2
mod_coupling::coupling
type(t_coupling), dimension(:), allocatable coupling
Definition mod_coupling.F:102
mod_diags
Definition mod_diags.F:2
mod_diags::diags
type(t_diags), dimension(:), allocatable diags
Definition mod_diags.F:100
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ncparam
Definition mod_ncparam.F:2
mod_ncparam::isvvel
integer isvvel
Definition mod_ncparam.F:507
mod_ncparam::isuvel
integer isuvel
Definition mod_ncparam.F:506
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
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_param::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
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::m3xvis
integer m3xvis
Definition mod_scalars.F:193
mod_scalars::gamma2
real(r8), dimension(:), allocatable gamma2
Definition mod_scalars.F:743
mod_scalars::m2hvis
integer m2hvis
Definition mod_scalars.F:168
mod_scalars::m3hvis
integer m3hvis
Definition mod_scalars.F:192
mod_scalars::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::m3yvis
integer m3yvis
Definition mod_scalars.F:194
mod_scalars::m2yvis
integer m2yvis
Definition mod_scalars.F:170
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::m2xvis
integer m2xvis
Definition mod_scalars.F:169
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
uv3dmix4_mod
Definition uv3dmix4_geo.h:1
uv3dmix4_mod::uv3dmix4
subroutine, public uv3dmix4(ng, tile)
Definition uv3dmix4_geo.h:43
uv3dmix4_mod::uv3dmix4_s_tile
subroutine uv3dmix4_s_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nnew, pmask, pmask_wet, hz, om_p, om_r, on_p, on_r, pm, pmon_p, pmon_r, pn, pnom_p, pnom_r, uvis3d_r, vvis3d_r, visc3d_r, visc4_p, visc4_r, diarufrc, diarvfrc, diau3wrk, diav3wrk, rufrc, rvfrc, u, v)
Definition uv3dmix4_s.h:147
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3