ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
uv3dmix4_mod Module Reference

Functions/Subroutines

subroutine, public uv3dmix4 (ng, tile)
 
subroutine uv3dmix4_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, uvis3d_r, vvis3d_r, visc3d_r, visc4_p, visc4_r, diarufrc, diarvfrc, diau3wrk, diav3wrk, u, v, rufrc, rvfrc)
 
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)
 

Function/Subroutine Documentation

◆ uv3dmix4()

subroutine public uv3dmix4_mod::uv3dmix4 ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 42 of file uv3dmix4_geo.h.

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, 33, __line__, myfile)
68#endif
69 CALL uv3dmix4_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# ifdef UV_U3ADV_SPLIT
99 & mixing(ng) % Uvis3d_r, &
100 & mixing(ng) % Vvis3d_r, &
101# else
102 & mixing(ng) % visc3d_r, &
103# endif
104#else
105 & mixing(ng) % visc4_p, &
106 & mixing(ng) % visc4_r, &
107#endif
108#ifdef DIAGNOSTICS_UV
109 & diags(ng) % DiaRUfrc, &
110 & diags(ng) % DiaRVfrc, &
111 & diags(ng) % DiaU3wrk, &
112 & diags(ng) % DiaV3wrk, &
113#endif
114 & ocean(ng) % u, &
115 & ocean(ng) % v, &
116 & coupling(ng) % rufrc, &
117 & coupling(ng) % rvfrc)
118#ifdef PROFILE
119 CALL wclock_off (ng, inlm, 33, __line__, myfile)
120#endif
121!
122 RETURN
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_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_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
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

References mod_coupling::coupling, mod_diags::diags, mod_grid::grid, mod_param::inlm, mod_mixing::mixing, mod_stepping::nnew, mod_stepping::nrhs, mod_ocean::ocean, uv3dmix4_geo_tile(), wclock_off(), and wclock_on().

Referenced by rhs3d_mod::rhs3d().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ uv3dmix4_geo_tile()

subroutine uv3dmix4_mod::uv3dmix4_geo_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,
integer, intent(in) nrhs,
integer, intent(in) nnew,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) rmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) rmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) umask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) vmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_p,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_r,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_u,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_p,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_r,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_u,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pm,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pn,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) uvis3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) vvis3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) visc3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) visc4_p,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) visc4_r,
real(r8), dimension(lbi:ubi,lbj:ubj,3,ndm2d-1), intent(inout) diarufrc,
real(r8), dimension(lbi:ubi,lbj:ubj,3,ndm2d-1), intent(inout) diarvfrc,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),ndm3d), intent(inout) diau3wrk,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),ndm3d), intent(inout) diav3wrk,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(inout) u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(inout) v,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(inout) rufrc,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(inout) rvfrc )
private

Definition at line 126 of file uv3dmix4_geo.h.

