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