ROMS
Loading...
Searching...
No Matches
uv3dmix2_geo.h
Go to the documentation of this file.
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 routine computes harmonic mixing of momentum, rotated along !
11! geopotentials, from the horizontal divergence of the stress !
12! tensor. A transverse isotropy is assumed so the stress tensor !
13! 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 uv3dmix2
38!
39 CONTAINS
40!
41!***********************************************************************
42 SUBROUTINE uv3dmix2 (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, 31, __line__, myfile)
68#endif
69 CALL uv3dmix2_geo_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 & grid(ng) % rmask, &
76 & grid(ng) % umask, &
77 & grid(ng) % vmask, &
78#endif
79#ifdef WET_DRY
80 & grid(ng) % pmask_wet, &
81 & grid(ng) % rmask_wet, &
82 & grid(ng) % umask_wet, &
83 & grid(ng) % vmask_wet, &
84#endif
85 & grid(ng) % om_p, &
86 & grid(ng) % om_r, &
87 & grid(ng) % om_u, &
88 & grid(ng) % om_v, &
89 & grid(ng) % on_p, &
90 & grid(ng) % on_r, &
91 & grid(ng) % on_u, &
92 & grid(ng) % on_v, &
93 & grid(ng) % pm, &
94 & grid(ng) % pn, &
95 & grid(ng) % Hz, &
96 & grid(ng) % z_r, &
97#ifdef VISC_3DCOEF
98 & mixing(ng) % visc3d_r, &
99#else
100 & mixing(ng) % visc2_p, &
101 & mixing(ng) % visc2_r, &
102#endif
103#ifdef DIAGNOSTICS_UV
104 & diags(ng) % DiaRUfrc, &
105 & diags(ng) % DiaRVfrc, &
106 & diags(ng) % DiaU3wrk, &
107 & diags(ng) % DiaV3wrk, &
108#endif
109 & ocean(ng) % u, &
110 & ocean(ng) % v, &
111 & coupling(ng) % rufrc, &
112 & coupling(ng) % rvfrc)
113#ifdef PROFILE
114 CALL wclock_off (ng, inlm, 31, __line__, myfile)
115#endif
116!
117 RETURN
118 END SUBROUTINE uv3dmix2
119!
120!***********************************************************************
121 SUBROUTINE uv3dmix2_geo_tile (ng, tile, &
122 & LBi, UBi, LBj, UBj, &
123 & IminS, ImaxS, JminS, JmaxS, &
124 & nrhs, nnew, &
125#ifdef MASKING
126 & pmask, rmask, umask, vmask, &
127#endif
128#ifdef WET_DRY
129 & pmask_wet, rmask_wet, &
130 & umask_wet, vmask_wet, &
131#endif
132 & om_p, om_r, om_u, om_v, &
133 & on_p, on_r, on_u, on_v, &
134 & pm, pn, &
135 & Hz, z_r, &
136#ifdef VISC_3DCOEF
137 & visc3d_r, &
138#else
139 & visc2_p, visc2_r, &
140#endif
141#ifdef DIAGNOSTICS_UV
142 & DiaRUfrc, DiaRVfrc, &
143 & DiaU3wrk, DiaV3wrk, &
144#endif
145 & u, v, &
146 & rufrc, rvfrc)
147!***********************************************************************
148!
149 USE mod_param
150 USE mod_scalars
151!
152! Imported variable declarations.
153!
154 integer, intent(in) :: ng, tile
155 integer, intent(in) :: LBi, UBi, LBj, UBj
156 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
157 integer, intent(in) :: nrhs, nnew
158
159#ifdef ASSUMED_SHAPE
160# ifdef MASKING
161 real(r8), intent(in) :: pmask(LBi:,LBj:)
162 real(r8), intent(in) :: rmask(LBi:,LBj:)
163 real(r8), intent(in) :: umask(LBi:,LBj:)
164 real(r8), intent(in) :: vmask(LBi:,LBj:)
165# endif
166# ifdef WET_DRY
167 real(r8), intent(in) :: pmask_wet(LBi:,LBj:)
168 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
169 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
170 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
171# endif
172 real(r8), intent(in) :: om_p(LBi:,LBj:)
173 real(r8), intent(in) :: om_r(LBi:,LBj:)
174 real(r8), intent(in) :: om_u(LBi:,LBj:)
175 real(r8), intent(in) :: om_v(LBi:,LBj:)
176 real(r8), intent(in) :: on_p(LBi:,LBj:)
177 real(r8), intent(in) :: on_r(LBi:,LBj:)
178 real(r8), intent(in) :: on_u(LBi:,LBj:)
179 real(r8), intent(in) :: on_v(LBi:,LBj:)
180 real(r8), intent(in) :: pm(LBi:,LBj:)
181 real(r8), intent(in) :: pn(LBi:,LBj:)
182 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
183 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
184# ifdef VISC_3DCOEF
185 real(r8), intent(in) :: visc3d_r(LBi:,LBj:,:)
186# else
187 real(r8), intent(in) :: visc2_p(LBi:,LBj:)
188 real(r8), intent(in) :: visc2_r(LBi:,LBj:)
189# endif
190# ifdef DIAGNOSTICS_UV
191 real(r8), intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
192 real(r8), intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
193 real(r8), intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
194 real(r8), intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
195# endif
196 real(r8), intent(inout) :: rufrc(LBi:,LBj:)
197 real(r8), intent(inout) :: rvfrc(LBi:,LBj:)
198 real(r8), intent(inout) :: u(LBi:,LBj:,:,:)
199 real(r8), intent(inout) :: v(LBi:,LBj:,:,:)
200#else
201# ifdef MASKING
202 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
203 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
204 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
205 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
206# endif
207# ifdef WET_DRY
208 real(r8), intent(in) :: pmask_wet(LBi:UBi,LBj:UBj)
209 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
210 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
211 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
212# endif
213 real(r8), intent(in) :: om_p(LBi:UBi,LBj:UBj)
214 real(r8), intent(in) :: om_r(LBi:UBi,LBj:UBj)
215 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
216 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
217 real(r8), intent(in) :: on_p(LBi:UBi,LBj:UBj)
218 real(r8), intent(in) :: on_r(LBi:UBi,LBj:UBj)
219 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
220 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
221 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
222 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
223 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
224 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
225# ifdef VISC_3DCOEF
226 real(r8), intent(in) :: visc3d_r(LBi:UBi,LBj:UBj,N(ng))
227# else
228 real(r8), intent(in) :: visc2_p(LBi:UBi,LBj:UBj)
229 real(r8), intent(in) :: visc2_r(LBi:UBi,LBj:UBj)
230# endif
231# ifdef DIAGNOSTICS_UV
232 real(r8), intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
233 real(r8), intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
234 real(r8), intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
235 real(r8), intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
236# endif
237 real(r8), intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
238 real(r8), intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
239 real(r8), intent(inout) :: u(LBi:UBi,LBj:UBj,N(ng),2)
240 real(r8), intent(inout) :: v(LBi:UBi,LBj:UBj,N(ng),2)
241#endif
242!
243! Local variable declarations.
244!
245 integer :: i, j, k, k1, k2
246
247 real(r8) :: cff, fac1, fac2, pm_p, pn_p
248 real(r8) :: cff1, cff2, cff3, cff4
249 real(r8) :: cff5, cff6, cff7, cff8
250 real(r8) :: dmUdz, dnUdz, dmVdz, dnVdz
251#ifdef VISC_3DCOEF
252 real(r8) :: visc_p
253#endif
254 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
255 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
256 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
257 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
258
259 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: UFse
260 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: UFsx
261 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: VFse
262 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: VFsx
263 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dmUde
264 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dmVde
265 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dnUdx
266 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dnVdx
267 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dUdz
268 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dVdz
269 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_p
270 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_r
271 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_p
272 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_r
273
274#include "set_bounds.h"
275!
276!-----------------------------------------------------------------------
277! Compute horizontal harmonic viscosity along geopotential surfaces.
278!-----------------------------------------------------------------------
279!
280! Compute horizontal and vertical gradients. Notice the recursive
281! blocking sequence. The vertical placement of the gradients is:
282!
283! dZdx_r, dZde_r, dnUdx, dmVde(:,:,k1) k rho-points
284! dZdx_r, dZde_r, dnUdx, dmVde(:,:,k2) k+1 rho-points
285! dZdx_p, dZde_p, dnVdx, dmUde(:,:,k1) k psi-points
286! dZdx_p, dZde_p, dnVdx, dmUde(:,:,k2) k+1 psi-points
287! UFse, UFsx, dUdz(:,:,k1) k-1/2 WU-points
288! UFse, UFsx, dUdz(:,:,k2) k+1/2 WU-points
289! VFse, VFsx, dVdz(:,:,k1) k-1/2 WV-points
290! VFse, VFsx, dVdz(:,:,k2) k+1/2 WV-points
291!
292 k2=1
293 k_loop : DO k=0,n(ng)
294 k1=k2
295 k2=3-k1
296 IF (k.lt.n(ng)) THEN
297!
298! Compute slopes (nondimensional) at RHO- and PSI-points.
299!
300 DO j=jstr-1,jend+1
301 DO i=istru-1,iend+1
302 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
303#ifdef MASKING
304 cff=cff*umask(i,j)
305#endif
306#ifdef WET_DRY
307 cff=cff*umask_wet(i,j)
308#endif
309 ufx(i,j)=cff*(z_r(i ,j,k+1)- &
310 & z_r(i-1,j,k+1))
311 END DO
312 END DO
313 DO j=jstrv-1,jend+1
314 DO i=istr-1,iend+1
315 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
316#ifdef MASKING
317 cff=cff*vmask(i,j)
318#endif
319#ifdef WET_DRY
320 cff=cff*vmask_wet(i,j)
321#endif
322 vfe(i,j)=cff*(z_r(i,j ,k+1)- &
323 & z_r(i,j-1,k+1))
324 END DO
325 END DO
326!
327 DO j=jstr,jend+1
328 DO i=istr,iend+1
329 dzdx_p(i,j,k2)=0.5_r8*(ufx(i,j-1)+ &
330 & ufx(i,j ))
331 dzde_p(i,j,k2)=0.5_r8*(vfe(i-1,j)+ &
332 & vfe(i ,j))
333 END DO
334 END DO
335 DO j=jstrv-1,jend
336 DO i=istru-1,iend
337 dzdx_r(i,j,k2)=0.5_r8*(ufx(i ,j)+ &
338 & ufx(i+1,j))
339 dzde_r(i,j,k2)=0.5_r8*(vfe(i,j )+ &
340 & vfe(i,j+1))
341 END DO
342 END DO
343!
344! Compute momentum horizontal (1/m/s) and vertical (1/s) gradients.
345!
346 DO j=jstrv-1,jend
347 DO i=istru-1,iend
348 cff=0.5_r8*pm(i,j)
349#ifdef MASKING
350 cff=cff*rmask(i,j)
351#endif
352#ifdef WET_DRY
353 cff=cff*rmask_wet(i,j)
354#endif
355 dnudx(i,j,k2)=cff*((pn(i ,j)+pn(i+1,j))* &
356 & u(i+1,j,k+1,nrhs)- &
357 & (pn(i-1,j)+pn(i ,j))* &
358 & u(i ,j,k+1,nrhs))
359 END DO
360 END DO
361
362 DO j=jstr,jend+1
363 DO i=istr,iend+1
364 cff=0.125_r8*(pn(i-1,j )+pn(i,j )+ &
365 & pn(i-1,j-1)+pn(i,j-1))
366#ifdef MASKING
367 cff=cff*pmask(i,j)
368#endif
369#ifdef WET_DRY
370 cff=cff*pmask_wet(i,j)
371#endif
372 dmude(i,j,k2)=cff*((pm(i-1,j )+pm(i,j ))* &
373 & u(i,j ,k+1,nrhs)- &
374 & (pm(i-1,j-1)+pm(i,j-1))* &
375 & u(i,j-1,k+1,nrhs))
376 END DO
377 END DO
378
379 DO j=jstr,jend+1
380 DO i=istr,iend+1
381 cff=0.125_r8*(pm(i-1,j )+pm(i,j )+ &
382 & pm(i-1,j-1)+pm(i,j-1))
383#ifdef MASKING
384 cff=cff*pmask(i,j)
385#endif
386#ifdef WET_DRY
387 cff=cff*pmask_wet(i,j)
388#endif
389 dnvdx(i,j,k2)=cff*((pn(i ,j-1)+pn(i ,j))* &
390 & v(i ,j,k+1,nrhs)- &
391 & (pn(i-1,j-1)+pn(i-1,j))* &
392 & v(i-1,j,k+1,nrhs))
393 END DO
394 END DO
395
396 DO j=jstrv-1,jend
397 DO i=istru-1,iend
398 cff=0.5_r8*pn(i,j)
399#ifdef MASKING
400 cff=cff*rmask(i,j)
401#endif
402#ifdef WET_DRY
403 cff=cff*rmask_wet(i,j)
404#endif
405 dmvde(i,j,k2)=cff*((pm(i,j )+pm(i,j+1))* &
406 & v(i,j+1,k+1,nrhs)- &
407 & (pm(i,j-1)+pm(i,j ))* &
408 & v(i,j ,k+1,nrhs))
409 END DO
410 END DO
411 END IF
412
413 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
414 DO j=jstr-1,jend+1
415 DO i=istru-1,iend+1
416 dudz(i,j,k2)=0.0_r8
417 END DO
418 END DO
419 DO j=jstrv-1,jend+1
420 DO i=istr-1,iend+1
421 dvdz(i,j,k2)=0.0_r8
422 END DO
423 END DO
424
425 DO j=jstr,jend
426 DO i=istru,iend
427 ufsx(i,j,k2)=0.0_r8
428 ufse(i,j,k2)=0.0_r8
429 END DO
430 END DO
431 DO j=jstrv,jend
432 DO i=istr,iend
433 vfsx(i,j,k2)=0.0_r8
434 vfse(i,j,k2)=0.0_r8
435 END DO
436 END DO
437 ELSE
438 DO j=jstr-1,jend+1
439 DO i=istru-1,iend+1
440 cff=1.0_r8/(0.5_r8*(z_r(i-1,j,k+1)-z_r(i-1,j,k)+ &
441 & z_r(i ,j,k+1)-z_r(i ,j,k)))
442 dudz(i,j,k2)=cff*(u(i,j,k+1,nrhs)- &
443 & u(i,j,k ,nrhs))
444 END DO
445 END DO
446
447 DO j=jstrv-1,jend+1
448 DO i=istr-1,iend+1
449 cff=1.0_r8/(0.5_r8*(z_r(i,j-1,k+1)-z_r(i,j-1,k)+ &
450 & z_r(i,j ,k+1)-z_r(i,j ,k)))
451 dvdz(i,j,k2)=cff*(v(i,j,k+1,nrhs)- &
452 & v(i,j,k ,nrhs))
453 END DO
454 END DO
455 END IF
456!
457! Compute components of the rotated viscous flux (m5/s2) along
458! geopotential surfaces in the XI- and ETA-directions.
459!
460 IF (k.gt.0) THEN
461 DO j=jstrv-1,jend
462 DO i=istru-1,iend
463 cff1=min(dzdx_r(i,j,k1),0.0_r8)
464 cff2=max(dzdx_r(i,j,k1),0.0_r8)
465 cff3=min(dzde_r(i,j,k1),0.0_r8)
466 cff4=max(dzde_r(i,j,k1),0.0_r8)
467 cff=hz(i,j,k)* &
468 & (on_r(i,j)*(dnudx(i,j,k1)- &
469 & 0.5_r8*pn(i,j)* &
470 & (cff1*(dudz(i ,j,k1)+ &
471 & dudz(i+1,j,k2))+ &
472 & cff2*(dudz(i ,j,k2)+ &
473 & dudz(i+1,j,k1))))- &
474 & om_r(i,j)*(dmvde(i,j,k1)- &
475 & 0.5_r8*pm(i,j)* &
476 & (cff3*(dvdz(i,j ,k1)+ &
477 & dvdz(i,j+1,k2))+ &
478 & cff4*(dvdz(i,j ,k2)+ &
479 & dvdz(i,j+1,k1)))))
480#ifdef MASKING
481 cff=cff*rmask(i,j)
482#endif
483#ifdef WET_DRY
484 cff=cff*rmask_wet(i,j)
485#endif
486#ifdef VISC_3DCOEF
487 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc3d_r(i,j,k)*cff
488 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc3d_r(i,j,k)*cff
489#else
490 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc2_r(i,j)*cff
491 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc2_r(i,j)*cff
492#endif
493 END DO
494 END DO
495
496 DO j=jstr,jend+1
497 DO i=istr,iend+1
498 pm_p=0.25_r8*(pm(i-1,j-1)+pm(i-1,j)+ &
499 & pm(i ,j-1)+pm(i ,j))
500 pn_p=0.25_r8*(pn(i-1,j-1)+pn(i-1,j)+ &
501 & pn(i ,j-1)+pn(i ,j))
502 cff1=min(dzdx_p(i,j,k1),0.0_r8)
503 cff2=max(dzdx_p(i,j,k1),0.0_r8)
504 cff3=min(dzde_p(i,j,k1),0.0_r8)
505 cff4=max(dzde_p(i,j,k1),0.0_r8)
506 cff=0.25_r8* &
507 & (hz(i-1,j ,k)+hz(i,j ,k)+ &
508 & hz(i-1,j-1,k)+hz(i,j-1,k))* &
509 & (on_p(i,j)*(dnvdx(i,j,k1)- &
510 & 0.5_r8*pn_p* &
511 & (cff1*(dvdz(i-1,j,k1)+ &
512 & dvdz(i ,j,k2))+ &
513 & cff2*(dvdz(i-1,j,k2)+ &
514 & dvdz(i ,j,k1))))+ &
515 & om_p(i,j)*(dmude(i,j,k1)- &
516 & 0.5_r8*pm_p* &
517 & (cff3*(dudz(i,j-1,k1)+ &
518 & dudz(i,j ,k2))+ &
519 & cff4*(dudz(i,j-1,k2)+ &
520 & dudz(i,j ,k1)))))
521#ifdef MASKING
522 cff=cff*pmask(i,j)
523#endif
524#ifdef WET_DRY
525 cff=cff*pmask_wet(i,j)
526#endif
527#ifdef VISC_3DCOEF
528 visc_p=0.25_r8* &
529 & (visc3d_r(i-1,j-1,k)+visc3d_r(i-1,j,k)+ &
530 & visc3d_r(i ,j-1,k)+visc3d_r(i ,j,k))
531 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc_p*cff
532 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc_p*cff
533#else
534 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc2_p(i,j)*cff
535 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc2_p(i,j)*cff
536#endif
537 END DO
538 END DO
539!
540! Compute vertical flux (m2/s2) due to sloping terrain-following
541! surfaces.
542!
543 IF (k.lt.n(ng)) THEN
544 DO j=jstr,jend
545 DO i=istru,iend
546#ifdef VISC_3DCOEF
547 cff=0.125_r8* &
548 & (visc3d_r(i-1,j,k )+visc3d_r(i,j,k )+ &
549 & visc3d_r(i-1,j,k+1)+visc3d_r(i,j,k+1))
550 fac1=cff*on_u(i,j)
551 fac2=cff*om_u(i,j)
552#else
553 cff=0.25_r8*(visc2_r(i-1,j)+visc2_r(i,j))
554 fac1=cff*on_u(i,j)
555 fac2=cff*om_u(i,j)
556#endif
557 cff=0.5_r8*(pn(i-1,j)+pn(i,j))
558 dnudz=cff*dudz(i,j,k2)
559 dnvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+ &
560 & dvdz(i ,j+1,k2)+ &
561 & dvdz(i-1,j ,k2)+ &
562 & dvdz(i ,j ,k2))
563 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
564 dmudz=cff*dudz(i,j,k2)
565 dmvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+ &
566 & dvdz(i ,j+1,k2)+ &
567 & dvdz(i-1,j ,k2)+ &
568 & dvdz(i ,j ,k2))
569
570 cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
571 cff2=min(dzdx_r(i ,j,k2),0.0_r8)
572 cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
573 cff4=max(dzdx_r(i ,j,k1),0.0_r8)
574 ufsx(i,j,k2)=fac1* &
575 & (cff1*(cff1*dnudz-dnudx(i-1,j,k1))+ &
576 & cff2*(cff2*dnudz-dnudx(i ,j,k2))+ &
577 & cff3*(cff3*dnudz-dnudx(i-1,j,k2))+ &
578 & cff4*(cff4*dnudz-dnudx(i ,j,k1)))
579
580 cff1=min(dzde_p(i,j ,k1),0.0_r8)
581 cff2=min(dzde_p(i,j+1,k2),0.0_r8)
582 cff3=max(dzde_p(i,j ,k2),0.0_r8)
583 cff4=max(dzde_p(i,j+1,k1),0.0_r8)
584 ufse(i,j,k2)=fac2* &
585 & (cff1*(cff1*dmudz-dmude(i,j ,k1))+ &
586 & cff2*(cff2*dmudz-dmude(i,j+1,k2))+ &
587 & cff3*(cff3*dmudz-dmude(i,j ,k2))+ &
588 & cff4*(cff4*dmudz-dmude(i,j+1,k1)))
589
590 cff1=min(dzde_p(i,j ,k1),0.0_r8)
591 cff2=min(dzde_p(i,j+1,k2),0.0_r8)
592 cff3=max(dzde_p(i,j ,k2),0.0_r8)
593 cff4=max(dzde_p(i,j+1,k1),0.0_r8)
594 cff5=min(dzdx_p(i,j ,k1),0.0_r8)
595 cff6=min(dzdx_p(i,j+1,k2),0.0_r8)
596 cff7=max(dzdx_p(i,j ,k2),0.0_r8)
597 cff8=max(dzdx_p(i,j+1,k1),0.0_r8)
598 ufsx(i,j,k2)=ufsx(i,j,k2)+ &
599 & fac1* &
600 & (cff1*(cff5*dnvdz-dnvdx(i,j ,k1))+ &
601 & cff2*(cff6*dnvdz-dnvdx(i,j+1,k2))+ &
602 & cff3*(cff7*dnvdz-dnvdx(i,j ,k2))+ &
603 & cff4*(cff8*dnvdz-dnvdx(i,j+1,k1)))
604
605 cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
606 cff2=min(dzdx_r(i ,j,k2),0.0_r8)
607 cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
608 cff4=max(dzdx_r(i ,j,k1),0.0_r8)
609 cff5=min(dzde_r(i-1,j,k1),0.0_r8)
610 cff6=min(dzde_r(i ,j,k2),0.0_r8)
611 cff7=max(dzde_r(i-1,j,k2),0.0_r8)
612 cff8=max(dzde_r(i ,j,k1),0.0_r8)
613 ufse(i,j,k2)=ufse(i,j,k2)- &
614 & fac2* &
615 & (cff1*(cff5*dmvdz-dmvde(i-1,j,k1))+ &
616 & cff2*(cff6*dmvdz-dmvde(i ,j,k2))+ &
617 & cff3*(cff7*dmvdz-dmvde(i-1,j,k2))+ &
618 & cff4*(cff8*dmvdz-dmvde(i ,j,k1)))
619 END DO
620 END DO
621!
622 DO j=jstrv,jend
623 DO i=istr,iend
624#ifdef VISC_3DCOEF
625 cff=0.125_r8* &
626 & (visc3d_r(i,j-1,k )+visc3d_r(i,j,k )+ &
627 & visc3d_r(i,j-1,k+1)+visc3d_r(i,j,k+1))
628 fac1=cff*on_v(i,j)
629 fac2=cff*om_v(i,j)
630#else
631 cff=0.25_r8*(visc2_r(i,j-1)+visc2_r(i,j))
632 fac1=cff*on_v(i,j)
633 fac2=cff*om_v(i,j)
634#endif
635 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
636 dnudz=cff*0.25_r8*(dudz(i ,j ,k2)+ &
637 & dudz(i+1,j ,k2)+ &
638 & dudz(i ,j-1,k2)+ &
639 & dudz(i+1,j-1,k2))
640 dnvdz=cff*dvdz(i,j,k2)
641 cff=0.5_r8*(pm(i,j-1)+pm(i,j))
642 dmudz=cff*0.25_r8*(dudz(i ,j ,k2)+ &
643 & dudz(i+1,j ,k2)+ &
644 & dudz(i ,j-1,k2)+ &
645 & dudz(i+1,j-1,k2))
646 dmvdz=cff*dvdz(i,j,k2)
647
648 cff1=min(dzdx_p(i ,j,k1),0.0_r8)
649 cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
650 cff3=max(dzdx_p(i ,j,k2),0.0_r8)
651 cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
652 vfsx(i,j,k2)=fac1* &
653 & (cff1*(cff1*dnvdz-dnvdx(i ,j,k1))+ &
654 & cff2*(cff2*dnvdz-dnvdx(i+1,j,k2))+ &
655 & cff3*(cff3*dnvdz-dnvdx(i ,j,k2))+ &
656 & cff4*(cff4*dnvdz-dnvdx(i+1,j,k1)))
657
658 cff1=min(dzde_r(i,j-1,k1),0.0_r8)
659 cff2=min(dzde_r(i,j ,k2),0.0_r8)
660 cff3=max(dzde_r(i,j-1,k2),0.0_r8)
661 cff4=max(dzde_r(i,j ,k1),0.0_r8)
662 vfse(i,j,k2)=fac2* &
663 & (cff1*(cff1*dmvdz-dmvde(i,j-1,k1))+ &
664 & cff2*(cff2*dmvdz-dmvde(i,j ,k2))+ &
665 & cff3*(cff3*dmvdz-dmvde(i,j-1,k2))+ &
666 & cff4*(cff4*dmvdz-dmvde(i,j ,k1)))
667
668 cff1=min(dzde_r(i,j-1,k1),0.0_r8)
669 cff2=min(dzde_r(i,j ,k2),0.0_r8)
670 cff3=max(dzde_r(i,j-1,k2),0.0_r8)
671 cff4=max(dzde_r(i,j ,k1),0.0_r8)
672 cff5=min(dzdx_r(i,j-1,k1),0.0_r8)
673 cff6=min(dzdx_r(i,j ,k2),0.0_r8)
674 cff7=max(dzdx_r(i,j-1,k2),0.0_r8)
675 cff8=max(dzdx_r(i,j ,k1),0.0_r8)
676 vfsx(i,j,k2)=vfsx(i,j,k2)- &
677 & fac1* &
678 & (cff1*(cff5*dnudz-dnudx(i,j-1,k1))+ &
679 & cff2*(cff6*dnudz-dnudx(i,j ,k2))+ &
680 & cff3*(cff7*dnudz-dnudx(i,j-1,k2))+ &
681 & cff4*(cff8*dnudz-dnudx(i,j ,k1)))
682
683 cff1=min(dzdx_p(i ,j,k1),0.0_r8)
684 cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
685 cff3=max(dzdx_p(i ,j,k2),0.0_r8)
686 cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
687 cff5=min(dzde_p(i ,j,k1),0.0_r8)
688 cff6=min(dzde_p(i+1,j,k2),0.0_r8)
689 cff7=max(dzde_p(i ,j,k2),0.0_r8)
690 cff8=max(dzde_p(i+1,j,k1),0.0_r8)
691 vfse(i,j,k2)=vfse(i,j,k2)+ &
692 & fac2* &
693 & (cff1*(cff5*dmudz-dmude(i ,j,k1))+ &
694 & cff2*(cff6*dmudz-dmude(i+1,j,k2))+ &
695 & cff3*(cff7*dmudz-dmude(i ,j,k2))+ &
696 & cff4*(cff8*dmudz-dmude(i+1,j,k1)))
697 END DO
698 END DO
699 END IF
700!
701! Time-step harmonic, geopotential viscosity term. Notice that
702! momentum at this stage is HzU and HzV and has m2/s units. Add
703! contribution for barotropic forcing terms.
704#ifdef DIAGNOSTICS_UV
705! The rotated vertical term cannot be split from the horizontal
706! terms because of the 2D/3D momentum coupling.
707#endif
708!
709 DO j=jstr,jend
710 DO i=istru,iend
711 cff=dt(ng)*0.25_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
712 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
713 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
714 cff3=ufsx(i,j,k2)-ufsx(i,j,k1)
715 cff4=ufse(i,j,k2)-ufse(i,j,k1)
716 cff5=cff*(cff1+cff2)
717 cff6=dt(ng)*(cff3+cff4)
718 rufrc(i,j)=rufrc(i,j)+cff1+cff2+cff3+cff4
719 u(i,j,k,nnew)=u(i,j,k,nnew)+cff5+cff6
720#ifdef DIAGNOSTICS_UV
721 diarufrc(i,j,3,m2hvis)=diarufrc(i,j,3,m2hvis)+cff1+cff2+ &
722 & cff3+cff4
723 diarufrc(i,j,3,m2xvis)=diarufrc(i,j,3,m2xvis)+cff1+cff3
724 diarufrc(i,j,3,m2yvis)=diarufrc(i,j,3,m2yvis)+cff2+cff4
725 diau3wrk(i,j,k,m3hvis)=cff5+cff6
726 diau3wrk(i,j,k,m3xvis)=cff*cff1+dt(ng)*cff3
727 diau3wrk(i,j,k,m3yvis)=cff*cff2+dt(ng)*cff4
728#endif
729 END DO
730 END DO
731
732 DO j=jstrv,jend
733 DO i=istr,iend
734 cff=dt(ng)*0.25_r8*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
735 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
736 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
737 cff3=vfsx(i,j,k2)-vfsx(i,j,k1)
738 cff4=vfse(i,j,k2)-vfse(i,j,k1)
739 cff5=cff*(cff1-cff2)
740 cff6=dt(ng)*(cff3+cff4)
741 rvfrc(i,j)=rvfrc(i,j)+cff1-cff2+cff3+cff4
742 v(i,j,k,nnew)=v(i,j,k,nnew)+cff5+cff6
743#ifdef DIAGNOSTICS_UV
744 diarvfrc(i,j,3,m2hvis)=diarvfrc(i,j,3,m2hvis)+cff1-cff2+ &
745 & cff3+cff4
746 diarvfrc(i,j,3,m2xvis)=diarvfrc(i,j,3,m2xvis)+cff1+cff3
747 diarvfrc(i,j,3,m2yvis)=diarvfrc(i,j,3,m2yvis)-cff2+cff4
748 diav3wrk(i,j,k,m3hvis)=cff5+cff6
749 diav3wrk(i,j,k,m3xvis)= cff*cff1+dt(ng)*cff3
750 diav3wrk(i,j,k,m3yvis)=-cff*cff2+dt(ng)*cff4
751#endif
752 END DO
753 END DO
754 END IF
755 END DO k_loop
756!
757 RETURN
758 END SUBROUTINE uv3dmix2_geo_tile
759
760 END MODULE uv3dmix2_mod
type(t_coupling), dimension(:), allocatable coupling
type(t_diags), dimension(:), allocatable diags
Definition mod_diags.F:100
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
integer, parameter inlm
Definition mod_param.F:662
real(dp), dimension(:), allocatable dt
integer m3xvis
integer m2hvis
integer m3hvis
integer m3yvis
integer m2yvis
integer m2xvis
integer, dimension(:), allocatable nrhs
integer, dimension(:), allocatable nnew
subroutine uv3dmix2_geo_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nnew, pmask, rmask, umask, vmask, pmask_wet, rmask_wet, umask_wet, vmask_wet, om_p, om_r, om_u, om_v, on_p, on_r, on_u, on_v, pm, pn, hz, z_r, visc3d_r, visc2_p, visc2_r, diarufrc, diarvfrc, diau3wrk, diav3wrk, u, v, rufrc, rvfrc)
subroutine, public uv3dmix2(ng, tile)
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3