156!***********************************************************************
157!
158 USE mod_param
159 USE mod_ncparam
160 USE mod_scalars
161!
162! Imported variable declarations.
163!
164 integer, intent(in) :: ng, tile
165 integer, intent(in) :: LBi, UBi, LBj, UBj
166 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
167 integer, intent(in) :: nrhs, nnew
168
169#ifdef ASSUMED_SHAPE
170# ifdef MASKING
171 real(r8), intent(in) :: pmask(LBi:,LBj:)
172 real(r8), intent(in) :: rmask(LBi:,LBj:)
173 real(r8), intent(in) :: umask(LBi:,LBj:)
174 real(r8), intent(in) :: vmask(LBi:,LBj:)
175# endif
176# ifdef WET_DRY
177 real(r8), intent(in) :: pmask_wet(LBi:,LBj:)
178 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
179 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
180 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
181# endif
182 real(r8), intent(in) :: om_p(LBi:,LBj:)
183 real(r8), intent(in) :: om_r(LBi:,LBj:)
184 real(r8), intent(in) :: om_u(LBi:,LBj:)
185 real(r8), intent(in) :: om_v(LBi:,LBj:)
186 real(r8), intent(in) :: on_p(LBi:,LBj:)
187 real(r8), intent(in) :: on_r(LBi:,LBj:)
188 real(r8), intent(in) :: on_u(LBi:,LBj:)
189 real(r8), intent(in) :: on_v(LBi:,LBj:)
190 real(r8), intent(in) :: pm(LBi:,LBj:)
191 real(r8), intent(in) :: pn(LBi:,LBj:)
192 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
193 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
194# ifdef VISC_3DCOEF
195# ifdef UV_U3ADV_SPLIT
196 real(r8), intent(in) :: Uvis3d_r(LBi:,LBj:,:)
197 real(r8), intent(in) :: Vvis3d_r(LBi:,LBj:,:)
198# else
199 real(r8), intent(in) :: visc3d_r(LBi:,LBj:,:)
200# endif
201# else
202 real(r8), intent(in) :: visc4_p(LBi:,LBj:)
203 real(r8), intent(in) :: visc4_r(LBi:,LBj:)
204# endif
205# ifdef DIAGNOSTICS_UV
206 real(r8), intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
207 real(r8), intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
208 real(r8), intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
209 real(r8), intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
210# endif
211 real(r8), intent(inout) :: rufrc(LBi:,LBj:)
212 real(r8), intent(inout) :: rvfrc(LBi:,LBj:)
213 real(r8), intent(inout) :: u(LBi:,LBj:,:,:)
214 real(r8), intent(inout) :: v(LBi:,LBj:,:,:)
215#else
216# ifdef MASKING
217 real(r8), intent(in) :: pmask(LBi:UBi,LBj:UBj)
218 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
219 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
220 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
221# endif
222# ifdef WET_DRY
223 real(r8), intent(in) :: pmask_wet(LBi:UBi,LBj:UBj)
224 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
225 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
226 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
227# endif
228 real(r8), intent(in) :: om_p(LBi:UBi,LBj:UBj)
229 real(r8), intent(in) :: om_r(LBi:UBi,LBj:UBj)
230 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
231 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
232 real(r8), intent(in) :: on_p(LBi:UBi,LBj:UBj)
233 real(r8), intent(in) :: on_r(LBi:UBi,LBj:UBj)
234 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
235 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
236 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
237 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
238 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
239 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
240# ifdef VISC_3DCOEF
241# ifdef UV_U3ADV_SPLIT
242 real(r8), intent(in) :: Uvis3d_r(LBi:UBi,LBj:UBj,N(ng))
243 real(r8), intent(in) :: Vvis3d_r(LBi:UBi,LBj:UBj,N(ng))
244# else
245 real(r8), intent(in) :: visc3d_r(LBi:UBi,LBj:UBj,N(ng))
246# endif
247# else
248 real(r8), intent(in) :: visc4_p(LBi:UBi,LBj:UBj)
249 real(r8), intent(in) :: visc4_r(LBi:UBi,LBj:UBj)
250# endif
251# ifdef DIAGNOSTICS_UV
252 real(r8), intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
253 real(r8), intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
254 real(r8), intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
255 real(r8), intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
256# endif
257 real(r8), intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
258 real(r8), intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
259 real(r8), intent(inout) :: u(LBi:UBi,LBj:UBj,N(ng),2)
260 real(r8), intent(inout) :: v(LBi:UBi,LBj:UBj,N(ng),2)
261#endif
262!
263! Local variable declarations.
264!
265 integer :: i, j, k, k1, k2
266
267 real(r8) :: cff, fac1, fac2, pm_p, pn_p
268 real(r8) :: cff1, cff2, cff3, cff4
269 real(r8) :: cff5, cff6, cff7, cff8
270 real(r8) :: dmUdz, dnUdz, dmVdz, dnVdz
271#ifdef VISC_3DCOEF
272 real(r8) :: Uvis_p, Vvis_p, visc_p
273#endif
274
275 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapU
276 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapV
277
278 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
279 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
280 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
281 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
282
283 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: UFse
284 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: UFsx
285 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: VFse
286 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: VFsx
287 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dmUde
288 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dmVde
289 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dnUdx
290 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dnVdx
291 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dUdz
292 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dVdz
293 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_p
294 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZde_r
295 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_p
296 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dZdx_r
297
298#include "set_bounds.h"
299!
300!-----------------------------------------------------------------------
301! Compute horizontal biharmonic viscosity along geopotential
302! surfaces. The biharmonic operator is computed by applying
303! the harmonic operator twice.
304!-----------------------------------------------------------------------
305!
306! Compute horizontal and vertical gradients. Notice the recursive
307! blocking sequence. It is assumed here that the mixing coefficients
308! are the squared root of the biharmonic viscosity coefficient. For
309! momentum balance purposes, the thickness "Hz" appears only when
310! computing the second harmonic operator. The vertical placement of
311! the gradients is:
312!
313! dZdx_r, dZde_r, dnUdx, dmVde(:,:,k1) k rho-points
314! dZdx_r, dZde_r, dnUdx, dmVde(:,:,k2) k+1 rho-points
315! dZdx_p, dZde_p, dnVdx, dmUde(:,:,k1) k psi-points
316! dZdx_p, dZde_p, dnVdx, dmUde(:,:,k2) k+1 psi-points
317! UFse, UFsx, dUdz(:,:,k1) k-1/2 WU-points
318! UFse, UFsx, dUdz(:,:,k2) k+1/2 WU-points
319! VFse, VFsx, dVdz(:,:,k1) k-1/2 WV-points
320! VFse, VFsx, dVdz(:,:,k2) k+1/2 WV-points
321!
322 k2=1
323 k_loop1 : DO k=0,n(ng)
324 k1=k2
325 k2=3-k1
326 IF (k.lt.n(ng)) THEN
327!
328! Compute slopes (nondimensional) at RHO- and PSI-points.
329!
330 DO j=jstrm2,jendp2
331 DO i=istrum2,iendp2
332 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
333#ifdef MASKING
334 cff=cff*umask(i,j)
335#endif
336#ifdef WET_DRY
337 cff=cff*umask_wet(i,j)
338#endif
339 ufx(i,j)=cff*(z_r(i ,j,k+1)- &
340 & z_r(i-1,j,k+1))
341 END DO
342 END DO
343 DO j=jstrvm2,jendp2
344 DO i=istrm2,iendp2
345 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
346#ifdef MASKING
347 cff=cff*vmask(i,j)
348#endif
349#ifdef WET_DRY
350 cff=cff*vmask_wet(i,j)
351#endif
352 vfe(i,j)=cff*(z_r(i,j ,k+1)- &
353 & z_r(i,j-1,k+1))
354 END DO
355 END DO
356!
357 DO j=jstrm1,jendp2
358 DO i=istrm1,iendp2
359 dzdx_p(i,j,k2)=0.5_r8*(ufx(i,j-1)+ &
360 & ufx(i,j ))
361 dzde_p(i,j,k2)=0.5_r8*(vfe(i-1,j)+ &
362 & vfe(i ,j))
363 END DO
364 END DO
365 DO j=jstrvm2,jendp1
366 DO i=istrum2,iendp1
367 dzdx_r(i,j,k2)=0.5_r8*(ufx(i ,j)+ &
368 & ufx(i+1,j))
369 dzde_r(i,j,k2)=0.5_r8*(vfe(i,j )+ &
370 & vfe(i,j+1))
371 END DO
372 END DO
373!
374! Compute momentum horizontal (1/m/s) and vertical (1/s) gradients.
375!
376 DO j=jstrvm2,jendp1
377 DO i=istrum2,iendp1
378 cff=0.5_r8*pm(i,j)
379#ifdef MASKING
380 cff=cff*rmask(i,j)
381#endif
382#ifdef WET_DRY
383 cff=cff*rmask_wet(i,j)
384#endif
385 dnudx(i,j,k2)=cff*((pn(i ,j)+pn(i+1,j))* &
386 & u(i+1,j,k+1,nrhs)- &
387 & (pn(i-1,j)+pn(i ,j))* &
388 & u(i ,j,k+1,nrhs))
389 END DO
390 END DO
391
392 DO j=jstrm1,jendp2
393 DO i=istrm1,iendp2
394 cff=0.125_r8*(pn(i-1,j )+pn(i,j )+ &
395 & pn(i-1,j-1)+pn(i,j-1))
396#ifdef MASKING
397 cff=cff*pmask(i,j)
398#endif
399#ifdef WET_DRY
400 cff=cff*pmask_wet(i,j)
401#endif
402 dmude(i,j,k2)=cff*((pm(i-1,j )+pm(i,j ))* &
403 & u(i,j ,k+1,nrhs)- &
404 & (pm(i-1,j-1)+pm(i,j-1))* &
405 & u(i,j-1,k+1,nrhs))
406 END DO
407 END DO
408
409 DO j=jstrm1,jendp2
410 DO i=istrm1,iendp2
411 cff=0.125_r8*(pm(i-1,j )+pm(i,j )+ &
412 & pm(i-1,j-1)+pm(i,j-1))
413#ifdef MASKING
414 cff=cff*pmask(i,j)
415#endif
416#ifdef WET_DRY
417 cff=cff*pmask_wet(i,j)
418#endif
419 dnvdx(i,j,k2)=cff*((pn(i ,j-1)+pn(i ,j))* &
420 & v(i ,j,k+1,nrhs)- &
421 & (pn(i-1,j-1)+pn(i-1,j))* &
422 & v(i-1,j,k+1,nrhs))
423 END DO
424 END DO
425
426 DO j=jstrvm2,jendp1
427 DO i=istrum2,iendp1
428 cff=0.5_r8*pn(i,j)
429#ifdef MASKING
430 cff=cff*rmask(i,j)
431#endif
432#ifdef WET_DRY
433 cff=cff*rmask_wet(i,j)
434#endif
435 dmvde(i,j,k2)=cff*((pm(i,j )+pm(i,j+1))* &
436 & v(i,j+1,k+1,nrhs)- &
437 & (pm(i,j-1)+pm(i,j ))* &
438 & v(i,j ,k+1,nrhs))
439 END DO
440 END DO
441 END IF
442
443 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
444 DO j=jstrm2,jendp2
445 DO i=istrum2,iendp2
446 dudz(i,j,k2)=0.0_r8
447 END DO
448 END DO
449 DO j=jstrvm2,jendp2
450 DO i=istrm2,iendp2
451 dvdz(i,j,k2)=0.0_r8
452 END DO
453 END DO
454
455 DO j=jstrm1,jendp1
456 DO i=istrum1,iendp1
457 ufsx(i,j,k2)=0.0_r8
458 ufse(i,j,k2)=0.0_r8
459 END DO
460 END DO
461 DO j=jstrvm1,jendp1
462 DO i=istrm1,iendp1
463 vfsx(i,j,k2)=0.0_r8
464 vfse(i,j,k2)=0.0_r8
465 END DO
466 END DO
467 ELSE
468 DO j=jstrm2,jendp2
469 DO i=istrum2,iendp2
470 cff=1.0_r8/(0.5_r8*(z_r(i-1,j,k+1)- &
471 & z_r(i-1,j,k )+ &
472 & z_r(i ,j,k+1)- &
473 & z_r(i ,j,k )))
474 dudz(i,j,k2)=cff*(u(i,j,k+1,nrhs)- &
475 & u(i,j,k ,nrhs))
476 END DO
477 END DO
478
479 DO j=jstrvm2,jendp2
480 DO i=istrm2,iendp2
481 cff=1.0_r8/(0.5_r8*(z_r(i,j-1,k+1)- &
482 & z_r(i,j-1,k )+ &
483 & z_r(i,j ,k+1)- &
484 & z_r(i,j ,k )))
485 dvdz(i,j,k2)=cff*(v(i,j,k+1,nrhs)- &
486 & v(i,j,k ,nrhs))
487 END DO
488 END DO
489 END IF
490!
491! Compute components of the rotated viscous flux (m^4 s-^3/2) along
492! geopotential surfaces in the XI- and ETA-directions.
493!
494 IF (k.gt.0) THEN
495 DO j=jstrvm2,jendp1
496 DO i=istrum2,iendp1
497 cff1=min(dzdx_r(i,j,k1),0.0_r8)
498 cff2=max(dzdx_r(i,j,k1),0.0_r8)
499 cff3=min(dzde_r(i,j,k1),0.0_r8)
500 cff4=max(dzde_r(i,j,k1),0.0_r8)
501 cff=on_r(i,j)*(dnudx(i,j,k1)- &
502 & 0.5_r8*pn(i,j)* &
503 & (cff1*(dudz(i ,j,k1)+ &
504 & dudz(i+1,j,k2))+ &
505 & cff2*(dudz(i ,j,k2)+ &
506 & dudz(i+1,j,k1))))- &
507 & om_r(i,j)*(dmvde(i,j,k1)- &
508 & 0.5_r8*pm(i,j)* &
509 & (cff3*(dvdz(i,j ,k1)+ &
510 & dvdz(i,j+1,k2))+ &
511 & cff4*(dvdz(i,j ,k2)+ &
512 & dvdz(i,j+1,k1))))
513#ifdef MASKING
514 cff=cff*rmask(i,j)
515#endif
516#ifdef WET_DRY
517 cff=cff*rmask_wet(i,j)
518#endif
519#ifdef VISC_3DCOEF
520# ifdef UV_U3ADV_SPLIT
521 ufx(i,j)=on_r(i,j)*on_r(i,j)*uvis3d_r(i,j,k)*cff
522 vfe(i,j)=om_r(i,j)*om_r(i,j)*vvis3d_r(i,j,k)*cff
523# else
524 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc3d_r(i,j,k)*cff
525 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc3d_r(i,j,k)*cff
526# endif
527#else
528 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc4_r(i,j)*cff
529 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc4_r(i,j)*cff
530#endif
531 END DO
532 END DO
533
534 DO j=jstrm1,jendp2
535 DO i=istrm1,iendp2
536 pm_p=0.25_r8*(pm(i-1,j-1)+pm(i-1,j)+ &
537 & pm(i ,j-1)+pm(i ,j))
538 pn_p=0.25_r8*(pn(i-1,j-1)+pn(i-1,j)+ &
539 & pn(i ,j-1)+pn(i ,j))
540 cff1=min(dzdx_p(i,j,k1),0.0_r8)
541 cff2=max(dzdx_p(i,j,k1),0.0_r8)
542 cff3=min(dzde_p(i,j,k1),0.0_r8)
543 cff4=max(dzde_p(i,j,k1),0.0_r8)
544 cff=on_p(i,j)*(dnvdx(i,j,k1)- &
545 & 0.5_r8*pn_p* &
546 & (cff1*(dvdz(i-1,j,k1)+ &
547 & dvdz(i ,j,k2))+ &
548 & cff2*(dvdz(i-1,j,k2)+ &
549 & dvdz(i ,j,k1))))+ &
550 & om_p(i,j)*(dmude(i,j,k1)- &
551 & 0.5_r8*pm_p* &
552 & (cff3*(dudz(i,j-1,k1)+ &
553 & dudz(i,j ,k2))+ &
554 & cff4*(dudz(i,j-1,k2)+ &
555 & dudz(i,j ,k1))))
556#ifdef MASKING
557 cff=cff*pmask(i,j)
558#endif
559#ifdef WET_DRY
560 cff=cff*pmask_wet(i,j)
561#endif
562#ifdef VISC_3DCOEF
563# ifdef UV_U3ADV_SPLIT
564 uvis_p=0.25_r8* &
565 & (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* &
568 & (vvis3d_r(i-1,j-1,k)+vvis3d_r(i-1,j,k)+ &
569 & vvis3d_r(i ,j-1,k)+vvis3d_r(i ,j,k))
570 ufe(i,j)=om_p(i,j)*om_p(i,j)*uvis_p*cff
571 vfx(i,j)=on_p(i,j)*on_p(i,j)*vvis_p*cff
572# else
573 visc_p=0.25_r8* &
574 & (visc3d_r(i-1,j-1,k)+visc3d_r(i-1,j,k)+ &
575 & visc3d_r(i ,j-1,k)+visc3d_r(i ,j,k))
576 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc_p*cff
577 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc_p*cff
578# endif
579#else
580 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc4_p(i,j)*cff
581 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc4_p(i,j)*cff
582#endif
583 END DO
584 END DO
585!
586! Compute vertical flux (m^2 s^-3/2) due to sloping terrain-following
587! surfaces.
588!
589 IF (k.lt.n(ng)) THEN
590 DO j=jstrm1,jendp1
591 DO i=istrum1,iendp1
592#ifdef VISC_3DCOEF
593# ifdef UV_U3ADV_SPLIT
594 cff=0.125_r8* &
595 & (uvis3d_r(i-1,j,k )+uvis3d_r(i,j,k )+ &
596 & uvis3d_r(i-1,j,k+1)+uvis3d_r(i,j,k+1))
597# else
598 cff=0.125_r8* &
599 & (visc3d_r(i-1,j,k )+visc3d_r(i,j,k )+ &
600 & visc3d_r(i-1,j,k+1)+visc3d_r(i,j,k+1))
601# endif
602 fac1=cff*on_u(i,j)
603 fac2=cff*om_u(i,j)
604#else
605 cff=0.25_r8*(visc4_r(i-1,j)+visc4_r(i,j))
606 fac1=cff*on_u(i,j)
607 fac2=cff*om_u(i,j)
608#endif
609 cff=0.5_r8*(pn(i-1,j)+pn(i,j))
610 dnudz=cff*dudz(i,j,k2)
611 dnvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+ &
612 & dvdz(i ,j+1,k2)+ &
613 & dvdz(i-1,j ,k2)+ &
614 & dvdz(i ,j ,k2))
615 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
616 dmudz=cff*dudz(i,j,k2)
617 dmvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+ &
618 & dvdz(i ,j+1,k2)+ &
619 & dvdz(i-1,j ,k2)+ &
620 & dvdz(i ,j ,k2))
621
622 cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
623 cff2=min(dzdx_r(i ,j,k2),0.0_r8)
624 cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
625 cff4=max(dzdx_r(i ,j,k1),0.0_r8)
626 ufsx(i,j,k2)=fac1* &
627 & (cff1*(cff1*dnudz-dnudx(i-1,j,k1))+ &
628 & cff2*(cff2*dnudz-dnudx(i ,j,k2))+ &
629 & cff3*(cff3*dnudz-dnudx(i-1,j,k2))+ &
630 & cff4*(cff4*dnudz-dnudx(i ,j,k1)))
631
632 cff1=min(dzde_p(i,j ,k1),0.0_r8)
633 cff2=min(dzde_p(i,j+1,k2),0.0_r8)
634 cff3=max(dzde_p(i,j ,k2),0.0_r8)
635 cff4=max(dzde_p(i,j+1,k1),0.0_r8)
636 ufse(i,j,k2)=fac2* &
637 & (cff1*(cff1*dmudz-dmude(i,j ,k1))+ &
638 & cff2*(cff2*dmudz-dmude(i,j+1,k2))+ &
639 & cff3*(cff3*dmudz-dmude(i,j ,k2))+ &
640 & cff4*(cff4*dmudz-dmude(i,j+1,k1)))
641
642 cff1=min(dzde_p(i,j ,k1),0.0_r8)
643 cff2=min(dzde_p(i,j+1,k2),0.0_r8)
644 cff3=max(dzde_p(i,j ,k2),0.0_r8)
645 cff4=max(dzde_p(i,j+1,k1),0.0_r8)
646 cff5=min(dzdx_p(i,j ,k1),0.0_r8)
647 cff6=min(dzdx_p(i,j+1,k2),0.0_r8)
648 cff7=max(dzdx_p(i,j ,k2),0.0_r8)
649 cff8=max(dzdx_p(i,j+1,k1),0.0_r8)
650 ufsx(i,j,k2)=ufsx(i,j,k2)+ &
651 & fac1* &
652 & (cff1*(cff5*dnvdz-dnvdx(i,j ,k1))+ &
653 & cff2*(cff6*dnvdz-dnvdx(i,j+1,k2))+ &
654 & cff3*(cff7*dnvdz-dnvdx(i,j ,k2))+ &
655 & cff4*(cff8*dnvdz-dnvdx(i,j+1,k1)))
656
657 cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
658 cff2=min(dzdx_r(i ,j,k2),0.0_r8)
659 cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
660 cff4=max(dzdx_r(i ,j,k1),0.0_r8)
661 cff5=min(dzde_r(i-1,j,k1),0.0_r8)
662 cff6=min(dzde_r(i ,j,k2),0.0_r8)
663 cff7=max(dzde_r(i-1,j,k2),0.0_r8)
664 cff8=max(dzde_r(i ,j,k1),0.0_r8)
665 ufse(i,j,k2)=ufse(i,j,k2)- &
666 & fac2* &
667 & (cff1*(cff5*dmvdz-dmvde(i-1,j,k1))+ &
668 & cff2*(cff6*dmvdz-dmvde(i ,j,k2))+ &
669 & cff3*(cff7*dmvdz-dmvde(i-1,j,k2))+ &
670 & cff4*(cff8*dmvdz-dmvde(i ,j,k1)))
671 END DO
672 END DO
673!
674 DO j=jstrvm1,jendp1
675 DO i=istrm1,iendp1
676#ifdef VISC_3DCOEF
677# ifdef UV_U3ADV_SPLIT
678 cff=0.125_r8* &
679 & (vvis3d_r(i,j-1,k )+vvis3d_r(i,j,k )+ &
680 & vvis3d_r(i,j-1,k+1)+vvis3d_r(i,j,k+1))
681# else
682 cff=0.125_r8* &
683 & (visc3d_r(i,j-1,k )+visc3d_r(i,j,k )+ &
684 & visc3d_r(i,j-1,k+1)+visc3d_r(i,j,k+1))
685# endif
686 fac1=cff*on_v(i,j)
687 fac2=cff*om_v(i,j)
688#else
689 cff=0.25_r8*(visc4_r(i,j-1)+visc4_r(i,j))
690 fac1=cff*on_v(i,j)
691 fac2=cff*om_v(i,j)
692#endif
693 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
694 dnudz=cff*0.25_r8*(dudz(i ,j ,k2)+ &
695 & dudz(i+1,j ,k2)+ &
696 & dudz(i ,j-1,k2)+ &
697 & dudz(i+1,j-1,k2))
698 dnvdz=cff*dvdz(i,j,k2)
699 cff=0.5_r8*(pm(i,j-1)+pm(i,j))
700 dmudz=cff*0.25_r8*(dudz(i ,j ,k2)+ &
701 & dudz(i+1,j ,k2)+ &
702 & dudz(i ,j-1,k2)+ &
703 & dudz(i+1,j-1,k2))
704 dmvdz=cff*dvdz(i,j,k2)
705
706 cff1=min(dzdx_p(i ,j,k1),0.0_r8)
707 cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
708 cff3=max(dzdx_p(i ,j,k2),0.0_r8)
709 cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
710 vfsx(i,j,k2)=fac1* &
711 & (cff1*(cff1*dnvdz-dnvdx(i ,j,k1))+ &
712 & cff2*(cff2*dnvdz-dnvdx(i+1,j,k2))+ &
713 & cff3*(cff3*dnvdz-dnvdx(i ,j,k2))+ &
714 & cff4*(cff4*dnvdz-dnvdx(i+1,j,k1)))
715
716 cff1=min(dzde_r(i,j-1,k1),0.0_r8)
717 cff2=min(dzde_r(i,j ,k2),0.0_r8)
718 cff3=max(dzde_r(i,j-1,k2),0.0_r8)
719 cff4=max(dzde_r(i,j ,k1),0.0_r8)
720 vfse(i,j,k2)=fac2* &
721 & (cff1*(cff1*dmvdz-dmvde(i,j-1,k1))+ &
722 & cff2*(cff2*dmvdz-dmvde(i,j ,k2))+ &
723 & cff3*(cff3*dmvdz-dmvde(i,j-1,k2))+ &
724 & cff4*(cff4*dmvdz-dmvde(i,j ,k1)))
725
726 cff1=min(dzde_r(i,j-1,k1),0.0_r8)
727 cff2=min(dzde_r(i,j ,k2),0.0_r8)
728 cff3=max(dzde_r(i,j-1,k2),0.0_r8)
729 cff4=max(dzde_r(i,j ,k1),0.0_r8)
730 cff5=min(dzdx_r(i,j-1,k1),0.0_r8)
731 cff6=min(dzdx_r(i,j ,k2),0.0_r8)
732 cff7=max(dzdx_r(i,j-1,k2),0.0_r8)
733 cff8=max(dzdx_r(i,j ,k1),0.0_r8)
734 vfsx(i,j,k2)=vfsx(i,j,k2)- &
735 & fac1* &
736 & (cff1*(cff5*dnudz-dnudx(i,j-1,k1))+ &
737 & cff2*(cff6*dnudz-dnudx(i,j ,k2))+ &
738 & cff3*(cff7*dnudz-dnudx(i,j-1,k2))+ &
739 & cff4*(cff8*dnudz-dnudx(i,j ,k1)))
740
741 cff1=min(dzdx_p(i ,j,k1),0.0_r8)
742 cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
743 cff3=max(dzdx_p(i ,j,k2),0.0_r8)
744 cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
745 cff5=min(dzde_p(i ,j,k1),0.0_r8)
746 cff6=min(dzde_p(i+1,j,k2),0.0_r8)
747 cff7=max(dzde_p(i ,j,k2),0.0_r8)
748 cff8=max(dzde_p(i+1,j,k1),0.0_r8)
749 vfse(i,j,k2)=vfse(i,j,k2)+ &
750 & fac2* &
751 & (cff1*(cff5*dmudz-dmude(i ,j,k1))+ &
752 & cff2*(cff6*dmudz-dmude(i+1,j,k2))+ &
753 & cff3*(cff7*dmudz-dmude(i ,j,k2))+ &
754 & cff4*(cff8*dmudz-dmude(i+1,j,k1)))
755 END DO
756 END DO
757 END IF
758!
759! Compute first harmonic operator (m s^-3/2).
760!
761 DO j=jstrm1,jendp1
762 DO i=istrum1,iendp1
763 cff=0.125_r8*(pm(i-1,j)+pm(i,j))* &
764 & (pn(i-1,j)+pn(i,j))
765 cff1=1.0_r8/(0.5_r8*(hz(i-1,j,k)+hz(i,j,k)))
766 lapu(i,j,k)=cff*((pn(i-1,j)+pn(i,j))* &
767 (ufx(i,j)-ufx(i-1,j))+ &
768 & (pm(i-1,j)+pm(i,j))* &
769 & (ufe(i,j+1)-ufe(i,j)))+ &
770 & cff1*((ufsx(i,j,k2)+ufse(i,j,k2))- &
771 & (ufsx(i,j,k1)+ufse(i,j,k1)))
772#ifdef MASKING
773 lapu(i,j,k)=lapu(i,j,k)*umask(i,j)
774#endif
775#ifdef WET_DRY
776 lapu(i,j,k)=lapu(i,j,k)*umask_wet(i,j)
777#endif
778 END DO
779 END DO
780
781 DO j=jstrvm1,jendp1
782 DO i=istrm1,iendp1
783 cff=0.125_r8*(pm(i,j)+pm(i,j-1))* &
784 & (pn(i,j)+pn(i,j-1))
785 cff1=1.0_r8/(0.5_r8*(hz(i,j-1,k)+hz(i,j,k)))
786 lapv(i,j,k)=cff*((pn(i,j-1)+pn(i,j))* &
787 & (vfx(i+1,j)-vfx(i,j))- &
788 & (pm(i,j-1)+pm(i,j))* &
789 & (vfe(i,j)-vfe(i,j-1)))+ &
790 & cff1*((vfsx(i,j,k2)+vfse(i,j,k2))- &
791 & (vfsx(i,j,k1)+vfse(i,j,k1)))
792#ifdef MASKING
793 lapv(i,j,k)=lapv(i,j,k)*vmask(i,j)
794#endif
795#ifdef WET_DRY
796 lapv(i,j,k)=lapv(i,j,k)*vmask_wet(i,j)
797#endif
798 END DO
799 END DO
800 END IF
801 END DO k_loop1
802!
803! Apply boundary conditions (closed or gradient; except periodic)
804! to the first harmonic operator.
805!
806 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
807 IF (domain(ng)%Western_Edge(tile)) THEN
808 IF (lbc(iwest,isuvel,ng)%closed) THEN
809 DO k=1,n(ng)
810 DO j=jstrm1,jendp1
811 lapu(istru-1,j,k)=0.0_r8
812 END DO
813 END DO
814 ELSE
815 DO k=1,n(ng)
816 DO j=jstrm1,jendp1
817 lapu(istru-1,j,k)=lapu(istru,j,k)
818 END DO
819 END DO
820 END IF
821 IF (lbc(iwest,isvvel,ng)%closed) THEN
822 DO k=1,n(ng)
823 DO j=jstrvm1,jendp1
824 lapv(istr-1,j,k)=gamma2(ng)*lapv(istr,j,k)
825 END DO
826 END DO
827 ELSE
828 DO k=1,n(ng)
829 DO j=jstrvm1,jendp1
830 lapv(istr-1,j,k)=0.0_r8
831 END DO
832 END DO
833 END IF
834 END IF
835 END IF
836!
837 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
838 IF (domain(ng)%Eastern_Edge(tile)) THEN
839 IF (lbc(ieast,isuvel,ng)%closed) THEN
840 DO k=1,n(ng)
841 DO j=jstrm1,jendp1
842 lapu(iend+1,j,k)=0.0_r8
843 END DO
844 END DO
845 ELSE
846 DO k=1,n(ng)
847 DO j=jstrm1,jendp1
848 lapu(iend+1,j,k)=lapu(iend,j,k)
849 END DO
850 END DO
851 END IF
852 IF (lbc(ieast,isvvel,ng)%closed) THEN
853 DO k=1,n(ng)
854 DO j=jstrvm1,jendp1
855 lapv(iend+1,j,k)=gamma2(ng)*lapv(iend,j,k)
856 END DO
857 END DO
858 ELSE
859 DO k=1,n(ng)
860 DO j=jstrvm1,jendp1
861 lapv(iend+1,j,k)=0.0_r8
862 END DO
863 END DO
864 END IF
865 END IF
866 END IF
867!
868 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
869 IF (domain(ng)%Southern_Edge(tile)) THEN
870 IF (lbc(isouth,isuvel,ng)%closed) THEN
871 DO k=1,n(ng)
872 DO i=istrum1,iendp1
873 lapu(i,jstr-1,k)=gamma2(ng)*lapu(i,jstr,k)
874 END DO
875 END DO
876 ELSE
877 DO k=1,n(ng)
878 DO i=istrum1,iendp1
879 lapu(i,jstr-1,k)=0.0_r8
880 END DO
881 END DO
882 END IF
883 IF (lbc(isouth,isvvel,ng)%closed) THEN
884 DO k=1,n(ng)
885 DO i=istrm1,iendp1
886 lapv(i,jstrv-1,k)=0.0_r8
887 END DO
888 END DO
889 ELSE
890 DO k=1,n(ng)
891 DO i=istrm1,iendp1
892 lapv(i,jstrv-1,k)=lapv(i,jstrv,k)
893 END DO
894 END DO
895 END IF
896 END IF
897 END IF
898!
899 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
900 IF (domain(ng)%Northern_Edge(tile)) THEN
901 IF (lbc(inorth,isuvel,ng)%closed) THEN
902 DO k=1,n(ng)
903 DO i=istrum1,iendp1
904 lapu(i,jend+1,k)=gamma2(ng)*lapu(i,jend,k)
905 END DO
906 END DO
907 ELSE
908 DO k=1,n(ng)
909 DO i=istrum1,iendp1
910 lapu(i,jend+1,k)=0.0_r8
911 END DO
912 END DO
913 END IF
914 IF (lbc(inorth,isvvel,ng)%closed) THEN
915 DO k=1,n(ng)
916 DO i=istrm1,iendp1
917 lapv(i,jend+1,k)=0.0_r8
918 END DO
919 END DO
920 ELSE
921 DO k=1,n(ng)
922 DO i=istrm1,iendp1
923 lapv(i,jend+1,k)=lapv(i,jend,k)
924 END DO
925 END DO
926 END IF
927 END IF
928 END IF
929!
930 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
931 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
932 IF (domain(ng)%SouthWest_Corner(tile)) THEN
933 DO k=1,n(ng)
934 lapu(istr ,jstr-1,k)=0.5_r8*(lapu(istr+1,jstr-1,k)+ &
935 & lapu(istr ,jstr ,k))
936 lapv(istr-1,jstr ,k)=0.5_r8*(lapv(istr-1,jstr+1,k)+ &
937 & lapv(istr ,jstr ,k))
938 END DO
939 END IF
940 END IF
941
942 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
943 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
944 IF (domain(ng)%SouthEast_Corner(tile)) THEN
945 DO k=1,n(ng)
946 lapu(iend+1,jstr-1,k)=0.5_r8*(lapu(iend ,jstr-1,k)+ &
947 & lapu(iend+1,jstr ,k))
948 lapv(iend+1,jstr ,k)=0.5_r8*(lapv(iend ,jstr ,k)+ &
949 & lapv(iend+1,jstr+1,k))
950 END DO
951 END IF
952 END IF
953
954 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
955 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
956 IF (domain(ng)%NorthWest_Corner(tile)) THEN
957 DO k=1,n(ng)
958 lapu(istr ,jend+1,k)=0.5_r8*(lapu(istr+1,jend+1,k)+ &
959 & lapu(istr ,jend ,k))
960 lapv(istr-1,jend+1,k)=0.5_r8*(lapv(istr ,jend+1,k)+ &
961 & lapv(istr-1,jend ,k))
962 END DO
963 END IF
964 END IF
965
966 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
967 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
968 IF (domain(ng)%NorthEast_Corner(tile)) THEN
969 DO k=1,n(ng)
970 lapu(iend+1,jend+1,k)=0.5_r8*(lapu(iend ,jend+1,k)+ &
971 & lapu(iend+1,jend ,k))
972 lapv(iend+1,jend+1,k)=0.5_r8*(lapv(iend ,jend+1,k)+ &
973 & lapv(iend+1,jend ,k))
974 END DO
975 END IF
976 END IF
977!
978! Compute horizontal and vertical gradients associated with the
979! second rotated harmonic operator.
980!
981 k2=1
982 k_loop2 : DO k=0,n(ng)
983 k1=k2
984 k2=3-k1
985 IF (k.lt.n(ng)) THEN
986!
987! Compute slopes (nondimensional) at RHO- and PSI-points.
988!
989 DO j=jstr-1,jend+1
990 DO i=istru-1,iend+1
991 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
992#ifdef MASKING
993 cff=cff*umask(i,j)
994#endif
995#ifdef WET_DRY
996 cff=cff*umask_wet(i,j)
997#endif
998 ufx(i,j)=cff*(z_r(i ,j,k+1)- &
999 & z_r(i-1,j,k+1))
1000 END DO
1001 END DO
1002 DO j=jstrv-1,jend+1
1003 DO i=istr-1,iend+1
1004 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
1005#ifdef MASKING
1006 cff=cff*vmask(i,j)
1007#endif
1008#ifdef WET_DRY
1009 cff=cff*vmask_wet(i,j)
1010#endif
1011 vfe(i,j)=cff*(z_r(i,j ,k+1)- &
1012 & z_r(i,j-1,k+1))
1013 END DO
1014 END DO
1015!
1016 DO j=jstr,jend+1
1017 DO i=istr,iend+1
1018 dzdx_p(i,j,k2)=0.5_r8*(ufx(i,j-1)+ &
1019 & ufx(i,j ))
1020 dzde_p(i,j,k2)=0.5_r8*(vfe(i-1,j)+ &
1021 & vfe(i ,j))
1022 END DO
1023 END DO
1024 DO j=jstrv-1,jend
1025 DO i=istru-1,iend
1026 dzdx_r(i,j,k2)=0.5_r8*(ufx(i ,j)+ &
1027 & ufx(i+1,j))
1028 dzde_r(i,j,k2)=0.5_r8*(vfe(i,j )+ &
1029 & vfe(i,j+1))
1030 END DO
1031 END DO
1032!
1033! Compute momentum horizontal (m^-1 s^-3/2) and vertical (s^-3/2)
1034! gradients.
1035!
1036 DO j=jstrv-1,jend
1037 DO i=istru-1,iend
1038 cff=0.5_r8*pm(i,j)
1039#ifdef MASKING
1040 cff=cff*rmask(i,j)
1041#endif
1042#ifdef WET_DRY
1043 cff=cff*rmask_wet(i,j)
1044#endif
1045 dnudx(i,j,k2)=cff*((pn(i ,j)+pn(i+1,j))* &
1046 & lapu(i+1,j,k+1)- &
1047 & (pn(i-1,j)+pn(i ,j))* &
1048 & lapu(i ,j,k+1))
1049 END DO
1050 END DO
1051
1052 DO j=jstr,jend+1
1053 DO i=istr,iend+1
1054 cff=0.125_r8*(pn(i-1,j )+pn(i,j )+ &
1055 & pn(i-1,j-1)+pn(i,j-1))
1056#ifdef MASKING
1057 cff=cff*pmask(i,j)
1058#endif
1059#ifdef WET_DRY
1060 cff=cff*pmask_wet(i,j)
1061#endif
1062 dmude(i,j,k2)=cff*((pm(i-1,j )+pm(i,j ))* &
1063 & lapu(i,j ,k+1)- &
1064 & (pm(i-1,j-1)+pm(i,j-1))* &
1065 & lapu(i,j-1,k+1))
1066 END DO
1067 END DO
1068
1069 DO j=jstr,jend+1
1070 DO i=istr,iend+1
1071 cff=0.125_r8*(pm(i-1,j )+pm(i,j )+ &
1072 & pm(i-1,j-1)+pm(i,j-1))
1073#ifdef MASKING
1074 cff=cff*pmask(i,j)
1075#endif
1076#ifdef WET_DRY
1077 cff=cff*pmask_wet(i,j)
1078#endif
1079 dnvdx(i,j,k2)=cff*((pn(i ,j-1)+pn(i ,j))* &
1080 & lapv(i ,j,k+1)- &
1081 & (pn(i-1,j-1)+pn(i-1,j))* &
1082 & lapv(i-1,j,k+1))
1083 END DO
1084 END DO
1085
1086 DO j=jstrv-1,jend
1087 DO i=istru-1,iend
1088 cff=0.5_r8*pn(i,j)
1089#ifdef MASKING
1090 cff=cff*rmask(i,j)
1091#endif
1092#ifdef WET_DRY
1093 cff=cff*rmask_wet(i,j)
1094#endif
1095 dmvde(i,j,k2)=cff*((pm(i,j )+pm(i,j+1))* &
1096 & lapv(i,j+1,k+1)- &
1097 & (pm(i,j-1)+pm(i,j ))* &
1098 & lapv(i,j ,k+1))
1099 END DO
1100 END DO
1101 END IF
1102
1103 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
1104 DO j=jstr-1,jend+1
1105 DO i=istru-1,iend+1
1106 dudz(i,j,k2)=0.0_r8
1107 END DO
1108 END DO
1109 DO j=jstrv-1,jend+1
1110 DO i=istr-1,iend+1
1111 dvdz(i,j,k2)=0.0_r8
1112 END DO
1113 END DO
1114
1115 DO j=jstr,jend
1116 DO i=istru,iend
1117 ufsx(i,j,k2)=0.0_r8
1118 ufse(i,j,k2)=0.0_r8
1119 END DO
1120 END DO
1121 DO j=jstrv,jend
1122 DO i=istr,iend
1123 vfsx(i,j,k2)=0.0_r8
1124 vfse(i,j,k2)=0.0_r8
1125 END DO
1126 END DO
1127 ELSE
1128 DO j=jstr-1,jend+1
1129 DO i=istru-1,iend+1
1130 cff=1.0_r8/(0.5_r8*(z_r(i-1,j,k+1)- &
1131 & z_r(i-1,j,k )+ &
1132 & z_r(i ,j,k+1)- &
1133 & z_r(i ,j,k )))
1134 dudz(i,j,k2)=cff*(lapu(i,j,k+1)- &
1135 & lapu(i,j,k ))
1136 END DO
1137 END DO
1138 DO j=jstrv-1,jend+1
1139 DO i=istr-1,iend+1
1140 cff=1.0_r8/(0.5_r8*(z_r(i,j-1,k+1)- &
1141 & z_r(i,j-1,k )+ &
1142 & z_r(i,j ,k+1)- &
1143 & z_r(i,j ,k )))
1144 dvdz(i,j,k2)=cff*(lapv(i,j,k+1)- &
1145 & lapv(i,j,k ))
1146 END DO
1147 END DO
1148 END IF
1149!
1150! Compute components of the rotated viscous flux (m5/s2) along
1151! geopotential surfaces in the XI- and ETA-directions.
1152!
1153 IF (k.gt.0) THEN
1154 DO j=jstrv-1,jend
1155 DO i=istru-1,iend
1156 cff1=min(dzdx_r(i,j,k1),0.0_r8)
1157 cff2=max(dzdx_r(i,j,k1),0.0_r8)
1158 cff3=min(dzde_r(i,j,k1),0.0_r8)
1159 cff4=max(dzde_r(i,j,k1),0.0_r8)
1160 cff=hz(i,j,k)* &
1161 & (on_r(i,j)*(dnudx(i,j,k1)- &
1162 & 0.5_r8*pn(i,j)* &
1163 & (cff1*(dudz(i ,j,k1)+ &
1164 & dudz(i+1,j,k2))+ &
1165 & cff2*(dudz(i ,j,k2)+ &
1166 & dudz(i+1,j,k1))))- &
1167 & om_r(i,j)*(dmvde(i,j,k1)- &
1168 & 0.5_r8*pm(i,j)* &
1169 & (cff3*(dvdz(i,j ,k1)+ &
1170 & dvdz(i,j+1,k2))+ &
1171 & cff4*(dvdz(i,j ,k2)+ &
1172 & dvdz(i,j+1,k1)))))
1173#ifdef MASKING
1174 cff=cff*rmask(i,j)
1175#endif
1176#ifdef WET_DRY
1177 cff=cff*rmask_wet(i,j)
1178#endif
1179#ifdef VISC_3DCOEF
1180# ifdef UV_U3ADV_SPLIT
1181 ufx(i,j)=on_r(i,j)*on_r(i,j)*uvis3d_r(i,j,k)*cff
1182 vfe(i,j)=om_r(i,j)*om_r(i,j)*vvis3d_r(i,j,k)*cff
1183# else
1184 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc3d_r(i,j,k)*cff
1185 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc3d_r(i,j,k)*cff
1186# endif
1187#else
1188 ufx(i,j)=on_r(i,j)*on_r(i,j)*visc4_r(i,j)*cff
1189 vfe(i,j)=om_r(i,j)*om_r(i,j)*visc4_r(i,j)*cff
1190#endif
1191 END DO
1192 END DO
1193
1194 DO j=jstr,jend+1
1195 DO i=istr,iend+1
1196 pm_p=0.25_r8*(pm(i-1,j-1)+pm(i-1,j)+ &
1197 & pm(i ,j-1)+pm(i ,j))
1198 pn_p=0.25_r8*(pn(i-1,j-1)+pn(i-1,j)+ &
1199 & pn(i ,j-1)+pn(i ,j))
1200 cff1=min(dzdx_p(i,j,k1),0.0_r8)
1201 cff2=max(dzdx_p(i,j,k1),0.0_r8)
1202 cff3=min(dzde_p(i,j,k1),0.0_r8)
1203 cff4=max(dzde_p(i,j,k1),0.0_r8)
1204 cff=0.25_r8* &
1205 & (hz(i-1,j ,k)+hz(i,j ,k)+ &
1206 & hz(i-1,j-1,k)+hz(i,j-1,k))* &
1207 & (on_p(i,j)*(dnvdx(i,j,k1)- &
1208 & 0.5_r8*pn_p* &
1209 & (cff1*(dvdz(i-1,j,k1)+ &
1210 & dvdz(i ,j,k2))+ &
1211 & cff2*(dvdz(i-1,j,k2)+ &
1212 & dvdz(i ,j,k1))))+ &
1213 & om_p(i,j)*(dmude(i,j,k1)- &
1214 & 0.5_r8*pm_p* &
1215 & (cff3*(dudz(i,j-1,k1)+ &
1216 & dudz(i,j ,k2))+ &
1217 & cff4*(dudz(i,j-1,k2)+ &
1218 & dudz(i,j ,k1)))))
1219#ifdef MASKING
1220 cff=cff*pmask(i,j)
1221#endif
1222#ifdef WET_DRY
1223 cff=cff*pmask_wet(i,j)
1224#endif
1225#ifdef VISC_3DCOEF
1226# ifdef UV_U3ADV_SPLIT
1227 uvis_p=0.25_r8* &
1228 & (uvis3d_r(i-1,j-1,k)+uvis3d_r(i-1,j,k)+ &
1229 & uvis3d_r(i ,j-1,k)+uvis3d_r(i ,j,k))
1230 vvis_p=0.25_r8* &
1231 & (vvis3d_r(i-1,j-1,k)+vvis3d_r(i-1,j,k)+ &
1232 & vvis3d_r(i ,j-1,k)+vvis3d_r(i ,j,k))
1233 ufe(i,j)=om_p(i,j)*om_p(i,j)*uvis_p*cff
1234 vfx(i,j)=on_p(i,j)*on_p(i,j)*vvis_p*cff
1235# else
1236 visc_p=0.25_r8* &
1237 & (visc3d_r(i-1,j-1,k)+visc3d_r(i-1,j,k)+ &
1238 & visc3d_r(i ,j-1,k)+visc3d_r(i ,j,k))
1239 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc_p*cff
1240 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc_p*cff
1241# endif
1242#else
1243 ufe(i,j)=om_p(i,j)*om_p(i,j)*visc4_p(i,j)*cff
1244 vfx(i,j)=on_p(i,j)*on_p(i,j)*visc4_p(i,j)*cff
1245#endif
1246 END DO
1247 END DO
1248!
1249! Compute vertical flux (m2/s2) due to sloping terrain-following
1250! surfaces.
1251!
1252 IF (k.lt.n(ng)) THEN
1253 DO j=jstr,jend
1254 DO i=istru,iend
1255#ifdef VISC_3DCOEF
1256# ifdef UV_U3ADV_SPLIT
1257 cff=0.125_r8* &
1258 & (uvis3d_r(i-1,j,k )+uvis3d_r(i,j,k )+ &
1259 & uvis3d_r(i-1,j,k+1)+uvis3d_r(i,j,k+1))
1260# else
1261 cff=0.125_r8* &
1262 & (visc3d_r(i-1,j,k )+visc3d_r(i,j,k )+ &
1263 & visc3d_r(i-1,j,k+1)+visc3d_r(i,j,k+1))
1264# endif
1265 fac1=cff*on_u(i,j)
1266 fac2=cff*om_u(i,j)
1267#else
1268 cff=0.25_r8*(visc4_r(i-1,j)+visc4_r(i,j))
1269 fac1=cff*on_u(i,j)
1270 fac2=cff*om_u(i,j)
1271#endif
1272 cff=0.5_r8*(pn(i-1,j)+pn(i,j))
1273 dnudz=cff*dudz(i,j,k2)
1274 dnvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+ &
1275 & dvdz(i ,j+1,k2)+ &
1276 & dvdz(i-1,j ,k2)+ &
1277 & dvdz(i ,j ,k2))
1278 cff=0.5_r8*(pm(i-1,j)+pm(i,j))
1279 dmudz=cff*dudz(i,j,k2)
1280 dmvdz=cff*0.25_r8*(dvdz(i-1,j+1,k2)+ &
1281 & dvdz(i ,j+1,k2)+ &
1282 & dvdz(i-1,j ,k2)+ &
1283 & dvdz(i ,j ,k2))
1284
1285 cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
1286 cff2=min(dzdx_r(i ,j,k2),0.0_r8)
1287 cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
1288 cff4=max(dzdx_r(i ,j,k1),0.0_r8)
1289 ufsx(i,j,k2)=fac1* &
1290 & (cff1*(cff1*dnudz-dnudx(i-1,j,k1))+ &
1291 & cff2*(cff2*dnudz-dnudx(i ,j,k2))+ &
1292 & cff3*(cff3*dnudz-dnudx(i-1,j,k2))+ &
1293 & cff4*(cff4*dnudz-dnudx(i ,j,k1)))
1294
1295 cff1=min(dzde_p(i,j ,k1),0.0_r8)
1296 cff2=min(dzde_p(i,j+1,k2),0.0_r8)
1297 cff3=max(dzde_p(i,j ,k2),0.0_r8)
1298 cff4=max(dzde_p(i,j+1,k1),0.0_r8)
1299 ufse(i,j,k2)=fac2* &
1300 & (cff1*(cff1*dmudz-dmude(i,j ,k1))+ &
1301 & cff2*(cff2*dmudz-dmude(i,j+1,k2))+ &
1302 & cff3*(cff3*dmudz-dmude(i,j ,k2))+ &
1303 & cff4*(cff4*dmudz-dmude(i,j+1,k1)))
1304
1305 cff1=min(dzde_p(i,j ,k1),0.0_r8)
1306 cff2=min(dzde_p(i,j+1,k2),0.0_r8)
1307 cff3=max(dzde_p(i,j ,k2),0.0_r8)
1308 cff4=max(dzde_p(i,j+1,k1),0.0_r8)
1309 cff5=min(dzdx_p(i,j ,k1),0.0_r8)
1310 cff6=min(dzdx_p(i,j+1,k2),0.0_r8)
1311 cff7=max(dzdx_p(i,j ,k2),0.0_r8)
1312 cff8=max(dzdx_p(i,j+1,k1),0.0_r8)
1313 ufsx(i,j,k2)=ufsx(i,j,k2)+ &
1314 & fac1* &
1315 & (cff1*(cff5*dnvdz-dnvdx(i,j ,k1))+ &
1316 & cff2*(cff6*dnvdz-dnvdx(i,j+1,k2))+ &
1317 & cff3*(cff7*dnvdz-dnvdx(i,j ,k2))+ &
1318 & cff4*(cff8*dnvdz-dnvdx(i,j+1,k1)))
1319
1320 cff1=min(dzdx_r(i-1,j,k1),0.0_r8)
1321 cff2=min(dzdx_r(i ,j,k2),0.0_r8)
1322 cff3=max(dzdx_r(i-1,j,k2),0.0_r8)
1323 cff4=max(dzdx_r(i ,j,k1),0.0_r8)
1324 cff5=min(dzde_r(i-1,j,k1),0.0_r8)
1325 cff6=min(dzde_r(i ,j,k2),0.0_r8)
1326 cff7=max(dzde_r(i-1,j,k2),0.0_r8)
1327 cff8=max(dzde_r(i ,j,k1),0.0_r8)
1328 ufse(i,j,k2)=ufse(i,j,k2)- &
1329 & fac2* &
1330 & (cff1*(cff5*dmvdz-dmvde(i-1,j,k1))+ &
1331 & cff2*(cff6*dmvdz-dmvde(i ,j,k2))+ &
1332 & cff3*(cff7*dmvdz-dmvde(i-1,j,k2))+ &
1333 & cff4*(cff8*dmvdz-dmvde(i ,j,k1)))
1334 END DO
1335 END DO
1336!
1337 DO j=jstrv,jend
1338 DO i=istr,iend
1339#ifdef VISC_3DCOEF
1340# ifdef UV_U3ADV_SPLIT
1341 cff=0.125_r8* &
1342 & (vvis3d_r(i,j-1,k )+vvis3d_r(i,j,k )+ &
1343 & vvis3d_r(i,j-1,k+1)+vvis3d_r(i,j,k+1))
1344# else
1345 cff=0.125_r8* &
1346 & (visc3d_r(i,j-1,k )+visc3d_r(i,j,k )+ &
1347 & visc3d_r(i,j-1,k+1)+visc3d_r(i,j,k+1))
1348# endif
1349 fac1=cff*on_v(i,j)
1350 fac2=cff*om_v(i,j)
1351#else
1352 cff=0.25_r8*(visc4_r(i,j-1)+visc4_r(i,j))
1353 fac1=cff*on_v(i,j)
1354 fac2=cff*om_v(i,j)
1355#endif
1356 cff=0.5_r8*(pn(i,j-1)+pn(i,j))
1357 dnudz=cff*0.25_r8*(dudz(i ,j ,k2)+ &
1358 & dudz(i+1,j ,k2)+ &
1359 & dudz(i ,j-1,k2)+ &
1360 & dudz(i+1,j-1,k2))
1361 dnvdz=cff*dvdz(i,j,k2)
1362 cff=0.5_r8*(pm(i,j-1)+pm(i,j))
1363 dmudz=cff*0.25_r8*(dudz(i ,j ,k2)+ &
1364 & dudz(i+1,j ,k2)+ &
1365 & dudz(i ,j-1,k2)+ &
1366 & dudz(i+1,j-1,k2))
1367 dmvdz=cff*dvdz(i,j,k2)
1368
1369 cff1=min(dzdx_p(i ,j,k1),0.0_r8)
1370 cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
1371 cff3=max(dzdx_p(i ,j,k2),0.0_r8)
1372 cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
1373 vfsx(i,j,k2)=fac1* &
1374 & (cff1*(cff1*dnvdz-dnvdx(i ,j,k1))+ &
1375 & cff2*(cff2*dnvdz-dnvdx(i+1,j,k2))+ &
1376 & cff3*(cff3*dnvdz-dnvdx(i ,j,k2))+ &
1377 & cff4*(cff4*dnvdz-dnvdx(i+1,j,k1)))
1378
1379 cff1=min(dzde_r(i,j-1,k1),0.0_r8)
1380 cff2=min(dzde_r(i,j ,k2),0.0_r8)
1381 cff3=max(dzde_r(i,j-1,k2),0.0_r8)
1382 cff4=max(dzde_r(i,j ,k1),0.0_r8)
1383 vfse(i,j,k2)=fac2* &
1384 & (cff1*(cff1*dmvdz-dmvde(i,j-1,k1))+ &
1385 & cff2*(cff2*dmvdz-dmvde(i,j ,k2))+ &
1386 & cff3*(cff3*dmvdz-dmvde(i,j-1,k2))+ &
1387 & cff4*(cff4*dmvdz-dmvde(i,j ,k1)))
1388
1389 cff1=min(dzde_r(i,j-1,k1),0.0_r8)
1390 cff2=min(dzde_r(i,j ,k2),0.0_r8)
1391 cff3=max(dzde_r(i,j-1,k2),0.0_r8)
1392 cff4=max(dzde_r(i,j ,k1),0.0_r8)
1393 cff5=min(dzdx_r(i,j-1,k1),0.0_r8)
1394 cff6=min(dzdx_r(i,j ,k2),0.0_r8)
1395 cff7=max(dzdx_r(i,j-1,k2),0.0_r8)
1396 cff8=max(dzdx_r(i,j ,k1),0.0_r8)
1397 vfsx(i,j,k2)=vfsx(i,j,k2)- &
1398 & fac1* &
1399 & (cff1*(cff5*dnudz-dnudx(i,j-1,k1))+ &
1400 & cff2*(cff6*dnudz-dnudx(i,j ,k2))+ &
1401 & cff3*(cff7*dnudz-dnudx(i,j-1,k2))+ &
1402 & cff4*(cff8*dnudz-dnudx(i,j ,k1)))
1403
1404 cff1=min(dzdx_p(i ,j,k1),0.0_r8)
1405 cff2=min(dzdx_p(i+1,j,k2),0.0_r8)
1406 cff3=max(dzdx_p(i ,j,k2),0.0_r8)
1407 cff4=max(dzdx_p(i+1,j,k1),0.0_r8)
1408 cff5=min(dzde_p(i ,j,k1),0.0_r8)
1409 cff6=min(dzde_p(i+1,j,k2),0.0_r8)
1410 cff7=max(dzde_p(i ,j,k2),0.0_r8)
1411 cff8=max(dzde_p(i+1,j,k1),0.0_r8)
1412 vfse(i,j,k2)=vfse(i,j,k2)+ &
1413 & fac2* &
1414 & (cff1*(cff5*dmudz-dmude(i ,j,k1))+ &
1415 & cff2*(cff6*dmudz-dmude(i+1,j,k2))+ &
1416 & cff3*(cff7*dmudz-dmude(i ,j,k2))+ &
1417 & cff4*(cff8*dmudz-dmude(i+1,j,k1)))
1418 END DO
1419 END DO
1420 END IF
1421!
1422! Time-step biharmonic, geopotential viscosity term. Notice that
1423! momentum at this stage is HzU and HzV and has m2/s units. Add
1424! contribution for barotropic forcing terms.
1425#ifdef DIAGNOSTICS_UV
1426! The rotated vertical term cannot be split from the horizontal
1427! terms because of the 2D/3D momentum coupling.
1428#endif
1429!
1430 DO j=jstr,jend
1431 DO i=istru,iend
1432 cff=dt(ng)*0.25_r8*(pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))
1433 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
1434 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
1435 cff3=ufsx(i,j,k2)-ufsx(i,j,k1)
1436 cff4=ufse(i,j,k2)-ufse(i,j,k1)
1437 cff5=cff*(cff1+cff2)
1438 cff6=dt(ng)*(cff3+cff4)
1439 rufrc(i,j)=rufrc(i,j)-cff1-cff2-cff3-cff4
1440 u(i,j,k,nnew)=u(i,j,k,nnew)-cff5-cff6
1441#ifdef DIAGNOSTICS_UV
1442 diarufrc(i,j,3,m2hvis)=diarufrc(i,j,3,m2hvis)-cff1-cff2- &
1443 & cff3-cff4
1444 diarufrc(i,j,3,m2xvis)=diarufrc(i,j,3,m2xvis)-cff1-cff3
1445 diarufrc(i,j,3,m2yvis)=diarufrc(i,j,3,m2yvis)-cff2-cff4
1446 diau3wrk(i,j,k,m3hvis)=-cff5-cff6
1447 diau3wrk(i,j,k,m3xvis)=-cff*cff1-dt(ng)*cff3
1448 diau3wrk(i,j,k,m3yvis)=-cff*cff2-dt(ng)*cff4
1449#endif
1450 END DO
1451 END DO
1452
1453 DO j=jstrv,jend
1454 DO i=istr,iend
1455 cff=dt(ng)*0.25_r8*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
1456 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
1457 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
1458 cff3=vfsx(i,j,k2)-vfsx(i,j,k1)
1459 cff4=vfse(i,j,k2)-vfse(i,j,k1)
1460 cff5=cff*(cff1-cff2)
1461 cff6=dt(ng)*(cff3+cff4)
1462 rvfrc(i,j)=rvfrc(i,j)-cff1+cff2-cff3-cff4
1463 v(i,j,k,nnew)=v(i,j,k,nnew)-cff5-cff6
1464#ifdef DIAGNOSTICS_UV
1465 diarvfrc(i,j,3,m2hvis)=diarvfrc(i,j,3,m2hvis)-cff1+cff2- &
1466 & cff3-cff4
1467 diarvfrc(i,j,3,m2xvis)=diarvfrc(i,j,3,m2xvis)-cff1-cff3
1468 diarvfrc(i,j,3,m2yvis)=diarvfrc(i,j,3,m2yvis)+cff2-cff4
1469 diav3wrk(i,j,k,m3hvis)=-cff5-cff6
1470 diav3wrk(i,j,k,m3xvis)=-cff*cff1-dt(ng)*cff3
1471 diav3wrk(i,j,k,m3yvis)= cff*cff2-dt(ng)*cff4
1472#endif
1473 END DO
1474 END DO
1475 END IF
1476 END DO k_loop2
1477!
1478 RETURN
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_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::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

