ROMS
Loading...
Searching...
No Matches
step3d_uv.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE step3d_uv_mod
3
4#if defined NONLINEAR && defined SOLVE3D
5!
6!git $Id$
7!=======================================================================
8! Copyright (c) 2002-2025 The ROMS Group !
9! Licensed under a MIT/X style license Hernan G. Arango !
10! See License_ROMS.md Alexander F. Shchepetkin !
11!==================================================== John C. Warner ===
12! !
13! This subroutine time-steps the nonlinear horizontal momentum !
14! equations. The vertical viscosity terms are time-stepped using !
15! an implicit algorithm. !
16! !
17!=======================================================================
18!
19 USE mod_param
20 USE mod_coupling
21# ifdef DIAGNOSTICS_UV
22 USE mod_diags
23# endif
24 USE mod_forces
25 USE mod_grid
26 USE mod_mixing
27 USE mod_ncparam
28 USE mod_ocean
29 USE mod_scalars
30 USE mod_sources
31!
32 USE exchange_2d_mod
33 USE exchange_3d_mod
34# ifdef DISTRIBUTE
35 USE mp_exchange_mod, ONLY : mp_exchange2d, mp_exchange3d
36# endif
37 USE uv_var_change_mod, ONLY : uv_c2a_grid
38 USE u3dbc_mod, ONLY : u3dbc_tile
39 USE v3dbc_mod, ONLY : v3dbc_tile
40!
41 implicit none
42!
43 PRIVATE
44 PUBLIC :: step3d_uv
45!
46 CONTAINS
47!
48!***********************************************************************
49 SUBROUTINE step3d_uv (ng, tile)
50!***********************************************************************
51!
52 USE mod_stepping
53!
54! Imported variable declarations.
55!
56 integer, intent(in) :: ng, tile
57!
58! Local variable declarations.
59!
60 character (len=*), parameter :: myfile = &
61 & __FILE__
62!
63# include "tile.h"
64!
65# ifdef PROFILE
66 CALL wclock_on (ng, inlm, 34, __line__, myfile)
67# endif
68 CALL step3d_uv_tile (ng, tile, &
69 & lbi, ubi, lbj, ubj, &
70 & imins, imaxs, jmins, jmaxs, &
71 & nrhs(ng), nstp(ng), nnew(ng), &
72# if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
73 & knew(ng), &
74# endif
75# ifdef MASKING
76 & grid(ng) % umask, &
77 & grid(ng) % vmask, &
78# endif
79# ifdef WET_DRY
80 & grid(ng) % umask_wet, &
81 & grid(ng) % vmask_wet, &
82# endif
83 & grid(ng) % om_v, &
84 & grid(ng) % on_u, &
85# ifdef OMEGA_IMPLICIT
86 & grid(ng) % om_u, &
87 & grid(ng) % on_v, &
88# endif
89 & grid(ng) % pm, &
90 & grid(ng) % pn, &
91 & grid(ng) % Hz, &
92 & grid(ng) % z_r, &
93 & grid(ng) % z_w, &
94 & mixing(ng) % Akv, &
95 & coupling(ng) % DU_avg1, &
96 & coupling(ng) % DV_avg1, &
97 & coupling(ng) % DU_avg2, &
98 & coupling(ng) % DV_avg2, &
99 & ocean(ng) % ru, &
100 & ocean(ng) % rv, &
101# ifdef DIAGNOSTICS_UV
102 & diags(ng) % DiaU2wrk, &
103 & diags(ng) % DiaV2wrk, &
104 & diags(ng) % DiaU2int, &
105 & diags(ng) % DiaV2int, &
106 & diags(ng) % DiaU3wrk, &
107 & diags(ng) % DiaV3wrk, &
108 & diags(ng) % DiaRU, &
109 & diags(ng) % DiaRV, &
110# endif
111 & ocean(ng) % u, &
112 & ocean(ng) % v, &
113 & ocean(ng) % ubar, &
114 & ocean(ng) % vbar, &
115# ifdef WEC
116 & ocean(ng) % ubar_stokes, &
117 & ocean(ng) % vbar_stokes, &
118 & ocean(ng) % u_stokes, &
119 & ocean(ng) % v_stokes, &
120# endif
121# ifdef OMEGA_IMPLICIT
122 & ocean(ng) % Wi, &
123# endif
124 & grid(ng) % Huon, &
125 & grid(ng) % Hvom)
126# ifdef PROFILE
127 CALL wclock_off (ng, inlm, 34, __line__, myfile)
128# endif
129!
130 RETURN
131 END SUBROUTINE step3d_uv
132!
133!***********************************************************************
134 SUBROUTINE step3d_uv_tile (ng, tile, &
135 & LBi, UBi, LBj, UBj, &
136 & IminS, ImaxS, JminS, JmaxS, &
137 & nrhs, nstp, nnew, &
138# if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
139 & knew, &
140# endif
141# ifdef MASKING
142 & umask, vmask, &
143# endif
144# ifdef WET_DRY
145 & umask_wet, vmask_wet, &
146# endif
147 & om_v, on_u, &
148# ifdef OMEGA_IMPLICIT
149 & om_u, on_v, &
150# endif
151 & pm, pn, &
152 & Hz, z_r, z_w, &
153 & Akv, &
154 & DU_avg1, DV_avg1, &
155 & DU_avg2, DV_avg2, &
156 & ru, rv, &
157# ifdef DIAGNOSTICS_UV
158 & DiaU2wrk, DiaV2wrk, &
159 & DiaU2int, DiaV2int, &
160 & DiaU3wrk, DiaV3wrk, &
161 & DiaRU, DiaRV, &
162# endif
163 & u, v, &
164 & ubar, vbar, &
165# ifdef WEC
166 & ubar_stokes, vbar_stokes, &
167 & u_stokes, v_stokes, &
168# endif
169# ifdef OMEGA_IMPLICIT
170 & Wi, &
171# endif
172 & Huon, Hvom)
173!***********************************************************************
174!
175! Imported variable declarations.
176!
177 integer, intent(in) :: ng, tile
178 integer, intent(in) :: LBi, UBi, LBj, UBj
179 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
180 integer, intent(in) :: nrhs, nstp, nnew
181# if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
182 integer, intent(in) :: knew
183# endif
184!
185# ifdef ASSUMED_SHAPE
186# ifdef MASKING
187 real(r8), intent(in) :: umask(LBi:,LBj:)
188 real(r8), intent(in) :: vmask(LBi:,LBj:)
189# endif
190# ifdef WET_DRY
191 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
192 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
193# endif
194 real(r8), intent(in) :: om_v(LBi:,LBj:)
195 real(r8), intent(in) :: on_u(LBi:,LBj:)
196# ifdef OMEGA_IMPLICIT
197 real(r8), intent(in) :: om_u(LBi:,LBj:)
198 real(r8), intent(in) :: on_v(LBi:,LBj:)
199# endif
200 real(r8), intent(in) :: pm(LBi:,LBj:)
201 real(r8), intent(in) :: pn(LBi:,LBj:)
202 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
203 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
204 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
205 real(r8), intent(in) :: Akv(LBi:,LBj:,0:)
206 real(r8), intent(in) :: DU_avg1(LBi:,LBj:)
207 real(r8), intent(in) :: DV_avg1(LBi:,LBj:)
208 real(r8), intent(in) :: DU_avg2(LBi:,LBj:)
209 real(r8), intent(in) :: DV_avg2(LBi:,LBj:)
210# ifdef WEC
211 real(r8), intent(in) :: ubar_stokes(LBi:,LBj:)
212 real(r8), intent(in) :: vbar_stokes(LBi:,LBj:)
213# endif
214# ifdef OMEGA_IMPLICIT
215 real(r8), intent(in) :: Wi(LBi:,LBj:,0:)
216# endif
217 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
218 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
219# ifdef DIAGNOSTICS_UV
220 real(r8), intent(inout) :: DiaU2wrk(LBi:,LBj:,:)
221 real(r8), intent(inout) :: DiaV2wrk(LBi:,LBj:,:)
222 real(r8), intent(inout) :: DiaU2int(LBi:,LBj:,:)
223 real(r8), intent(inout) :: DiaV2int(LBi:,LBj:,:)
224 real(r8), intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
225 real(r8), intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
226 real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
227 real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
228# endif
229 real(r8), intent(inout) :: u(LBi:,LBj:,:,:)
230 real(r8), intent(inout) :: v(LBi:,LBj:,:,:)
231# ifdef WEC
232 real(r8), intent(inout) :: u_stokes(LBi:,LBj:,:)
233 real(r8), intent(inout) :: v_stokes(LBi:,LBj:,:)
234# endif
235 real(r8), intent(out) :: ubar(LBi:,LBj:,:)
236 real(r8), intent(out) :: vbar(LBi:,LBj:,:)
237 real(r8), intent(out) :: Huon(LBi:,LBj:,:)
238 real(r8), intent(out) :: Hvom(LBi:,LBj:,:)
239
240# else
241
242# ifdef MASKING
243 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
244 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
245# endif
246# ifdef WET_DRY
247 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
248 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
249# endif
250 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
251 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
252# ifdef OMEGA_IMPLICIT
253 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
254 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
255# endif
256 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
257 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
258 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
259 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
260 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
261 real(r8), intent(in) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
262 real(r8), intent(in) :: DU_avg1(LBi:UBi,LBj:UBj)
263 real(r8), intent(in) :: DV_avg1(LBi:UBi,LBj:UBj)
264 real(r8), intent(in) :: DU_avg2(LBi:UBi,LBj:UBj)
265 real(r8), intent(in) :: DV_avg2(LBi:UBi,LBj:UBj)
266# ifdef WEC
267 real(r8), intent(in) :: ubar_stokes(LBi:UBi,LBj:UBj)
268 real(r8), intent(in) :: vbar_stokes(LBi:UBi,LBj:UBj)
269# endif
270# ifdef OMEGA_IMPLICIT
271 real(r8), intent(in) :: Wi(LBi:UBi,LBj:UBj,0:N(ng))
272# endif
273 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:N(ng),2)
274 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:N(ng),2)
275# ifdef DIAGNOSTICS_UV
276 real(r8), intent(inout) :: DiaU2wrk(LBi:UBi,LBj:UBj,NDM2d)
277 real(r8), intent(inout) :: DiaV2wrk(LBi:UBi,LBj:UBj,NDM2d)
278 real(r8), intent(inout) :: DiaU2int(LBi:UBi,LBj:UBj,NDM2d)
279 real(r8), intent(inout) :: DiaV2int(LBi:UBi,LBj:UBj,NDM2d)
280 real(r8), intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
281 real(r8), intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
282 real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
283 real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
284# endif
285 real(r8), intent(inout) :: u(LBi:UBi,LBj:UBj,N(ng),2)
286 real(r8), intent(inout) :: v(LBi:UBi,LBj:UBj,N(ng),2)
287# ifdef WEC
288 real(r8), intent(inout) :: u_stokes(LBi:UBi,LBj:UBj,N(ng))
289 real(r8), intent(inout) :: v_stokes(LBi:UBi,LBj:UBj,N(ng))
290# endif
291 real(r8), intent(out) :: ubar(LBi:UBi,LBj:UBj,:)
292 real(r8), intent(out) :: vbar(LBi:UBi,LBj:UBj,:)
293 real(r8), intent(out) :: Huon(LBi:UBi,LBj:UBj,N(ng))
294 real(r8), intent(out) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
295# endif
296!
297! Local variable declarations.
298!
299 integer :: i, idiag, is, j, k
300!
301 real(r8) :: cff, cff1, cff2
302!
303 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: AK
304 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
305 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
306 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
307 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
308# ifdef OMEGA_IMPLICIT
309 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCmin
310 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCmax
311 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: WK
312# endif
313# ifdef WEC
314 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CFs
315 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DCs
316# endif
317 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hzk
318 real(r8), dimension(IminS:ImaxS,N(ng)) :: oHz
319# ifdef DIAGNOSTICS_UV
320 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: wrk
321 real(r8), dimension(IminS:ImaxS,1:NDM2d-1) :: Dwrk
322# endif
323
324# include "set_bounds.h"
325!
326!-----------------------------------------------------------------------
327! Time step momentum equation in the XI-direction.
328!-----------------------------------------------------------------------
329!
330 DO j=jstr,jend
331 DO i=istru,iend
332 ak(i,0)=0.5_r8*(akv(i-1,j,0)+ &
333 & akv(i ,j,0))
334 DO k=1,n(ng)
335 ak(i,k)=0.5_r8*(akv(i-1,j,k)+ &
336 & akv(i ,j,k))
337 hzk(i,k)=0.5_r8*(hz(i-1,j,k)+ &
338 & hz(i ,j,k))
339# if defined SPLINES_VVISC || defined DIAGNOSTICS_UV
340 ohz(i,k)=1.0_r8/hzk(i,k)
341# endif
342 END DO
343 END DO
344!
345! Time step right-hand-side terms.
346!
347 IF (iic(ng).eq.ntfirst(ng)) THEN
348 cff=0.25_r8*dt(ng)
349 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
350 cff=0.25_r8*dt(ng)*3.0_r8/2.0_r8
351 ELSE
352 cff=0.25_r8*dt(ng)*23.0_r8/12.0_r8
353 END IF
354 DO i=istru,iend
355 dc(i,0)=cff*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
356 END DO
357 DO k=1,n(ng)
358 DO i=istru,iend
359 u(i,j,k,nnew)=u(i,j,k,nnew)+ &
360 & dc(i,0)*ru(i,j,k,nrhs)
361# ifdef SPLINES_VVISC
362 u(i,j,k,nnew)=u(i,j,k,nnew)*ohz(i,k)
363# endif
364# ifdef DIAGNOSTICS_UV
365 DO idiag=1,m3pgrd
366 diau3wrk(i,j,k,idiag)=(diau3wrk(i,j,k,idiag)+ &
367 & dc(i,0)*diaru(i,j,k,nrhs,idiag))* &
368 & ohz(i,k)
369 END DO
370# if defined UV_VIS2 || defined UV_VIS4
371 diau3wrk(i,j,k,m3xvis)=diau3wrk(i,j,k,m3xvis)*ohz(i,k)
372 diau3wrk(i,j,k,m3yvis)=diau3wrk(i,j,k,m3yvis)*ohz(i,k)
373 diau3wrk(i,j,k,m3hvis)=diau3wrk(i,j,k,m3hvis)*ohz(i,k)
374# endif
375 diau3wrk(i,j,k,m3vvis)=diau3wrk(i,j,k,m3vvis)*ohz(i,k)
376# ifdef BODYFORCE
377 diau3wrk(i,j,k,m3vvis)=diau3wrk(i,j,k,m3vvis)+ &
378 & dc(i,0)*diaru(i,j,k,nrhs,m3vvis)* &
379 & ohz(i,k)
380# endif
381 diau3wrk(i,j,k,m3rate)=diau3wrk(i,j,k,m3rate)*ohz(i,k)
382# endif
383 END DO
384 END DO
385
386# ifdef SPLINES_VVISC
387!
388! Use conservative, parabolic spline reconstruction of vertical
389! viscosity derivatives. Then, time step vertical viscosity term
390! implicitly by solving a tridiagonal system.
391!
392 cff1=1.0_r8/6.0_r8
393 DO k=1,n(ng)-1
394 DO i=istru,iend
395 fc(i,k)=cff1*hzk(i,k )-dt(ng)*ak(i,k-1)*ohz(i,k )
396 cf(i,k)=cff1*hzk(i,k+1)-dt(ng)*ak(i,k+1)*ohz(i,k+1)
397 END DO
398 END DO
399 DO i=istru,iend
400 cf(i,0)=0.0_r8
401 dc(i,0)=0.0_r8
402 END DO
403!
404! LU decomposition and forward substitution.
405!
406 cff1=1.0_r8/3.0_r8
407 DO k=1,n(ng)-1
408 DO i=istru,iend
409 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
410 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
411 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
412 cf(i,k)=cff*cf(i,k)
413 dc(i,k)=cff*(u(i,j,k+1,nnew)-u(i,j,k,nnew)- &
414 & fc(i,k)*dc(i,k-1))
415 END DO
416 END DO
417!
418! Backward substitution.
419!
420 DO i=istru,iend
421 dc(i,n(ng))=0.0_r8
422 END DO
423 DO k=n(ng)-1,1,-1
424 DO i=istru,iend
425 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
426 END DO
427 END DO
428!
429 DO k=1,n(ng)
430 DO i=istru,iend
431 dc(i,k)=dc(i,k)*ak(i,k)
432 cff=dt(ng)*ohz(i,k)*(dc(i,k)-dc(i,k-1))
433 u(i,j,k,nnew)=u(i,j,k,nnew)+cff
434# ifdef DIAGNOSTICS_UV
435 diau3wrk(i,j,k,m3vvis)=diau3wrk(i,j,k,m3vvis)+cff
436# endif
437 END DO
438 END DO
439# else
440!
441! Compute off-diagonal coefficients [lambda*dt*Akv/Hz] for the
442! implicit vertical viscosity term at horizontal U-points and
443! vertical W-points.
444!
445 cff=-lambda*dt(ng)/0.5_r8
446 DO k=1,n(ng)-1
447 DO i=istru,iend
448 cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
449 & z_r(i,j,k )-z_r(i-1,j,k ))
450 fc(i,k)=cff*cff1*ak(i,k)
451 END DO
452 END DO
453 DO i=istru,iend
454 fc(i,0)=0.0_r8
455 fc(i,n(ng))=0.0_r8
456 END DO
457!
458! Solve the tridiagonal system.
459!
460 DO k=1,n(ng)
461 DO i=istru,iend
462 dc(i,k)=u(i,j,k,nnew)
463 bc(i,k)=hzk(i,k)-fc(i,k)-fc(i,k-1)
464 END DO
465 END DO
466 DO i=istru,iend
467 cff=1.0_r8/bc(i,1)
468 cf(i,1)=cff*fc(i,1)
469 dc(i,1)=cff*dc(i,1)
470 END DO
471 DO k=2,n(ng)-1
472 DO i=istru,iend
473 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
474 cf(i,k)=cff*fc(i,k)
475 dc(i,k)=cff*(dc(i,k)-fc(i,k-1)*dc(i,k-1))
476 END DO
477 END DO
478!
479! Compute new solution by back substitution.
480!
481 DO i=istru,iend
482# ifdef DIAGNOSTICS_UV
483 wrk(i,n(ng))=u(i,j,n(ng),nnew)*ohz(i,n(ng))
484# endif
485 dc(i,n(ng))=(dc(i,n(ng))-fc(i,n(ng)-1)*dc(i,n(ng)-1))/ &
486 & (bc(i,n(ng))-fc(i,n(ng)-1)*cf(i,n(ng)-1))
487 u(i,j,n(ng),nnew)=dc(i,n(ng))
488# ifdef DIAGNOSTICS_UV
489 diau3wrk(i,j,n(ng),m3vvis)=diau3wrk(i,j,n(ng),m3vvis)+ &
490 & u(i,j,n(ng),nnew)-wrk(i,n(ng))
491# endif
492 END DO
493 DO k=n(ng)-1,1,-1
494 DO i=istru,iend
495# ifdef DIAGNOSTICS_UV
496 wrk(i,k)=u(i,j,k,nnew)*ohz(i,k)
497# endif
498 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
499 u(i,j,k,nnew)=dc(i,k)
500# ifdef DIAGNOSTICS_UV
501 diau3wrk(i,j,k,m3vvis)=diau3wrk(i,j,k,m3vvis)+ &
502 & u(i,j,k,nnew)-wrk(i,k)
503# endif
504 END DO
505 END DO
506# endif
507# ifdef OMEGA_IMPLICIT
508!
509! Adaptive, Courant-number based implicit vertical advection
510! contribution for u-momentum.
511!
512 DO i=istru,iend
513 wk(i,0)=0.5_r8*(wi(i-1,j,0)+ &
514 & wi(i ,j,0))
515 DO k=1,n(ng)
516 wk(i,k)=0.5_r8*(wi(i-1,j,k)+ &
517 & wi(i ,j,k))
518 hzk(i,k)=0.5_r8*(hz(i-1,j,k)+ &
519 & hz(i ,j,k))
520 END DO
521 END DO
522!
523! Compute off-diagonal coefficients [dt*Wi*pm*pn] for the
524! implicit vertical viscosity term at horizontal U-points and
525! vertical W-points.
526!
527 cff=dt(ng)
528 DO k=1,n(ng)-1
529 DO i=istru,iend
530 cff1=cff/(on_u(i,j)*om_u(i,j))
531 fcmax(i,k)=max(wk(i,k),0.0_r8)*cff1
532 fcmin(i,k)=min(wk(i,k),0.0_r8)*cff1
533 END DO
534 END DO
535 DO i=istru,iend
536 fcmax(i,0)=0.0_r8
537 fcmin(i,0)=0.0_r8
538 fcmax(i,n(ng))=0.0_r8
539 fcmin(i,n(ng))=0.0_r8
540 END DO
541!
542! Solve the tridiagonal system.
543!
544 DO k=1,n(ng)
545 DO i=istru,iend
546 bc(i,k)=hzk(i,k)+fcmax(i,k)-fcmin(i,k-1)
547 dc(i,k)=u(i,j,k,nnew)*hzk(i,k)
548 END DO
549 END DO
550 DO i=istru,iend
551 cff=1.0_r8/bc(i,1)
552 cf(i,1)=cff*fcmin(i,1)
553 dc(i,1)=cff*dc(i,1)
554 END DO
555 DO k=2,n(ng)-1
556 DO i=istru,iend
557 cff=1.0_r8/(bc(i,k)+fcmax(i,k-1)*cf(i,k-1))
558 cf(i,k)=cff*fcmin(i,k)
559 dc(i,k)=cff*(dc(i,k)+fcmax(i,k-1)*dc(i,k-1))
560 END DO
561 END DO
562!
563! Compute new solution by back substitution.
564!
565 DO i=istru,iend
566# ifdef DIAGNOSTICS_UV
567 cff1=u(i,j,n(ng),nnew)
568# endif
569 cff=1.0_r8/(bc(i,n(ng))+fcmax(i,n(ng)-1)*cf(i,n(ng)-1))
570 dc(i,n(ng))=cff*(dc(i,n(ng))+ &
571 & fcmax(i,n(ng)-1)*dc(i,n(ng)-1))
572 u(i,j,n(ng),nnew)=dc(i,n(ng))
573# ifdef DIAGNOSTICS_UV
574 diaru(i,j,n(ng),nrhs,m3vadv)=diaru(i,j,n(ng),nrhs,m3vadv)+ &
575 & u(i,j,n(ng),nnew)-cff1
576# endif
577 END DO
578!
579 DO k=n(ng)-1,1,-1
580 DO i=istru,iend
581# ifdef DIAGNOSTICS_UV
582 cff1=u(i,j,k,nnew)
583# endif
584 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
585 u(i,j,k,nnew)=dc(i,k)
586# ifdef DIAGNOSTICS_UV
587 diaru(i,j,k,nrhs,m3vadv)=diaru(i,j,k,nrhs,m3vadv)+ &
588 & u(i,j,k,nnew)-cff1
589# endif
590 END DO
591 END DO
592# endif
593!
594! Replace INTERIOR POINTS incorrect vertical mean with more accurate
595! barotropic component, ubar=DU_avg1/(D*on_u). Recall that, D=CF(:,0).
596!
597 DO i=istru,iend
598 cf(i,0)=hzk(i,1)
599 dc(i,0)=u(i,j,1,nnew)*hzk(i,1)
600# ifdef WEC
601 dcs(i,0)=u_stokes(i,j,1)*hzk(i,1)
602# endif
603# ifdef DIAGNOSTICS_UV
604 dwrk(i,m2pgrd)=diau3wrk(i,j,1,m3pgrd)*hzk(i,1)
605 dwrk(i,m2bstr)=diau3wrk(i,j,1,m3vvis)*hzk(i,1)
606# ifdef UV_COR
607 dwrk(i,m2fcor)=diau3wrk(i,j,1,m3fcor)*hzk(i,1)
608# endif
609# if defined UV_VIS2 || defined UV_VIS4
610 dwrk(i,m2xvis)=diau3wrk(i,j,1,m3xvis)*hzk(i,1)
611 dwrk(i,m2yvis)=diau3wrk(i,j,1,m3yvis)*hzk(i,1)
612 dwrk(i,m2hvis)=diau3wrk(i,j,1,m3hvis)*hzk(i,1)
613# endif
614# ifdef UV_ADV
615 dwrk(i,m2xadv)=diau3wrk(i,j,1,m3xadv)*hzk(i,1)
616 dwrk(i,m2yadv)=diau3wrk(i,j,1,m3yadv)*hzk(i,1)
617 dwrk(i,m2hadv)=diau3wrk(i,j,1,m3hadv)*hzk(i,1)
618# endif
619# ifdef WEC_VF
620 dwrk(i,m2hjvf)=diau3wrk(i,j,1,m3hjvf)*hzk(i,1)
621 dwrk(i,m2kvrf)=diau3wrk(i,j,1,m3kvrf)*hzk(i,1)
622# ifdef UV_COR
623 dwrk(i,m2fsco)=diau3wrk(i,j,1,m3fsco)*hzk(i,1)
624# endif
625# ifdef BOTTOM_STREAMING
626 dwrk(i,m2bstm)=diau3wrk(i,j,1,m3bstm)*hzk(i,1)
627# endif
628# ifdef SURFACE_STREAMING
629 dwrk(i,m2sstm)=diau3wrk(i,j,1,m3sstm)*hzk(i,1)
630# endif
631 dwrk(i,m2wrol)=diau3wrk(i,j,1,m3wrol)*hzk(i,1)
632 dwrk(i,m2wbrk)=diau3wrk(i,j,1,m3wbrk)*hzk(i,1)
633# endif
634# endif
635 END DO
636 DO k=2,n(ng)
637 DO i=istru,iend
638 cf(i,0)=cf(i,0)+hzk(i,k)
639 dc(i,0)=dc(i,0)+u(i,j,k,nnew)*hzk(i,k)
640# ifdef WEC
641 dcs(i,0)=dcs(i,0)+u_stokes(i,j,k)*hzk(i,k)
642# endif
643# ifdef DIAGNOSTICS_UV
644 dwrk(i,m2pgrd)=dwrk(i,m2pgrd)+ &
645 & diau3wrk(i,j,k,m3pgrd)*hzk(i,k)
646 dwrk(i,m2bstr)=dwrk(i,m2bstr)+ &
647 & diau3wrk(i,j,k,m3vvis)*hzk(i,k)
648# ifdef UV_COR
649 dwrk(i,m2fcor)=dwrk(i,m2fcor)+ &
650 & diau3wrk(i,j,k,m3fcor)*hzk(i,k)
651# endif
652# if defined UV_VIS2 || defined UV_VIS4
653 dwrk(i,m2xvis)=dwrk(i,m2xvis)+ &
654 & diau3wrk(i,j,k,m3xvis)*hzk(i,k)
655 dwrk(i,m2yvis)=dwrk(i,m2yvis)+ &
656 & diau3wrk(i,j,k,m3yvis)*hzk(i,k)
657 dwrk(i,m2hvis)=dwrk(i,m2hvis)+ &
658 & diau3wrk(i,j,k,m3hvis)*hzk(i,k)
659# endif
660# ifdef UV_ADV
661 dwrk(i,m2xadv)=dwrk(i,m2xadv)+ &
662 & diau3wrk(i,j,k,m3xadv)*hzk(i,k)
663 dwrk(i,m2yadv)=dwrk(i,m2yadv)+ &
664 & diau3wrk(i,j,k,m3yadv)*hzk(i,k)
665 dwrk(i,m2hadv)=dwrk(i,m2hadv)+ &
666 & diau3wrk(i,j,k,m3hadv)*hzk(i,k)
667# endif
668# ifdef WEC_VF
669 dwrk(i,m2hjvf)=dwrk(i,m2hjvf)+ &
670 & diau3wrk(i,j,k,m3hjvf)*hzk(i,k)
671 dwrk(i,m2kvrf)=dwrk(i,m2kvrf)+ &
672 & diau3wrk(i,j,k,m3kvrf)*hzk(i,k)
673# ifdef UV_COR
674 dwrk(i,m2fsco)=dwrk(i,m2fsco)+ &
675 & diau3wrk(i,j,k,m3fsco)*hzk(i,k)
676# endif
677# ifdef BOTTOM_STREAMING
678 dwrk(i,m2bstm)=dwrk(i,m2bstm)+ &
679 & diau3wrk(i,j,k,m3bstm)*hzk(i,k)
680# endif
681# ifdef SURFACE_STREAMING
682 dwrk(i,m2sstm)=dwrk(i,m2sstm)+ &
683 & diau3wrk(i,j,k,m3sstm)*hzk(i,k)
684# endif
685 dwrk(i,m2wrol)=dwrk(i,m2wrol)+ &
686 & diau3wrk(i,j,k,m3wrol)*hzk(i,k)
687 dwrk(i,m2wbrk)=dwrk(i,m2wbrk)+ &
688 & diau3wrk(i,j,k,m3wbrk)*hzk(i,k)
689# endif
690# endif
691 END DO
692 END DO
693 DO i=istru,iend
694 cff1=1.0_r8/(cf(i,0)*on_u(i,j))
695 dc(i,0)=(dc(i,0)*on_u(i,j)-du_avg1(i,j))*cff1 ! recursive
696# ifdef WEC
697 cff2=1.0_r8/cf(i,0)
698 dcs(i,0)=dcs(i,0)*cff2-ubar_stokes(i,j) ! recursive
699# endif
700# ifdef DIAGNOSTICS_UV
701 DO idiag=1,m2pgrd
702 dwrk(i,idiag)=(dwrk(i,idiag)*on_u(i,j)- &
703 & diau2wrk(i,j,idiag))*cff1
704 END DO
705 dwrk(i,m2bstr)=(dwrk(i,m2bstr)*on_u(i,j)- &
706 & diau2wrk(i,j,m2bstr)-diau2wrk(i,j,m2sstr))* &
707 & cff1
708# endif
709 END DO
710!
711! Couple and update new solution.
712!
713 DO k=1,n(ng)
714 DO i=istru,iend
715 u(i,j,k,nnew)=u(i,j,k,nnew)-dc(i,0)
716# ifdef MASKING
717 u(i,j,k,nnew)=u(i,j,k,nnew)*umask(i,j)
718# endif
719# ifdef WET_DRY
720 u(i,j,k,nnew)=u(i,j,k,nnew)*umask_wet(i,j)
721 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask_wet(i,j)
722# endif
723# ifdef WEC
724 u_stokes(i,j,k)=u_stokes(i,j,k)-dcs(i,0)
725# ifdef MASKING
726 u_stokes(i,j,k)=u_stokes(i,j,k)*umask(i,j)
727# endif
728# ifdef WET_DRY
729 u_stokes(i,j,k)=u_stokes(i,j,k)*umask_wet(i,j)
730# endif
731# endif
732# ifdef DIAGNOSTICS_UV
733 diau3wrk(i,j,k,m3pgrd)=diau3wrk(i,j,k,m3pgrd)- &
734 & dwrk(i,m2pgrd)
735 diau3wrk(i,j,k,m3vvis)=diau3wrk(i,j,k,m3vvis)- &
736 & dwrk(i,m2bstr)
737# ifdef UV_COR
738 diau3wrk(i,j,k,m3fcor)=diau3wrk(i,j,k,m3fcor)- &
739 & dwrk(i,m2fcor)
740# endif
741# if defined UV_VIS2 || defined UV_VIS4
742 diau3wrk(i,j,k,m3xvis)=diau3wrk(i,j,k,m3xvis)- &
743 & dwrk(i,m2xvis)
744 diau3wrk(i,j,k,m3yvis)=diau3wrk(i,j,k,m3yvis)- &
745 & dwrk(i,m2yvis)
746 diau3wrk(i,j,k,m3hvis)=diau3wrk(i,j,k,m3hvis)- &
747 & dwrk(i,m2hvis)
748# endif
749# ifdef UV_ADV
750 diau3wrk(i,j,k,m3xadv)=diau3wrk(i,j,k,m3xadv)- &
751 & dwrk(i,m2xadv)
752 diau3wrk(i,j,k,m3yadv)=diau3wrk(i,j,k,m3yadv)- &
753 & dwrk(i,m2yadv)
754 diau3wrk(i,j,k,m3hadv)=diau3wrk(i,j,k,m3hadv)- &
755 & dwrk(i,m2hadv)
756# endif
757# ifdef WEC_VF
758 diau3wrk(i,j,k,m3hjvf)=diau3wrk(i,j,k,m3hjvf)- &
759 & dwrk(i,m2hjvf)
760 diau3wrk(i,j,k,m3kvrf)=diau3wrk(i,j,k,m3kvrf)- &
761 & dwrk(i,m2kvrf)
762# ifdef UV_COR
763 diau3wrk(i,j,k,m3fsco)=diau3wrk(i,j,k,m3fsco)- &
764 & dwrk(i,m2fsco)
765# endif
766# ifdef BOTTOM_STREAMING
767 diau3wrk(i,j,k,m3bstm)=diau3wrk(i,j,k,m3bstm)- &
768 & dwrk(i,m2bstm)
769# endif
770# ifdef SURFACE_STREAMING
771 diau3wrk(i,j,k,m3sstm)=diau3wrk(i,j,k,m3sstm)- &
772 & dwrk(i,m2sstm)
773# endif
774 diau3wrk(i,j,k,m3wrol)=diau3wrk(i,j,k,m3wrol)- &
775 & dwrk(i,m2wrol)
776 diau3wrk(i,j,k,m3wbrk)=diau3wrk(i,j,k,m3wbrk)- &
777 & dwrk(i,m2wbrk)
778# endif
779# endif
780 END DO
781 END DO
782
783# if defined DIAGNOSTICS_UV && defined MASKING
784 DO k=1,n(ng)
785 DO i=istru,iend
786 DO idiag=1,ndm3d
787 diau3wrk(i,j,k,idiag)=diau3wrk(i,j,k,idiag)*umask(i,j)
788 END DO
789 END DO
790 END DO
791# endif
792!
793!-----------------------------------------------------------------------
794! Time step momentum equation in the ETA-direction.
795!-----------------------------------------------------------------------
796!
797 IF (j.ge.jstrv) THEN
798 DO i=istr,iend
799 ak(i,0)=0.5_r8*(akv(i,j-1,0)+ &
800 & akv(i,j ,0))
801 DO k=1,n(ng)
802 ak(i,k)=0.5_r8*(akv(i,j-1,k)+ &
803 & akv(i,j ,k))
804 hzk(i,k)=0.5_r8*(hz(i,j-1,k)+ &
805 & hz(i,j ,k))
806# if defined SPLINES_VVISC || defined DIAGNOSTICS_UV
807 ohz(i,k)=1.0_r8/hzk(i,k)
808# endif
809 END DO
810 END DO
811!
812! Time step right-hand-side terms.
813!
814 IF (iic(ng).eq.ntfirst(ng)) THEN
815 cff=0.25_r8*dt(ng)
816 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
817 cff=0.25_r8*dt(ng)*3.0_r8/2.0_r8
818 ELSE
819 cff=0.25_r8*dt(ng)*23.0_r8/12.0_r8
820 END IF
821 DO i=istr,iend
822 dc(i,0)=cff*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
823 END DO
824 DO k=1,n(ng)
825 DO i=istr,iend
826 v(i,j,k,nnew)=v(i,j,k,nnew)+dc(i,0)*rv(i,j,k,nrhs)
827# ifdef SPLINES_VVISC
828 v(i,j,k,nnew)=v(i,j,k,nnew)*ohz(i,k)
829# endif
830# ifdef DIAGNOSTICS_UV
831 DO idiag=1,m3pgrd
832 diav3wrk(i,j,k,idiag)=(diav3wrk(i,j,k,idiag)+ &
833 & dc(i,0)* &
834 & diarv(i,j,k,nrhs,idiag))* &
835 & ohz(i,k)
836 END DO
837# if defined UV_VIS2 || defined UV_VIS4
838 diav3wrk(i,j,k,m3xvis)=diav3wrk(i,j,k,m3xvis)*ohz(i,k)
839 diav3wrk(i,j,k,m3yvis)=diav3wrk(i,j,k,m3yvis)*ohz(i,k)
840 diav3wrk(i,j,k,m3hvis)=diav3wrk(i,j,k,m3hvis)*ohz(i,k)
841# endif
842 diav3wrk(i,j,k,m3vvis)=diav3wrk(i,j,k,m3vvis)*ohz(i,k)
843# ifdef BODYFORCE
844 diav3wrk(i,j,k,m3vvis)=diav3wrk(i,j,k,m3vvis)+ &
845 & dc(i,0)*diarv(i,j,k,nrhs,m3vvis)* &
846 & ohz(i,k)
847# endif
848 diav3wrk(i,j,k,m3rate)=diav3wrk(i,j,k,m3rate)*ohz(i,k)
849# endif
850 END DO
851 END DO
852
853# ifdef SPLINES_VVISC
854!
855! Use conservative, parabolic spline reconstruction of vertical
856! viscosity derivatives. Then, time step vertical viscosity term
857! implicitly by solving a tridiagonal system.
858!
859 cff1=1.0_r8/6.0_r8
860 DO k=1,n(ng)-1
861 DO i=istr,iend
862 fc(i,k)=cff1*hzk(i,k )-dt(ng)*ak(i,k-1)*ohz(i,k )
863 cf(i,k)=cff1*hzk(i,k+1)-dt(ng)*ak(i,k+1)*ohz(i,k+1)
864 END DO
865 END DO
866 DO i=istr,iend
867 cf(i,0)=0.0_r8
868 dc(i,0)=0.0_r8
869 END DO
870!
871! LU decomposition and forward substitution.
872!
873 cff1=1.0_r8/3.0_r8
874 DO k=1,n(ng)-1
875 DO i=istr,iend
876 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
877 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
878 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
879 cf(i,k)=cff*cf(i,k)
880 dc(i,k)=cff*(v(i,j,k+1,nnew)-v(i,j,k,nnew)- &
881 & fc(i,k)*dc(i,k-1))
882 END DO
883 END DO
884!
885! Backward substitution.
886!
887 DO i=istr,iend
888 dc(i,n(ng))=0.0_r8
889 END DO
890 DO k=n(ng)-1,1,-1
891 DO i=istr,iend
892 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
893 END DO
894 END DO
895!
896 DO k=1,n(ng)
897 DO i=istr,iend
898 dc(i,k)=dc(i,k)*ak(i,k)
899 cff=dt(ng)*ohz(i,k)*(dc(i,k)-dc(i,k-1))
900 v(i,j,k,nnew)=v(i,j,k,nnew)+cff
901# ifdef DIAGNOSTICS_UV
902 diav3wrk(i,j,k,m3vvis)=diav3wrk(i,j,k,m3vvis)+cff
903# endif
904 END DO
905 END DO
906# else
907!
908! Compute off-diagonal coefficients [lambda*dt*Akv/Hz] for the
909! implicit vertical viscosity term at horizontal V-points and
910! vertical W-points.
911!
912 cff=-lambda*dt(ng)/0.5_r8
913 DO k=1,n(ng)-1
914 DO i=istr,iend
915 cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
916 & z_r(i,j,k )-z_r(i,j-1,k ))
917 fc(i,k)=cff*cff1*ak(i,k)
918 END DO
919 END DO
920 DO i=istr,iend
921 fc(i,0)=0.0_r8
922 fc(i,n(ng))=0.0_r8
923 END DO
924!
925! Solve the tridiagonal system.
926!
927 DO k=1,n(ng)
928 DO i=istr,iend
929 dc(i,k)=v(i,j,k,nnew)
930 bc(i,k)=hzk(i,k)-fc(i,k)-fc(i,k-1)
931 END DO
932 END DO
933 DO i=istr,iend
934 cff=1.0_r8/bc(i,1)
935 cf(i,1)=cff*fc(i,1)
936 dc(i,1)=cff*dc(i,1)
937 END DO
938 DO k=2,n(ng)-1
939 DO i=istr,iend
940 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
941 cf(i,k)=cff*fc(i,k)
942 dc(i,k)=cff*(dc(i,k)-fc(i,k-1)*dc(i,k-1))
943 END DO
944 END DO
945!
946! Compute new solution by back substitution.
947!
948 DO i=istr,iend
949# ifdef DIAGNOSTICS_UV
950 wrk(i,n(ng))=v(i,j,n(ng),nnew)*ohz(i,n(ng))
951# endif
952 dc(i,n(ng))=(dc(i,n(ng))-fc(i,n(ng)-1)*dc(i,n(ng)-1))/ &
953 & (bc(i,n(ng))-fc(i,n(ng)-1)*cf(i,n(ng)-1))
954 v(i,j,n(ng),nnew)=dc(i,n(ng))
955# ifdef DIAGNOSTICS_UV
956 diav3wrk(i,j,n(ng),m3vvis)=diav3wrk(i,j,n(ng),m3vvis)+ &
957 & v(i,j,n(ng),nnew)-wrk(i,n(ng))
958# endif
959 END DO
960 DO k=n(ng)-1,1,-1
961 DO i=istr,iend
962# ifdef DIAGNOSTICS_UV
963 wrk(i,k)=v(i,j,k,nnew)*ohz(i,k)
964# endif
965 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
966 v(i,j,k,nnew)=dc(i,k)
967# ifdef DIAGNOSTICS_UV
968 diav3wrk(i,j,k,m3vvis)=diav3wrk(i,j,k,m3vvis)+ &
969 & v(i,j,k,nnew)-wrk(i,k)
970# endif
971 END DO
972 END DO
973# endif
974# ifdef OMEGA_IMPLICIT
975!
976! Adaptive, Courant-number based implicit vertical advection
977! contribution for v-momentum.
978!
979 DO i=istr,iend
980 wk(i,0)=0.5_r8*(wi(i,j-1,0)+ &
981 & wi(i,j ,0))
982 DO k=1,n(ng)
983 wk(i,k)=0.5_r8*(wi(i,j-1,k)+ &
984 & wi(i,j ,k))
985 hzk(i,k)=0.5_r8*(hz(i,j-1,k)+ &
986 & hz(i,j ,k))
987 END DO
988 END DO
989!
990! Compute off-diagonal coefficients [dt*Wi*pm*pn] for the
991! implicit vertical viscosity term at horizontal V-points and
992! vertical W-points.
993!
994 cff=dt(ng)
995 DO k=1,n(ng)-1
996 DO i=istr,iend
997 cff1=cff/(on_v(i,j)*om_v(i,j))
998 fcmax(i,k)=max(wk(i,k),0.0_r8)*cff1
999 fcmin(i,k)=min(wk(i,k),0.0_r8)*cff1
1000 END DO
1001 END DO
1002 DO i=istr,iend
1003 fcmax(i,0)=0.0_r8
1004 fcmin(i,0)=0.0_r8
1005 fcmax(i,n(ng))=0.0_r8
1006 fcmin(i,n(ng))=0.0_r8
1007 END DO
1008!
1009! Solve the tridiagonal system.
1010!
1011 DO k=1,n(ng)
1012 DO i=istr,iend
1013 bc(i,k)=hzk(i,k)+fcmax(i,k)-fcmin(i,k-1)
1014 dc(i,k)=v(i,j,k,nnew)*hzk(i,k)
1015 END DO
1016 END DO
1017 DO i=istr,iend
1018 cff=1.0_r8/bc(i,1)
1019 cf(i,1)=cff*fcmin(i,1)
1020 dc(i,1)=cff*dc(i,1)
1021 END DO
1022 DO k=2,n(ng)-1
1023 DO i=istr,iend
1024 cff=1.0_r8/(bc(i,k)+fcmax(i,k-1)*cf(i,k-1))
1025 cf(i,k)=cff*fcmin(i,k)
1026 dc(i,k)=cff*(dc(i,k)+fcmax(i,k-1)*dc(i,k-1))
1027 END DO
1028 END DO
1029!
1030! Compute new solution by back substitution.
1031!
1032 DO i=istr,iend
1033# ifdef DIAGNOSTICS_UV
1034 cff1=v(i,j,n(ng),nnew)
1035# endif
1036 cff=1.0_r8/(bc(i,n(ng))+fcmax(i,n(ng)-1)*cf(i,n(ng)-1))
1037 dc(i,n(ng))=cff*(dc(i,n(ng))+ &
1038 & fcmax(i,n(ng)-1)*dc(i,n(ng)-1))
1039 v(i,j,n(ng),nnew)=dc(i,n(ng))
1040# ifdef DIAGNOSTICS_UV
1041 diarv(i,j,n(ng),nrhs,m3vadv)=diarv(i,j,n(ng),nrhs,m3vadv)+ &
1042 & v(i,j,n(ng),nnew)-cff1
1043# endif
1044 END DO
1045!
1046 DO k=n(ng)-1,1,-1
1047 DO i=istr,iend
1048# ifdef DIAGNOSTICS_UV
1049 cff1=v(i,j,k,nnew)
1050# endif
1051 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1052 v(i,j,k,nnew)=dc(i,k)
1053# ifdef DIAGNOSTICS_UV
1054 diarv(i,j,k,nrhs,m3vadv)=diarv(i,j,k,nrhs,m3vadv)+ &
1055 & v(i,j,k,nnew)-cff1
1056# endif
1057 END DO
1058 END DO
1059# endif
1060!
1061! Replace INTERIOR POINTS incorrect vertical mean with more accurate
1062! barotropic component, vbar=DV_avg1/(D*om_v). Recall that, D=CF(:,0).
1063!
1064 DO i=istr,iend
1065 cf(i,0)=hzk(i,1)
1066 dc(i,0)=v(i,j,1,nnew)*hzk(i,1)
1067# ifdef WEC
1068 dcs(i,0)=v_stokes(i,j,1)*hzk(i,1)
1069# endif
1070# ifdef DIAGNOSTICS_UV
1071 dwrk(i,m2pgrd)=diav3wrk(i,j,1,m3pgrd)*hzk(i,1)
1072 dwrk(i,m2bstr)=diav3wrk(i,j,1,m3vvis)*hzk(i,1)
1073# ifdef UV_COR
1074 dwrk(i,m2fcor)=diav3wrk(i,j,1,m3fcor)*hzk(i,1)
1075# endif
1076# if defined UV_VIS2 || defined UV_VIS4
1077 dwrk(i,m2xvis)=diav3wrk(i,j,1,m3xvis)*hzk(i,1)
1078 dwrk(i,m2yvis)=diav3wrk(i,j,1,m3yvis)*hzk(i,1)
1079 dwrk(i,m2hvis)=diav3wrk(i,j,1,m3hvis)*hzk(i,1)
1080# endif
1081# ifdef UV_ADV
1082 dwrk(i,m2xadv)=diav3wrk(i,j,1,m3xadv)*hzk(i,1)
1083 dwrk(i,m2yadv)=diav3wrk(i,j,1,m3yadv)*hzk(i,1)
1084 dwrk(i,m2hadv)=diav3wrk(i,j,1,m3hadv)*hzk(i,1)
1085# endif
1086# ifdef WEC_VF
1087 dwrk(i,m2hjvf)=diav3wrk(i,j,1,m3hjvf)*hzk(i,1)
1088 dwrk(i,m2kvrf)=diav3wrk(i,j,1,m3kvrf)*hzk(i,1)
1089# ifdef UV_COR
1090 dwrk(i,m2fsco)=diav3wrk(i,j,1,m3fsco)*hzk(i,1)
1091# endif
1092# ifdef BOTTOM_STREAMING
1093 dwrk(i,m2bstm)=diav3wrk(i,j,1,m3bstm)*hzk(i,1)
1094# endif
1095# ifdef SURFACE_STREAMING
1096 dwrk(i,m2sstm)=diav3wrk(i,j,1,m3sstm)*hzk(i,1)
1097# endif
1098 dwrk(i,m2wrol)=diav3wrk(i,j,1,m3wrol)*hzk(i,1)
1099 dwrk(i,m2wbrk)=diav3wrk(i,j,1,m3wbrk)*hzk(i,1)
1100# endif
1101# endif
1102 END DO
1103 DO k=2,n(ng)
1104 DO i=istr,iend
1105 cf(i,0)=cf(i,0)+hzk(i,k)
1106 dc(i,0)=dc(i,0)+v(i,j,k,nnew)*hzk(i,k)
1107# ifdef WEC
1108 dcs(i,0)=dcs(i,0)+v_stokes(i,j,k)*hzk(i,k)
1109# endif
1110# ifdef DIAGNOSTICS_UV
1111 dwrk(i,m2pgrd)=dwrk(i,m2pgrd)+ &
1112 & diav3wrk(i,j,k,m3pgrd)*hzk(i,k)
1113 dwrk(i,m2bstr)=dwrk(i,m2bstr)+ &
1114 & diav3wrk(i,j,k,m3vvis)*hzk(i,k)
1115# ifdef UV_COR
1116 dwrk(i,m2fcor)=dwrk(i,m2fcor)+ &
1117 & diav3wrk(i,j,k,m3fcor)*hzk(i,k)
1118# endif
1119# if defined UV_VIS2 || defined UV_VIS4
1120 dwrk(i,m2xvis)=dwrk(i,m2xvis)+ &
1121 & diav3wrk(i,j,k,m3xvis)*hzk(i,k)
1122 dwrk(i,m2yvis)=dwrk(i,m2yvis)+ &
1123 & diav3wrk(i,j,k,m3yvis)*hzk(i,k)
1124 dwrk(i,m2hvis)=dwrk(i,m2hvis)+ &
1125 & diav3wrk(i,j,k,m3hvis)*hzk(i,k)
1126# endif
1127# ifdef UV_ADV
1128 dwrk(i,m2xadv)=dwrk(i,m2xadv)+ &
1129 & diav3wrk(i,j,k,m3xadv)*hzk(i,k)
1130 dwrk(i,m2yadv)=dwrk(i,m2yadv)+ &
1131 & diav3wrk(i,j,k,m3yadv)*hzk(i,k)
1132 dwrk(i,m2hadv)=dwrk(i,m2hadv)+ &
1133 & diav3wrk(i,j,k,m3hadv)*hzk(i,k)
1134# endif
1135# ifdef WEC_VF
1136 dwrk(i,m2hjvf)=dwrk(i,m2hjvf)+ &
1137 & diav3wrk(i,j,k,m3hjvf)*hzk(i,k)
1138 dwrk(i,m2kvrf)=dwrk(i,m2kvrf)+ &
1139 & diav3wrk(i,j,k,m3kvrf)*hzk(i,k)
1140# ifdef UV_COR
1141 dwrk(i,m2fsco)=dwrk(i,m2fsco)+ &
1142 & diav3wrk(i,j,k,m3fsco)*hzk(i,k)
1143# endif
1144# ifdef BOTTOM_STREAMING
1145 dwrk(i,m2bstm)=dwrk(i,m2bstm)+ &
1146 & diav3wrk(i,j,k,m3bstm)*hzk(i,k)
1147# endif
1148# ifdef SURFACE_STREAMING
1149 dwrk(i,m2sstm)=dwrk(i,m2sstm)+ &
1150 & diav3wrk(i,j,k,m3sstm)*hzk(i,k)
1151# endif
1152 dwrk(i,m2wrol)=dwrk(i,m2wrol)+ &
1153 & diav3wrk(i,j,k,m3wrol)*hzk(i,k)
1154 dwrk(i,m2wbrk)=dwrk(i,m2wbrk)+ &
1155 & diav3wrk(i,j,k,m3wbrk)*hzk(i,k)
1156# endif
1157# endif
1158 END DO
1159 END DO
1160 DO i=istr,iend
1161 cff1=1.0_r8/(cf(i,0)*om_v(i,j))
1162 dc(i,0)=(dc(i,0)*om_v(i,j)-dv_avg1(i,j))*cff1 ! recursive
1163# ifdef WEC
1164 cff2=1.0_r8/cf(i,0)
1165 dcs(i,0)=dcs(i,0)*cff2-vbar_stokes(i,j) ! recursive
1166# endif
1167# ifdef DIAGNOSTICS_UV
1168 DO idiag=1,m2pgrd
1169 dwrk(i,idiag)=(dwrk(i,idiag)*om_v(i,j)- &
1170 & diav2wrk(i,j,idiag))*cff1
1171 END DO
1172 dwrk(i,m2bstr)=(dwrk(i,m2bstr)*om_v(i,j)- &
1173 & diav2wrk(i,j,m2bstr)-diav2wrk(i,j,m2sstr))* &
1174 & cff1
1175# endif
1176 END DO
1177!
1178! Couple and update new solution.
1179!
1180 DO k=1,n(ng)
1181 DO i=istr,iend
1182 v(i,j,k,nnew)=v(i,j,k,nnew)-dc(i,0)
1183# ifdef MASKING
1184 v(i,j,k,nnew)=v(i,j,k,nnew)*vmask(i,j)
1185# endif
1186# ifdef WET_DRY
1187 v(i,j,k,nnew)=v(i,j,k,nnew)*vmask_wet(i,j)
1188 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask_wet(i,j)
1189# endif
1190# ifdef WEC
1191 v_stokes(i,j,k)=v_stokes(i,j,k)-dcs(i,0)
1192# ifdef MASKING
1193 v_stokes(i,j,k)=v_stokes(i,j,k)*vmask(i,j)
1194# endif
1195# ifdef WET_DRY
1196 v_stokes(i,j,k)=v_stokes(i,j,k)*vmask_wet(i,j)
1197# endif
1198# endif
1199# ifdef DIAGNOSTICS_UV
1200 diav3wrk(i,j,k,m3pgrd)=diav3wrk(i,j,k,m3pgrd)- &
1201 & dwrk(i,m2pgrd)
1202 diav3wrk(i,j,k,m3vvis)=diav3wrk(i,j,k,m3vvis)- &
1203 & dwrk(i,m2bstr)
1204# ifdef UV_COR
1205 diav3wrk(i,j,k,m3fcor)=diav3wrk(i,j,k,m3fcor)- &
1206 & dwrk(i,m2fcor)
1207# endif
1208# if defined UV_VIS2 || defined UV_VIS4
1209 diav3wrk(i,j,k,m3xvis)=diav3wrk(i,j,k,m3xvis)- &
1210 & dwrk(i,m2xvis)
1211 diav3wrk(i,j,k,m3yvis)=diav3wrk(i,j,k,m3yvis)- &
1212 & dwrk(i,m2yvis)
1213 diav3wrk(i,j,k,m3hvis)=diav3wrk(i,j,k,m3hvis)- &
1214 & dwrk(i,m2hvis)
1215# endif
1216# ifdef UV_ADV
1217 diav3wrk(i,j,k,m3xadv)=diav3wrk(i,j,k,m3xadv)- &
1218 & dwrk(i,m2xadv)
1219 diav3wrk(i,j,k,m3yadv)=diav3wrk(i,j,k,m3yadv)- &
1220 & dwrk(i,m2yadv)
1221 diav3wrk(i,j,k,m3hadv)=diav3wrk(i,j,k,m3hadv)- &
1222 & dwrk(i,m2hadv)
1223# endif
1224# ifdef WEC_VF
1225 diav3wrk(i,j,k,m3hjvf)=diav3wrk(i,j,k,m3hjvf)- &
1226 & dwrk(i,m2hjvf)
1227 diav3wrk(i,j,k,m3kvrf)=diav3wrk(i,j,k,m3kvrf)- &
1228 & dwrk(i,m2kvrf)
1229# ifdef UV_COR
1230 diav3wrk(i,j,k,m3fsco)=diav3wrk(i,j,k,m3fsco)- &
1231 & dwrk(i,m2fsco)
1232# endif
1233# ifdef BOTTOM_STREAMING
1234 diav3wrk(i,j,k,m3bstm)=diav3wrk(i,j,k,m3bstm)- &
1235 & dwrk(i,m2bstm)
1236# endif
1237# ifdef SURFACE_STREAMING
1238 diav3wrk(i,j,k,m3sstm)=diav3wrk(i,j,k,m3sstm)- &
1239 & dwrk(i,m2sstm)
1240# endif
1241 diav3wrk(i,j,k,m3wrol)=diav3wrk(i,j,k,m3wrol)- &
1242 & dwrk(i,m2wrol)
1243 diav3wrk(i,j,k,m3wbrk)=diav3wrk(i,j,k,m3wbrk)- &
1244 & dwrk(i,m2wbrk)
1245# endif
1246# endif
1247 END DO
1248 END DO
1249
1250# if defined DIAGNOSTICS_UV && defined MASKING
1251 DO k=1,n(ng)
1252 DO i=istr,iend
1253 DO idiag=1,ndm3d
1254 diav3wrk(i,j,k,idiag)=diav3wrk(i,j,k,idiag)*vmask(i,j)
1255 END DO
1256 END DO
1257 END DO
1258# endif
1259 END IF
1260 END DO
1261!
1262!-----------------------------------------------------------------------
1263! Set lateral boundary conditions.
1264!-----------------------------------------------------------------------
1265!
1266 CALL u3dbc_tile (ng, tile, &
1267 & lbi, ubi, lbj, ubj, n(ng), &
1268 & imins, imaxs, jmins, jmaxs, &
1269 & nstp, nnew, &
1270 & u)
1271 CALL v3dbc_tile (ng, tile, &
1272 & lbi, ubi, lbj, ubj, n(ng), &
1273 & imins, imaxs, jmins, jmaxs, &
1274 & nstp, nnew, &
1275 & v)
1276!
1277!-----------------------------------------------------------------------
1278! Apply momentum transport point sources (like river runoff), if any.
1279!-----------------------------------------------------------------------
1280!
1281 IF (luvsrc(ng)) THEN
1282 DO is=1,nsrc(ng)
1283 i=sources(ng)%Isrc(is)
1284 j=sources(ng)%Jsrc(is)
1285 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1286 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1287 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
1288 DO k=1,n(ng)
1289 cff1=1.0_r8/(on_u(i,j)* &
1290 & 0.5_r8*(z_w(i-1,j,k)-z_w(i-1,j,k-1)+ &
1291 & z_w(i ,j,k)-z_w(i ,j,k-1)))
1292 u(i,j,k,nnew)=sources(ng)%Qsrc(is,k)*cff1
1293 END DO
1294 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
1295 DO k=1,n(ng)
1296 cff1=1.0_r8/(om_v(i,j)* &
1297 & 0.5_r8*(z_w(i,j-1,k)-z_w(i,j-1,k-1)+ &
1298 & z_w(i,j ,k)-z_w(i,j ,k-1)))
1299 v(i,j,k,nnew)=sources(ng)%Qsrc(is,k)*cff1
1300 END DO
1301 END IF
1302 END IF
1303 END DO
1304 END IF
1305!
1306!-----------------------------------------------------------------------
1307! Couple 2D and 3D momentum equations.
1308!-----------------------------------------------------------------------
1309!
1310! Couple velocity component in the XI-direction.
1311!
1312 DO j=jstrt,jendt
1313 DO i=istrp,iendt
1314 dc(i,0)=0.0_r8
1315 cf(i,0)=0.0_r8
1316# ifdef WEC
1317 cfs(i,0)=0.0_r8
1318# endif
1319 fc(i,0)=0.0_r8
1320 END DO
1321!
1322! Compute thicknesses of U-boxes DC(i,1:N), total depth of the water
1323! column DC(i,0), and incorrect vertical mean CF(i,0). Notice that
1324! barotropic component is replaced with its fast-time averaged
1325! values.
1326!
1327 DO k=1,n(ng)
1328 DO i=istrp,iendt
1329 cff=0.5_r8*on_u(i,j)
1330 dc(i,k)=cff*(hz(i,j,k)+hz(i-1,j,k))
1331 dc(i,0)=dc(i,0)+dc(i,k)
1332 cf(i,0)=cf(i,0)+ &
1333 & dc(i,k)*u(i,j,k,nnew)
1334# ifdef WEC
1335 cfs(i,0)=cfs(i,0)+ &
1336 & dc(i,k)*u_stokes(i,j,k)
1337# endif
1338 END DO
1339 END DO
1340 DO i=istrp,iendt
1341 cff1=dc(i,0) ! intermediate
1342# ifdef WEC
1343!! cff2=on_u(i,j)*DC(i,0)*ubar_stokes(i,j)
1344 cff2=dc(i,0)*ubar_stokes(i,j)
1345# endif
1346 dc(i,0)=1.0_r8/dc(i,0) ! recursive
1347 cf(i,0)=dc(i,0)*(cf(i,0)-du_avg1(i,j)) ! recursive
1348# ifdef WEC
1349 cfs(i,0)=dc(i,0)*(cfs(i,0)-cff2) ! recursive
1350# endif
1351# if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
1352 ubar(i,j,knew)=dc(i,0)*du_avg1(i,j)
1353# ifdef WET_DRY
1354 ubar(i,j,knew)=ubar(i,j,knew)*umask_wet(i,j)
1355# endif
1356# else
1357 ubar(i,j,1)=dc(i,0)*du_avg1(i,j)
1358# ifdef WET_DRY
1359 ubar(i,j,1)=ubar(i,j,1)*umask_wet(i,j)
1360# endif
1361 ubar(i,j,2)=ubar(i,j,1)
1362# endif
1363# ifdef DIAGNOSTICS_UV
1364 diau2wrk(i,j,m2rate)=ubar(i,j,1)-diau2int(i,j,m2rate)*dc(i,0)
1365 diau2int(i,j,m2rate)=ubar(i,j,1)*cff1
1366# endif
1367 END DO
1368# ifdef DIAGNOSTICS_UV
1369!
1370! Convert the units of the fast-time integrated change in mass flux
1371! diagnostic values to change in velocity (m s-1).
1372!
1373 DO idiag=1,ndm2d-1
1374 DO i=istrp,iendt
1375 diau2wrk(i,j,idiag)=dc(i,0)*diau2wrk(i,j,idiag)
1376# ifdef MASKING
1377 diau2wrk(i,j,idiag)=diau2wrk(i,j,idiag)*umask(i,j)
1378# endif
1379 END DO
1380 END DO
1381# endif
1382!
1383! Replace only BOUNDARY POINTS incorrect vertical mean with more
1384! accurate barotropic component, ubar=DU_avg1/(D*on_u). Recall that,
1385! D=CF(:,0).
1386!
1387! NOTE: Only the BOUNDARY POINTS need to be replaced. Avoid redundant
1388! update in the interior again for computational purposes which
1389! will not affect the nonlinear code. However, the adjoint
1390! code is wrong because the interior solution is corrected
1391! twice. The replacement is avoided if the boundary edge is
1392! periodic. The J-loop is pipelined, so we need to do a special
1393! test for the southern and northern domain edges.
1394!
1395 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1396 IF (domain(ng)%Western_Edge(tile)) THEN
1397 DO k=1,n(ng)
1398 u(istr,j,k,nnew)=u(istr,j,k,nnew)-cf(istr,0)
1399# ifdef MASKING
1400 u(istr,j,k,nnew)=u(istr,j,k,nnew)* &
1401 & umask(istr,j)
1402# endif
1403# ifdef WET_DRY
1404 u(istr,j,k,nnew)=u(istr,j,k,nnew)* &
1405 & umask_wet(istr,j)
1406# endif
1407# ifdef WEC
1408 u_stokes(istr,j,k)=u_stokes(istr,j,k)-cfs(istr,0)
1409# ifdef MASKING
1410 u_stokes(istr,j,k)=u_stokes(istr,j,k)* &
1411 & umask(istr,j)
1412# endif
1413# ifdef WET_DRY
1414 u_stokes(istr,j,k)=u_stokes(istr,j,k)* &
1415 & umask_wet(istr,j)
1416# endif
1417# endif
1418 END DO
1419 END IF
1420 END IF
1421!
1422 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1423 IF (domain(ng)%Eastern_Edge(tile)) THEN
1424 DO k=1,n(ng)
1425 u(iend+1,j,k,nnew)=u(iend+1,j,k,nnew)-cf(iend+1,0)
1426# ifdef MASKING
1427 u(iend+1,j,k,nnew)=u(iend+1,j,k,nnew)* &
1428 & umask(iend+1,j)
1429# endif
1430# ifdef WET_DRY
1431 u(iend+1,j,k,nnew)=u(iend+1,j,k,nnew)* &
1432 & umask_wet(iend+1,j)
1433# endif
1434# ifdef WEC
1435 u_stokes(iend+1,j,k)=u_stokes(iend+1,j,k)-cfs(iend+1,0)
1436# ifdef MASKING
1437 u_stokes(iend+1,j,k)=u_stokes(iend+1,j,k)* &
1438 & umask(iend+1,j)
1439# endif
1440# ifdef WET_DRY
1441 u_stokes(iend+1,j,k)=u_stokes(iend+1,j,k)* &
1442 & umask_wet(iend+1,j)
1443# endif
1444# endif
1445 END DO
1446 END IF
1447 END IF
1448!
1449 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1450 IF (j.eq.0) THEN ! southern boundary
1451 DO k=1,n(ng) ! J-loop pipelined
1452 DO i=istru,iend
1453 u(i,j,k,nnew)=u(i,j,k,nnew)-cf(i,0)
1454# ifdef MASKING
1455 u(i,j,k,nnew)=u(i,j,k,nnew)* &
1456 & umask(i,j)
1457# endif
1458# ifdef WET_DRY
1459 u(i,j,k,nnew)=u(i,j,k,nnew)* &
1460 & umask_wet(i,j)
1461# endif
1462# ifdef WEC
1463 u_stokes(i,j,k)=u_stokes(i,j,k)-cfs(i,0)
1464# ifdef MASKING
1465 u_stokes(i,j,k)=u_stokes(i,j,k)* &
1466 & umask(i,j)
1467# endif
1468# ifdef WET_DRY
1469 u_stokes(i,j,k)=u_stokes(i,j,k)* &
1470 & umask_wet(i,j)
1471# endif
1472# endif
1473 END DO
1474 END DO
1475 END IF
1476 END IF
1477!
1478 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1479 IF (j.eq.mm(ng)+1) THEN ! northern boundary
1480 DO k=1,n(ng) ! J-loop pipelined
1481 DO i=istru,iend
1482 u(i,j,k,nnew)=u(i,j,k,nnew)-cf(i,0)
1483# ifdef MASKING
1484 u(i,j,k,nnew)=u(i,j,k,nnew)* &
1485 & umask(i,j)
1486# endif
1487# ifdef WET_DRY
1488 u(i,j,k,nnew)=u(i,j,k,nnew)* &
1489 & umask_wet(i,j)
1490# endif
1491# ifdef WEC
1492 u_stokes(i,j,k)=u_stokes(i,j,k)-cfs(i,0)
1493# ifdef MASKING
1494 u_stokes(i,j,k)=u_stokes(i,j,k)* &
1495 & umask(i,j)
1496# endif
1497# ifdef WET_DRY
1498 u_stokes(i,j,k)=u_stokes(i,j,k)* &
1499 & umask_wet(i,j)
1500# endif
1501# endif
1502 END DO
1503 END DO
1504 END IF
1505 END IF
1506!
1507! Compute correct mass flux, Hz*u/n.
1508!
1509 DO k=n(ng),1,-1
1510 DO i=istrp,iendt
1511 huon(i,j,k)=0.5_r8*(huon(i,j,k)+u(i,j,k,nnew)*dc(i,k))
1512# ifdef WEC
1513 huon(i,j,k)=huon(i,j,k)+0.5_r8*u_stokes(i,j,k)*dc(i,k)
1514# endif
1515 fc(i,0)=fc(i,0)+huon(i,j,k)
1516# ifdef DIAGNOSTICS_UV
1517 diau3wrk(i,j,k,m3rate)=u(i,j,k,nnew)-diau3wrk(i,j,k,m3rate)
1518# endif
1519 END DO
1520 END DO
1521 DO i=istrp,iendt
1522 fc(i,0)=dc(i,0)*(fc(i,0)-du_avg2(i,j)) ! recursive
1523 END DO
1524 DO k=1,n(ng)
1525 DO i=istrp,iendt
1526 huon(i,j,k)=huon(i,j,k)-dc(i,k)*fc(i,0)
1527 END DO
1528 END DO
1529!
1530! Couple velocity component in the ETA-direction.
1531!
1532 IF (j.ge.jstr) THEN
1533 DO i=istrt,iendt
1534 dc(i,0)=0.0_r8
1535 cf(i,0)=0.0_r8
1536# ifdef WEC
1537 cfs(i,0)=0.0_r8
1538# endif
1539 fc(i,0)=0.0_r8
1540 END DO
1541!
1542! Compute thicknesses of V-boxes DC(i,1:N), total depth of the water
1543! column DC(i,0), and incorrect vertical mean CF(i,0). Notice that
1544! barotropic component is replaced with its fast-time averaged
1545! values.
1546!
1547 DO k=1,n(ng)
1548 DO i=istrt,iendt
1549 cff=0.5_r8*om_v(i,j)
1550 dc(i,k)=cff*(hz(i,j,k)+hz(i,j-1,k))
1551 dc(i,0)=dc(i,0)+dc(i,k)
1552 cf(i,0)=cf(i,0)+ &
1553 & dc(i,k)*v(i,j,k,nnew)
1554# ifdef WEC
1555 cfs(i,0)=cfs(i,0)+ &
1556 & dc(i,k)*v_stokes(i,j,k)
1557# endif
1558 END DO
1559 END DO
1560 DO i=istrt,iendt
1561 cff1=dc(i,0) ! Intermediate
1562# ifdef WEC
1563!! cff2=om_v(i,j)*DC(i,0)*vbar_stokes(i,j)
1564 cff2=dc(i,0)*vbar_stokes(i,j)
1565# endif
1566 dc(i,0)=1.0_r8/dc(i,0) ! recursive
1567 cf(i,0)=dc(i,0)*(cf(i,0)-dv_avg1(i,j)) ! recursive
1568# ifdef WEC
1569 cfs(i,0)=dc(i,0)*(cfs(i,0)-cff2) ! recursive
1570# endif
1571# if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
1572 vbar(i,j,knew)=dc(i,0)*dv_avg1(i,j)
1573# ifdef WET_DRY
1574 vbar(i,j,knew)=vbar(i,j,knew)*vmask_wet(i,j)
1575# endif
1576# else
1577 vbar(i,j,1)=dc(i,0)*dv_avg1(i,j)
1578# ifdef WET_DRY
1579 vbar(i,j,1)=vbar(i,j,1)*vmask_wet(i,j)
1580# endif
1581 vbar(i,j,2)=vbar(i,j,1)
1582# endif
1583# ifdef DIAGNOSTICS_UV
1584 diav2wrk(i,j,m2rate)=vbar(i,j,1)- &
1585 & diav2int(i,j,m2rate)*dc(i,0)
1586 diav2int(i,j,m2rate)=vbar(i,j,1)*cff1
1587!! DiaV2wrk(i,j,M2rate)=vbar(i,j,1)-DiaV2int(i,j,M2rate)
1588!! DiaV2int(i,j,M2rate)=vbar(i,j,1)
1589# endif
1590 END DO
1591# ifdef DIAGNOSTICS_UV
1592!
1593! Convert the units of the fast-time integrated change in mass flux
1594! diagnostic values to change in velocity (m s-1).
1595!
1596 DO idiag=1,ndm2d-1
1597 DO i=istrt,iendt
1598 diav2wrk(i,j,idiag)=dc(i,0)*diav2wrk(i,j,idiag)
1599# ifdef MASKING
1600 diav2wrk(i,j,idiag)=diav2wrk(i,j,idiag)*vmask(i,j)
1601# endif
1602 END DO
1603 END DO
1604# endif
1605!
1606! Replace only BOUNDARY POINTS incorrect vertical mean with more
1607! accurate barotropic component, vbar=DV_avg1/(D*om_v). Recall that,
1608! D=CF(:,0).
1609!
1610! NOTE: Only the BOUNDARY POINTS need to be replaced. Avoid redundant
1611! update in the interior again for computational purposes which
1612! will not affect the nonlinear code. However, the adjoint
1613! code is wrong because the interior solution is corrected
1614! twice. The replacement is avoided if the boundary edge is
1615! periodic. The J-loop is pipelined, so we need to do a special
1616! test for the southern and northern domain edges.
1617!
1618 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1619 IF (domain(ng)%Western_Edge(tile)) THEN
1620 DO k=1,n(ng)
1621 v(istr-1,j,k,nnew)=v(istr-1,j,k,nnew)-cf(istr-1,0)
1622# ifdef MASKING
1623 v(istr-1,j,k,nnew)=v(istr-1,j,k,nnew)* &
1624 & vmask(istr-1,j)
1625# endif
1626# ifdef WET_DRY
1627 v(istr-1,j,k,nnew)=v(istr-1,j,k,nnew)* &
1628 & vmask_wet(istr-1,j)
1629# endif
1630# ifdef WEC
1631 v_stokes(istr-1,j,k)=v_stokes(istr-1,j,k)- &
1632 & cfs(istr-1,0)
1633# ifdef MASKING
1634 v_stokes(istr-1,j,k)=v_stokes(istr-1,j,k)* &
1635 & vmask(istr-1,j)
1636# endif
1637# ifdef WET_DRY
1638 v_stokes(istr-1,j,k)=v_stokes(istr-1,j,k)* &
1639 & vmask_wet(istr-1,j)
1640# endif
1641# endif
1642 END DO
1643 END IF
1644 END IF
1645!
1646 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1647 IF (domain(ng)%Eastern_Edge(tile)) THEN
1648 DO k=1,n(ng)
1649 v(iend+1,j,k,nnew)=v(iend+1,j,k,nnew)-cf(iend+1,0)
1650# ifdef MASKING
1651 v(iend+1,j,k,nnew)=v(iend+1,j,k,nnew)* &
1652 & vmask(iend+1,j)
1653# endif
1654# ifdef WET_DRY
1655 v(iend+1,j,k,nnew)=v(iend+1,j,k,nnew)* &
1656 & vmask_wet(iend+1,j)
1657# endif
1658# ifdef WEC
1659 v_stokes(iend+1,j,k)=v_stokes(iend+1,j,k)- &
1660 & cfs(iend+1,0)
1661# ifdef MASKING
1662 v_stokes(iend+1,j,k)=v_stokes(iend+1,j,k)* &
1663 & vmask(iend+1,j)
1664# endif
1665# ifdef WET_DRY
1666 v_stokes(iend+1,j,k)=v_stokes(iend+1,j,k)* &
1667 & vmask_wet(iend+1,j)
1668# endif
1669# endif
1670 END DO
1671 END IF
1672 END IF
1673!
1674 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1675 IF (j.eq.1) THEN ! southern boundary
1676 DO k=1,n(ng) ! J-loop pipelined
1677 DO i=istr,iend
1678 v(i,j,k,nnew)=v(i,j,k,nnew)-cf(i,0)
1679# ifdef MASKING
1680 v(i,j,k,nnew)=v(i,j,k,nnew)* &
1681 & vmask(i,j)
1682# endif
1683# ifdef WET_DRY
1684 v(i,j,k,nnew)=v(i,j,k,nnew)* &
1685 & vmask_wet(i,j)
1686# endif
1687# ifdef WEC
1688 v_stokes(i,j,k)=v_stokes(i,j,k)-cfs(i,0)
1689# ifdef MASKING
1690 v_stokes(i,j,k)=v_stokes(i,j,k)* &
1691 & vmask(i,j)
1692# endif
1693# ifdef WET_DRY
1694 v_stokes(i,j,k)=v_stokes(i,j,k)* &
1695 & vmask_wet(i,j)
1696# endif
1697# endif
1698 END DO
1699 END DO
1700 END IF
1701 END IF
1702!
1703 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1704 IF (j.eq.mm(ng)+1) THEN ! northern boundary
1705 DO k=1,n(ng) ! J-loop pipelined
1706 DO i=istr,iend
1707 v(i,j,k,nnew)=v(i,j,k,nnew)-cf(i,0)
1708# ifdef MASKING
1709 v(i,j,k,nnew)=v(i,j,k,nnew)* &
1710 & vmask(i,j)
1711# endif
1712# ifdef WET_DRY
1713 v(i,j,k,nnew)=v(i,j,k,nnew)* &
1714 & vmask_wet(i,j)
1715# endif
1716# ifdef WEC
1717 v_stokes(i,j,k)=v_stokes(i,j,k)-cfs(i,0)
1718# ifdef MASKING
1719 v_stokes(i,j,k)=v_stokes(i,j,k)* &
1720 & vmask(i,j)
1721# endif
1722# ifdef WET_DRY
1723 v_stokes(i,j,k)=v_stokes(i,j,k)* &
1724 & vmask_wet(i,j)
1725# endif
1726# endif
1727 END DO
1728 END DO
1729 END IF
1730 END IF
1731!
1732! Compute correct mass flux, Hz*v/m.
1733!
1734 DO k=n(ng),1,-1
1735 DO i=istrt,iendt
1736 hvom(i,j,k)=0.5_r8*(hvom(i,j,k)+v(i,j,k,nnew)*dc(i,k))
1737# ifdef WEC
1738 hvom(i,j,k)=hvom(i,j,k)+0.5_r8*v_stokes(i,j,k)*dc(i,k)
1739# endif
1740 fc(i,0)=fc(i,0)+hvom(i,j,k)
1741# ifdef DIAGNOSTICS_UV
1742 diav3wrk(i,j,k,m3rate)=v(i,j,k,nnew)- &
1743 & diav3wrk(i,j,k,m3rate)
1744# endif
1745 END DO
1746 END DO
1747 DO i=istrt,iendt
1748 fc(i,0)=dc(i,0)*(fc(i,0)-dv_avg2(i,j)) ! recursive
1749 END DO
1750 DO k=1,n(ng)
1751 DO i=istrt,iendt
1752 hvom(i,j,k)=hvom(i,j,k)-dc(i,k)*fc(i,0)
1753 END DO
1754 END DO
1755 END IF
1756 END DO
1757!
1758!-----------------------------------------------------------------------
1759! Exchange boundary data.
1760!-----------------------------------------------------------------------
1761!
1762 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1763 CALL exchange_u3d_tile (ng, tile, &
1764 & lbi, ubi, lbj, ubj, 1, n(ng), &
1765 & u(:,:,:,nnew))
1766 CALL exchange_v3d_tile (ng, tile, &
1767 & lbi, ubi, lbj, ubj, 1, n(ng), &
1768 & v(:,:,:,nnew))
1769# ifdef WEC
1770 CALL exchange_u3d_tile (ng, tile, &
1771 & lbi, ubi, lbj, ubj, 1, n(ng), &
1772 & u_stokes(:,:,:))
1773 CALL exchange_v3d_tile (ng, tile, &
1774 & lbi, ubi, lbj, ubj, 1, n(ng), &
1775 & v_stokes(:,:,:))
1776# endif
1777 CALL exchange_u3d_tile (ng, tile, &
1778 & lbi, ubi, lbj, ubj, 1, n(ng), &
1779 & huon)
1780 CALL exchange_v3d_tile (ng, tile, &
1781 & lbi, ubi, lbj, ubj, 1, n(ng), &
1782 & hvom)
1783!
1784# if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
1785 CALL exchange_u2d_tile (ng, tile, &
1786 & lbi, ubi, lbj, ubj, &
1787 & ubar(:,:,knew))
1788 CALL exchange_v2d_tile (ng, tile, &
1789 & lbi, ubi, lbj, ubj, &
1790 & vbar(:,:,knew))
1791# else
1792 DO k=1,2
1793 CALL exchange_u2d_tile (ng, tile, &
1794 & lbi, ubi, lbj, ubj, &
1795 & ubar(:,:,k))
1796 CALL exchange_v2d_tile (ng, tile, &
1797 & lbi, ubi, lbj, ubj, &
1798 & vbar(:,:,k))
1799 END DO
1800# endif
1801 END IF
1802
1803# ifdef DISTRIBUTE
1804!
1805 CALL mp_exchange3d (ng, tile, inlm, 4, &
1806 & lbi, ubi, lbj, ubj, 1, n(ng), &
1807 & nghostpoints, &
1808 & ewperiodic(ng), nsperiodic(ng), &
1809 & u(:,:,:,nnew), v(:,:,:,nnew), &
1810 & huon, hvom)
1811!
1812# if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
1813 CALL mp_exchange2d (ng, tile, inlm, 2, &
1814 & lbi, ubi, lbj, ubj, &
1815 & nghostpoints, &
1816 & ewperiodic(ng), nsperiodic(ng), &
1817 & ubar(:,:,knew), vbar(:,:,knew))
1818# else
1819 CALL mp_exchange2d (ng, tile, inlm, 4, &
1820 & lbi, ubi, lbj, ubj, &
1821 & nghostpoints, &
1822 & ewperiodic(ng), nsperiodic(ng), &
1823 & ubar(:,:,1), vbar(:,:,1), &
1824 & ubar(:,:,2), vbar(:,:,2))
1825# endif
1826# ifdef WEC
1827 CALL mp_exchange3d (ng, tile, inlm, 2, &
1828 & lbi, ubi, lbj, ubj, 1, n(ng), &
1829 & nghostpoints, &
1830 & ewperiodic(ng), nsperiodic(ng), &
1831 & u_stokes(:,:,:), v_stokes(:,:,:))
1832# endif
1833# endif
1834!
1835!-----------------------------------------------------------------------
1836! Compute 3D momentum (ua, va) at RHO-points (A-Grid) for output
1837! purposes and data assimilation where the observations and state
1838! vector is located at the cell-center.
1839!-----------------------------------------------------------------------
1840!
1841 CALL uv_c2a_grid (ng, tile, inlm, nnew)
1842!
1843 RETURN
1844 END SUBROUTINE step3d_uv_tile
1845#endif
1846 END MODULE step3d_uv_mod
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_u2d_tile
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:440
exchange_2d_mod::exchange_v2d_tile
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:627
exchange_3d_mod
Definition exchange_3d.F:2
exchange_3d_mod::exchange_u3d_tile
subroutine exchange_u3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:495
exchange_3d_mod::exchange_v3d_tile
subroutine exchange_v3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:707
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_forces
Definition mod_forces.F:2
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_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
mod_scalars
Definition mod_scalars.F:2
mod_scalars::m3vvis
integer m3vvis
Definition mod_scalars.F:191
mod_scalars::luvsrc
logical, dimension(:), allocatable luvsrc
Definition mod_scalars.F:542
mod_scalars::m2fcor
integer m2fcor
Definition mod_scalars.F:150
mod_scalars::m3hadv
integer m3hadv
Definition mod_scalars.F:177
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
mod_scalars::m3xadv
integer m3xadv
Definition mod_scalars.F:178
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::lambda
real(r8) lambda
Definition mod_scalars.F:750
mod_scalars::m2pgrd
integer m2pgrd
Definition mod_scalars.F:167
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::m2xadv
integer m2xadv
Definition mod_scalars.F:152
mod_scalars::m2yadv
integer m2yadv
Definition mod_scalars.F:153
mod_scalars::m3vadv
integer m3vadv
Definition mod_scalars.F:176
mod_scalars::m3xvis
integer m3xvis
Definition mod_scalars.F:193
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::m3rate
integer m3rate
Definition mod_scalars.F:195
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::m3yadv
integer m3yadv
Definition mod_scalars.F:179
mod_scalars::m3yvis
integer m3yvis
Definition mod_scalars.F:194
mod_scalars::m2rate
integer m2rate
Definition mod_scalars.F:173
mod_scalars::m2yvis
integer m2yvis
Definition mod_scalars.F:170
mod_scalars::ntfirst
integer, dimension(:), allocatable ntfirst
Definition mod_scalars.F:377
mod_scalars::m3fcor
integer m3fcor
Definition mod_scalars.F:175
mod_scalars::m2sstr
integer m2sstr
Definition mod_scalars.F:171
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::m2hadv
integer m2hadv
Definition mod_scalars.F:151
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::m3pgrd
integer m3pgrd
Definition mod_scalars.F:190
mod_scalars::m2xvis
integer m2xvis
Definition mod_scalars.F:169
mod_scalars::m2bstr
integer m2bstr
Definition mod_scalars.F:172
mod_sources
Definition mod_sources.F:2
mod_sources::sources
type(t_sources), dimension(:), allocatable sources
Definition mod_sources.F:90
mod_sources::nsrc
integer, dimension(:), allocatable nsrc
Definition mod_sources.F:97
mod_stepping
Definition mod_stepping.F:2
mod_stepping::knew
integer, dimension(:), allocatable knew
Definition mod_stepping.F:64
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
step3d_uv_mod
Definition step3d_uv.F:2
step3d_uv_mod::step3d_uv_tile
subroutine step3d_uv_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nstp, nnew, umask, vmask, umask_wet, vmask_wet, om_v, on_u, pm, pn, hz, z_r, z_w, akv, du_avg1, dv_avg1, du_avg2, dv_avg2, ru, rv, diau2wrk, diav2wrk, diau2int, diav2int, diau3wrk, diav3wrk, diaru, diarv, u, v, ubar, vbar, huon, hvom)
Definition step3d_uv.F:173
step3d_uv_mod::step3d_uv
subroutine, public step3d_uv(ng, tile)
Definition step3d_uv.F:50
u3dbc_mod
Definition u3dbc_im.F:2
u3dbc_mod::u3dbc_tile
subroutine, public u3dbc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nstp, nout, u)
Definition u3dbc_im.F:55
uv_var_change_mod
Definition uv_var_change.F:2
uv_var_change_mod::uv_c2a_grid
subroutine, public uv_c2a_grid(ng, tile, model, ninp)
Definition uv_var_change.F:75
v3dbc_mod
Definition v3dbc_im.F:2
v3dbc_mod::v3dbc_tile
subroutine, public v3dbc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nstp, nout, v)
Definition v3dbc_im.F:55
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