ROMS
Loading...
Searching...
No Matches
rhs3d.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE rhs3d_mod
3#if defined NONLINEAR && defined SOLVE3D
4!
5!git $Id$
6!=======================================================================
7! Copyright (c) 2002-2025 The ROMS Group !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md Hernan G. Arango !
10!========================================== Alexander F. Shchepetkin ===
11! !
12! This subroutine evaluates right-hand-side terms for 3D momentum !
13! and tracers equations. !
14! !
15!=======================================================================
16!
17 implicit none
18!
19 PRIVATE
20 PUBLIC :: rhs3d
21!
22 CONTAINS
23!
24!***********************************************************************
25 SUBROUTINE rhs3d (ng, tile)
26!***********************************************************************
27!
28 USE mod_param
29 USE mod_coupling
30# ifdef DIAGNOSTICS_UV
31 USE mod_diags
32# endif
33 USE mod_forces
34 USE mod_grid
35# ifdef WEC
36 USE mod_mixing
37# endif
38 USE mod_ocean
39 USE mod_stepping
40!
41 USE pre_step3d_mod, ONLY : pre_step3d
42 USE prsgrd_mod, ONLY : prsgrd
43# ifndef TS_FIXED
44# ifdef TS_DIF2
45 USE t3dmix2_mod, ONLY : t3dmix2
46# endif
47# ifdef TS_DIF4
48 USE t3dmix4_mod, ONLY : t3dmix4
49# endif
50# endif
51# ifdef UV_VIS2
52 USE uv3dmix2_mod, ONLY : uv3dmix2
53# endif
54# ifdef UV_VIS4
55 USE uv3dmix4_mod, ONLY : uv3dmix4
56# endif
57# ifdef WEC_VF
58 USE wec_vf_mod, ONLY : wec_vf
59# if defined BOTTOM_STREAMING || defined SURFACE_STREAMING
60 USE wec_streaming_mod, ONLY : wec_streaming
61# endif
62# endif
63!
64! Imported variable declarations.
65!
66 integer, intent(in) :: ng, tile
67!
68! Local variable declarations.
69!
70 character (len=*), parameter :: myfile = &
71 & __FILE__
72!
73# include "tile.h"
74!
75!-----------------------------------------------------------------------
76! Initialize computations for new time step of the 3D primitive
77! variables.
78!-----------------------------------------------------------------------
79!
80 CALL pre_step3d (ng, tile)
81!
82!-----------------------------------------------------------------------
83! Compute baroclinic pressure gradient.
84!-----------------------------------------------------------------------
85!
86 CALL prsgrd (ng, tile)
87
88# ifdef WEC_VF
89!
90!-----------------------------------------------------------------------
91! Compute Waves Effect on Currents.
92!-----------------------------------------------------------------------
93!
94 CALL wec_vf (ng, tile)
95# if defined BOTTOM_STREAMING || defined SURFACE_STREAMING
96 CALL wec_streaming (ng, tile)
97# endif
98# endif
99
100# ifndef TS_FIXED
101# ifdef TS_DIF2
102!
103!-----------------------------------------------------------------------
104! Compute horizontal harmonic mixing of tracer type variables.
105!-----------------------------------------------------------------------
106!
107 CALL t3dmix2 (ng, tile)
108# endif
109# ifdef TS_DIF4
110!
111!-----------------------------------------------------------------------
112! Compute horizontal biharmonic mixing of tracer type variables.
113!-----------------------------------------------------------------------
114!
115 CALL t3dmix4 (ng, tile)
116# endif
117# endif
118!
119!-----------------------------------------------------------------------
120! Compute right-hand-side terms for the 3D momentum equations.
121!-----------------------------------------------------------------------
122!
123# ifdef PROFILE
124 CALL wclock_on (ng, inlm, 21, __line__, myfile)
125# endif
126 CALL rhs3d_tile (ng, tile, &
127 & lbi, ubi, lbj, ubj, &
128 & imins, imaxs, jmins, jmaxs, &
129 & nrhs(ng), &
130 & grid(ng) % Hz, &
131 & grid(ng) % Huon, &
132 & grid(ng) % Hvom, &
133# if defined CURVGRID && defined UV_ADV
134 & grid(ng) % dmde, &
135 & grid(ng) % dndx, &
136# endif
137 & grid(ng) % fomn, &
138 & grid(ng) % om_u, &
139 & grid(ng) % om_v, &
140 & grid(ng) % on_u, &
141 & grid(ng) % on_v, &
142 & grid(ng) % pm, &
143 & grid(ng) % pn, &
144# ifdef WET_DRY
145 & grid(ng)%umask_wet, &
146 & grid(ng)%vmask_wet, &
147# endif
148 & forces(ng) % bustr, &
149 & forces(ng) % bvstr, &
150 & forces(ng) % sustr, &
151 & forces(ng) % svstr, &
152 & ocean(ng) % u, &
153 & ocean(ng) % v, &
154 & ocean(ng) % W, &
155# ifdef WEC
156 & ocean(ng) % u_stokes, &
157 & ocean(ng) % v_stokes, &
158 & ocean(ng) % W_stokes, &
159 & mixing(ng) % rustr3d, &
160 & mixing(ng) % rvstr3d, &
161# endif
162 & coupling(ng) % rufrc, &
163 & coupling(ng) % rvfrc, &
164# ifdef DIAGNOSTICS_UV
165 & diags(ng) % DiaRUfrc, &
166 & diags(ng) % DiaRVfrc, &
167 & diags(ng) % DiaRU, &
168 & diags(ng) % DiaRV, &
169# endif
170 & ocean(ng) % ru, &
171 & ocean(ng) % rv)
172# ifdef PROFILE
173 CALL wclock_off (ng, inlm, 21, __line__, myfile)
174# endif
175# ifdef UV_VIS2
176!
177!-----------------------------------------------------------------------
178! Compute horizontal, harmonic mixing of momentum.
179!-----------------------------------------------------------------------
180!
181 CALL uv3dmix2 (ng, tile)
182# endif
183# ifdef UV_VIS4
184!
185!-----------------------------------------------------------------------
186! Compute horizontal, biharmonic mixing of momentum.
187!-----------------------------------------------------------------------
188!
189 CALL uv3dmix4 (ng, tile)
190# endif
191!
192 RETURN
193 END SUBROUTINE rhs3d
194!
195!***********************************************************************
196 SUBROUTINE rhs3d_tile (ng, tile, &
197 & LBi, UBi, LBj, UBj, &
198 & IminS, ImaxS, JminS, JmaxS, &
199 & nrhs, &
200 & Hz, Huon, Hvom, &
201# if defined CURVGRID && defined UV_ADV
202 & dmde, dndx, &
203# endif
204 & fomn, &
205 & om_u, om_v, on_u, on_v, pm, pn, &
206# ifdef WET_DRY
207 & umask_wet, vmask_wet, &
208# endif
209 & bustr, bvstr, &
210 & sustr, svstr, &
211 & u, v, W, &
212# ifdef WEC
213 & u_stokes, v_stokes, W_stokes, &
214 & rustr3d, rvstr3d, &
215# endif
216 & rufrc, rvfrc, &
217# ifdef DIAGNOSTICS_UV
218 & DiaRUfrc, DiaRVfrc, &
219 & DiaRU, DiaRV, &
220# endif
221 & ru, rv)
222!***********************************************************************
223!
224 USE mod_param
225 USE mod_clima
226 USE mod_scalars
227!
228! Imported variable declarations.
229!
230 integer, intent(in) :: ng, tile
231 integer, intent(in) :: LBi, UBi, LBj, UBj
232 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
233 integer, intent(in) :: nrhs
234!
235# ifdef ASSUMED_SHAPE
236 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
237 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
238 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
239# if defined CURVGRID && defined UV_ADV
240 real(r8), intent(in) :: dmde(LBi:,LBj:)
241 real(r8), intent(in) :: dndx(LBi:,LBj:)
242# endif
243 real(r8), intent(in) :: fomn(LBi:,LBj:)
244 real(r8), intent(in) :: om_u(LBi:,LBj:)
245 real(r8), intent(in) :: om_v(LBi:,LBj:)
246 real(r8), intent(in) :: on_u(LBi:,LBj:)
247 real(r8), intent(in) :: on_v(LBi:,LBj:)
248 real(r8), intent(in) :: pm(LBi:,LBj:)
249 real(r8), intent(in) :: pn(LBi:,LBj:)
250# ifdef WET_DRY
251 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
252 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
253# endif
254 real(r8), intent(in) :: bustr(LBi:,LBj:)
255 real(r8), intent(in) :: bvstr(LBi:,LBj:)
256 real(r8), intent(in) :: sustr(LBi:,LBj:)
257 real(r8), intent(in) :: svstr(LBi:,LBj:)
258 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
259 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
260 real(r8), intent(in) :: W(LBi:,LBj:,0:)
261# ifdef WEC
262 real(r8), intent(in) :: u_stokes(LBi:,LBj:,:)
263 real(r8), intent(in) :: v_stokes(LBi:,LBj:,:)
264 real(r8), intent(in) :: W_stokes(LBi:,LBj:,0:)
265 real(r8), intent(in) :: rustr3d(LBi:,LBj:,:)
266 real(r8), intent(in) :: rvstr3d(LBi:,LBj:,:)
267# endif
268 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
269 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
270# ifdef DIAGNOSTICS_UV
271 real(r8), intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
272 real(r8), intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
273 real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
274 real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
275# endif
276
277 real(r8), intent(out) :: rufrc(LBi:,LBj:)
278 real(r8), intent(out) :: rvfrc(LBi:,LBj:)
279# else
280 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
281 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
282 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
283# if defined CURVGRID && defined UV_ADV
284 real(r8), intent(in) :: dmde(LBi:UBi,LBj:UBj)
285 real(r8), intent(in) :: dndx(LBi:UBi,LBj:UBj)
286# endif
287 real(r8), intent(in) :: fomn(LBi:UBi,LBj:UBj)
288 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
289 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
290 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
291 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
292 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
293 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
294# ifdef WET_DRY
295 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
296 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
297# endif
298 real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
299 real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
300 real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj)
301 real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj)
302 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
303 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
304 real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
305# ifdef WEC
306 real(r8), intent(in) :: u_stokes(LBi:UBi,LBj:UBj,N(ng))
307 real(r8), intent(in) :: v_stokes(LBi:UBi,LBj:UBj,N(ng))
308 real(r8), intent(in) :: W_stokes(LBi:UBi,LBj:UBj,0:N(ng))
309 real(r8), intent(in) :: rustr3d(LBi:UBi,LBj:UBj,N(ng))
310 real(r8), intent(in) :: rvstr3d(LBi:UBi,LBj:UBj,N(ng))
311# endif
312 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:N(ng),2)
313 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:N(ng),2)
314# ifdef DIAGNOSTICS_UV
315 real(r8), intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
316 real(r8), intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
317 real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
318 real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
319# endif
320
321 real(r8), intent(out) :: rufrc(LBi:UBi,LBj:UBj)
322 real(r8), intent(out) :: rvfrc(LBi:UBi,LBj:UBj)
323# endif
324!
325! Local variable declarations.
326!
327 integer :: i, j, k
328
329 real(r8), parameter :: Gadv = -0.25_r8
330
331 real(r8) :: cff, cff1, cff2, cff3, cff4
332 real(r8) :: fac, fac1, fac2
333
334 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
335 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
336 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
337# ifdef WEC_VF
338 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCs
339# endif
340
341 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Huee
342 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Huxx
343 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hvee
344 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Hvxx
345 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
346 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
347 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Uwrk
348 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
349 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
350 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vwrk
351 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: uee
352 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: uxx
353 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: vee
354 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: vxx
355 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: wrk
356
357# include "set_bounds.h"
358
359# ifdef BODYFORCE
360!
361!-----------------------------------------------------------------------
362! Apply surface stress as a bodyforce: determine the thickness (m)
363! of the surface layer; then add in surface stress as a bodyfoce.
364!-----------------------------------------------------------------------
365!
366# ifdef DIAGNOSTICS_UV
367 DO k=1,n(ng)
368 DO j=jstr,jend
369 DO i=istr,iend
370 diaru(i,j,k,nrhs,m3vvis)=0.0_r8
371 diarv(i,j,k,nrhs,m3vvis)=0.0_r8
372 END DO
373 END DO
374 END DO
375 DO j=jstr,jend
376 DO i=istru,iend
377 diarufrc(i,j,3,m2sstr)=0.0_r8
378 diarufrc(i,j,3,m2bstr)=0.0_r8
379 END DO
380 END DO
381 DO j=jstrv,jend
382 DO i=istr,iend
383 diarvfrc(i,j,3,m2sstr)=0.0_r8
384 diarvfrc(i,j,3,m2bstr)=0.0_r8
385 END DO
386 END DO
387# endif
388 DO j=jstrv-1,jend
389 DO i=istru-1,iend
390 wrk(i,j)=0.0_r8
391 END DO
392 END DO
393 DO k=n(ng),levsfrc(ng),-1
394 DO j=jstrv-1,jend
395 DO i=istru-1,iend
396 wrk(i,j)=wrk(i,j)+hz(i,j,k)
397 END DO
398 END DO
399 END DO
400 DO j=jstr,jend
401 DO i=istru,iend
402 cff=0.25_r8*(pm(i-1,j)+pm(i,j))* &
403 & (pn(i-1,j)+pn(i,j))
404 cff1=1.0_r8/(cff*(wrk(i-1,j)+wrk(i,j)))
405 uwrk(i,j)=sustr(i,j)*cff1
406 END DO
407 END DO
408 DO j=jstrv,jend
409 DO i=istr,iend
410 cff=0.25*(pm(i,j-1)+pm(i,j))* &
411 & (pn(i,j-1)+pn(i,j))
412 cff1=1.0_r8/(cff*(wrk(i,j-1)+wrk(i,j)))
413 vwrk(i,j)=svstr(i,j)*cff1
414 END DO
415 END DO
416 DO k=levsfrc(ng),n(ng)
417 DO j=jstr,jend
418 DO i=istru,iend
419 cff=uwrk(i,j)*(hz(i ,j,k)+ &
420 & hz(i-1,j,k))
421 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)+cff
422# ifdef DIAGNOSTICS_UV
423 diaru(i,j,k,nrhs,m3vvis)=diaru(i,j,k,nrhs,m3vvis)+cff
424 diarufrc(i,j,3,m2sstr)=diarufrc(i,j,3,m2sstr)+cff
425# endif
426 END DO
427 END DO
428 DO j=jstrv,jend
429 DO i=istr,iend
430 cff=vwrk(i,j)*(hz(i,j ,k)+ &
431 & hz(i,j-1,k))
432 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)+cff
433# ifdef DIAGNOSTICS_UV
434 diarv(i,j,k,nrhs,m3vvis)=diarv(i,j,k,nrhs,m3vvis)+cff
435 diarvfrc(i,j,3,m2sstr)=diarvfrc(i,j,3,m2sstr)+cff
436# endif
437 END DO
438 END DO
439 END DO
440!
441! Apply bottom stress as a bodyforce: determine the thickness (m)
442! of the bottom layer; then add in bottom stress as a bodyfoce.
443!
444 DO j=jstrv-1,jend
445 DO i=istru-1,iend
446 wrk(i,j)=0.0_r8
447 END DO
448 END DO
449 DO k=1,levbfrc(ng)
450 DO j=jstrv-1,jend
451 DO i=istru-1,iend
452 wrk(i,j)=wrk(i,j)+hz(i,j,k)
453 END DO
454 END DO
455 END DO
456 DO j=jstr,jend
457 DO i=istru,iend
458 cff=0.25_r8*(pm(i-1,j)+pm(i,j))* &
459 & (pn(i-1,j)+pn(i,j))
460 cff1=1.0_r8/(cff*(wrk(i-1,j)+wrk(i,j)))
461 uwrk(i,j)=bustr(i,j)*cff1
462 END DO
463 END DO
464 DO j=jstrv,jend
465 DO i=istr,iend
466 cff=0.25_r8*(pm(i,j-1)+pm(i,j))* &
467 & (pn(i,j-1)+pn(i,j))
468 cff1=1.0_r8/(cff*(wrk(i,j-1)+wrk(i,j)))
469 vwrk(i,j)=bvstr(i,j)*cff1
470 END DO
471 END DO
472 DO k=1,levbfrc(ng)
473 DO j=jstr,jend
474 DO i=istru,iend
475 cff=uwrk(i,j)*(hz(i ,j,k)+ &
476 & hz(i-1,j,k))
477 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)-cff
478# ifdef DIAGNOSTICS_UV
479 diaru(i,j,k,nrhs,m3vvis)=diaru(i,j,k,nrhs,m3vvis)-cff
480 diarufrc(i,j,3,m2bstr)=diarufrc(i,j,3,m2bstr)-cff
481# endif
482 END DO
483 END DO
484 DO j=jstrv,jend
485 DO i=istr,iend
486 cff=vwrk(i,j)*(hz(i,j ,k)+ &
487 & hz(i,j-1,k))
488 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)-cff
489# ifdef DIAGNOSTICS_UV
490 diarv(i,j,k,nrhs,m3vvis)=diarv(i,j,k,nrhs,m3vvis)-cff
491 diarvfrc(i,j,3,m2bstr)=diarvfrc(i,j,3,m2bstr)-cff
492# endif
493 END DO
494 END DO
495 END DO
496# endif
497!
498 k_loop : DO k=1,n(ng)
499
500# ifdef UV_COR
501!
502!-----------------------------------------------------------------------
503! Add in Coriolis terms.
504!-----------------------------------------------------------------------
505!
506 DO j=jstrv-1,jend
507 DO i=istru-1,iend
508 cff=0.5_r8*hz(i,j,k)*fomn(i,j)
509 ufx(i,j)=cff*(v(i,j ,k,nrhs)+ &
510 & v(i,j+1,k,nrhs))
511 vfe(i,j)=cff*(u(i ,j,k,nrhs)+ &
512 & u(i+1,j,k,nrhs))
513 END DO
514 END DO
515 DO j=jstr,jend
516 DO i=istru,iend
517 cff1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
518 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)+cff1
519# ifdef DIAGNOSTICS_UV
520 diaru(i,j,k,nrhs,m3fcor)=cff1
521# endif
522 END DO
523 END DO
524 DO j=jstrv,jend
525 DO i=istr,iend
526 cff1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
527 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)-cff1
528# ifdef DIAGNOSTICS_UV
529 diarv(i,j,k,nrhs,m3fcor)=-cff1
530# endif
531 END DO
532 END DO
533!
534# ifdef WEC
535 DO j=jstrv-1,jend
536 DO i=istru-1,iend
537 cff=0.5_r8*hz(i,j,k)*fomn(i,j)
538 ufx(i,j)=cff*(v_stokes(i,j ,k)+ &
539 & v_stokes(i,j+1,k))
540 vfe(i,j)=cff*(u_stokes(i ,j,k)+ &
541 & u_stokes(i+1,j,k))
542 END DO
543 END DO
544 DO j=jstr,jend
545 DO i=istru,iend
546 cff1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
547 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)+cff1
548# ifdef DIAGNOSTICS_UV
549 diaru(i,j,k,nrhs,m3fsco)=cff1
550# endif
551 END DO
552 END DO
553 DO j=jstrv,jend
554 DO i=istr,iend
555 cff1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
556 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)-cff1
557# ifdef DIAGNOSTICS_UV
558 diarv(i,j,k,nrhs,m3fsco)=-cff1
559# endif
560 END DO
561 END DO
562# endif
563# endif
564# if defined CURVGRID && defined UV_ADV
565!
566!-----------------------------------------------------------------------
567! Add in curvilinear transformation terms.
568!-----------------------------------------------------------------------
569!
570 DO j=jstrv-1,jend
571 DO i=istru-1,iend
572 cff1=0.5_r8*(v(i,j ,k,nrhs)+ &
573# ifdef WEC
574 & v_stokes(i,j ,k)+ &
575 & v_stokes(i,j+1,k)+ &
576# endif
577 & v(i,j+1,k,nrhs))
578 cff2=0.5_r8*(u(i ,j,k,nrhs)+ &
579# ifdef WEC
580 & u_stokes(i ,j,k)+ &
581 & u_stokes(i+1,j,k)+ &
582# endif
583 & u(i+1,j,k,nrhs))
584 cff3=cff1*dndx(i,j)
585 cff4=cff2*dmde(i,j)
586# ifdef WEC_VF
587 cff5=0.5_r8*(v_stokes(i,j ,k)+ &
588 & v_stokes(i,j+1,k))
589 cff6=0.5_r8*(u_stokes(i ,j,k)+ &
590 & u_stokes(i+1,j,k))
591 cff7=cff5*dndx(i,j)
592 cff8=cff6*dmde(i,j)
593# endif
594 cff=hz(i,j,k)*(cff3-cff4)
595 ufx(i,j)=cff*cff1
596 vfe(i,j)=cff*cff2
597# ifdef WEC_VF
598 ufx(i,j)=ufx(i,j)-(cff5*hz(i,j,k)*(cff7-cff8))
599 vfe(i,j)=vfe(i,j)+(cff6*hz(i,j,k)*(cff7-cff8))
600# endif
601# if defined DIAGNOSTICS_UV
602 cff=hz(i,j,k)*cff4
603 uwrk(i,j)=-cff*cff1 ! u equation, ETA-term
604 vwrk(i,j)=-cff*cff2 ! v equation, ETA-term
605# ifdef WEC_VF
606 uwrk(i,j)=uwrk(i,j)+hz(i,j,k)*cff5*cff8
607 vwrk(i,j)=vwrk(i,j)-hz(i,j,k)*cff6*cff8
608# endif
609# endif
610 END DO
611 END DO
612 DO j=jstr,jend
613 DO i=istru,iend
614 cff1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
615 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)+cff1
616# ifdef DIAGNOSTICS_UV
617 cff2=0.5_r8*(uwrk(i,j)+uwrk(i-1,j))
618# ifdef WEC_VF
619 diaru(i,j,k,nrhs,m3xadv)=diaru(i,j,k,nrhs,m3xadv)+cff1-cff2
620 diaru(i,j,k,nrhs,m3yadv)=diaru(i,j,k,nrhs,m3yadv)+cff2
621 diaru(i,j,k,nrhs,m3hadv)=diaru(i,j,k,nrhs,m3hadv)+cff1
622# else
623 diaru(i,j,k,nrhs,m3xadv)=cff1-cff2
624 diaru(i,j,k,nrhs,m3yadv)=cff2
625 diaru(i,j,k,nrhs,m3hadv)=cff1
626# endif
627# endif
628 END DO
629 END DO
630 DO j=jstrv,jend
631 DO i=istr,iend
632 cff1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
633 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)-cff1
634# ifdef DIAGNOSTICS_UV
635 cff2=0.5_r8*(vwrk(i,j)+vwrk(i,j-1))
636# ifdef WEC_VF
637 diarv(i,j,k,nrhs,m3xadv)=diarv(i,j,k,nrhs,m3xadv)-cff1+cff2
638 diarv(i,j,k,nrhs,m3yadv)=diarv(i,j,k,nrhs,m3yadv)-cff2
639 diarv(i,j,k,nrhs,m3hadv)=diarv(i,j,k,nrhs,m3hadv)-cff1
640# else
641 diarv(i,j,k,nrhs,m3xadv)=-cff1+cff2
642 diarv(i,j,k,nrhs,m3yadv)=-cff2
643 diarv(i,j,k,nrhs,m3hadv)=-cff1
644# endif
645# endif
646 END DO
647 END DO
648# endif
649!
650!-----------------------------------------------------------------------
651! Add in nudging of 3D momentum climatology.
652!-----------------------------------------------------------------------
653!
654 IF (lnudgem3clm(ng)) THEN
655 DO j=jstr,jend
656 DO i=istru,iend
657 cff=0.25_r8*(clima(ng)%M3nudgcof(i-1,j,k)+ &
658 & clima(ng)%M3nudgcof(i ,j,k))* &
659 & om_u(i,j)*on_u(i,j)
660 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)+ &
661 & cff*(hz(i-1,j,k)+hz(i,j,k))* &
662 & (clima(ng)%uclm(i,j,k)- &
663 & u(i,j,k,nrhs))
664 END DO
665 END DO
666 DO j=jstrv,jend
667 DO i=istr,iend
668 cff=0.25_r8*(clima(ng)%M3nudgcof(i,j-1,k)+ &
669 & clima(ng)%M3nudgcof(i,j ,k))* &
670 & om_v(i,j)*on_v(i,j)
671 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)+ &
672 & cff*(hz(i,j-1,k)+hz(i,j,k))* &
673 & (clima(ng)%vclm(i,j,k)- &
674 & v(i,j,k,nrhs))
675 END DO
676 END DO
677 END IF
678
679# ifdef UV_ADV
680!
681!-----------------------------------------------------------------------
682! Add in horizontal advection of momentum.
683!-----------------------------------------------------------------------
684!
685! Compute diagonal [UFx,VFe] and off-diagonal [UFe,VFx] components
686! of tensor of momentum flux due to horizontal advection.
687!
688# ifdef UV_C2ADVECTION
689!
690! Second-order, centered differences advection.
691!
692 DO j=jstr,jend
693 DO i=istru-1,iend
694 ufx(i,j)=0.25_r8*(u(i ,j,k,nrhs)+ &
695 & u(i+1,j,k,nrhs))* &
696 & (huon(i ,j,k)+ &
697 & huon(i+1,j,k))
698 END DO
699 END DO
700 DO j=jstr,jend+1
701 DO i=istru,iend
702 ufe(i,j)=0.25_r8*(u(i,j-1,k,nrhs)+ &
703 & u(i,j ,k,nrhs))* &
704 & (hvom(i-1,j,k)+ &
705 & hvom(i ,j,k))
706 END DO
707 END DO
708 DO j=jstrv,jend
709 DO i=istr,iend+1
710 vfx(i,j)=0.25_r8*(v(i-1,j,k,nrhs)+ &
711 & v(i ,j,k,nrhs))* &
712 & (huon(i,j-1,k)+ &
713 & huon(i,j ,k))
714 END DO
715 END DO
716 DO j=jstrv-1,jend
717 DO i=istr,iend
718 vfe(i,j)=0.25_r8*(v(i,j ,k,nrhs)+ &
719 & v(i,j+1,k,nrhs))* &
720 & (hvom(i,j ,k)+ &
721 & hvom(i,j+1,k))
722 END DO
723 END DO
724# else
725 DO j=jstr,jend
726 DO i=istrum1,iendp1
727 uxx(i,j)=u(i-1,j,k,nrhs)-2.0_r8*u(i,j,k,nrhs)+ &
728 & u(i+1,j,k,nrhs)
729 huxx(i,j)=huon(i-1,j,k)-2.0_r8*huon(i,j,k)+huon(i+1,j,k)
730 END DO
731 END DO
732 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
733 IF (domain(ng)%Western_Edge(tile)) THEN
734 DO j=jstr,jend
735 uxx(istr,j)=uxx(istr+1,j)
736 huxx(istr,j)=huxx(istr+1,j)
737 END DO
738 END IF
739 END IF
740 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
741 IF (domain(ng)%Eastern_Edge(tile)) THEN
742 DO j=jstr,jend
743 uxx(iend+1,j)=uxx(iend,j)
744 huxx(iend+1,j)=huxx(iend,j)
745 END DO
746 END IF
747 END IF
748# ifdef UV_C4ADVECTION
749!
750! Fourth-order, centered differences u-momentum horizontal advection.
751!
752 cff=1.0_r8/6.0_r8
753 DO j=jstr,jend
754 DO i=istru-1,iend
755 ufx(i,j)=0.25_r8*(u(i ,j,k,nrhs)+ &
756 & u(i+1,j,k,nrhs)- &
757 & cff*(uxx(i ,j)+ &
758 & uxx(i+1,j)))* &
759 & (huon(i ,j,k)+ &
760 & huon(i+1,j,k)- &
761 & cff*(huxx(i ,j)+ &
762 & huxx(i+1,j)))
763 END DO
764 END DO
765# else
766!
767! Third-order, upstream bias u-momentum advection with velocity
768! dependent hyperdiffusion.
769!
770 DO j=jstr,jend
771 DO i=istru-1,iend
772 cff1=u(i ,j,k,nrhs)+ &
773 & u(i+1,j,k,nrhs)
774 IF (cff1.gt.0.0_r8) THEN
775 cff=uxx(i,j)
776 ELSE
777 cff=uxx(i+1,j)
778 END IF
779 ufx(i,j)=0.25_r8*(cff1+gadv*cff)* &
780 & (huon(i ,j,k)+ &
781 & huon(i+1,j,k)+ &
782 & gadv*0.5_r8*(huxx(i ,j)+ &
783 & huxx(i+1,j)))
784 END DO
785 END DO
786# endif
787 DO j=jstrm1,jendp1
788 DO i=istru,iend
789 uee(i,j)=u(i,j-1,k,nrhs)-2.0_r8*u(i,j,k,nrhs)+ &
790 & u(i,j+1,k,nrhs)
791 END DO
792 END DO
793 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
794 IF (domain(ng)%Southern_Edge(tile)) THEN
795 DO i=istru,iend
796 uee(i,jstr-1)=uee(i,jstr)
797 END DO
798 END IF
799 END IF
800 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
801 IF (domain(ng)%Northern_Edge(tile)) THEN
802 DO i=istru,iend
803 uee(i,jend+1)=uee(i,jend)
804 END DO
805 END IF
806 END IF
807 DO j=jstr,jend+1
808 DO i=istru-1,iend
809 hvxx(i,j)=hvom(i-1,j,k)-2.0_r8*hvom(i,j,k)+hvom(i+1,j,k)
810 END DO
811 END DO
812# ifdef UV_C4ADVECTION
813 cff=1.0_r8/6.0_r8
814 DO j=jstr,jend+1
815 DO i=istru,iend
816 ufe(i,j)=0.25_r8*(u(i,j ,k,nrhs)+ &
817 & u(i,j-1,k,nrhs)- &
818 & cff*(uee(i,j )+ &
819 & uee(i,j-1)))* &
820 & (hvom(i ,j,k)+ &
821 & hvom(i-1,j,k)- &
822 & cff*(hvxx(i ,j)+ &
823 & hvxx(i-1,j)))
824 END DO
825 END DO
826# else
827 DO j=jstr,jend+1
828 DO i=istru,iend
829 cff1=u(i,j ,k,nrhs)+ &
830 & u(i,j-1,k,nrhs)
831 cff2=hvom(i,j,k)+hvom(i-1,j,k)
832 IF (cff2.gt.0.0_r8) THEN
833 cff=uee(i,j-1)
834 ELSE
835 cff=uee(i,j)
836 END IF
837 ufe(i,j)=0.25_r8*(cff1+gadv*cff)* &
838 & (cff2+gadv*0.5_r8*(hvxx(i ,j)+ &
839 & hvxx(i-1,j)))
840 END DO
841 END DO
842# endif
843 DO j=jstrv,jend
844 DO i=istrm1,iendp1
845 vxx(i,j)=v(i-1,j,k,nrhs)-2.0_r8*v(i,j,k,nrhs)+ &
846 & v(i+1,j,k,nrhs)
847 END DO
848 END DO
849 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
850 IF (domain(ng)%Western_Edge(tile)) THEN
851 DO j=jstrv,jend
852 vxx(istr-1,j)=vxx(istr,j)
853 END DO
854 END IF
855 END IF
856 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
857 IF (domain(ng)%Eastern_Edge(tile)) THEN
858 DO j=jstrv,jend
859 vxx(iend+1,j)=vxx(iend,j)
860 END DO
861 END IF
862 END IF
863 DO j=jstrv-1,jend
864 DO i=istr,iend+1
865 huee(i,j)=huon(i,j-1,k)-2.0_r8*huon(i,j,k)+huon(i,j+1,k)
866 END DO
867 END DO
868# ifdef UV_C4ADVECTION
869!
870! Fourth-order, centered differences v-momentum horizontal advection.
871!
872 cff=1.0_r8/6.0_r8
873 DO j=jstrv,jend
874 DO i=istr,iend+1
875 vfx(i,j)=0.25_r8*(v(i ,j,k,nrhs)+ &
876 & v(i-1,j,k,nrhs)- &
877 & cff*(vxx(i ,j)+ &
878 & vxx(i-1,j)))* &
879 & (huon(i,j ,k)+ &
880 & huon(i,j-1,k)- &
881 & cff*(huee(i,j )+ &
882 & huee(i,j-1)))
883 END DO
884 END DO
885# else
886!
887! Third-order, upstream bias v-momentum advection with velocity
888! dependent hyperdiffusion.
889!
890 DO j=jstrv,jend
891 DO i=istr,iend+1
892 cff1=v(i ,j,k,nrhs)+ &
893 & v(i-1,j,k,nrhs)
894 cff2=huon(i,j,k)+huon(i,j-1,k)
895 IF (cff2.gt.0.0_r8) THEN
896 cff=vxx(i-1,j)
897 ELSE
898 cff=vxx(i,j)
899 END IF
900 vfx(i,j)=0.25_r8*(cff1+gadv*cff)* &
901 & (cff2+gadv*0.5_r8*(huee(i,j )+ &
902 & huee(i,j-1)))
903 END DO
904 END DO
905# endif
906 DO j=jstrvm1,jendp1
907 DO i=istr,iend
908 vee(i,j)=v(i,j-1,k,nrhs)-2.0_r8*v(i,j,k,nrhs)+ &
909 & v(i,j+1,k,nrhs)
910 hvee(i,j)=hvom(i,j-1,k)-2.0_r8*hvom(i,j,k)+hvom(i,j+1,k)
911 END DO
912 END DO
913 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
914 IF (domain(ng)%Southern_Edge(tile)) THEN
915 DO i=istr,iend
916 vee(i,jstr)=vee(i,jstr+1)
917 hvee(i,jstr)=hvee(i,jstr+1)
918 END DO
919 END IF
920 END IF
921 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
922 IF (domain(ng)%Northern_Edge(tile)) THEN
923 DO i=istr,iend
924 vee(i,jend+1)=vee(i,jend)
925 hvee(i,jend+1)=hvee(i,jend)
926 END DO
927 END IF
928 END IF
929# ifdef UV_C4ADVECTION
930 cff=1.0_r8/6.0_r8
931 DO j=jstrv-1,jend
932 DO i=istr,iend
933 vfe(i,j)=0.25_r8*(v(i,j ,k,nrhs)+ &
934 & v(i,j+1,k,nrhs)- &
935 & cff*(vee(i,j )+ &
936 & vee(i,j+1)))* &
937 & (hvom(i,j ,k)+ &
938 & hvom(i,j+1,k)- &
939 & cff*(hvee(i,j )+ &
940 & hvee(i,j+1)))
941 END DO
942 END DO
943# else
944 DO j=jstrv-1,jend
945 DO i=istr,iend
946 cff1=v(i,j ,k,nrhs)+ &
947 & v(i,j+1,k,nrhs)
948 IF (cff1.gt.0.0_r8) THEN
949 cff=vee(i,j)
950 ELSE
951 cff=vee(i,j+1)
952 END IF
953 vfe(i,j)=0.25_r8*(cff1+gadv*cff)* &
954 & (hvom(i,j ,k)+ &
955 & hvom(i,j+1,k)+ &
956 & gadv*0.5_r8*(hvee(i,j )+ &
957 & hvee(i,j+1)))
958 END DO
959 END DO
960# endif
961# endif
962!
963! Add in horizontal advection.
964!
965 DO j=jstr,jend
966 DO i=istru,iend
967 cff1=ufx(i,j)-ufx(i-1,j)
968 cff2=ufe(i,j+1)-ufe(i,j)
969 cff=cff1+cff2
970 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)-cff
971# ifdef DIAGNOSTICS_UV
972# if defined CURVGRID || defined WEC_VF
973 diaru(i,j,k,nrhs,m3xadv)=diaru(i,j,k,nrhs,m3xadv)-cff1
974 diaru(i,j,k,nrhs,m3yadv)=diaru(i,j,k,nrhs,m3yadv)-cff2
975 diaru(i,j,k,nrhs,m3hadv)=diaru(i,j,k,nrhs,m3hadv)-cff
976# else
977 diaru(i,j,k,nrhs,m3xadv)=-cff1
978 diaru(i,j,k,nrhs,m3yadv)=-cff2
979 diaru(i,j,k,nrhs,m3hadv)=-cff
980# endif
981# endif
982 END DO
983 END DO
984 DO j=jstrv,jend
985 DO i=istr,iend
986 cff1=vfx(i+1,j)-vfx(i,j)
987 cff2=vfe(i,j)-vfe(i,j-1)
988 cff=cff1+cff2
989 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)-cff
990# ifdef DIAGNOSTICS_UV
991# if defined CURVGRID || defined WEC_VF
992 diarv(i,j,k,nrhs,m3xadv)=diarv(i,j,k,nrhs,m3xadv)-cff1
993 diarv(i,j,k,nrhs,m3yadv)=diarv(i,j,k,nrhs,m3yadv)-cff2
994 diarv(i,j,k,nrhs,m3hadv)=diarv(i,j,k,nrhs,m3hadv)-cff
995# else
996 diarv(i,j,k,nrhs,m3xadv)=-cff1
997 diarv(i,j,k,nrhs,m3yadv)=-cff2
998 diarv(i,j,k,nrhs,m3hadv)=-cff
999# endif
1000# endif
1001 END DO
1002 END DO
1003# endif
1004# ifdef WEC
1005!
1006!-----------------------------------------------------------------------
1007! Add in radiation stress or Vortex Force terms.
1008! Convert stresses to m4/s2.
1009!-----------------------------------------------------------------------
1010!
1011 DO j=jstr,jend
1012 DO i=istru,iend
1013 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)- &
1014 & rustr3d(i,j,k)*om_u(i,j)*on_u(i,j)
1015 END DO
1016 END DO
1017 DO j=jstrv,jend
1018 DO i=istr,iend
1019 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)- &
1020 & rvstr3d(i,j,k)*om_v(i,j)*on_v(i,j)
1021 END DO
1022 END DO
1023# endif
1024
1025 END DO k_loop
1026!
1027 j_loop : DO j=jstr,jend
1028# ifdef UV_ADV
1029!
1030!-----------------------------------------------------------------------
1031! Add in vertical advection.
1032!-----------------------------------------------------------------------
1033!
1034# ifdef UV_SADVECTION
1035!
1036! Construct conservative parabolic splines for the vertical
1037! derivatives "CF" of u-momentum.
1038!
1039 cff1=9.0_r8/16.0_r8
1040 cff2=1.0_r8/16.0_r8
1041 DO k=1,n(ng)
1042 DO i=istru,iend
1043 dc(i,k)=cff1*(hz(i ,j,k)+ &
1044 & hz(i-1,j,k))- &
1045 & cff2*(hz(i+1,j,k)+ &
1046 & hz(i-2,j,k))
1047 END DO
1048 END DO
1049 DO i=istru,iend
1050 fc(i,0)=0.0_r8
1051 cf(i,0)=0.0_r8
1052 END DO
1053 DO k=1,n(ng)-1
1054 DO i=istru,iend
1055 cff=1.0_r8/(2.0_r8*dc(i,k+1)+dc(i,k)*(2.0_r8-fc(i,k-1)))
1056 fc(i,k)=cff*dc(i,k+1)
1057 cf(i,k)=cff*(6.0_r8*(u(i,j,k+1,nrhs)- &
1058 & u(i,j,k ,nrhs))- &
1059 & dc(i,k)*cf(i,k-1))
1060 END DO
1061 END DO
1062 DO i=istru,iend
1063 cf(i,n(ng))=0.0_r8
1064 END DO
1065 DO k=n(ng)-1,1,-1
1066 DO i=istru,iend
1067 cf(i,k)=cf(i,k)-fc(i,k)*cf(i,k+1)
1068 END DO
1069 END DO
1070!
1071! Compute spline-interpolated, vertical advective u-momentum flux.
1072!
1073 cff3=1.0_r8/3.0_r8
1074 cff4=1.0_r8/6.0_r8
1075 DO k=1,n(ng)-1
1076 DO i=istru,iend
1077 fc(i,k)=(cff1*(w(i ,j,k)+ &
1078 & w(i-1,j,k))- &
1079 & cff2*(w(i+1,j,k)+ &
1080 & w(i-2,j,k)))* &
1081 & (u(i,j,k,nrhs)+ &
1082 & dc(i,k)*(cff3*cf(i,k )+ &
1083 & cff4*cf(i,k-1)))
1084 END DO
1085 END DO
1086 DO i=istru,iend
1087 fc(i,n(ng))=0.0_r8
1088 fc(i,0)=0.0_r8
1089 END DO
1090# elif defined UV_C2ADVECTION
1091 DO k=1,n(ng)-1
1092 DO i=istru,iend
1093 fc(i,k)=0.25_r8*(u(i,j,k ,nrhs)+ &
1094 & u(i,j,k+1,nrhs))* &
1095 & (w(i ,j,k)+ &
1096 & w(i-1,j,k))
1097 END DO
1098 END DO
1099 DO i=istru,iend
1100 fc(i,0)=0.0_r8
1101 fc(i,n(ng))=0.0_r8
1102 END DO
1103# elif defined UV_C4ADVECTION
1104 cff1=9.0_r8/32.0_r8
1105 cff2=1.0_r8/32.0_r8
1106 DO k=2,n(ng)-2
1107 DO i=istru,iend
1108 fc(i,k)=(cff1*(u(i,j,k ,nrhs)+ &
1109 & u(i,j,k+1,nrhs))- &
1110 & cff2*(u(i,j,k-1,nrhs)+ &
1111 & u(i,j,k+2,nrhs)))* &
1112 & (w(i ,j,k)+ &
1113 & w(i-1,j,k))
1114 END DO
1115 END DO
1116 DO i=istru,iend
1117 fc(i,n(ng))=0.0_r8
1118 fc(i,n(ng)-1)=(cff1*(u(i,j,n(ng)-1,nrhs)+ &
1119 & u(i,j,n(ng) ,nrhs))- &
1120 & cff2*(u(i,j,n(ng)-2,nrhs)+ &
1121 & u(i,j,n(ng) ,nrhs)))* &
1122 & (w(i ,j,n(ng)-1)+ &
1123 & w(i-1,j,n(ng)-1))
1124 fc(i,1)=(cff1*(u(i,j,1,nrhs)+ &
1125 & u(i,j,2,nrhs))- &
1126 & cff2*(u(i,j,1,nrhs)+ &
1127 & u(i,j,3,nrhs)))* &
1128 & (w(i ,j,1)+ &
1129 & w(i-1,j,1))
1130 fc(i,0)=0.0_r8
1131 END DO
1132# else
1133 cff1=9.0_r8/16.0_r8
1134 cff2=1.0_r8/16.0_r8
1135 DO k=2,n(ng)-2
1136 DO i=istru,iend
1137 fc(i,k)=(cff1*(u(i,j,k ,nrhs)+ &
1138 & u(i,j,k+1,nrhs))- &
1139 & cff2*(u(i,j,k-1,nrhs)+ &
1140 & u(i,j,k+2,nrhs)))* &
1141 & (cff1*(w(i ,j,k)+ &
1142 & w(i-1,j,k))- &
1143 & cff2*(w(i+1,j,k)+ &
1144 & w(i-2,j,k)))
1145 END DO
1146 END DO
1147 DO i=istru,iend
1148 fc(i,n(ng))=0.0_r8
1149 fc(i,n(ng)-1)=(cff1*(u(i,j,n(ng)-1,nrhs)+ &
1150 & u(i,j,n(ng) ,nrhs))- &
1151 & cff2*(u(i,j,n(ng)-2,nrhs)+ &
1152 & u(i,j,n(ng) ,nrhs)))* &
1153 & (cff1*(w(i ,j,n(ng)-1)+ &
1154 & w(i-1,j,n(ng)-1))- &
1155 & cff2*(w(i+1,j,n(ng)-1)+ &
1156 & w(i-2,j,n(ng)-1)))
1157 fc(i,1)=(cff1*(u(i,j,1,nrhs)+ &
1158 & u(i,j,2,nrhs))- &
1159 & cff2*(u(i,j,1,nrhs)+ &
1160 & u(i,j,3,nrhs)))* &
1161 & (cff1*(w(i ,j,1)+ &
1162 & w(i-1,j,1))- &
1163 & cff2*(w(i+1,j,1)+ &
1164 & w(i-2,j,1)))
1165 fc(i,0)=0.0_r8
1166 END DO
1167# endif
1168 DO k=1,n(ng)
1169 DO i=istru,iend
1170 cff=fc(i,k)-fc(i,k-1)
1171 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)-cff
1172# ifdef DIAGNOSTICS_UV
1173 diaru(i,j,k,nrhs,m3vadv)=-cff
1174# endif
1175 END DO
1176 END DO
1177 IF (j.ge.jstrv) THEN
1178# ifdef UV_SADVECTION
1179!
1180! Construct conservative parabolic splines for the vertical
1181! derivatives "CF" of v-momentum.
1182!
1183 cff1=9.0_r8/16.0_r8
1184 cff2=1.0_r8/16.0_r8
1185 DO k=1,n(ng)
1186 DO i=istr,iend
1187 dc(i,k)=(cff1*(hz(i,j ,k)+ &
1188 & hz(i,j-1,k))- &
1189 & cff2*(hz(i,j+1,k)+ &
1190 & hz(i,j-2,k)))
1191 END DO
1192 END DO
1193 DO i=istr,iend
1194 fc(i,0)=0.0_r8
1195 cf(i,0)=0.0_r8
1196 END DO
1197 DO k=1,n(ng)-1
1198 DO i=istr,iend
1199 cff=1.0_r8/(2.0_r8*dc(i,k+1)+dc(i,k)*(2.0_r8-fc(i,k-1)))
1200 fc(i,k)=cff*dc(i,k+1)
1201 cf(i,k)=cff*(6.0_r8*(v(i,j,k+1,nrhs)- &
1202 & v(i,j,k ,nrhs))- &
1203 & dc(i,k)*cf(i,k-1))
1204 END DO
1205 END DO
1206 DO i=istr,iend
1207 cf(i,n(ng))=0.0_r8
1208 END DO
1209 DO k=n(ng)-1,1,-1
1210 DO i=istr,iend
1211 cf(i,k)=cf(i,k)-fc(i,k)*cf(i,k+1)
1212 END DO
1213 END DO
1214!
1215! Compute spline-interpolated, vertical advective v-momentum flux.
1216!
1217 cff3=1.0_r8/3.0_r8
1218 cff4=1.0_r8/6.0_r8
1219 DO k=1,n(ng)-1
1220 DO i=istr,iend
1221 fc(i,k)=(cff1*(w(i,j ,k)+ &
1222 & w(i,j-1,k))- &
1223 & cff2*(w(i,j+1,k)+ &
1224 & w(i,j-2,k)))* &
1225 & (v(i,j,k,nrhs)+ &
1226 & dc(i,k)*(cff3*cf(i,k )+ &
1227 & cff4*cf(i,k-1)))
1228 END DO
1229 END DO
1230 DO i=istr,iend
1231 fc(i,n(ng))=0.0_r8
1232 fc(i,0)=0.0_r8
1233 END DO
1234# elif defined UV_C2ADVECTION
1235!
1236! Second-order, centered differences vertical advection.
1237!
1238 DO k=1,n(ng)-1
1239 DO i=istr,iend
1240 fc(i,k)=0.25_r8*(v(i,j,k ,nrhs)+ &
1241 & v(i,j,k+1,nrhs))* &
1242 & (w(i,j ,k)+ &
1243 & w(i,j-1,k))
1244 END DO
1245 END DO
1246 DO i=istr,iend
1247 fc(i,0)=0.0_r8
1248 fc(i,n(ng))=0.0_r8
1249 END DO
1250# elif defined UV_C4ADVECTION
1251!
1252! Forth-order, centered differences vertical advection.
1253!
1254 cff1=9.0_r8/32.0_r8
1255 cff2=1.0_r8/32.0_r8
1256 DO k=2,n(ng)-2
1257 DO i=istr,iend
1258 fc(i,k)=(cff1*(v(i,j,k ,nrhs)+ &
1259 & v(i,j,k+1,nrhs))- &
1260 & cff2*(v(i,j,k-1,nrhs)+ &
1261 & v(i,j,k+2,nrhs)))* &
1262 & (w(i,j ,k)+ &
1263 & w(i,j-1,k))
1264 END DO
1265 END DO
1266 DO i=istr,iend
1267 fc(i,n(ng))=0.0_r8
1268 fc(i,n(ng)-1)=(cff1*(v(i,j,n(ng)-1,nrhs)+ &
1269 & v(i,j,n(ng) ,nrhs))- &
1270 & cff2*(v(i,j,n(ng)-2,nrhs)+ &
1271 & v(i,j,n(ng) ,nrhs)))* &
1272 & (w(i,j ,n(ng)-1)+ &
1273 & w(i,j-1,n(ng)-1))
1274 fc(i,1)=(cff1*(v(i,j,1,nrhs)+ &
1275 & v(i,j,2,nrhs))- &
1276 & cff2*(v(i,j,1,nrhs)+ &
1277 & v(i,j,3,nrhs)))* &
1278 & (w(i,j ,1)+ &
1279 & w(i,j-1,1))
1280 fc(i,0)=0.0_r8
1281 END DO
1282# else
1283 cff1=9.0_r8/16.0_r8
1284 cff2=1.0_r8/16.0_r8
1285 DO k=2,n(ng)-2
1286 DO i=istr,iend
1287 fc(i,k)=(cff1*(v(i,j,k ,nrhs)+ &
1288 & v(i,j,k+1,nrhs))- &
1289 & cff2*(v(i,j,k-1,nrhs)+ &
1290 & v(i,j,k+2,nrhs)))* &
1291 & (cff1*(w(i,j ,k)+ &
1292 & w(i,j-1,k))- &
1293 & cff2*(w(i,j+1,k)+ &
1294 & w(i,j-2,k)))
1295 END DO
1296 END DO
1297 DO i=istr,iend
1298 fc(i,n(ng))=0.0_r8
1299 fc(i,n(ng)-1)=(cff1*(v(i,j,n(ng)-1,nrhs)+ &
1300 & v(i,j,n(ng) ,nrhs))- &
1301 & cff2*(v(i,j,n(ng)-2,nrhs)+ &
1302 & v(i,j,n(ng) ,nrhs)))* &
1303 & (cff1*(w(i,j ,n(ng)-1)+ &
1304 & w(i,j-1,n(ng)-1))- &
1305 & cff2*(w(i,j+1,n(ng)-1)+ &
1306 & w(i,j-2,n(ng)-1)))
1307 fc(i,1)=(cff1*(v(i,j,1,nrhs)+ &
1308 & v(i,j,2,nrhs))- &
1309 & cff2*(v(i,j,1,nrhs)+ &
1310 & v(i,j,3,nrhs)))* &
1311 & (cff1*(w(i,j ,1)+ &
1312 & w(i,j-1,1))- &
1313 & cff2*(w(i,j+1,1)+ &
1314 & w(i,j-2,1)))
1315 fc(i,0)=0.0_r8
1316 END DO
1317# endif
1318 DO k=1,n(ng)
1319 DO i=istr,iend
1320 cff=fc(i,k)-fc(i,k-1)
1321 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)-cff
1322# ifdef DIAGNOSTICS_UV
1323 diarv(i,j,k,nrhs,m3vadv)=-cff
1324# endif
1325 END DO
1326 END DO
1327 END IF
1328# ifdef WEC_VF
1329!
1330!-----------------------------------------------------------------------
1331! Add in vertical advection from W_stokes.
1332!-----------------------------------------------------------------------
1333!
1334# ifdef UV_SADVECTION
1335!
1336! Construct conservative parabolic splines for the vertical
1337! derivatives "CF" of u-momentum.
1338!
1339 cff1=9.0_r8/16.0_r8
1340 cff2=1.0_r8/16.0_r8
1341 DO k=1,n(ng)
1342 DO i=istru,iend
1343 dc(i,k)=cff1*(hz(i ,j,k)+ &
1344 & hz(i-1,j,k))- &
1345 & cff2*(hz(i+1,j,k)+ &
1346 & hz(i-2,j,k))
1347 END DO
1348 END DO
1349 DO i=istru,iend
1350 fc(i,0)=0.0_r8
1351 cf(i,0)=0.0_r8
1352 END DO
1353 DO k=1,n(ng)-1
1354 DO i=istru,iend
1355 cff=1.0_r8/(2.0_r8*dc(i,k+1)+dc(i,k)*(2.0_r8-fc(i,k-1)))
1356 fc(i,k)=cff*dc(i,k+1)
1357 cf(i,k)=cff*(6.0_r8*(u(i,j,k+1,nrhs)- &
1358 & u(i,j,k ,nrhs))- &
1359 & dc(i,k)*cf(i,k-1))
1360 END DO
1361 END DO
1362 DO i=istru,iend
1363 cf(i,n(ng))=0.0_r8
1364 END DO
1365 DO k=n(ng)-1,1,-1
1366 DO i=istru,iend
1367 cf(i,k)=cf(i,k)-fc(i,k)*cf(i,k+1)
1368 END DO
1369 END DO
1370!
1371! Compute spline-interpolated, vertical advective u-momentum flux.
1372!
1373 cff3=1.0_r8/3.0_r8
1374 cff4=1.0_r8/6.0_r8
1375 DO k=1,n(ng)-1
1376 DO i=istru,iend
1377 fc(i,k)=(cff1*(w_stokes(i ,j,k)+ &
1378 & w_stokes(i-1,j,k))- &
1379 & cff2*(w_stokes(i+1,j,k)+ &
1380 & w_stokes(i-2,j,k)))* &
1381 & (u(i,j,k,nrhs)+ &
1382 & dc(i,k)*(cff3*cf(i,k )+ &
1383 & cff4*cf(i,k-1)))
1384 fcs(i,k)=(cff1*(w_stokes(i ,j,k)+ &
1385 & w_stokes(i-1,j,k))- &
1386 & cff2*(w_stokes(i+1,j,k)+ &
1387 & w_stokes(i-2,j,k)))
1388 END DO
1389 END DO
1390 DO i=istru,iend
1391 fc(i,n(ng))=0.0_r8
1392 fc(i,0)=0.0_r8
1393 fcs(i,n(ng))=0.0_r8
1394 fcs(i,0)=0.0_r8
1395 END DO
1396# elif defined UV_C2ADVECTION
1397 DO k=1,n(ng)-1
1398 DO i=istru,iend
1399 fc(i,k)=0.25_r8*(u(i,j,k ,nrhs)+ &
1400 & u(i,j,k+1,nrhs))* &
1401 & (w_stokes(i ,j,k)+ &
1402 & w_stokes(i-1,j,k))
1403 fcs(i,k)=0.5_r8*(w_stokes(i ,j,k)+ &
1404 & w_stokes(i-1,j,k))
1405 END DO
1406 END DO
1407 DO i=istru,iend
1408 fc(i,0)=0.0_r8
1409 fc(i,n(ng))=0.0_r8
1410 fcs(i,0)=0.0_r8
1411 fcs(i,n(ng))=0.0_r8
1412 END DO
1413# elif defined UV_C4ADVECTION
1414 cff1=9.0_r8/32.0_r8
1415 cff2=1.0_r8/32.0_r8
1416 DO k=2,n(ng)-2
1417 DO i=istru,iend
1418 fc(i,k)=(cff1*(u(i,j,k ,nrhs)+ &
1419 & u(i,j,k+1,nrhs))- &
1420 & cff2*(u(i,j,k-1,nrhs)+ &
1421 & u(i,j,k+2,nrhs)))* &
1422 & (w_stokes(i ,j,k)+ &
1423 & w_stokes(i-1,j,k))
1424 fcs(i,k)=0.5_r8*(w_stokes(i ,j,k)+ &
1425 & w_stokes(i-1,j,k))
1426 END DO
1427 END DO
1428 DO i=istru,iend
1429 fc(i,n(ng))=0.0_r8
1430 fcs(i,n(ng))=0.0_r8
1431 fc(i,n(ng)-1)=(cff1*(u(i,j,n(ng)-1,nrhs)+ &
1432 & u(i,j,n(ng) ,nrhs))- &
1433 & cff2*(u(i,j,n(ng)-2,nrhs)+ &
1434 & u(i,j,n(ng) ,nrhs)))* &
1435 & (w_stokes(i ,j,n(ng)-1)+ &
1436 & w_stokes(i-1,j,n(ng)-1))
1437 fcs(i,n(ng)-1)=0.5_r8*(w_stokes(i ,j,n(ng)-1)+ &
1438 & w_stokes(i-1,j,n(ng)-1))
1439 fc(i,1)=(cff1*(u(i,j,1,nrhs)+ &
1440 & u(i,j,2,nrhs))- &
1441 & cff2*(u(i,j,1,nrhs)+ &
1442 & u(i,j,3,nrhs)))* &
1443 & (w_stokes(i ,j,1)+ &
1444 & w_stokes(i-1,j,1))
1445 fcs(i,1)=0.5_r8*(w_stokes(i ,j,1)+ &
1446 & w_stokes(i-1,j,1))
1447 fc(i,0)=0.0_r8
1448 fcs(i,0)=0.0_r8
1449 END DO
1450# else
1451 cff1=9.0_r8/16.0_r8
1452 cff2=1.0_r8/16.0_r8
1453 DO k=2,n(ng)-2
1454 DO i=istru,iend
1455 fc(i,k)=(cff1*(u(i,j,k ,nrhs)+ &
1456 & u(i,j,k+1,nrhs))- &
1457 & cff2*(u(i,j,k-1,nrhs)+ &
1458 & u(i,j,k+2,nrhs)))* &
1459 & (cff1*(w_stokes(i ,j,k)+ &
1460 & w_stokes(i-1,j,k))- &
1461 & cff2*(w_stokes(i+1,j,k)+ &
1462 & w_stokes(i-2,j,k)))
1463 fcs(i,k)=cff1*(w_stokes(i ,j,k)+ &
1464 & w_stokes(i-1,j,k))- &
1465 & cff2*(w_stokes(i+1,j,k)+ &
1466 & w_stokes(i-2,j,k))
1467 END DO
1468 END DO
1469 DO i=istru,iend
1470 fc(i,n(ng))=0.0_r8
1471 fcs(i,n(ng))=0.0_r8
1472 fc(i,n(ng)-1)=(cff1*(u(i,j,n(ng)-1,nrhs)+ &
1473 & u(i,j,n(ng) ,nrhs))- &
1474 & cff2*(u(i,j,n(ng)-2,nrhs)+ &
1475 & u(i,j,n(ng) ,nrhs)))* &
1476 & (cff1*(w_stokes(i ,j,n(ng)-1)+ &
1477 & w_stokes(i-1,j,n(ng)-1))- &
1478 & cff2*(w_stokes(i+1,j,n(ng)-1)+ &
1479 & w_stokes(i-2,j,n(ng)-1)))
1480 fcs(i,n(ng)-1)=(cff1*(w_stokes(i ,j,n(ng)-1)+ &
1481 & w_stokes(i-1,j,n(ng)-1))- &
1482 & cff2*(w_stokes(i+1,j,n(ng)-1)+ &
1483 & w_stokes(i-2,j,n(ng)-1)))
1484 fc(i,1)=(cff1*(u(i,j,1,nrhs)+ &
1485 & u(i,j,2,nrhs))- &
1486 & cff2*(u(i,j,1,nrhs)+ &
1487 & u(i,j,3,nrhs)))* &
1488 & (cff1*(w_stokes(i ,j,1)+ &
1489 & w_stokes(i-1,j,1))- &
1490 & cff2*(w_stokes(i+1,j,1)+ &
1491 & w_stokes(i-2,j,1)))
1492 fcs(i,1)=(cff1*(w_stokes(i ,j,1)+ &
1493 & w_stokes(i-1,j,1))- &
1494 & cff2*(w_stokes(i+1,j,1)+ &
1495 & w_stokes(i-2,j,1)))
1496 fc(i,0)=0.0_r8
1497 fcs(i,0)=0.0_r8
1498 END DO
1499# endif
1500 DO k=1,n(ng)
1501 DO i=istru,iend
1502 cff=fc(i,k)-fc(i,k-1)
1503 cff1=u(i,j,k,nrhs)*(fcs(i,k)-fcs(i,k-1))
1504 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)-cff
1505# ifdef DIAGNOSTICS_UV
1506 diaru(i,j,k,nrhs,m3vadv)=diaru(i,j,k,nrhs,m3vadv)-cff1
1507 diaru(i,j,k,nrhs,m3vjvf)=-(cff-cff1)
1508# endif
1509 END DO
1510 END DO
1511 IF (j.ge.jstrv) THEN
1512# ifdef UV_SADVECTION
1513!
1514! Construct conservative parabolic splines for the vertical
1515! derivatives "CF" of v-momentum.
1516!
1517 cff1=9.0_r8/16.0_r8
1518 cff2=1.0_r8/16.0_r8
1519 DO k=1,n(ng)
1520 DO i=istr,iend
1521 dc(i,k)=(cff1*(hz(i,j ,k)+ &
1522 & hz(i,j-1,k))- &
1523 & cff2*(hz(i,j+1,k)+ &
1524 & hz(i,j-2,k)))
1525 END DO
1526 END DO
1527 DO i=istr,iend
1528 fc(i,0)=0.0_r8
1529 cf(i,0)=0.0_r8
1530 END DO
1531 DO k=1,n(ng)-1
1532 DO i=istr,iend
1533 cff=1.0_r8/(2.0_r8*dc(i,k+1)+dc(i,k)*(2.0_r8-fc(i,k-1)))
1534 fc(i,k)=cff*dc(i,k+1)
1535 cf(i,k)=cff*(6.0_r8*(v(i,j,k+1,nrhs)- &
1536 & v(i,j,k ,nrhs))- &
1537 & dc(i,k)*cf(i,k-1))
1538 END DO
1539 END DO
1540 DO i=istr,iend
1541 cf(i,n(ng))=0.0_r8
1542 END DO
1543 DO k=n(ng)-1,1,-1
1544 DO i=istr,iend
1545 cf(i,k)=cf(i,k)-fc(i,k)*cf(i,k+1)
1546 END DO
1547 END DO
1548!
1549! Compute spline-interpolated, vertical advective v-momentum flux.
1550!
1551 cff3=1.0_r8/3.0_r8
1552 cff4=1.0_r8/6.0_r8
1553 DO k=1,n(ng)-1
1554 DO i=istr,iend
1555 fc(i,k)=(cff1*(w_stokes(i,j ,k)+ &
1556 & w_stokes(i,j-1,k))- &
1557 & cff2*(w_stokes(i,j+1,k)+ &
1558 & w_stokes(i,j-2,k)))* &
1559 & (v(i,j,k,nrhs)+ &
1560 & dc(i,k)*(cff3*cf(i,k )+ &
1561 & cff4*cf(i,k-1)))
1562 fcs(i,k)=(cff1*(w_stokes(i,j ,k)+ &
1563 & w_stokes(i,j-1,k))- &
1564 & cff2*(w_stokes(i,j+1,k)+ &
1565 & w_stokes(i,j-2,k)))
1566 END DO
1567 END DO
1568 DO i=istr,iend
1569 fc(i,n(ng))=0.0_r8
1570 fc(i,0)=0.0_r8
1571 fcs(i,n(ng))=0.0_r8
1572 fcs(i,0)=0.0_r8
1573 END DO
1574# elif defined UV_C2ADVECTION
1575!
1576! Second-order, centered differences vertical advection.
1577!
1578 DO k=1,n(ng)-1
1579 DO i=istr,iend
1580 fc(i,k)=0.25_r8*(v(i,j,k ,nrhs)+ &
1581 & v(i,j,k+1,nrhs))* &
1582 & (w_stokes(i,j ,k)+ &
1583 & w_stokes(i,j-1,k))
1584 fcs(i,k)=0.5_r8*(w_stokes(i,j ,k)+ &
1585 & w_stokes(i,j-1,k))
1586 END DO
1587 END DO
1588 DO i=istr,iend
1589 fc(i,0)=0.0_r8
1590 fc(i,n(ng))=0.0_r8
1591 fcs(i,0)=0.0_r8
1592 fcs(i,n(ng))=0.0_r8
1593 END DO
1594# elif defined UV_C4ADVECTION
1595!
1596! Forth-order, centered differences vertical advection.
1597!
1598 cff1=9.0_r8/32.0_r8
1599 cff2=1.0_r8/32.0_r8
1600 DO k=2,n(ng)-2
1601 DO i=istr,iend
1602 fc(i,k)=(cff1*(v(i,j,k ,nrhs)+ &
1603 & v(i,j,k+1,nrhs))- &
1604 & cff2*(v(i,j,k-1,nrhs)+ &
1605 & v(i,j,k+2,nrhs)))* &
1606 & (w_stokes(i,j ,k)+ &
1607 & w_stokes(i,j-1,k))
1608 fcs(i,k)=0.5_r8*(w_stokes(i,j ,k)+ &
1609 & w_stokes(i,j-1,k))
1610 END DO
1611 END DO
1612 DO i=istr,iend
1613 fc(i,n(ng))=0.0_r8
1614 fcs(i,n(ng))=0.0_r8
1615 fc(i,n(ng)-1)=(cff1*(v(i,j,n(ng)-1,nrhs)+ &
1616 & v(i,j,n(ng) ,nrhs))- &
1617 & cff2*(v(i,j,n(ng)-2,nrhs)+ &
1618 & v(i,j,n(ng) ,nrhs)))* &
1619 & (w_stokes(i,j ,n(ng)-1)+ &
1620 & w_stokes(i,j-1,n(ng)-1))
1621 fcs(i,n(ng)-1)=0.5_r8*(w_stokes(i,j ,n(ng)-1)+ &
1622 & w_stokes(i,j-1,n(ng)-1))
1623 fc(i,1)=(cff1*(v(i,j,1,nrhs)+ &
1624 & v(i,j,2,nrhs))- &
1625 & cff2*(v(i,j,1,nrhs)+ &
1626 & v(i,j,3,nrhs)))* &
1627 & (w_stokes(i,j ,1)+ &
1628 & w_stokes(i,j-1,1))
1629 fcs(i,1)=0.5_r8*(w_stokes(i,j ,1)+ &
1630 & w_stokes(i,j-1,1))
1631 fc(i,0)=0.0_r8
1632 fcs(i,0)=0.0_r8
1633 END DO
1634# else
1635 cff1=9.0_r8/16.0_r8
1636 cff2=1.0_r8/16.0_r8
1637 DO k=2,n(ng)-2
1638 DO i=istr,iend
1639 fc(i,k)=(cff1*(v(i,j,k ,nrhs)+ &
1640 & v(i,j,k+1,nrhs))- &
1641 & cff2*(v(i,j,k-1,nrhs)+ &
1642 & v(i,j,k+2,nrhs)))* &
1643 & (cff1*(w_stokes(i,j ,k)+ &
1644 & w_stokes(i,j-1,k))- &
1645 & cff2*(w_stokes(i,j+1,k)+ &
1646 & w_stokes(i,j-2,k)))
1647 fcs(i,k)=(cff1*(w_stokes(i,j ,k)+ &
1648 & w_stokes(i,j-1,k))- &
1649 & cff2*(w_stokes(i,j+1,k)+ &
1650 & w_stokes(i,j-2,k)))
1651 END DO
1652 END DO
1653 DO i=istr,iend
1654 fc(i,n(ng))=0.0_r8
1655 fcs(i,n(ng))=0.0_r8
1656 fc(i,n(ng)-1)=(cff1*(v(i,j,n(ng)-1,nrhs)+ &
1657 & v(i,j,n(ng) ,nrhs))- &
1658 & cff2*(v(i,j,n(ng)-2,nrhs)+ &
1659 & v(i,j,n(ng) ,nrhs)))* &
1660 & (cff1*(w_stokes(i,j ,n(ng)-1)+ &
1661 & w_stokes(i,j-1,n(ng)-1))- &
1662 & cff2*(w_stokes(i,j+1,n(ng)-1)+ &
1663 & w_stokes(i,j-2,n(ng)-1)))
1664 fcs(i,n(ng)-1)=(cff1*(w_stokes(i,j ,n(ng)-1)+ &
1665 & w_stokes(i,j-1,n(ng)-1))- &
1666 & cff2*(w_stokes(i,j+1,n(ng)-1)+ &
1667 & w_stokes(i,j-2,n(ng)-1)))
1668 fc(i,1)=(cff1*(v(i,j,1,nrhs)+ &
1669 & v(i,j,2,nrhs))- &
1670 & cff2*(v(i,j,1,nrhs)+ &
1671 & v(i,j,3,nrhs)))* &
1672 & (cff1*(w_stokes(i,j ,1)+ &
1673 & w_stokes(i,j-1,1))- &
1674 & cff2*(w_stokes(i,j+1,1)+ &
1675 & w_stokes(i,j-2,1)))
1676 fcs(i,1)=(cff1*(w_stokes(i,j ,1)+ &
1677 & w_stokes(i,j-1,1))- &
1678 & cff2*(w_stokes(i,j+1,1)+ &
1679 & w_stokes(i,j-2,1)))
1680 fc(i,0)=0.0_r8
1681 fcs(i,0)=0.0_r8
1682 END DO
1683# endif
1684 DO k=1,n(ng)
1685 DO i=istr,iend
1686 cff=fc(i,k)-fc(i,k-1)
1687 cff1=v(i,j,k,nrhs)*(fcs(i,k)-fcs(i,k-1))
1688 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)-cff
1689# ifdef DIAGNOSTICS_UV
1690 diarv(i,j,k,nrhs,m3vadv)=diarv(i,j,k,nrhs,m3vadv)-cff1
1691 diarv(i,j,k,nrhs,m3vjvf)=-(cff-cff1)
1692# endif
1693 END DO
1694 END DO
1695 END IF
1696# endif
1697# endif
1698!
1699!-----------------------------------------------------------------------
1700! Compute forcing term for the 2D momentum equations.
1701!-----------------------------------------------------------------------
1702!
1703! Vertically integrate baroclinic right-hand-side terms. If not
1704! body force stresses, add in the difference between surface and
1705! bottom stresses.
1706!
1707 DO i=istru,iend
1708# ifdef WET_DRY
1709 ru(i,j,1,nrhs)=ru(i,j,1,nrhs)*umask_wet(i,j)
1710# endif
1711 rufrc(i,j)=ru(i,j,1,nrhs)
1712# ifdef DIAGNOSTICS_UV
1713 diarufrc(i,j,3,m2pgrd)=diaru(i,j,1,nrhs,m3pgrd)
1714# ifdef UV_COR
1715 diarufrc(i,j,3,m2fcor)=diaru(i,j,1,nrhs,m3fcor)
1716# endif
1717# ifdef UV_ADV
1718 diarufrc(i,j,3,m2xadv)=diaru(i,j,1,nrhs,m3xadv)
1719 diarufrc(i,j,3,m2yadv)=diaru(i,j,1,nrhs,m3yadv)
1720 diarufrc(i,j,3,m2hadv)=diaru(i,j,1,nrhs,m3hadv)
1721# endif
1722# ifdef WEC_VF
1723 diarufrc(i,j,3,m2hjvf)=diaru(i,j,1,nrhs,m3hjvf)
1724 diarufrc(i,j,3,m2kvrf)=diaru(i,j,1,nrhs,m3kvrf)
1725# ifdef UV_COR
1726 diarufrc(i,j,3,m2fsco)=diaru(i,j,1,nrhs,m3fsco)
1727# endif
1728# ifdef BOTTOM_STREAMING
1729 diarufrc(i,j,3,m2bstm)=diaru(i,j,1,nrhs,m3bstm)
1730# endif
1731# ifdef SURFACE_STREAMING
1732 diarufrc(i,j,3,m2sstm)=diaru(i,j,1,nrhs,m3sstm)
1733# endif
1734 diarufrc(i,j,3,m2wrol)=diaru(i,j,1,nrhs,m3wrol)
1735 diarufrc(i,j,3,m2wbrk)=diaru(i,j,1,nrhs,m3wbrk)
1736# endif
1737# if defined UV_VIS2 || defined UV_VIS4
1738 diarufrc(i,j,3,m2xvis)=0.0_r8
1739 diarufrc(i,j,3,m2yvis)=0.0_r8
1740 diarufrc(i,j,3,m2hvis)=0.0_r8
1741# endif
1742# ifdef BODYFORCE
1743!! DiaRUfrc(i,j,3,M2strs)=DiaRU(i,j,1,nrhs,M3vvis)
1744# endif
1745# endif
1746 END DO
1747 DO k=2,n(ng)
1748 DO i=istru,iend
1749# ifdef WET_DRY
1750 ru(i,j,k,nrhs)=ru(i,j,k,nrhs)*umask_wet(i,j)
1751# endif
1752 rufrc(i,j)=rufrc(i,j)+ru(i,j,k,nrhs)
1753# ifdef DIAGNOSTICS_UV
1754 diarufrc(i,j,3,m2pgrd)=diarufrc(i,j,3,m2pgrd)+ &
1755 & diaru(i,j,k,nrhs,m3pgrd)
1756# ifdef UV_COR
1757 diarufrc(i,j,3,m2fcor)=diarufrc(i,j,3,m2fcor)+ &
1758 & diaru(i,j,k,nrhs,m3fcor)
1759# endif
1760# ifdef UV_ADV
1761 diarufrc(i,j,3,m2xadv)=diarufrc(i,j,3,m2xadv)+ &
1762 & diaru(i,j,k,nrhs,m3xadv)
1763 diarufrc(i,j,3,m2yadv)=diarufrc(i,j,3,m2yadv)+ &
1764 & diaru(i,j,k,nrhs,m3yadv)
1765 diarufrc(i,j,3,m2hadv)=diarufrc(i,j,3,m2hadv)+ &
1766 & diaru(i,j,k,nrhs,m3hadv)
1767# endif
1768# ifdef WEC_VF
1769 diarufrc(i,j,3,m2hjvf)=diarufrc(i,j,3,m2hjvf)+ &
1770 & diaru(i,j,k,nrhs,m3hjvf)
1771 diarufrc(i,j,3,m2kvrf)=diarufrc(i,j,3,m2kvrf)+ &
1772 & diaru(i,j,k,nrhs,m3kvrf)
1773# ifdef UV_COR
1774 diarufrc(i,j,3,m2fsco)=diarufrc(i,j,3,m2fsco)+ &
1775 & diaru(i,j,k,nrhs,m3fsco)
1776# endif
1777# ifdef BOTTOM_STREAMING
1778 diarufrc(i,j,3,m2bstm)=diarufrc(i,j,3,m2bstm)+ &
1779 & diaru(i,j,k,nrhs,m3bstm)
1780# endif
1781# ifdef SURFACE_STREAMING
1782 diarufrc(i,j,3,m2sstm)=diarufrc(i,j,3,m2sstm)+ &
1783 & diaru(i,j,k,nrhs,m3sstm)
1784# endif
1785 diarufrc(i,j,3,m2wrol)=diarufrc(i,j,3,m2wrol)+ &
1786 & diaru(i,j,k,nrhs,m3wrol)
1787 diarufrc(i,j,3,m2wbrk)=diarufrc(i,j,3,m2wbrk)+ &
1788 & diaru(i,j,k,nrhs,m3wbrk)
1789# endif
1790# ifdef BODYFORCE
1791!! DiaRUfrc(i,j,3,M2strs)=DiaRUfrc(i,j,3,M2strs)+ &
1792!! & DiaRU(i,j,k,nrhs,M3vvis)
1793# endif
1794# endif
1795 END DO
1796 END DO
1797# ifndef BODYFORCE
1798 DO i=istru,iend
1799 cff=om_u(i,j)*on_u(i,j)
1800 cff1= sustr(i,j)*cff
1801 cff2=-bustr(i,j)*cff
1802 rufrc(i,j)=rufrc(i,j)+cff1+cff2
1803# ifdef WET_DRY
1804 rufrc(i,j)=rufrc(i,j)*umask_wet(i,j)
1805# endif
1806# ifdef DIAGNOSTICS_UV
1807 diarufrc(i,j,3,m2sstr)=cff1
1808 diarufrc(i,j,3,m2bstr)=cff2
1809# endif
1810 END DO
1811# endif
1812 IF (j.ge.jstrv) THEN
1813 DO i=istr,iend
1814# ifdef WET_DRY
1815 rv(i,j,1,nrhs)=rv(i,j,1,nrhs)*vmask_wet(i,j)
1816# endif
1817 rvfrc(i,j)=rv(i,j,1,nrhs)
1818# ifdef DIAGNOSTICS_UV
1819 diarvfrc(i,j,3,m2pgrd)=diarv(i,j,1,nrhs,m3pgrd)
1820# ifdef UV_COR
1821 diarvfrc(i,j,3,m2fcor)=diarv(i,j,1,nrhs,m3fcor)
1822# endif
1823# ifdef UV_ADV
1824 diarvfrc(i,j,3,m2xadv)=diarv(i,j,1,nrhs,m3xadv)
1825 diarvfrc(i,j,3,m2yadv)=diarv(i,j,1,nrhs,m3yadv)
1826 diarvfrc(i,j,3,m2hadv)=diarv(i,j,1,nrhs,m3hadv)
1827# endif
1828# ifdef WEC_VF
1829 diarvfrc(i,j,3,m2hjvf)=diarv(i,j,1,nrhs,m3hjvf)
1830 diarvfrc(i,j,3,m2kvrf)=diarv(i,j,1,nrhs,m3kvrf)
1831# ifdef UV_COR
1832 diarvfrc(i,j,3,m2fsco)=diarv(i,j,1,nrhs,m3fsco)
1833# endif
1834# ifdef BOTTOM_STREAMING
1835 diarvfrc(i,j,3,m2bstm)=diarv(i,j,1,nrhs,m3bstm)
1836# endif
1837# ifdef SURFACE_STREAMING
1838 diarvfrc(i,j,3,m2sstm)=diarv(i,j,1,nrhs,m3sstm)
1839# endif
1840 diarvfrc(i,j,3,m2wrol)=diarv(i,j,1,nrhs,m3wrol)
1841 diarvfrc(i,j,3,m2wbrk)=diarv(i,j,1,nrhs,m3wbrk)
1842# endif
1843# if defined UV_VIS2 || defined UV_VIS4
1844 diarvfrc(i,j,3,m2hvis)=0.0_r8
1845 diarvfrc(i,j,3,m2xvis)=0.0_r8
1846 diarvfrc(i,j,3,m2yvis)=0.0_r8
1847# endif
1848# ifdef BODYFORCE
1849! DiaRVfrc(i,j,3,M2strs)=DiaRV(i,j,1,nrhs,M3vvis)
1850# endif
1851# endif
1852 END DO
1853 DO k=2,n(ng)
1854 DO i=istr,iend
1855# ifdef WET_DRY
1856 rv(i,j,k,nrhs)=rv(i,j,k,nrhs)*vmask_wet(i,j)
1857# endif
1858 rvfrc(i,j)=rvfrc(i,j)+rv(i,j,k,nrhs)
1859# ifdef DIAGNOSTICS_UV
1860 diarvfrc(i,j,3,m2pgrd)=diarvfrc(i,j,3,m2pgrd)+ &
1861 & diarv(i,j,k,nrhs,m3pgrd)
1862# ifdef UV_COR
1863 diarvfrc(i,j,3,m2fcor)=diarvfrc(i,j,3,m2fcor)+ &
1864 & diarv(i,j,k,nrhs,m3fcor)
1865# endif
1866# ifdef UV_ADV
1867 diarvfrc(i,j,3,m2xadv)=diarvfrc(i,j,3,m2xadv)+ &
1868 & diarv(i,j,k,nrhs,m3xadv)
1869 diarvfrc(i,j,3,m2yadv)=diarvfrc(i,j,3,m2yadv)+ &
1870 & diarv(i,j,k,nrhs,m3yadv)
1871 diarvfrc(i,j,3,m2hadv)=diarvfrc(i,j,3,m2hadv)+ &
1872 & diarv(i,j,k,nrhs,m3hadv)
1873# endif
1874# ifdef WEC_VF
1875 diarvfrc(i,j,3,m2hjvf)=diarvfrc(i,j,3,m2hjvf)+ &
1876 & diarv(i,j,k,nrhs,m3hjvf)
1877 diarvfrc(i,j,3,m2kvrf)=diarvfrc(i,j,3,m2kvrf)+ &
1878 & diarv(i,j,k,nrhs,m3kvrf)
1879# ifdef UV_COR
1880 diarvfrc(i,j,3,m2fsco)=diarvfrc(i,j,3,m2fsco)+ &
1881 & diarv(i,j,k,nrhs,m3fsco)
1882# endif
1883# ifdef BOTTOM_STREAMING
1884 diarvfrc(i,j,3,m2bstm)=diarvfrc(i,j,3,m2bstm)+ &
1885 & diarv(i,j,k,nrhs,m3bstm)
1886# endif
1887# ifdef SURFACE_STREAMING
1888 diarvfrc(i,j,3,m2sstm)=diarvfrc(i,j,3,m2sstm)+ &
1889 & diarv(i,j,k,nrhs,m3sstm)
1890# endif
1891 diarvfrc(i,j,3,m2wrol)=diarvfrc(i,j,3,m2wrol)+ &
1892 & diarv(i,j,k,nrhs,m3wrol)
1893 diarvfrc(i,j,3,m2wbrk)=diarvfrc(i,j,3,m2wbrk)+ &
1894 & diarv(i,j,k,nrhs,m3wbrk)
1895# endif
1896# ifdef BODYFORCE
1897!! DiaRVfrc(i,j,3,M2strs)=DiaRVfrc(i,j,3,M2strs)+ &
1898!! & DiaRV(i,j,k,nrhs,M3vvis)
1899# endif
1900# endif
1901 END DO
1902 END DO
1903# ifndef BODYFORCE
1904 DO i=istr,iend
1905 cff=om_v(i,j)*on_v(i,j)
1906 cff1= svstr(i,j)*cff
1907 cff2=-bvstr(i,j)*cff
1908 rvfrc(i,j)=rvfrc(i,j)+cff1+cff2
1909# ifdef WET_DRY
1910 rvfrc(i,j)=rvfrc(i,j)*vmask_wet(i,j)
1911# endif
1912# ifdef DIAGNOSTICS_UV
1913 diarvfrc(i,j,3,m2sstr)=cff1
1914 diarvfrc(i,j,3,m2bstr)=cff2
1915# endif
1916 END DO
1917# endif
1918 END IF
1919 END DO j_loop
1920!
1921 RETURN
1922 END SUBROUTINE rhs3d_tile
1923#endif
1924 END MODULE rhs3d_mod
mod_clima
Definition mod_clima.F:2
mod_clima::clima
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
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_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
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_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_scalars
Definition mod_scalars.F:2
mod_scalars::m3vvis
integer m3vvis
Definition mod_scalars.F:191
mod_scalars::m2fcor
integer m2fcor
Definition mod_scalars.F:150
mod_scalars::m3hadv
integer m3hadv
Definition mod_scalars.F:177
mod_scalars::m3xadv
integer m3xadv
Definition mod_scalars.F:178
mod_scalars::levbfrc
integer, dimension(:), allocatable levbfrc
Definition mod_scalars.F:587
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::lnudgem3clm
logical, dimension(:), allocatable lnudgem3clm
Definition mod_scalars.F:534
mod_scalars::m2hvis
integer m2hvis
Definition mod_scalars.F:168
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::m3yadv
integer m3yadv
Definition mod_scalars.F:179
mod_scalars::m2yvis
integer m2yvis
Definition mod_scalars.F:170
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::levsfrc
integer, dimension(:), allocatable levsfrc
Definition mod_scalars.F:586
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_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
pre_step3d_mod
Definition pre_step3d.F:2
pre_step3d_mod::pre_step3d
subroutine, public pre_step3d(ng, tile)
Definition pre_step3d.F:40
prsgrd_mod
Definition prsgrd31.h:1
prsgrd_mod::prsgrd
subroutine, public prsgrd(ng, tile)
Definition prsgrd31.h:36
rhs3d_mod
Definition rhs3d.F:2
rhs3d_mod::rhs3d_tile
subroutine rhs3d_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, hz, huon, hvom, dmde, dndx, fomn, om_u, om_v, on_u, on_v, pm, pn, umask_wet, vmask_wet, bustr, bvstr, sustr, svstr, u, v, w, rufrc, rvfrc, diarufrc, diarvfrc, diaru, diarv, ru, rv)
Definition rhs3d.F:222
rhs3d_mod::rhs3d
subroutine, public rhs3d(ng, tile)
Definition rhs3d.F:26
t3dmix2_mod
Definition t3dmix2_geo.h:1
t3dmix2_mod::t3dmix2
subroutine, public t3dmix2(ng, tile)
Definition t3dmix2_geo.h:24
t3dmix4_mod
Definition t3dmix4_geo.h:1
t3dmix4_mod::t3dmix4
subroutine, public t3dmix4(ng, tile)
Definition t3dmix4_geo.h:24
uv3dmix2_mod
Definition uv3dmix2_geo.h:1
uv3dmix2_mod::uv3dmix2
subroutine, public uv3dmix2(ng, tile)
Definition uv3dmix2_geo.h:43
uv3dmix4_mod
Definition uv3dmix4_geo.h:1
uv3dmix4_mod::uv3dmix4
subroutine, public uv3dmix4(ng, tile)
Definition uv3dmix4_geo.h:43
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