References mod_scalars::compositegrid, mod_param::domain, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::gamma2, mod_scalars::ieast, mod_scalars::inorth, mod_scalars::isouth, mod_ncparam::isuvel, mod_ncparam::isvvel, mod_scalars::iwest, mod_param::lbc, mod_scalars::m2hvis, mod_scalars::m2xvis, mod_scalars::m2yvis, mod_scalars::m3hvis, mod_scalars::m3xvis, mod_scalars::m3yvis, and mod_scalars::nsperiodic.

Referenced by uv3dmix4().

Here is the caller graph for this function:

◆ uv3dmix4_s_tile()

subroutine uv3dmix4_mod::uv3dmix4_s_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,
integer, intent(in) nrhs,
integer, intent(in) nnew,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pmask,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pmask_wet,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) hz,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_p,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) om_r,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_p,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) on_r,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pm,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pmon_p,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pmon_r,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pn,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pnom_p,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pnom_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) uvis3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) vvis3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) visc3d_r,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) visc4_p,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) visc4_r,
real(r8), dimension(lbi:ubi,lbj:ubj,3,ndm2d-1), intent(inout) diarufrc,
real(r8), dimension(lbi:ubi,lbj:ubj,3,ndm2d-1), intent(inout) diarvfrc,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),ndm3d), intent(inout) diau3wrk,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),ndm3d), intent(inout) diav3wrk,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(inout) rufrc,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(inout) rvfrc,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(inout) u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(inout) v )
private

Definition at line 119 of file uv3dmix4_s.h.

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
mod_param::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456

References mod_scalars::compositegrid, mod_param::domain, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::gamma2, mod_scalars::ieast, mod_scalars::inorth, mod_scalars::isouth, mod_ncparam::isuvel, mod_ncparam::isvvel, mod_scalars::iwest, mod_param::lbc, mod_param::lm, mod_scalars::m2hvis, mod_scalars::m2xvis, mod_scalars::m2yvis, mod_scalars::m3hvis, mod_scalars::m3xvis, mod_scalars::m3yvis, mod_param::mm, and mod_scalars::nsperiodic.