ROMS
Loading...
Searching...
No Matches
step2d_LF_AM3.h
Go to the documentation of this file.
1 MODULE step2d_mod
2!
3!git $Id$
4!=======================================================================
5! !
6! Nonlinear shallow-water primitive equations predictor (Leap-frog) !
7! and corrector (Adams-Moulton) time-stepping engine. !
8! !
9!=======================================================================
10!
11 implicit none
12!
13 PRIVATE
14 PUBLIC :: step2d
15!
16 CONTAINS
17!
18 SUBROUTINE step2d (ng, tile)
19
20 USE mod_param
21#ifdef SOLVE3D
22 USE mod_coupling
23#endif
24#ifdef DIAGNOSTICS_UV
25 USE mod_diags
26#endif
27 USE mod_forces
28 USE mod_grid
29#if defined UV_VIS2 || defined UV_VIS4
30 USE mod_mixing
31#endif
32 USE mod_ocean
33#if defined SEDIMENT && defined SED_MORPH && defined SOLVE3D
34 USE mod_sedbed
35#endif
36 USE mod_stepping
37!
38! Imported variable declarations.
39!
40 integer, intent(in) :: ng, tile
41!
42! Local variable declarations.
43!
44 character (len=*), parameter :: MyFile = &
45 & __FILE__
46!
47#include "tile.h"
48!
49#ifdef PROFILE
50 CALL wclock_on (ng, inlm, 9, __line__, myfile)
51#endif
52 CALL step2d_tile (ng, tile, &
53 & lbi, ubi, lbj, ubj, n(ng), &
54 & imins, imaxs, jmins, jmaxs, &
55 & krhs(ng), kstp(ng), knew(ng), &
56#ifdef SOLVE3D
57 & nstp(ng), nnew(ng), &
58#endif
59#ifdef MASKING
60 & grid(ng) % pmask, grid(ng) % rmask, &
61 & grid(ng) % umask, grid(ng) % vmask, &
62#endif
63#ifdef WET_DRY
64 & grid(ng) % pmask_wet, grid(ng) % pmask_full, &
65 & grid(ng) % rmask_wet, grid(ng) % rmask_full, &
66 & grid(ng) % umask_wet, grid(ng) % umask_full, &
67 & grid(ng) % vmask_wet, grid(ng) % vmask_full, &
68# ifdef SOLVE3D
69 & grid(ng) % rmask_wet_avg, &
70# endif
71#endif
72 & grid(ng) % fomn, grid(ng) % h, &
73 & grid(ng) % om_u, grid(ng) % om_v, &
74 & grid(ng) % on_u, grid(ng) % on_v, &
75 & grid(ng) % omn, &
76 & grid(ng) % pm, grid(ng) % pn, &
77#if defined CURVGRID && defined UV_ADV
78 & grid(ng) % dndx, grid(ng) % dmde, &
79#endif
80#if defined UV_VIS2 || defined UV_VIS4
81 & grid(ng) % pmon_r, grid(ng) % pnom_r, &
82 & grid(ng) % pmon_p, grid(ng) % pnom_p, &
83 & grid(ng) % om_r, grid(ng) % on_r, &
84 & grid(ng) % om_p, grid(ng) % on_p, &
85# ifdef UV_VIS2
86 & mixing(ng) % visc2_p, mixing(ng) % visc2_r, &
87# endif
88# ifdef UV_VIS4
89 & mixing(ng) % visc4_p, mixing(ng) % visc4_r, &
90# endif
91#endif
92#if defined SEDIMENT && defined SED_MORPH
93 & sedbed(ng) % bed_thick, &
94#endif
95
96#ifdef WEC
97# ifdef WEC_VF
98# ifdef WEC_ROLLER
99 & mixing(ng) % rurol2d, &
100 & mixing(ng) % rvrol2d, &
101# endif
102# ifdef BOTTOM_STREAMING
103 & mixing(ng) % rubst2d, &
104 & mixing(ng) % rvbst2d, &
105# endif
106# ifdef SURFACE_STREAMING
107 & mixing(ng) % russt2d, &
108 & mixing(ng) % rvsst2d, &
109# endif
110 & mixing(ng) % rubrk2d, &
111 & mixing(ng) % rvbrk2d, &
112 & mixing(ng) % rukvf2d, &
113 & mixing(ng) % rvkvf2d, &
114 & ocean(ng) % bh, &
115 & ocean(ng) % qsp, &
116 & ocean(ng) % zetaw, &
117# endif
118 & ocean(ng) % ubar_stokes, &
119 & ocean(ng) % vbar_stokes, &
120#endif
121#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
122 & ocean(ng) % eq_tide, &
123#endif
124#ifndef SOLVE3D
125 & forces(ng) % sustr, forces(ng) % svstr, &
126 & forces(ng) % bustr, forces(ng) % bvstr, &
127# ifdef ATM_PRESS
128 & forces(ng) % Pair, &
129# endif
130#else
131# ifdef VAR_RHO_2D
132 & coupling(ng) % rhoA, coupling(ng) % rhoS, &
133# endif
134 & coupling(ng) % DU_avg1, coupling(ng) % DU_avg2, &
135 & coupling(ng) % DV_avg1, coupling(ng) % DV_avg2, &
136 & coupling(ng) % Zt_avg1, &
137 & coupling(ng) % rufrc, coupling(ng) % rvfrc, &
138 & ocean(ng) % ru, ocean(ng) % rv, &
139#endif
140#ifdef DIAGNOSTICS_UV
141 & diags(ng) % DiaU2wrk, diags(ng) % DiaV2wrk, &
142 & diags(ng) % DiaRUbar, diags(ng) % DiaRVbar, &
143# ifdef SOLVE3D
144 & diags(ng) % DiaU2int, diags(ng) % DiaV2int, &
145 & diags(ng) % DiaRUfrc, diags(ng) % DiaRVfrc, &
146# endif
147#endif
148#if defined NESTING && !defined SOLVE3D
149 & ocean(ng) % DU_flux, ocean(ng) % DV_flux, &
150#endif
151 & ocean(ng) % rubar, ocean(ng) % rvbar, &
152 & ocean(ng) % rzeta, &
153 & ocean(ng) % ubar, ocean(ng) % vbar, &
154 & ocean(ng) % zeta)
155#ifdef PROFILE
156 CALL wclock_off (ng, inlm, 9, __line__, myfile)
157#endif
158!
159 RETURN
160 END SUBROUTINE step2d
161!
162!***********************************************************************
163 SUBROUTINE step2d_tile (ng, tile, &
164 & LBi, UBi, LBj, UBj, UBk, &
165 & IminS, ImaxS, JminS, JmaxS, &
166 & krhs, kstp, knew, &
167#ifdef SOLVE3D
168 & nstp, nnew, &
169#endif
170#ifdef MASKING
171 & pmask, rmask, umask, vmask, &
172#endif
173#ifdef WET_DRY
174 & pmask_wet, pmask_full, &
175 & rmask_wet, rmask_full, &
176 & umask_wet, umask_full, &
177 & vmask_wet, vmask_full, &
178# ifdef SOLVE3D
179 & rmask_wet_avg, &
180# endif
181#endif
182 & fomn, h, &
183 & om_u, om_v, on_u, on_v, omn, pm, pn, &
184#if defined CURVGRID && defined UV_ADV
185 & dndx, dmde, &
186#endif
187#if defined UV_VIS2 || defined UV_VIS4
188 & pmon_r, pnom_r, pmon_p, pnom_p, &
189 & om_r, on_r, om_p, on_p, &
190# ifdef UV_VIS2
191 & visc2_p, visc2_r, &
192# endif
193# ifdef UV_VIS4
194 & visc4_p, visc4_r, &
195# endif
196#endif
197#if defined SEDIMENT && defined SED_MORPH
198 & bed_thick, &
199#endif
200#ifdef WEC
201# ifdef WEC_VF
202# ifdef WEC_ROLLER
203 & rurol2d, rvrol2d, &
204# endif
205# ifdef BOTTOM_STREAMING
206 & rubst2d, rvbst2d, &
207# endif
208# ifdef SURFACE_STREAMING
209 & russt2d, rvsst2d, &
210# endif
211 & rubrk2d, rvbrk2d, &
212 & rukvf2d, rvkvf2d, &
213 & bh, qsp, zetaw, &
214# endif
215 & ubar_stokes, vbar_stokes, &
216#endif
217#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
218 & eq_tide, &
219#endif
220#ifndef SOLVE3D
221 & sustr, svstr, bustr, bvstr, &
222# ifdef ATM_PRESS
223 & Pair, &
224# endif
225#else
226# ifdef VAR_RHO_2D
227 & rhoA, rhoS, &
228# endif
229 & DU_avg1, DU_avg2, &
230 & DV_avg1, DV_avg2, &
231 & Zt_avg1, &
232 & rufrc, rvfrc, ru, rv, &
233#endif
234#ifdef DIAGNOSTICS_UV
235 & DiaU2wrk, DiaV2wrk, &
236 & DiaRUbar, DiaRVbar, &
237# ifdef SOLVE3D
238 & DiaU2int, DiaV2int, &
239 & DiaRUfrc, DiaRVfrc, &
240# endif
241#endif
242#if defined NESTING && !defined SOLVE3D
243 & DU_flux, DV_flux, &
244#endif
245 & rubar, rvbar, rzeta, &
246 & ubar, vbar, zeta)
247!***********************************************************************
248!
249 USE mod_param
250 USE mod_clima
251 USE mod_ncparam
252 USE mod_scalars
253#if defined SEDIMENT && defined SED_MORPH
254 USE mod_sedbed
255#endif
256 USE mod_sources
257!
258 USE exchange_2d_mod
259#ifdef DISTRIBUTE
260 USE mp_exchange_mod, ONLY : mp_exchange2d
261#endif
262 USE obc_volcons_mod, ONLY : obc_flux_tile, set_duv_bc_tile
263 USE u2dbc_mod, ONLY : u2dbc_tile
264 USE v2dbc_mod, ONLY : v2dbc_tile
265 USE zetabc_mod, ONLY : zetabc_tile
266#ifdef WET_DRY
267 USE wetdry_mod, ONLY : wetdry_tile
268#endif
269!
270! Imported variable declarations.
271!
272 integer, intent(in ) :: ng, tile
273 integer, intent(in ) :: LBi, UBi, LBj, UBj, UBk
274 integer, intent(in ) :: IminS, ImaxS, JminS, JmaxS
275 integer, intent(in ) :: krhs, kstp, knew
276#ifdef SOLVE3D
277 integer, intent(in ) :: nstp, nnew
278#endif
279!
280#ifdef ASSUMED_SHAPE
281# ifdef MASKING
282 real(r8), intent(in ) :: pmask(LBi:,LBj:)
283 real(r8), intent(in ) :: rmask(LBi:,LBj:)
284 real(r8), intent(in ) :: umask(LBi:,LBj:)
285 real(r8), intent(in ) :: vmask(LBi:,LBj:)
286# endif
287 real(r8), intent(in ) :: fomn(LBi:,LBj:)
288# if defined SEDIMENT && defined SED_MORPH
289 real(r8), intent(inout) :: h(LBi:,LBj:)
290# else
291 real(r8), intent(in ) :: h(LBi:,LBj:)
292# endif
293 real(r8), intent(in ) :: om_u(LBi:,LBj:)
294 real(r8), intent(in ) :: om_v(LBi:,LBj:)
295 real(r8), intent(in ) :: on_u(LBi:,LBj:)
296 real(r8), intent(in ) :: on_v(LBi:,LBj:)
297 real(r8), intent(in ) :: omn(LBi:,LBj:)
298 real(r8), intent(in ) :: pm(LBi:,LBj:)
299 real(r8), intent(in ) :: pn(LBi:,LBj:)
300# if defined CURVGRID && defined UV_ADV
301 real(r8), intent(in ) :: dndx(LBi:,LBj:)
302 real(r8), intent(in ) :: dmde(LBi:,LBj:)
303# endif
304# if defined UV_VIS2 || defined UV_VIS4
305 real(r8), intent(in ) :: pmon_r(LBi:,LBj:)
306 real(r8), intent(in ) :: pnom_r(LBi:,LBj:)
307 real(r8), intent(in ) :: pmon_p(LBi:,LBj:)
308 real(r8), intent(in ) :: pnom_p(LBi:,LBj:)
309 real(r8), intent(in ) :: om_r(LBi:,LBj:)
310 real(r8), intent(in ) :: on_r(LBi:,LBj:)
311 real(r8), intent(in ) :: om_p(LBi:,LBj:)
312 real(r8), intent(in ) :: on_p(LBi:,LBj:)
313# ifdef UV_VIS2
314 real(r8), intent(in ) :: visc2_p(LBi:,LBj:)
315 real(r8), intent(in ) :: visc2_r(LBi:,LBj:)
316# endif
317# ifdef UV_VIS4
318 real(r8), intent(in ) :: visc4_p(LBi:,LBj:)
319 real(r8), intent(in ) :: visc4_r(LBi:,LBj:)
320# endif
321# endif
322# if defined SEDIMENT && defined SED_MORPH
323 real(r8), intent(in ) :: bed_thick(LBi:,LBj:,:)
324# endif
325# ifdef WEC
326# ifdef WEC_VF
327# ifdef WEC_ROLLER
328 real(r8), intent(in) :: rurol2d(LBi:,LBj:)
329 real(r8), intent(in) :: rvrol2d(LBi:,LBj:)
330# endif
331# ifdef BOTTOM_STREAMING
332 real(r8), intent(in) :: rubst2d(LBi:,LBj:)
333 real(r8), intent(in) :: rvbst2d(LBi:,LBj:)
334# endif
335# ifdef SURFACE_STREAMING
336 real(r8), intent(in) :: russt2d(LBi:,LBj:)
337 real(r8), intent(in) :: rvsst2d(LBi:,LBj:)
338# endif
339 real(r8), intent(in) :: rubrk2d(LBi:,LBj:)
340 real(r8), intent(in) :: rvbrk2d(LBi:,LBj:)
341 real(r8), intent(in) :: rukvf2d(LBi:,LBj:)
342 real(r8), intent(in) :: rvkvf2d(LBi:,LBj:)
343 real(r8), intent(in) :: bh(LBi:,LBj:)
344 real(r8), intent(in) :: qsp(LBi:,LBj:)
345 real(r8), intent(in) :: zetaw(LBi:,LBj:)
346# endif
347 real(r8), intent(in) :: ubar_stokes(LBi:,LBj:)
348 real(r8), intent(in) :: vbar_stokes(LBi:,LBj:)
349# endif
350# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
351 real(r8), intent(in ) :: eq_tide(LBi:,LBj:)
352# endif
353# ifndef SOLVE3D
354 real(r8), intent(in ) :: sustr(LBi:,LBj:)
355 real(r8), intent(in ) :: svstr(LBi:,LBj:)
356 real(r8), intent(in ) :: bustr(LBi:,LBj:)
357 real(r8), intent(in ) :: bvstr(LBi:,LBj:)
358# ifdef ATM_PRESS
359 real(r8), intent(in ) :: Pair(LBi:,LBj:)
360# endif
361# else
362# ifdef VAR_RHO_2D
363 real(r8), intent(in ) :: rhoA(LBi:,LBj:)
364 real(r8), intent(in ) :: rhoS(LBi:,LBj:)
365# endif
366 real(r8), intent(inout) :: DU_avg1(LBi:,LBj:)
367 real(r8), intent(inout) :: DU_avg2(LBi:,LBj:)
368 real(r8), intent(inout) :: DV_avg1(LBi:,LBj:)
369 real(r8), intent(inout) :: DV_avg2(LBi:,LBj:)
370 real(r8), intent(inout) :: Zt_avg1(LBi:,LBj:)
371 real(r8), intent(inout) :: rufrc(LBi:,LBj:)
372 real(r8), intent(inout) :: rvfrc(LBi:,LBj:)
373 real(r8), intent(inout) :: ru(LBi:,LBj:,0:,:)
374 real(r8), intent(inout) :: rv(LBi:,LBj:,0:,:)
375# endif
376# ifdef WET_DRY
377 real(r8), intent(inout) :: pmask_full(LBi:,LBj:)
378 real(r8), intent(inout) :: rmask_full(LBi:,LBj:)
379 real(r8), intent(inout) :: umask_full(LBi:,LBj:)
380 real(r8), intent(inout) :: vmask_full(LBi:,LBj:)
381
382 real(r8), intent(inout) :: pmask_wet(LBi:,LBj:)
383 real(r8), intent(inout) :: rmask_wet(LBi:,LBj:)
384 real(r8), intent(inout) :: umask_wet(LBi:,LBj:)
385 real(r8), intent(inout) :: vmask_wet(LBi:,LBj:)
386# ifdef SOLVE3D
387 real(r8), intent(inout) :: rmask_wet_avg(LBi:,LBj:)
388# endif
389# endif
390# ifdef DIAGNOSTICS_UV
391 real(r8), intent(inout) :: DiaU2wrk(LBi:,LBj:,:)
392 real(r8), intent(inout) :: DiaV2wrk(LBi:,LBj:,:)
393 real(r8), intent(inout) :: DiaRUbar(LBi:,LBj:,:,:)
394 real(r8), intent(inout) :: DiaRVbar(LBi:,LBj:,:,:)
395# ifdef SOLVE3D
396 real(r8), intent(inout) :: DiaU2int(LBi:,LBj:,:)
397 real(r8), intent(inout) :: DiaV2int(LBi:,LBj:,:)
398 real(r8), intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
399 real(r8), intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
400# endif
401# endif
402 real(r8), intent(inout) :: rubar(LBi:,LBj:,:)
403 real(r8), intent(inout) :: rvbar(LBi:,LBj:,:)
404 real(r8), intent(inout) :: rzeta(LBi:,LBj:,:)
405 real(r8), intent(inout) :: ubar(LBi:,LBj:,:)
406 real(r8), intent(inout) :: vbar(LBi:,LBj:,:)
407 real(r8), intent(inout) :: zeta(LBi:,LBj:,:)
408# if defined NESTING && !defined SOLVE3D
409 real(r8), intent(out ) :: DU_flux(LBi:,LBj:)
410 real(r8), intent(out ) :: DV_flux(LBi:,LBj:)
411# endif
412
413#else
414
415# ifdef MASKING
416 real(r8), intent(in ) :: pmask(LBi:UBi,LBj:UBj)
417 real(r8), intent(in ) :: rmask(LBi:UBi,LBj:UBj)
418 real(r8), intent(in ) :: umask(LBi:UBi,LBj:UBj)
419 real(r8), intent(in ) :: vmask(LBi:UBi,LBj:UBj)
420# endif
421 real(r8), intent(in ) :: fomn(LBi:UBi,LBj:UBj)
422# if defined SEDIMENT && defined SED_MORPH
423 real(r8), intent(inout) :: h(LBi:UBi,LBj:UBj)
424# else
425 real(r8), intent(in ) :: h(LBi:UBi,LBj:UBj)
426# endif
427 real(r8), intent(in ) :: om_u(LBi:UBi,LBj:UBj)
428 real(r8), intent(in ) :: om_v(LBi:UBi,LBj:UBj)
429 real(r8), intent(in ) :: on_u(LBi:UBi,LBj:UBj)
430 real(r8), intent(in ) :: on_v(LBi:UBi,LBj:UBj)
431 real(r8), intent(in ) :: omn(LBi:UBi,LBj:UBj)
432 real(r8), intent(in ) :: pm(LBi:UBi,LBj:UBj)
433 real(r8), intent(in ) :: pn(LBi:UBi,LBj:UBj)
434# if defined CURVGRID && defined UV_ADV
435 real(r8), intent(in ) :: dndx(LBi:UBi,LBj:UBj)
436 real(r8), intent(in ) :: dmde(LBi:UBi,LBj:UBj)
437# endif
438# if defined UV_VIS2 || defined UV_VIS4
439 real(r8), intent(in ) :: pmon_r(LBi:UBi,LBj:UBj)
440 real(r8), intent(in ) :: pnom_r(LBi:UBi,LBj:UBj)
441 real(r8), intent(in ) :: pmon_p(LBi:UBi,LBj:UBj)
442 real(r8), intent(in ) :: pnom_p(LBi:UBi,LBj:UBj)
443 real(r8), intent(in ) :: om_r(LBi:UBi,LBj:UBj)
444 real(r8), intent(in ) :: on_r(LBi:UBi,LBj:UBj)
445 real(r8), intent(in ) :: om_p(LBi:UBi,LBj:UBj)
446 real(r8), intent(in ) :: on_p(LBi:UBi,LBj:UBj)
447# ifdef UV_VIS2
448 real(r8), intent(in ) :: visc2_p(LBi:UBi,LBj:UBj)
449 real(r8), intent(in ) :: visc2_r(LBi:UBi,LBj:UBj)
450# endif
451# ifdef UV_VIS4
452 real(r8), intent(in ) :: visc4_p(LBi:UBi,LBj:UBj)
453 real(r8), intent(in ) :: visc4_r(LBi:UBi,LBj:UBj)
454# endif
455# endif
456# if defined SEDIMENT && defined SED_MORPH
457 real(r8), intent(in ) :: bed_thick(LBi:UBi,LBj:UBj,1:3)
458# endif
459# ifdef WEC
460# ifdef WEC_VF
461# ifdef WEC_ROLLER
462 real(r8), intent(in) :: rurol2d(LBi:UBi,LBj:UBj)
463 real(r8), intent(in) :: rvrol2d(LBi:UBi,LBj:UBj)
464# endif
465# ifdef BOTTOM_STREAMING
466 real(r8), intent(in) :: rubst2d(LBi:UBi,LBj:UBj)
467 real(r8), intent(in) :: rvbst2d(LBi:UBi,LBj:UBj)
468# endif
469# ifdef SURFACE_STREAMING
470 real(r8), intent(in) :: russt2d(LBi:UBi,LBj:UBj)
471 real(r8), intent(in) :: rvsst2d(LBi:UBi,LBj:UBj)
472# endif
473 real(r8), intent(in) :: rubrk2d(LBi:UBi,LBj:UBj)
474 real(r8), intent(in) :: rvbrk2d(LBi:UBi,LBj:UBj)
475 real(r8), intent(in) :: rukvf2d(LBi:UBi,LBj:UBj)
476 real(r8), intent(in) :: rvkvf2d(LBi:UBi,LBj:UBj)
477 real(r8), intent(in) :: bh(LBi:UBi,LBj:UBj)
478 real(r8), intent(in) :: qsp(LBi:UBi,LBj:UBj)
479 real(r8), intent(in) :: zetaw(LBi:UBi,LBj:UBj)
480# endif
481 real(r8), intent(in) :: ubar_stokes(LBi:UBi,LBj:UBj)
482 real(r8), intent(in) :: vbar_stokes(LBi:UBi,LBj:UBj)
483# endif
484# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
485 real(r8), intent(in ) :: eq_tide(LBi:UBi,LBj:UBj)
486# endif
487# ifndef SOLVE3D
488 real(r8), intent(in ) :: sustr(LBi:UBi,LBj:UBj)
489 real(r8), intent(in ) :: svstr(LBi:UBi,LBj:UBj)
490 real(r8), intent(in ) :: bustr(LBi:UBi,LBj:UBj)
491 real(r8), intent(in ) :: bvstr(LBi:UBi,LBj:UBj)
492# ifdef ATM_PRESS
493 real(r8), intent(in ) :: Pair(LBi:UBi,LBj:UBj)
494# endif
495# else
496# ifdef VAR_RHO_2D
497 real(r8), intent(in ) :: rhoA(LBi:UBi,LBj:UBj)
498 real(r8), intent(in ) :: rhoS(LBi:UBi,LBj:UBj)
499# endif
500 real(r8), intent(inout) :: DU_avg1(LBi:UBi,LBj:UBj)
501 real(r8), intent(inout) :: DU_avg2(LBi:UBi,LBj:UBj)
502 real(r8), intent(inout) :: DV_avg1(LBi:UBi,LBj:UBj)
503 real(r8), intent(inout) :: DV_avg2(LBi:UBi,LBj:UBj)
504 real(r8), intent(inout) :: Zt_avg1(LBi:UBi,LBj:UBj)
505 real(r8), intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
506 real(r8), intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
507 real(r8), intent(inout) :: ru(LBi:UBi,LBj:UBj,0:UBk,2)
508 real(r8), intent(inout) :: rv(LBi:UBi,LBj:UBj,0:UBk,2)
509# endif
510# ifdef WET_DRY
511 real(r8), intent(inout) :: pmask_full(LBi:UBi,LBj:UBj)
512 real(r8), intent(inout) :: rmask_full(LBi:UBi,LBj:UBj)
513 real(r8), intent(inout) :: umask_full(LBi:UBi,LBj:UBj)
514 real(r8), intent(inout) :: vmask_full(LBi:UBi,LBj:UBj)
515
516 real(r8), intent(inout) :: pmask_wet(LBi:UBi,LBj:UBj)
517 real(r8), intent(inout) :: rmask_wet(LBi:UBi,LBj:UBj)
518 real(r8), intent(inout) :: umask_wet(LBi:UBi,LBj:UBj)
519 real(r8), intent(inout) :: vmask_wet(LBi:UBi,LBj:UBj)
520# ifdef SOLVE3D
521 real(r8), intent(inout) :: rmask_wet_avg(LBi:UBi,LBj:UBj)
522# endif
523# endif
524# ifdef DIAGNOSTICS_UV
525 real(r8), intent(inout) :: DiaU2wrk(LBi:UBi,LBj:UBj,NDM2d)
526 real(r8), intent(inout) :: DiaV2wrk(LBi:UBi,LBj:UBj,NDM2d)
527 real(r8), intent(inout) :: DiaRUbar(LBi:UBi,LBj:UBj,2,NDM2d-1)
528 real(r8), intent(inout) :: DiaRVbar(LBi:UBi,LBj:UBj,2,NDM2d-1)
529# ifdef SOLVE3D
530 real(r8), intent(inout) :: DiaU2int(LBi:UBi,LBj:UBj,NDM2d)
531 real(r8), intent(inout) :: DiaV2int(LBi:UBi,LBj:UBj,NDM2d)
532 real(r8), intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
533 real(r8), intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
534# endif
535# endif
536 real(r8), intent(inout) :: rubar(LBi:UBi,LBj:UBj,2)
537 real(r8), intent(inout) :: rvbar(LBi:UBi,LBj:UBj,2)
538 real(r8), intent(inout) :: rzeta(LBi:UBi,LBj:UBj,2)
539 real(r8), intent(inout) :: ubar(LBi:UBi,LBj:UBj,:)
540 real(r8), intent(inout) :: vbar(LBi:UBi,LBj:UBj,:)
541 real(r8), intent(inout) :: zeta(LBi:UBi,LBj:UBj,:)
542# if defined NESTING && !defined SOLVE3D
543 real(r8), intent(out ) :: DU_flux(LBi:UBi,LBj:UBj)
544 real(r8), intent(out ) :: DV_flux(LBi:UBi,LBj:UBj)
545# endif
546#endif
547!
548! Local variable declarations.
549!
550 logical :: CORRECTOR_2D_STEP
551!
552 integer :: i, is, j, ptsk
553#ifdef DIAGNOSTICS_UV
554 integer :: idiag
555#endif
556!
557 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5, cff6, cff7
558 real(r8) :: fac, fac1, fac2, fac3
559!
560 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dgrad
561 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dnew
562 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs
563 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs_p
564 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dstp
565 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUon
566 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVom
567#ifdef WEC
568 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUSon
569 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVSom
570#endif
571#ifdef WEC_VF
572 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUSom
573 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVSon
574#endif
575#ifdef UV_VIS4
576 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LapU
577 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LapV
578#endif
579 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
580 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
581 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
582 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
583 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
584 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gzeta
585 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gzeta2
586#if defined VAR_RHO_2D && defined SOLVE3D
587 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gzetaSA
588#endif
589 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rhs_ubar
590 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rhs_vbar
591 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rhs_zeta
592 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zeta_new
593 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zwrk
594#ifdef WET_DRY
595 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: wetdry
596#endif
597#ifdef DIAGNOSTICS_UV
598 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Uwrk
599 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vwrk
600 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaU2rhs
601 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaV2rhs
602#endif
603
604#include "set_bounds.h"
605!
606 ptsk=3-kstp
607 corrector_2d_step=.not.predictor_2d_step(ng)
608!
609!-----------------------------------------------------------------------
610! Compute total depth (m) and vertically integrated mass fluxes.
611!-----------------------------------------------------------------------
612!
613#if defined DISTRIBUTE && !defined NESTING
614
615! In distributed-memory, the I- and J-ranges are different and a
616! special exchange is done to avoid having three ghost points for
617! high order numerical stencils. Notice that a private array is
618! passed below to the exchange routine. It also applies periodic
619! boundary conditions, if appropriate and no partitions in I- or
620! J-directions.
621!
622 DO j=jstrv-2,jendp2
623 DO i=istru-2,iendp2
624 drhs(i,j)=zeta(i,j,krhs)+h(i,j)
625 END DO
626 END DO
627 DO j=jstrv-2,jendp2
628 DO i=istru-1,iendp2
629 cff=0.5_r8*on_u(i,j)
630 cff1=cff*(drhs(i,j)+drhs(i-1,j))
631 duon(i,j)=ubar(i,j,krhs)*cff1
632# ifdef WEC
633# ifdef WET_DRY
634 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
635 cff6=0.5_r8+dsign(0.5_r8,ubar_stokes(i,j))*umask_wet(i,j)
636 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
637 cff1=cff1*cff7
638# endif
639 duson(i,j)=ubar_stokes(i,j)*cff1
640 duon(i,j)=duon(i,j)+duson(i,j)
641# endif
642 END DO
643 END DO
644 DO j=jstrv-1,jendp2
645 DO i=istru-2,iendp2
646 cff=0.5_r8*om_v(i,j)
647 cff1=cff*(drhs(i,j)+drhs(i,j-1))
648 dvom(i,j)=vbar(i,j,krhs)*cff1
649# ifdef WEC
650# ifdef WET_DRY
651 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
652 cff6=0.5_r8+dsign(0.5_r8,vbar_stokes(i,j))*vmask_wet(i,j)
653 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
654 cff1=cff1*cff7
655# endif
656 dvsom(i,j)=vbar_stokes(i,j)*cff1
657 dvom(i,j)=dvom(i,j)+dvsom(i,j)
658# endif
659 END DO
660 END DO
661
662#else
663
664 DO j=jstrvm2-1,jendp2
665 DO i=istrum2-1,iendp2
666 drhs(i,j)=zeta(i,j,krhs)+h(i,j)
667 END DO
668 END DO
669 DO j=jstrvm2-1,jendp2
670 DO i=istrum2,iendp2
671 cff=0.5_r8*on_u(i,j)
672 cff1=cff*(drhs(i,j)+drhs(i-1,j))
673 duon(i,j)=ubar(i,j,krhs)*cff1
674# ifdef WEC
675# ifdef WET_DRY
676 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
677 cff6=0.5_r8+dsign(0.5_r8,ubar_stokes(i,j))*umask_wet(i,j)
678 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
679 cff1=cff1*cff7
680# endif
681 duson(i,j)=ubar_stokes(i,j)*cff1
682 duon(i,j)=duon(i,j)+duson(i,j)
683# endif
684 END DO
685 END DO
686 DO j=jstrvm2,jendp2
687 DO i=istrum2-1,iendp2
688 cff=0.5_r8*om_v(i,j)
689 cff1=cff*(drhs(i,j)+drhs(i,j-1))
690 dvom(i,j)=vbar(i,j,krhs)*cff1
691# ifdef WEC
692# ifdef WET_DRY
693 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
694 cff6=0.5_r8+dsign(0.5_r8,vbar_stokes(i,j))*vmask_wet(i,j)
695 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
696 cff1=cff1*cff7
697# endif
698 dvsom(i,j)=vbar_stokes(i,j)*cff1
699 dvom(i,j)=dvom(i,j)+dvsom(i,j)
700# endif
701 END DO
702 END DO
703#endif
704#ifdef DISTRIBUTE
705!
706 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
707 CALL exchange_u2d_tile (ng, tile, &
708 & imins, imaxs, jmins, jmaxs, &
709 & duon)
710 CALL exchange_v2d_tile (ng, tile, &
711 & imins, imaxs, jmins, jmaxs, &
712 & dvom)
713 END IF
714 CALL mp_exchange2d (ng, tile, inlm, 2, &
715 & imins, imaxs, jmins, jmaxs, &
716 & nghostpoints, &
717 & ewperiodic(ng), nsperiodic(ng), &
718 & duon, dvom)
719#endif
720!
721! Set vertically integrated mass fluxes DUon and DVom along the open
722! boundaries in such a way that the integral volume is conserved.
723!
724 IF (any(volcons(:,ng))) THEN
725 CALL set_duv_bc_tile (ng, tile, &
726 & lbi, ubi, lbj, ubj, &
727 & imins, imaxs, jmins, jmaxs, &
728 & krhs, &
729#ifdef MASKING
730 & umask, vmask, &
731#endif
732 & om_v, on_u, &
733 & ubar, vbar, &
734 & drhs, duon, dvom)
735 END IF
736#ifdef SOLVE3D
737!
738!-----------------------------------------------------------------------
739! Compute time averaged fields over all short time-steps.
740!-----------------------------------------------------------------------
741!
742 IF (predictor_2d_step(ng)) THEN
743 IF (first_2d_step) THEN
744!
745! Reset arrays for 2D fields averaged within the short time-steps.
746!
747 cff2=(-1.0_r8/12.0_r8)*weight(2,iif(ng)+1,ng)
748 DO j=jstrr,jendr
749 DO i=istrr,iendr
750 zt_avg1(i,j)=0.0_r8
751 END DO
752 DO i=istr,iendr
753 du_avg1(i,j)=0.0_r8
754 du_avg2(i,j)=cff2*duon(i,j)
755 END DO
756 END DO
757 DO j=jstr,jendr
758 DO i=istrr,iendr
759 dv_avg1(i,j)=0.0_r8
760 dv_avg2(i,j)=cff2*dvom(i,j)
761 END DO
762 END DO
763 ELSE
764!
765! Accumulate field averages of previous time-step after they are
766! computed in the previous corrector step, updated their boundaries,
767! and synchronized.
768!
769 cff1=weight(1,iif(ng)-1,ng)
770 cff2=(8.0_r8/12.0_r8)*weight(2,iif(ng) ,ng)- &
771 & (1.0_r8/12.0_r8)*weight(2,iif(ng)+1,ng)
772 DO j=jstrr,jendr
773 DO i=istrr,iendr
774 zt_avg1(i,j)=zt_avg1(i,j)+cff1*zeta(i,j,krhs)
775 END DO
776 DO i=istr,iendr
777 du_avg1(i,j)=du_avg1(i,j)+cff1*duon(i,j)
778# ifdef WEC
779 du_avg1(i,j)=du_avg1(i,j)-cff1*duson(i,j)
780# endif
781 du_avg2(i,j)=du_avg2(i,j)+cff2*duon(i,j)
782 END DO
783 END DO
784 DO j=jstr,jendr
785 DO i=istrr,iendr
786 dv_avg1(i,j)=dv_avg1(i,j)+cff1*dvom(i,j)
787# ifdef WEC
788 dv_avg1(i,j)=dv_avg1(i,j)-cff1*dvsom(i,j)
789# endif
790 dv_avg2(i,j)=dv_avg2(i,j)+cff2*dvom(i,j)
791 END DO
792 END DO
793 END IF
794 ELSE
795 IF (first_2d_step) THEN
796 cff2=weight(2,iif(ng),ng)
797 ELSE
798 cff2=(5.0_r8/12.0_r8)*weight(2,iif(ng),ng)
799 END IF
800 DO j=jstrr,jendr
801 DO i=istr,iendr
802 du_avg2(i,j)=du_avg2(i,j)+cff2*duon(i,j)
803 END DO
804 END DO
805 DO j=jstr,jendr
806 DO i=istrr,iendr
807 dv_avg2(i,j)=dv_avg2(i,j)+cff2*dvom(i,j)
808 END DO
809 END DO
810 END IF
811!
812! After all fast time steps are completed, apply boundary conditions
813! to time averaged fields.
814# ifdef NESTING
815! In nesting applications with refinement grids, we need to exchange
816! the DU_avg2 and DV_avg2 fluxes boundary information for the case
817! that a contact point is at a tile partition. Notice that in such
818! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
819# endif
820!
821 IF ((iif(ng).eq.(nfast(ng)+1)).and.predictor_2d_step(ng)) THEN
822 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
823 CALL exchange_r2d_tile (ng, tile, &
824 & lbi, ubi, lbj, ubj, &
825 & zt_avg1)
826 CALL exchange_u2d_tile (ng, tile, &
827 & lbi, ubi, lbj, ubj, &
828 & du_avg1)
829 CALL exchange_v2d_tile (ng, tile, &
830 & lbi, ubi, lbj, ubj, &
831 & dv_avg1)
832# ifdef NESTING
833 CALL exchange_u2d_tile (ng, tile, &
834 & lbi, ubi, lbj, ubj, &
835 & du_avg2)
836 CALL exchange_v2d_tile (ng, tile, &
837 & lbi, ubi, lbj, ubj, &
838 & dv_avg2)
839# endif
840 END IF
841# ifdef DISTRIBUTE
842 CALL mp_exchange2d (ng, tile, inlm, 3, &
843 & lbi, ubi, lbj, ubj, &
844 & nghostpoints, &
845 & ewperiodic(ng), nsperiodic(ng), &
846 & zt_avg1, du_avg1, dv_avg1)
847# ifdef NESTING
848 CALL mp_exchange2d (ng, tile, inlm, 2, &
849 & lbi, ubi, lbj, ubj, &
850 & nghostpoints, &
851 & ewperiodic(ng), nsperiodic(ng), &
852 & du_avg2, dv_avg2)
853# endif
854# endif
855 END IF
856#endif
857#ifdef WET_DRY
858!
859!-----------------------------------------------------------------------
860! Compute new wet/dry masks.
861!-----------------------------------------------------------------------
862!
863 CALL wetdry_tile (ng, tile, &
864 & lbi, ubi, lbj, ubj, &
865 & imins, imaxs, jmins, jmaxs, &
866# ifdef MASKING
867 & pmask, rmask, umask, vmask, &
868# endif
869 & h, zeta(:,:,kstp), &
870# ifdef SOLVE3D
871 & du_avg1, dv_avg1, &
872 & rmask_wet_avg, &
873# endif
874 & pmask_wet, pmask_full, &
875 & rmask_wet, rmask_full, &
876 & umask_wet, umask_full, &
877 & vmask_wet, vmask_full)
878#endif
879!
880! Do not perform the actual time stepping during the auxiliary
881! (nfast(ng)+1) time step.
882!
883 IF (iif(ng).gt.nfast(ng)) RETURN
884!
885!=======================================================================
886! Time step free-surface equation.
887!=======================================================================
888!
889! During the first time-step, the predictor step is Forward-Euler
890! and the corrector step is Backward-Euler. Otherwise, the predictor
891! step is Leap-frog and the corrector step is Adams-Moulton.
892#if defined VAR_RHO_2D && defined SOLVE3D
893! Recall that the vertical averaged density (rhoA) and density
894! pertubation (rhoS) are nondimensional quantities.
895!
896 fac=1000.0_r8/rho0 ! nondimensional
897#endif
898!
899 IF (first_2d_step) THEN
900 cff1=dtfast(ng)
901 DO j=jstrv-1,jend
902 DO i=istru-1,iend
903 rhs_zeta(i,j)=(duon(i,j)-duon(i+1,j))+ &
904 & (dvom(i,j)-dvom(i,j+1))
905 zeta_new(i,j)=zeta(i,j,kstp)+ &
906 & pm(i,j)*pn(i,j)*cff1*rhs_zeta(i,j)
907#ifdef MASKING
908 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
909#endif
910 dnew(i,j)=zeta_new(i,j)+h(i,j)
911!
912 zwrk(i,j)=0.5_r8*(zeta(i,j,kstp)+zeta_new(i,j))
913#if defined VAR_RHO_2D && defined SOLVE3D
914 gzeta(i,j)=(fac+rhos(i,j))*zwrk(i,j)
915 gzeta2(i,j)=gzeta(i,j)*zwrk(i,j)
916 gzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
917#else
918 gzeta(i,j)=zwrk(i,j)
919 gzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
920#endif
921 END DO
922 END DO
923 ELSE IF (predictor_2d_step(ng)) THEN
924 cff1=2.0_r8*dtfast(ng)
925 cff4=4.0_r8/25.0_r8
926 cff5=1.0_r8-2.0_r8*cff4
927 DO j=jstrv-1,jend
928 DO i=istru-1,iend
929 rhs_zeta(i,j)=(duon(i,j)-duon(i+1,j))+ &
930 & (dvom(i,j)-dvom(i,j+1))
931 zeta_new(i,j)=zeta(i,j,kstp)+ &
932 & pm(i,j)*pn(i,j)*cff1*rhs_zeta(i,j)
933#ifdef MASKING
934 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
935#endif
936 dnew(i,j)=zeta_new(i,j)+h(i,j)
937!
938 zwrk(i,j)=cff5*zeta(i,j,krhs)+ &
939 & cff4*(zeta(i,j,kstp)+zeta_new(i,j))
940#if defined VAR_RHO_2D && defined SOLVE3D
941 gzeta(i,j)=(fac+rhos(i,j))*zwrk(i,j)
942 gzeta2(i,j)=gzeta(i,j)*zwrk(i,j)
943 gzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
944#else
945 gzeta(i,j)=zwrk(i,j)
946 gzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
947#endif
948 END DO
949 END DO
950 ELSE IF (corrector_2d_step) THEN
951 cff1=dtfast(ng)*5.0_r8/12.0_r8
952 cff2=dtfast(ng)*8.0_r8/12.0_r8
953 cff3=dtfast(ng)*1.0_r8/12.0_r8
954 cff4=2.0_r8/5.0_r8
955 cff5=1.0_r8-cff4
956 DO j=jstrv-1,jend
957 DO i=istru-1,iend
958 cff=cff1*((duon(i,j)-duon(i+1,j))+ &
959 & (dvom(i,j)-dvom(i,j+1)))
960 zeta_new(i,j)=zeta(i,j,kstp)+ &
961 & pm(i,j)*pn(i,j)*(cff+ &
962 & cff2*rzeta(i,j,kstp)- &
963 & cff3*rzeta(i,j,ptsk))
964#ifdef MASKING
965 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
966#endif
967 dnew(i,j)=zeta_new(i,j)+h(i,j)
968!
969 zwrk(i,j)=cff5*zeta_new(i,j)+cff4*zeta(i,j,krhs)
970#if defined VAR_RHO_2D && defined SOLVE3D
971 gzeta(i,j)=(fac+rhos(i,j))*zwrk(i,j)
972 gzeta2(i,j)=gzeta(i,j)*zwrk(i,j)
973 gzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
974#else
975 gzeta(i,j)=zwrk(i,j)
976 gzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
977#endif
978 END DO
979 END DO
980 END IF
981!
982! Load new free-surface values into shared array at both predictor
983! and corrector steps.
984#ifdef WET_DRY
985! Modify new free-surface to Ensure that depth is > Dcrit for masked
986! cells.
987#endif
988!
989 DO j=jstr,jend
990 DO i=istr,iend
991 zeta(i,j,knew)=zeta_new(i,j)
992#if defined WET_DRY && defined MASKING
993 zeta(i,j,knew)=zeta(i,j,knew)+ &
994 & (dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
995#endif
996 END DO
997 END DO
998!
999! If predictor step, load right-side-term into shared array.
1000!
1001 IF (predictor_2d_step(ng)) THEN
1002 DO j=jstr,jend
1003 DO i=istr,iend
1004 rzeta(i,j,krhs)=rhs_zeta(i,j)
1005 END DO
1006 END DO
1007 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1008 CALL exchange_r2d_tile (ng, tile, &
1009 & lbi, ubi, lbj, ubj, &
1010 & rzeta(:,:,krhs))
1011 END IF
1012#ifdef DISTRIBUTE
1013 CALL mp_exchange2d (ng, tile, inlm, 1, &
1014 & lbi, ubi, lbj, ubj, &
1015 & nghostpoints, &
1016 & ewperiodic(ng), nsperiodic(ng), &
1017 & rzeta(:,:,krhs))
1018#endif
1019 END IF
1020!
1021! Apply mass point sources (volume vertical influx), if any.
1022!
1023! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
1024!
1025 IF (lwsrc(ng)) THEN
1026 DO is=1,nsrc(ng)
1027 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
1028 i=sources(ng)%Isrc(is)
1029 j=sources(ng)%Jsrc(is)
1030 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1031 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1032 zeta(i,j,knew)=zeta(i,j,knew)+ &
1033 & sources(ng)%Qbar(is)* &
1034 & pm(i,j)*pn(i,j)*dtfast(ng)
1035 END IF
1036 END IF
1037 END DO
1038 END IF
1039
1040#if defined SEDIMENT && defined SED_MORPH
1041!
1042! Scale the bed change with the fast time stepping. The half is
1043! becasue we do predictor and corrector. The "ndtfast/nfast" is
1044! becasue we do "nfast" steps to here.
1045!
1046 fac=0.5_r8*dtfast(ng)*ndtfast(ng)/(nfast(ng)*dt(ng))
1047 DO j=jstr,jend
1048 DO i=istr,iend
1049 cff=fac*(bed_thick(i,j,nstp)-bed_thick(i,j,nnew))
1050 h(i,j)=h(i,j)-cff
1051 END DO
1052 END DO
1053#endif
1054!
1055! Set free-surface lateral boundary conditions.
1056!
1057 CALL zetabc_tile (ng, tile, &
1058 & lbi, ubi, lbj, ubj, &
1059 & imins, imaxs, jmins, jmaxs, &
1060 & krhs, kstp, knew, &
1061 & zeta)
1062 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1063 CALL exchange_r2d_tile (ng, tile, &
1064 & lbi, ubi, lbj, ubj, &
1065 & zeta(:,:,knew))
1066 END IF
1067#ifdef DISTRIBUTE
1068 CALL mp_exchange2d (ng, tile, inlm, 1, &
1069 & lbi, ubi, lbj, ubj, &
1070 & nghostpoints, &
1071 & ewperiodic(ng), nsperiodic(ng), &
1072 & zeta(:,:,knew))
1073#endif
1074!
1075!=======================================================================
1076! Compute right-hand-side for the 2D momentum equations.
1077!=======================================================================
1078!
1079!-----------------------------------------------------------------------
1080! Compute pressure gradient terms.
1081!-----------------------------------------------------------------------
1082!
1083 cff1=0.5_r8*g
1084 cff2=1.0_r8/3.0_r8
1085#if !defined SOLVE3D && defined ATM_PRESS
1086 fac3=0.5_r8*100.0_r8/rho0
1087#endif
1088 DO j=jstr,jend
1089 DO i=istru,iend
1090 rhs_ubar(i,j)=cff1*on_u(i,j)* &
1091 & ((h(i-1,j)+ &
1092 & h(i ,j))* &
1093 & (gzeta(i-1,j)- &
1094 & gzeta(i ,j))+ &
1095#if defined VAR_RHO_2D && defined SOLVE3D
1096 & (h(i-1,j)- &
1097 & h(i ,j))* &
1098 & (gzetasa(i-1,j)+ &
1099 & gzetasa(i ,j)+ &
1100 & cff2*(rhoa(i-1,j)- &
1101 & rhoa(i ,j))* &
1102 & (zwrk(i-1,j)- &
1103 & zwrk(i ,j)))+ &
1104#endif
1105 & (gzeta2(i-1,j)- &
1106 & gzeta2(i ,j)))
1107#if defined ATM_PRESS && !defined SOLVE3D
1108 rhs_ubar(i,j)=rhs_ubar(i,j)- &
1109 & fac3*on_u(i,j)* &
1110 & (h(i-1,j)+h(i,j)+ &
1111 & gzeta(i-1,j)+gzeta(i,j))* &
1112 & (pair(i,j)-pair(i-1,j))
1113#endif
1114#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
1115 rhs_ubar(i,j)=rhs_ubar(i,j)- &
1116 & cff1*on_u(i,j)* &
1117 & (h(i-1,j)+h(i,j)+ &
1118 & gzeta(i-1,j)+gzeta(i,j))* &
1119 & (eq_tide(i,j)-eq_tide(i-1,j))
1120#endif
1121#ifdef DIAGNOSTICS_UV
1122 diau2rhs(i,j,m2pgrd)=rhs_ubar(i,j)
1123#endif
1124#if defined WEC_VF
1125 cff3=0.5_r8*on_u(i,j)* &
1126 & (h(i-1,j)+h(i,j)+ &
1127 & gzeta(i-1,j)+gzeta(i,j))
1128 cff4=cff3*g*(zetaw(i-1,j)-zetaw(i,j))
1129 cff5=cff3*g*(qsp(i-1,j)-qsp(i,j))
1130 cff6=cff3*(bh(i-1,j)-bh(i,j))
1131 cff7=rukvf2d(i,j)
1132 rhs_ubar(i,j)=rhs_ubar(i,j)-cff4-cff5+cff6+cff7
1133# ifdef DIAGNOSTICS_UV
1134 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2pgrd)
1135 diau2rhs(i,j,m2pgrd)=diau2rhs(i,j,m2pgrd)-cff4-cff5+cff6
1136 diau2rhs(i,j,m2zetw)=-cff4
1137 diau2rhs(i,j,m2zqsp)=-cff5
1138 diau2rhs(i,j,m2zbeh)=cff6
1139 diau2rhs(i,j,m2kvrf)=cff7
1140# ifndef UV_ADV
1141 diau2rhs(i,j,m2hjvf)=0.0_r8
1142# endif
1143# endif
1144#endif
1145 END DO
1146 IF (j.ge.jstrv) THEN
1147 DO i=istr,iend
1148 rhs_vbar(i,j)=cff1*om_v(i,j)* &
1149 & ((h(i,j-1)+ &
1150 & h(i,j ))* &
1151 & (gzeta(i,j-1)- &
1152 & gzeta(i,j ))+ &
1153#if defined VAR_RHO_2D && defined SOLVE3D
1154 & (h(i,j-1)- &
1155 & h(i,j ))* &
1156 & (gzetasa(i,j-1)+ &
1157 & gzetasa(i,j )+ &
1158 & cff2*(rhoa(i,j-1)- &
1159 & rhoa(i,j ))* &
1160 & (zwrk(i,j-1)- &
1161 & zwrk(i,j )))+ &
1162#endif
1163 & (gzeta2(i,j-1)- &
1164 & gzeta2(i,j )))
1165#if defined ATM_PRESS && !defined SOLVE3D
1166 rhs_vbar(i,j)=rhs_vbar(i,j)- &
1167 & fac3*om_v(i,j)* &
1168 & (h(i,j-1)+h(i,j)+ &
1169 & gzeta(i,j-1)+gzeta(i,j))* &
1170 & (pair(i,j)-pair(i,j-1))
1171#endif
1172#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
1173 rhs_vbar(i,j)=rhs_vbar(i,j)- &
1174 & cff1*om_v(i,j)* &
1175 & (h(i,j-1)+h(i,j)+ &
1176 & gzeta(i,j-1)+gzeta(i,j))* &
1177 & (eq_tide(i,j)-eq_tide(i,j-1))
1178#endif
1179#ifdef DIAGNOSTICS_UV
1180 diav2rhs(i,j,m2pgrd)=rhs_vbar(i,j)
1181#endif
1182#if defined WEC_VF
1183 cff3=0.5_r8*om_v(i,j)* &
1184 & (h(i,j-1)+h(i,j)+ &
1185 & gzeta(i,j-1)+gzeta(i,j))
1186 cff4=cff3*g*(zetaw(i,j-1)-zetaw(i,j))
1187 cff5=cff3*g*(qsp(i,j-1)-qsp(i,j))
1188 cff6=cff3*(bh(i,j-1)-bh(i,j))
1189 cff7=rvkvf2d(i,j)
1190 rhs_vbar(i,j)=rhs_vbar(i,j)-cff4-cff5+cff6+cff7
1191# ifdef DIAGNOSTICS_UV
1192 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2pgrd)
1193 diav2rhs(i,j,m2pgrd)=diav2rhs(i,j,m2pgrd)-cff4-cff5+cff6
1194 diav2rhs(i,j,m2zetw)=-cff4
1195 diav2rhs(i,j,m2zqsp)=-cff5
1196 diav2rhs(i,j,m2zbeh)=cff6
1197 diav2rhs(i,j,m2kvrf)=cff7
1198# ifndef UV_ADV
1199 diav2rhs(i,j,m2hjvf)=0.0_r8
1200# endif
1201# endif
1202#endif
1203 END DO
1204 END IF
1205 END DO
1206#ifdef UV_ADV
1207!
1208!-----------------------------------------------------------------------
1209! Add in horizontal advection of momentum.
1210!-----------------------------------------------------------------------
1211
1212# ifdef UV_C2ADVECTION
1213!
1214! Second-order, centered differences advection fluxes.
1215!
1216 DO j=jstr,jend
1217 DO i=istru-1,iend
1218 ufx(i,j)=0.25_r8*(duon(i,j)+duon(i+1,j))* &
1219 & (ubar(i ,j,krhs)+ &
1220 & ubar(i+1,j,krhs))
1221 END DO
1222 END DO
1223!
1224 DO j=jstr,jend+1
1225 DO i=istru,iend
1226 ufe(i,j)=0.25_r8*(dvom(i,j)+dvom(i-1,j))* &
1227 & (ubar(i,j ,krhs)+ &
1228 & ubar(i,j-1,krhs))
1229 END DO
1230 END DO
1231!
1232 DO j=jstrv,jend
1233 DO i=istr,iend+1
1234 vfx(i,j)=0.25_r8*(duon(i,j)+duon(i,j-1))* &
1235 & (vbar(i ,j,krhs)+ &
1236 & vbar(i-1,j,krhs))
1237 END DO
1238 END DO
1239!
1240 DO j=jstrv-1,jend
1241 DO i=istr,iend
1242 vfe(i,j)=0.25_r8*(dvom(i,j)+dvom(i,j+1))* &
1243 & (vbar(i,j ,krhs)+ &
1244 & vbar(i,j+1,krhs))
1245 END DO
1246 END DO
1247# else
1248!
1249! Fourth-order, centered differences advection fluxes.
1250!
1251 DO j=jstr,jend
1252 DO i=istrum1,iendp1
1253 grad(i,j)=ubar(i-1,j,krhs)-2.0_r8*ubar(i,j,krhs)+ &
1254 & ubar(i+1,j,krhs)
1255 dgrad(i,j)=duon(i-1,j)-2.0_r8*duon(i,j)+duon(i+1,j)
1256 END DO
1257 END DO
1258 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1259 IF (domain(ng)%Western_Edge(tile)) THEN
1260 DO j=jstr,jend
1261 grad(istr,j)=grad(istr+1,j)
1262 dgrad(istr,j)=dgrad(istr+1,j)
1263 END DO
1264 END IF
1265 END IF
1266 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1267 IF (domain(ng)%Eastern_Edge(tile)) THEN
1268 DO j=jstr,jend
1269 grad(iend+1,j)=grad(iend,j)
1270 dgrad(iend+1,j)=dgrad(iend,j)
1271 END DO
1272 END IF
1273 END IF
1274
1275 cff=1.0_r8/6.0_r8
1276 DO j=jstr,jend
1277 DO i=istru-1,iend
1278 ufx(i,j)=0.25_r8*(ubar(i ,j,krhs)+ &
1279 & ubar(i+1,j,krhs)- &
1280 & cff*(grad(i,j)+grad(i+1,j)))* &
1281 & (duon(i,j)+duon(i+1,j)- &
1282 & cff*(dgrad(i,j)+dgrad(i+1,j)))
1283 END DO
1284 END DO
1285!
1286 DO j=jstrm1,jendp1
1287 DO i=istru,iend
1288 grad(i,j)=ubar(i,j-1,krhs)-2.0_r8*ubar(i,j,krhs)+ &
1289 & ubar(i,j+1,krhs)
1290 END DO
1291 END DO
1292 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1293 IF (domain(ng)%Southern_Edge(tile)) THEN
1294 DO i=istru,iend
1295 grad(i,jstr-1)=grad(i,jstr)
1296 END DO
1297 END IF
1298 END IF
1299 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1300 IF (domain(ng)%Northern_Edge(tile)) THEN
1301 DO i=istru,iend
1302 grad(i,jend+1)=grad(i,jend)
1303 END DO
1304 END IF
1305 END IF
1306 DO j=jstr,jend+1
1307 DO i=istru-1,iend
1308 dgrad(i,j)=dvom(i-1,j)-2.0_r8*dvom(i,j)+dvom(i+1,j)
1309 END DO
1310 END DO
1311
1312 cff=1.0_r8/6.0_r8
1313 DO j=jstr,jend+1
1314 DO i=istru,iend
1315 ufe(i,j)=0.25_r8*(ubar(i,j ,krhs)+ &
1316 & ubar(i,j-1,krhs)- &
1317 & cff*(grad(i,j)+grad(i,j-1)))* &
1318 & (dvom(i,j)+dvom(i-1,j)- &
1319 & cff*(dgrad(i,j)+dgrad(i-1,j)))
1320 END DO
1321 END DO
1322!
1323 DO j=jstrv,jend
1324 DO i=istrm1,iendp1
1325 grad(i,j)=vbar(i-1,j,krhs)-2.0_r8*vbar(i,j,krhs)+ &
1326 & vbar(i+1,j,krhs)
1327 END DO
1328 END DO
1329 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1330 IF (domain(ng)%Western_Edge(tile)) THEN
1331 DO j=jstrv,jend
1332 grad(istr-1,j)=grad(istr,j)
1333 END DO
1334 END IF
1335 END IF
1336 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1337 IF (domain(ng)%Eastern_Edge(tile)) THEN
1338 DO j=jstrv,jend
1339 grad(iend+1,j)=grad(iend,j)
1340 END DO
1341 END IF
1342 END IF
1343 DO j=jstrv-1,jend
1344 DO i=istr,iend+1
1345 dgrad(i,j)=duon(i,j-1)-2.0_r8*duon(i,j)+duon(i,j+1)
1346 END DO
1347 END DO
1348
1349 cff=1.0_r8/6.0_r8
1350 DO j=jstrv,jend
1351 DO i=istr,iend+1
1352 vfx(i,j)=0.25_r8*(vbar(i ,j,krhs)+ &
1353 & vbar(i-1,j,krhs)- &
1354 & cff*(grad(i,j)+grad(i-1,j)))* &
1355 & (duon(i,j)+duon(i,j-1)- &
1356 & cff*(dgrad(i,j)+dgrad(i,j-1)))
1357 END DO
1358 END DO
1359!
1360 DO j=jstrvm1,jendp1
1361 DO i=istr,iend
1362 grad(i,j)=vbar(i,j-1,krhs)-2.0_r8*vbar(i,j,krhs)+ &
1363 & vbar(i,j+1,krhs)
1364 dgrad(i,j)=dvom(i,j-1)-2.0_r8*dvom(i,j)+dvom(i,j+1)
1365 END DO
1366 END DO
1367 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1368 IF (domain(ng)%Southern_Edge(tile)) THEN
1369 DO i=istr,iend
1370 grad(i,jstr)=grad(i,jstr+1)
1371 dgrad(i,jstr)=dgrad(i,jstr+1)
1372 END DO
1373 END IF
1374 END IF
1375 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1376 IF (domain(ng)%Northern_Edge(tile)) THEN
1377 DO i=istr,iend
1378 grad(i,jend+1)=grad(i,jend)
1379 dgrad(i,jend+1)=dgrad(i,jend)
1380 END DO
1381 END IF
1382 END IF
1383
1384 cff=1.0_r8/6.0_r8
1385 DO j=jstrv-1,jend
1386 DO i=istr,iend
1387 vfe(i,j)=0.25_r8*(vbar(i,j ,krhs)+ &
1388 & vbar(i,j+1,krhs)- &
1389 & cff*(grad(i,j)+grad(i,j+1)))* &
1390 & (dvom(i,j)+dvom(i,j+1)- &
1391 & cff*(dgrad(i,j)+dgrad(i,j+1)))
1392 END DO
1393 END DO
1394# endif
1395!
1396! Add advection to RHS terms.
1397!
1398 DO j=jstr,jend
1399 DO i=istru,iend
1400 cff1=ufx(i,j)-ufx(i-1,j)
1401 cff2=ufe(i,j+1)-ufe(i,j)
1402 fac=cff1+cff2
1403 rhs_ubar(i,j)=rhs_ubar(i,j)-fac
1404# if defined DIAGNOSTICS_UV
1405 diau2rhs(i,j,m2xadv)=-cff1
1406 diau2rhs(i,j,m2yadv)=-cff2
1407 diau2rhs(i,j,m2hadv)=-fac
1408# endif
1409 END DO
1410 END DO
1411 DO j=jstrv,jend
1412 DO i=istr,iend
1413 cff1=vfx(i+1,j)-vfx(i,j)
1414 cff2=vfe(i,j)-vfe(i,j-1)
1415 fac=cff1+cff2
1416 rhs_vbar(i,j)=rhs_vbar(i,j)-fac
1417# if defined DIAGNOSTICS_UV
1418 diav2rhs(i,j,m2xadv)=-cff1
1419 diav2rhs(i,j,m2yadv)=-cff2
1420 diav2rhs(i,j,m2hadv)=-fac
1421# endif
1422 END DO
1423 END DO
1424#endif
1425
1426#ifdef UV_COR
1427!
1428!-----------------------------------------------------------------------
1429! Add in Coriolis term.
1430!-----------------------------------------------------------------------
1431!
1432 DO j=jstrv-1,jend
1433 DO i=istru-1,iend
1434 cff=0.5_r8*drhs(i,j)*fomn(i,j)
1435 ufx(i,j)=cff*(vbar(i,j ,krhs)+ &
1436 & vbar(i,j+1,krhs))
1437 vfe(i,j)=cff*(ubar(i ,j,krhs)+ &
1438 & ubar(i+1,j,krhs))
1439 END DO
1440 END DO
1441 DO j=jstr,jend
1442 DO i=istru,iend
1443 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1444 rhs_ubar(i,j)=rhs_ubar(i,j)+fac1
1445# if defined DIAGNOSTICS_UV
1446 diau2rhs(i,j,m2fcor)=fac1
1447# endif
1448 END DO
1449 END DO
1450 DO j=jstrv,jend
1451 DO i=istr,iend
1452 fac1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1453 rhs_vbar(i,j)=rhs_vbar(i,j)-fac1
1454# if defined DIAGNOSTICS_UV
1455 diav2rhs(i,j,m2fcor)=-fac1
1456# endif
1457 END DO
1458 END DO
1459!
1460# ifdef WEC
1461 DO j=jstrv-1,jend
1462 DO i=istru-1,iend
1463 cff=0.5_r8*drhs(i,j)*fomn(i,j)
1464 ufx(i,j)=cff*(vbar_stokes(i,j )+ &
1465 & vbar_stokes(i,j+1))
1466 vfe(i,j)=cff*(ubar_stokes(i ,j)+ &
1467 & ubar_stokes(i+1,j))
1468 END DO
1469 END DO
1470 DO j=jstr,jend
1471 DO i=istru,iend
1472 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1473 rhs_ubar(i,j)=rhs_ubar(i,j)+fac1
1474# if defined DIAGNOSTICS_UV
1475 diau2rhs(i,j,m2fsco)=fac1
1476# endif
1477 END DO
1478 END DO
1479 DO j=jstrv,jend
1480 DO i=istr,iend
1481 fac1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1482 rhs_vbar(i,j)=rhs_vbar(i,j)-fac1
1483# if defined DIAGNOSTICS_UV
1484 diav2rhs(i,j,m2fsco)=-fac1
1485# endif
1486 END DO
1487 END DO
1488# endif
1489#endif
1490
1491#if defined CURVGRID && defined UV_ADV
1492!
1493!-----------------------------------------------------------------------
1494! Add in curvilinear transformation terms.
1495!-----------------------------------------------------------------------
1496!
1497 DO j=jstrv-1,jend
1498 DO i=istru-1,iend
1499 cff1=0.5_r8*(vbar(i,j ,krhs)+ &
1500# ifdef WEC
1501 & vbar_stokes(i,j )+ &
1502 & vbar_stokes(i,j+1)+ &
1503# endif
1504 & vbar(i,j+1,krhs))
1505 cff2=0.5_r8*(ubar(i ,j,krhs)+ &
1506# ifdef WEC
1507 & ubar_stokes(i ,j)+ &
1508 & ubar_stokes(i+1,j)+ &
1509# endif
1510 & ubar(i+1,j,krhs))
1511 cff3=cff1*dndx(i,j)
1512 cff4=cff2*dmde(i,j)
1513# ifdef WEC_VF
1514 cff5=0.5_r8*(vbar_stokes(i,j )+ &
1515 & vbar_stokes(i,j+1))
1516 cff6=0.5_r8*(ubar_stokes(i ,j)+ &
1517 & ubar_stokes(i+1,j))
1518 cff7=cff5*dndx(i,j)
1519 cff8=cff6*dmde(i,j)
1520# endif
1521 cff=drhs(i,j)*(cff3-cff4)
1522 ufx(i,j)=cff*cff1
1523 vfe(i,j)=cff*cff2
1524# ifdef WEC_VF
1525 ufx(i,j)=ufx(i,j)-(cff5*drhs(i,j)*(cff7-cff8))
1526 vfe(i,j)=vfe(i,j)-(cff6*drhs(i,j)*(cff7-cff8))
1527# endif
1528# if defined DIAGNOSTICS_UV
1529 cff=drhs(i,j)*cff4
1530 uwrk(i,j)=-cff*cff1 ! ubar equation, ETA-term
1531 vwrk(i,j)=-cff*cff2 ! vbar equation, ETA-term
1532# ifdef WEC_VF
1533 uwrk(i,j)=uwrk(i,j)+drhs(i,j)*cff5*cff8
1534 vwrk(i,j)=vwrk(i,j)-drhs(i,j)*cff6*cff8
1535# endif
1536# endif
1537 END DO
1538 END DO
1539 DO j=jstr,jend
1540 DO i=istru,iend
1541 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1542 rhs_ubar(i,j)=rhs_ubar(i,j)+fac1
1543# if defined DIAGNOSTICS_UV
1544 fac2=0.5_r8*(uwrk(i,j)+uwrk(i-1,j))
1545 diau2rhs(i,j,m2xadv)=diau2rhs(i,j,m2xadv)+fac1-fac2
1546 diau2rhs(i,j,m2yadv)=diau2rhs(i,j,m2yadv)+fac2
1547 diau2rhs(i,j,m2hadv)=diau2rhs(i,j,m2hadv)+fac1
1548# endif
1549 END DO
1550 END DO
1551 DO j=jstrv,jend
1552 DO i=istr,iend
1553 fac1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1554 rhs_vbar(i,j)=rhs_vbar(i,j)-fac1
1555# if defined DIAGNOSTICS_UV
1556 fac2=0.5_r8*(vwrk(i,j)+vwrk(i,j-1))
1557 diav2rhs(i,j,m2xadv)=diav2rhs(i,j,m2xadv)-fac1+fac2
1558 diav2rhs(i,j,m2yadv)=diav2rhs(i,j,m2yadv)-fac2
1559 diav2rhs(i,j,m2hadv)=diav2rhs(i,j,m2hadv)-fac1
1560# endif
1561 END DO
1562 END DO
1563#endif
1564#if defined UV_VIS2 || defined UV_VIS4
1565!
1566!-----------------------------------------------------------------------
1567! If horizontal mixing, compute total depth at PSI-points.
1568!-----------------------------------------------------------------------
1569!
1570# ifdef UV_VIS4
1571 DO j=jstrm1,jendp2
1572 DO i=istrm1,iendp2
1573# else
1574 DO j=jstr,jend+1
1575 DO i=istr,iend+1
1576# endif
1577 drhs_p(i,j)=0.25_r8*(drhs(i,j )+drhs(i-1,j )+ &
1578 & drhs(i,j-1)+drhs(i-1,j-1))
1579 END DO
1580 END DO
1581#endif
1582#ifdef UV_VIS2
1583!
1584!-----------------------------------------------------------------------
1585! Add in horizontal harmonic viscosity.
1586!-----------------------------------------------------------------------
1587!
1588! Compute flux-components of the horizontal divergence of the stress
1589! tensor (m5/s2) in XI- and ETA-directions.
1590!
1591 DO j=jstrv-1,jend
1592 DO i=istru-1,iend
1593 cff=visc2_r(i,j)*drhs(i,j)*0.5_r8* &
1594 & (pmon_r(i,j)* &
1595 & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,krhs)- &
1596 & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,krhs))- &
1597 & pnom_r(i,j)* &
1598 & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,krhs)- &
1599 & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,krhs)))
1600 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
1601 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
1602 END DO
1603 END DO
1604 DO j=jstr,jend+1
1605 DO i=istr,iend+1
1606 cff=visc2_p(i,j)*drhs_p(i,j)*0.5_r8* &
1607 & (pmon_p(i,j)* &
1608 & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,krhs)- &
1609 & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,krhs))+ &
1610 & pnom_p(i,j)* &
1611 & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,krhs)- &
1612 & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,krhs)))
1613# ifdef MASKING
1614 cff=cff*pmask(i,j)
1615# endif
1616# ifdef WET_DRY
1617 cff=cff*pmask_wet(i,j)
1618# endif
1619 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
1620 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
1621 END DO
1622 END DO
1623!
1624! Add in harmonic viscosity.
1625!
1626 DO j=jstr,jend
1627 DO i=istru,iend
1628 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
1629 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
1630 fac=cff1+cff2
1631 rhs_ubar(i,j)=rhs_ubar(i,j)+fac
1632# if defined DIAGNOSTICS_UV
1633 diau2rhs(i,j,m2hvis)=fac
1634 diau2rhs(i,j,m2xvis)=cff1
1635 diau2rhs(i,j,m2yvis)=cff2
1636# endif
1637 END DO
1638 END DO
1639 DO j=jstrv,jend
1640 DO i=istr,iend
1641 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
1642 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
1643 fac=cff1-cff2
1644 rhs_vbar(i,j)=rhs_vbar(i,j)+fac
1645# if defined DIAGNOSTICS_UV
1646 diav2rhs(i,j,m2hvis)=fac
1647 diav2rhs(i,j,m2xvis)= cff1
1648 diav2rhs(i,j,m2yvis)=-cff2
1649# endif
1650 END DO
1651 END DO
1652#endif
1653#ifdef UV_VIS4
1654!
1655!-----------------------------------------------------------------------
1656! Add in horizontal biharmonic viscosity. The biharmonic operator
1657! is computed by applying the harmonic operator twice.
1658!-----------------------------------------------------------------------
1659!
1660! Compute flux-components of the horizontal divergence of the stress
1661! tensor (m4 s^-3/2) in XI- and ETA-directions. It is assumed here
1662! that "visc4_r" and "visc4_p" are the squared root of the biharmonic
1663! viscosity coefficient. For momentum balance purposes, the total
1664! thickness "D" appears only when computing the second harmonic
1665! operator.
1666!
1667 DO j=jstrvm2,jendp1
1668 DO i=istrum2,iendp1
1669 cff=visc4_r(i,j)*0.5_r8* &
1670 & (pmon_r(i,j)* &
1671 & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,krhs)- &
1672 & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,krhs))- &
1673 & pnom_r(i,j)* &
1674 & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,krhs)- &
1675 & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,krhs)))
1676 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
1677 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
1678 END DO
1679 END DO
1680 DO j=jstrm1,jendp2
1681 DO i=istrm1,iendp2
1682 cff=visc4_p(i,j)*0.5_r8* &
1683 & (pmon_p(i,j)* &
1684 & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,krhs)- &
1685 & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,krhs))+ &
1686 & pnom_p(i,j)* &
1687 & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,krhs)- &
1688 & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,krhs)))
1689# ifdef MASKING
1690 cff=cff*pmask(i,j)
1691# endif
1692# ifdef WET_DRY
1693 cff=cff*pmask_wet(i,j)
1694# endif
1695 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
1696 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
1697 END DO
1698 END DO
1699!
1700! Compute first harmonic operator (m s^-3/2).
1701!
1702 DO j=jstrm1,jendp1
1703 DO i=istrum1,iendp1
1704 lapu(i,j)=0.125_r8* &
1705 & (pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))* &
1706 & ((pn(i-1,j)+pn(i,j))* &
1707 & (ufx(i,j )-ufx(i-1,j))+ &
1708 & (pm(i-1,j)+pm(i,j))* &
1709 & (ufe(i,j+1)-ufe(i ,j)))
1710 END DO
1711 END DO
1712 DO j=jstrvm1,jendp1
1713 DO i=istrm1,iendp1
1714 lapv(i,j)=0.125_r8* &
1715 & (pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))* &
1716 & ((pn(i,j-1)+pn(i,j))* &
1717 & (vfx(i+1,j)-vfx(i,j ))- &
1718 & (pm(i,j-1)+pm(i,j))* &
1719 & (vfe(i ,j)-vfe(i,j-1)))
1720 END DO
1721 END DO
1722!
1723! Apply boundary conditions (other than periodic) to the first
1724! harmonic operator. These are gradient or closed (free slip or
1725! no slip) boundary conditions.
1726!
1727 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1728 IF (domain(ng)%Western_Edge(tile)) THEN
1729 IF (lbc(iwest,isubar,ng)%closed) THEN
1730 DO j=jstrm1,jendp1
1731 lapu(istru-1,j)=0.0_r8
1732 END DO
1733 ELSE
1734 DO j=jstrm1,jendp1
1735 lapu(istru-1,j)=lapu(istru,j)
1736 END DO
1737 END IF
1738 IF (lbc(iwest,isvbar,ng)%closed) THEN
1739 DO j=jstrvm1,jendp1
1740 lapv(istr-1,j)=gamma2(ng)*lapv(istr,j)
1741 END DO
1742 ELSE
1743 DO j=jstrvm1,jendp1
1744 lapv(istr-1,j)=0.0_r8
1745 END DO
1746 END IF
1747 END IF
1748 END IF
1749!
1750 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1751 IF (domain(ng)%Eastern_Edge(tile)) THEN
1752 IF (lbc(ieast,isubar,ng)%closed) THEN
1753 DO j=jstrm1,jendp1
1754 lapu(iend+1,j)=0.0_r8
1755 END DO
1756 ELSE
1757 DO j=jstrm1,jendp1
1758 lapu(iend+1,j)=lapu(iend,j)
1759 END DO
1760 END IF
1761 IF (lbc(ieast,isvbar,ng)%closed) THEN
1762 DO j=jstrvm1,jendp1
1763 lapv(iend+1,j)=gamma2(ng)*lapv(iend,j)
1764 END DO
1765 ELSE
1766 DO j=jstrvm1,jendp1
1767 lapv(iend+1,j)=0.0_r8
1768 END DO
1769 END IF
1770 END IF
1771 END IF
1772!
1773 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1774 IF (domain(ng)%Southern_Edge(tile)) THEN
1775 IF (lbc(isouth,isubar,ng)%closed) THEN
1776 DO i=istrum1,iendp1
1777 lapu(i,jstr-1)=gamma2(ng)*lapu(i,jstr)
1778 END DO
1779 ELSE
1780 DO i=istrum1,iendp1
1781 lapu(i,jstr-1)=0.0_r8
1782 END DO
1783 END IF
1784 IF (lbc(isouth,isvbar,ng)%closed) THEN
1785 DO i=istrm1,iendp1
1786 lapv(i,jstrv-1)=0.0_r8
1787 END DO
1788 ELSE
1789 DO i=istrm1,iendp1
1790 lapv(i,jstrv-1)=lapv(i,jstrv)
1791 END DO
1792 END IF
1793 END IF
1794 END IF
1795!
1796 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1797 IF (domain(ng)%Northern_Edge(tile)) THEN
1798 IF (lbc(inorth,isubar,ng)%closed) THEN
1799 DO i=istrum1,iendp1
1800 lapu(i,jend+1)=gamma2(ng)*lapu(i,jend)
1801 END DO
1802 ELSE
1803 DO i=istrum1,iendp1
1804 lapu(i,jend+1)=0.0_r8
1805 END DO
1806 END IF
1807 IF (lbc(inorth,isvbar,ng)%closed) THEN
1808 DO i=istrm1,iendp1
1809 lapv(i,jend+1)=0.0_r8
1810 END DO
1811 ELSE
1812 DO i=istrm1,iendp1
1813 lapv(i,jend+1)=lapv(i,jend)
1814 END DO
1815 END IF
1816 END IF
1817 END IF
1818!
1819 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
1820 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
1821 IF (domain(ng)%SouthWest_Corner(tile)) THEN
1822 lapu(istr ,jstr-1)=0.5_r8*(lapu(istr+1,jstr-1)+ &
1823 & lapu(istr ,jstr ))
1824 lapv(istr-1,jstr )=0.5_r8*(lapv(istr-1,jstr+1)+ &
1825 & lapv(istr ,jstr ))
1826 END IF
1827 END IF
1828
1829 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
1830 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
1831 IF (domain(ng)%SouthEast_Corner(tile)) THEN
1832 lapu(iend+1,jstr-1)=0.5_r8*(lapu(iend ,jstr-1)+ &
1833 & lapu(iend+1,jstr ))
1834 lapv(iend+1,jstr )=0.5_r8*(lapv(iend ,jstr )+ &
1835 & lapv(iend+1,jstr+1))
1836 END IF
1837 END IF
1838
1839 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
1840 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
1841 IF (domain(ng)%NorthWest_Corner(tile)) THEN
1842 lapu(istr ,jend+1)=0.5_r8*(lapu(istr+1,jend+1)+ &
1843 & lapu(istr ,jend ))
1844 lapv(istr-1,jend+1)=0.5_r8*(lapv(istr ,jend+1)+ &
1845 & lapv(istr-1,jend ))
1846 END IF
1847 END IF
1848
1849 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
1850 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
1851 IF (domain(ng)%NorthEast_Corner(tile)) THEN
1852 lapu(iend+1,jend+1)=0.5_r8*(lapu(iend ,jend+1)+ &
1853 & lapu(iend+1,jend ))
1854 lapv(iend+1,jend+1)=0.5_r8*(lapv(iend ,jend+1)+ &
1855 & lapv(iend+1,jend ))
1856 END IF
1857 END IF
1858!
1859! Compute flux-components of the horizontal divergence of the
1860! biharmonic stress tensor (m4/s2) in XI- and ETA-directions.
1861!
1862 DO j=jstrv-1,jend
1863 DO i=istru-1,iend
1864 cff=visc4_r(i,j)*drhs(i,j)*0.5_r8* &
1865 & (pmon_r(i,j)* &
1866 & ((pn(i ,j)+pn(i+1,j))*lapu(i+1,j)- &
1867 & (pn(i-1,j)+pn(i ,j))*lapu(i ,j))- &
1868 & pnom_r(i,j)* &
1869 & ((pm(i,j )+pm(i,j+1))*lapv(i,j+1)- &
1870 & (pm(i,j-1)+pm(i,j ))*lapv(i,j )))
1871 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
1872 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
1873 END DO
1874 END DO
1875 DO j=jstr,jend+1
1876 DO i=istr,iend+1
1877 cff=visc4_p(i,j)*drhs_p(i,j)*0.5_r8* &
1878 & (pmon_p(i,j)* &
1879 & ((pn(i ,j-1)+pn(i ,j))*lapv(i ,j)- &
1880 & (pn(i-1,j-1)+pn(i-1,j))*lapv(i-1,j))+ &
1881 & pnom_p(i,j)* &
1882 & ((pm(i-1,j )+pm(i,j ))*lapu(i,j )- &
1883 & (pm(i-1,j-1)+pm(i,j-1))*lapu(i,j-1)))
1884# ifdef MASKING
1885 cff=cff*pmask(i,j)
1886# endif
1887# ifdef WET_DRY
1888 cff=cff*pmask_wet(i,j)
1889# endif
1890 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
1891 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
1892 END DO
1893 END DO
1894!
1895! Add in biharmonic viscosity.
1896!
1897 DO j=jstr,jend
1898 DO i=istru,iend
1899 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
1900 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
1901 fac=cff1+cff2
1902 rhs_ubar(i,j)=rhs_ubar(i,j)-fac
1903# if defined DIAGNOSTICS_UV
1904 diau2rhs(i,j,m2hvis)=-fac
1905 diau2rhs(i,j,m2xvis)=-cff1
1906 diau2rhs(i,j,m2yvis)=-cff2
1907# endif
1908 END DO
1909 END DO
1910 DO j=jstrv,jend
1911 DO i=istr,iend
1912 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
1913 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
1914 fac=cff1-cff2
1915 rhs_vbar(i,j)=rhs_vbar(i,j)-fac
1916# if defined DIAGNOSTICS_UV
1917 diav2rhs(i,j,m2hvis)=-fac
1918 diav2rhs(i,j,m2xvis)=-cff1
1919 diav2rhs(i,j,m2yvis)= cff2
1920# endif
1921 END DO
1922 END DO
1923#endif
1924#if defined WEC_VF && defined SOLVE3D
1925!
1926!-----------------------------------------------------------------------
1927! Add in non-conservative roller terms.
1928!-----------------------------------------------------------------------
1929!
1930 DO j=jstr,jend
1931 DO i=istru,iend
1932 cff1=rubrk2d(i,j)
1933 rhs_ubar(i,j)=rhs_ubar(i,j)+cff1
1934# ifdef DIAGNOSTICS_UV
1935 diau2rhs(i,j,m2wbrk)=cff1
1936# endif
1937# ifdef WEC_ROLLER
1938 cff1=rurol2d(i,j)
1939 rhs_ubar(i,j)=rhs_ubar(i,j)+cff1
1940# ifdef DIAGNOSTICS_UV
1941 diau2rhs(i,j,m2wrol)=cff1
1942# endif
1943# else
1944# ifdef DIAGNOSTICS_UV
1945 diau2rhs(i,j,m2wrol)=0.0_r8
1946# endif
1947# endif
1948# ifdef BOTTOM_STREAMING
1949 cff1=rubst2d(i,j)
1950 rhs_ubar(i,j)=rhs_ubar(i,j)+cff1
1951# ifdef DIAGNOSTICS_UV
1952 diau2rhs(i,j,m2bstm)=cff1
1953# endif
1954# endif
1955# ifdef SURFACE_STREAMING
1956 cff1=russt2d(i,j)
1957 rhs_ubar(i,j)=rhs_ubar(i,j)+cff1
1958# ifdef DIAGNOSTICS_UV
1959 diau2rhs(i,j,m2sstm)=cff1
1960# endif
1961# endif
1962 END DO
1963 IF (j.ge.jstrv) THEN
1964 DO i=istr,iend
1965 cff1=rvbrk2d(i,j)
1966 rhs_vbar(i,j)=rhs_vbar(i,j)+cff1
1967# ifdef DIAGNOSTICS_UV
1968 diav2rhs(i,j,m2wbrk)=cff1
1969# endif
1970# ifdef WEC_ROLLER
1971 cff1=rvrol2d(i,j)
1972 rhs_vbar(i,j)=rhs_vbar(i,j)+cff1
1973# ifdef DIAGNOSTICS_UV
1974 diav2rhs(i,j,m2wrol)=cff1
1975# endif
1976# else
1977# ifdef DIAGNOSTICS_UV
1978 diav2rhs(i,j,m2wrol)=0.0_r8
1979# endif
1980# endif
1981# ifdef BOTTOM_STREAMING
1982 cff1=rvbst2d(i,j)
1983 rhs_vbar(i,j)=rhs_vbar(i,j)+cff1
1984# ifdef DIAGNOSTICS_UV
1985 diav2rhs(i,j,m2bstm)=cff1
1986# endif
1987# endif
1988# ifdef SURFACE_STREAMING
1989 cff1=rvsst2d(i,j)
1990 rhs_vbar(i,j)=rhs_vbar(i,j)+cff1
1991# ifdef DIAGNOSTICS_UV
1992 diav2rhs(i,j,m2sstm)=cff1
1993# endif
1994# endif
1995 END DO
1996 END IF
1997 END DO
1998# ifdef UV_ADV
1999# ifdef DIAGNOSTICS_UV
2000!
2001!---------------------------------------------------------------------------
2002! To obtain the full horizotal 'J' vortex force term:
2003! Compute term for diagnostics only. Subtract from hadv and add to vorf.
2004!---------------------------------------------------------------------------
2005!
2006 DO j=jstr,jend
2007 DO i=istru,iend
2008 cff=0.5_r8*(drhs(i-1,j)+drhs(i,j))
2009 dvsom(i,j)=0.25_r8*cff*om_u(i,j)* &
2010 & (vbar_stokes(i ,j )+ &
2011 & vbar_stokes(i ,j+1)+ &
2012 & vbar_stokes(i-1,j )+ &
2013 & vbar_stokes(i-1,j+1))
2014 END DO
2015 END DO
2016 DO j=jstr,jend+1
2017 DO i=istru,iend
2018 ufx(i,j)=0.5_r8*(ubar(i ,j-1,krhs)+ &
2019 ubar(i ,j ,krhs))
2020 END DO
2021 END DO
2022 DO j=jstr,jend
2023 DO i=istru,iend
2024 cff1=ufx(i,j+1)-ufx(i,j)
2025 cff=cff1*dvsom(i,j)
2026 diau2rhs(i,j,m2xadv)=diau2rhs(i,j,m2xadv)+cff
2027 diau2rhs(i,j,m2hadv)=diau2rhs(i,j,m2hadv)+cff
2028 diau2rhs(i,j,m2hjvf)=-cff
2029 END DO
2030 END DO
2031 DO j=jstrv,jend
2032 DO i=istr,iend
2033 cff=0.5_r8*(drhs(i,j)+drhs(i,j-1))
2034 duson(i,j)=cff*0.25_r8*on_v(i,j)* &
2035 & (ubar_stokes(i ,j )+ &
2036 & ubar_stokes(i+1,j )+ &
2037 & ubar_stokes(i ,j-1)+ &
2038 & ubar_stokes(i+1,j-1))
2039 END DO
2040 END DO
2041 DO j=jstrv,jend
2042 DO i=istr,iend+1
2043 vfe(i,j)=0.5_r8*(vbar(i-1,j ,krhs)+ &
2044 & vbar(i ,j ,krhs))
2045 END DO
2046 END DO
2047 DO i=istr,iend
2048 DO j=jstrv,jend
2049 cff2=vfe(i+1,j)-vfe(i,j)
2050 cff=cff2*duson(i,j)
2051 diav2rhs(i,j,m2yadv)=diav2rhs(i,j,m2yadv)+cff
2052 diav2rhs(i,j,m2hadv)=diav2rhs(i,j,m2hadv)+cff
2053 diav2rhs(i,j,m2hjvf)=-cff
2054 END DO
2055 END DO
2056# endif
2057!
2058!---------------------------------------------------------------------------
2059! Contribution of a term corresponding to product of
2060! Stokes and Eulerian Velocity Eqn. 26 and 27.
2061! This removes terms that were unneccessarily added in flux form.
2062!---------------------------------------------------------------------------
2063!
2064 DO j=jstr,jend
2065 DO i=istru,iend
2066 cff=0.5_r8*(drhs(i-1,j)+drhs(i,j))
2067 duson(i,j)=cff*on_u(i,j)*ubar_stokes(i,j)
2068 dvson(i,j)=0.25_r8*cff*on_u(i,j)* &
2069 & (vbar_stokes(i ,j )+ &
2070 & vbar_stokes(i ,j+1)+ &
2071 & vbar_stokes(i-1,j )+ &
2072 & vbar_stokes(i-1,j+1))
2073 END DO
2074 DO i=istru-1,iend
2075 ufx(i,j)=0.5_r8*(ubar(i ,j ,krhs)+ &
2076 ubar(i+1,j ,krhs))
2077 vfx(i,j)=0.5_r8*(vbar(i ,j ,krhs)+ &
2078 & vbar(i ,j+1,krhs))
2079 END DO
2080 END DO
2081 DO j=jstrv,jend
2082 DO i=istr,iend
2083 cff=0.5_r8*(drhs(i,j)+drhs(i,j-1))
2084 dusom(i,j)=cff*0.25_r8*om_v(i,j)* &
2085 & (ubar_stokes(i ,j )+ &
2086 & ubar_stokes(i+1,j )+ &
2087 & ubar_stokes(i ,j-1)+ &
2088 & ubar_stokes(i+1,j-1))
2089 dvsom(i,j)=cff*om_v(i,j)*vbar_stokes(i,j)
2090 END DO
2091 END DO
2092 DO j=jstrv-1,jend
2093 DO i=istr,iend
2094 cff=0.5_r8*(drhs(i,j)+drhs(i,j-1))
2095 ufe(i,j)=0.5_r8*(ubar(i+1,j ,krhs)+ &
2096 & ubar(i ,j ,krhs))
2097 vfe(i,j)=0.5_r8*(vbar(i ,j ,krhs)+ &
2098 & vbar(i ,j+1,krhs))
2099 END DO
2100 END DO
2101 DO j=jstr,jend
2102 DO i=istru,iend
2103 cff1=ufx(i,j)-ufx(i-1,j)
2104 cff2=vfx(i,j)-vfx(i-1,j)
2105 cff3=duson(i,j)*cff1
2106 cff4=dvson(i,j)*cff2
2107 rhs_ubar(i,j)=rhs_ubar(i,j)+cff3+cff4
2108! rustr2d(i,j)=rustr2d(i,j)-cff3-cff4
2109# ifdef DIAGNOSTICS_UV
2110 diau2rhs(i,j,m2xadv)=diau2rhs(i,j,m2xadv)+cff3
2111 diau2rhs(i,j,m2hadv)=diau2rhs(i,j,m2hadv)+cff3
2112 diau2rhs(i,j,m2hjvf)=diau2rhs(i,j,m2hjvf)+cff4
2113# endif
2114 END DO
2115 END DO
2116 DO i=istr,iend
2117 DO j=jstrv,jend
2118 cff1=ufe(i,j)-ufe(i,j-1)
2119 cff2=vfe(i,j)-vfe(i,j-1)
2120 cff3=dusom(i,j)*cff1
2121 cff4=dvsom(i,j)*cff2
2122 rhs_vbar(i,j)=rhs_vbar(i,j)+cff3+cff4
2123! rvstr2d(i,j)=rvstr2d(i,j,k)-cff3-cff4
2124# ifdef DIAGNOSTICS_UV
2125 diav2rhs(i,j,m2yadv)=diav2rhs(i,j,m2yadv)+cff4
2126 diav2rhs(i,j,m2hadv)=diav2rhs(i,j,m2hadv)+cff4
2127 diav2rhs(i,j,m2hjvf)=diav2rhs(i,j,m2hjvf)+cff3
2128# endif
2129 END DO
2130 END DO
2131# endif
2132#endif
2133#ifndef SOLVE3D
2134!
2135!-----------------------------------------------------------------------
2136! Add in bottom stress.
2137!-----------------------------------------------------------------------
2138!
2139 DO j=jstr,jend
2140 DO i=istru,iend
2141 fac=bustr(i,j)*om_u(i,j)*on_u(i,j)
2142 rhs_ubar(i,j)=rhs_ubar(i,j)-fac
2143# ifdef DIAGNOSTICS_UV
2144 diau2rhs(i,j,m2bstr)=-fac
2145# endif
2146 END DO
2147 END DO
2148 DO j=jstrv,jend
2149 DO i=istr,iend
2150 fac=bvstr(i,j)*om_v(i,j)*on_v(i,j)
2151 rhs_vbar(i,j)=rhs_vbar(i,j)-fac
2152# ifdef DIAGNOSTICS_UV
2153 diav2rhs(i,j,m2bstr)=-fac
2154# endif
2155 END DO
2156 END DO
2157#else
2158# ifdef DIAGNOSTICS_UV
2159!
2160! Initialize the stress term if no bottom friction is defined.
2161!
2162 DO j=jstr,jend
2163 DO i=istru,iend
2164 diau2rhs(i,j,m2bstr)=0.0_r8
2165 END DO
2166 END DO
2167 DO j=jstrv,jend
2168 DO i=istr,iend
2169 diav2rhs(i,j,m2bstr)=0.0_r8
2170 END DO
2171 END DO
2172# endif
2173#endif
2174!
2175!-----------------------------------------------------------------------
2176! Add in nudging of 2D momentum climatology.
2177!-----------------------------------------------------------------------
2178!
2179 IF (lnudgem2clm(ng)) THEN
2180 DO j=jstr,jend
2181 DO i=istru,iend
2182 cff=0.25_r8*(clima(ng)%M2nudgcof(i-1,j)+ &
2183 & clima(ng)%M2nudgcof(i ,j))* &
2184 & om_u(i,j)*on_u(i,j)
2185 rhs_ubar(i,j)=rhs_ubar(i,j)+ &
2186 & cff*(drhs(i-1,j)+drhs(i,j))* &
2187 & (clima(ng)%ubarclm(i,j)- &
2188 & ubar(i,j,krhs))
2189 END DO
2190 END DO
2191 DO j=jstrv,jend
2192 DO i=istr,iend
2193 cff=0.25_r8*(clima(ng)%M2nudgcof(i,j-1)+ &
2194 & clima(ng)%M2nudgcof(i,j ))* &
2195 & om_v(i,j)*on_v(i,j)
2196 rhs_vbar(i,j)=rhs_vbar(i,j)+ &
2197 & cff*(drhs(i,j-1)+drhs(i,j))* &
2198 & (clima(ng)%vbarclm(i,j)- &
2199 & vbar(i,j,krhs))
2200 END DO
2201 END DO
2202 END IF
2203
2204#ifdef SOLVE3D
2205# ifdef WET_DRY
2206 DO j=jstr,jend
2207 DO i=istru,iend
2208 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
2209 cff6=0.5_r8+dsign(0.5_r8,rhs_ubar(i,j))*umask_wet(i,j)
2210 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2211 rhs_ubar(i,j)=rhs_ubar(i,j)*cff7
2212 END DO
2213 END DO
2214 DO j=jstrv,jend
2215 DO i=istr,iend
2216 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
2217 cff6=0.5_r8+dsign(0.5_r8,rhs_vbar(i,j))*vmask_wet(i,j)
2218 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2219 rhs_vbar(i,j)=rhs_vbar(i,j)*cff7
2220 END DO
2221 END DO
2222# endif
2223!
2224!-----------------------------------------------------------------------
2225! Coupling between 2D and 3D equations.
2226!-----------------------------------------------------------------------
2227!
2228! Before the predictor step of the first barotropic time-step,
2229! arrays "rufrc" and "rvfrc" contain the vertical integrals of
2230! the 3D right-hand-side terms for momentum equations (including
2231! surface and bottom stresses, if so prescribed).
2232!
2233! Convert them into forcing terms by subtracting the fast time
2234! "rhs_ubar" and "rhs_vbar" from them; Also, immediately apply
2235! these forcing terms "rhs_ubar" and "rhs_vbar".
2236!
2237! From now on, these newly computed forcing terms will remain
2238! constant during the fast time stepping and will added to
2239! "rhs_ubar" and "rhs_vbar" during all subsequent time steps.
2240!
2241 IF (first_2d_step.and.predictor_2d_step(ng)) THEN
2242 IF (iic(ng).eq.ntfirst(ng)) THEN
2243 DO j=jstr,jend
2244 DO i=istru,iend
2245 rufrc(i,j)=rufrc(i,j)-rhs_ubar(i,j)
2246 rhs_ubar(i,j)=rhs_ubar(i,j)+rufrc(i,j)
2247 ru(i,j,0,nstp)=rufrc(i,j)
2248# ifdef DIAGNOSTICS_UV
2249 DO idiag=1,m2pgrd
2250 diarufrc(i,j,3,idiag)=diarufrc(i,j,3,idiag)- &
2251 & diau2rhs(i,j,idiag)
2252 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)+ &
2253 & diarufrc(i,j,3,idiag)
2254 diarufrc(i,j,nstp,idiag)=diarufrc(i,j,3,idiag)
2255 END DO
2256 diau2rhs(i,j,m2sstr)=diarufrc(i,j,3,m2sstr)
2257 diarufrc(i,j,nstp,m2sstr)=diarufrc(i,j,3,m2sstr)
2258 diau2rhs(i,j,m2bstr)=diarufrc(i,j,3,m2bstr)
2259 diarufrc(i,j,nstp,m2bstr)=diarufrc(i,j,3,m2bstr)
2260# ifdef WEC_VF
2261 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2zeta)+ &
2262 & diarufrc(i,j,3,m2pgrd)
2263# endif
2264# endif
2265 END DO
2266 END DO
2267 DO j=jstrv,jend
2268 DO i=istr,iend
2269 rvfrc(i,j)=rvfrc(i,j)-rhs_vbar(i,j)
2270 rhs_vbar(i,j)=rhs_vbar(i,j)+rvfrc(i,j)
2271 rv(i,j,0,nstp)=rvfrc(i,j)
2272# ifdef DIAGNOSTICS_UV
2273 DO idiag=1,m2pgrd
2274 diarvfrc(i,j,3,idiag)=diarvfrc(i,j,3,idiag)- &
2275 & diav2rhs(i,j,idiag)
2276 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)+ &
2277 & diarvfrc(i,j,3,idiag)
2278 diarvfrc(i,j,nstp,idiag)=diarvfrc(i,j,3,idiag)
2279 END DO
2280 diav2rhs(i,j,m2sstr)=diarvfrc(i,j,3,m2sstr)
2281 diarvfrc(i,j,nstp,m2sstr)=diarvfrc(i,j,3,m2sstr)
2282 diav2rhs(i,j,m2bstr)=diarvfrc(i,j,3,m2bstr)
2283 diarvfrc(i,j,nstp,m2bstr)=diarvfrc(i,j,3,m2bstr)
2284# ifdef WEC_VF
2285 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2zeta)+ &
2286 & diarvfrc(i,j,3,m2pgrd)
2287# endif
2288# endif
2289 END DO
2290 END DO
2291 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
2292 DO j=jstr,jend
2293 DO i=istru,iend
2294 rufrc(i,j)=rufrc(i,j)-rhs_ubar(i,j)
2295 rhs_ubar(i,j)=rhs_ubar(i,j)+ &
2296 & 1.5_r8*rufrc(i,j)-0.5_r8*ru(i,j,0,nnew)
2297 ru(i,j,0,nstp)=rufrc(i,j)
2298# ifdef DIAGNOSTICS_UV
2299 DO idiag=1,m2pgrd
2300 diarufrc(i,j,3,idiag)=diarufrc(i,j,3,idiag)- &
2301 & diau2rhs(i,j,idiag)
2302 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)+ &
2303 & 1.5_r8*diarufrc(i,j,3,idiag)- &
2304 & 0.5_r8*diarufrc(i,j,nnew,idiag)
2305 diarufrc(i,j,nstp,idiag)=diarufrc(i,j,3,idiag)
2306 END DO
2307 diau2rhs(i,j,m2sstr)=1.5_r8*diarufrc(i,j,3,m2sstr)- &
2308 & 0.5_r8*diarufrc(i,j,nnew,m2sstr)
2309 diarufrc(i,j,nstp,m2sstr)=diarufrc(i,j,3,m2sstr)
2310 diau2rhs(i,j,m2bstr)=1.5_r8*diarufrc(i,j,3,m2bstr)- &
2311 & 0.5_r8*diarufrc(i,j,nnew,m2bstr)
2312 diarufrc(i,j,nstp,m2bstr)=diarufrc(i,j,3,m2bstr)
2313# ifdef WEC_VF
2314 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2zeta)+ &
2315 & 1.5_r8*diarufrc(i,j,3,m2pgrd)- &
2316 & 0.5_r8*diarufrc(i,j,nnew,m2pgrd)
2317# endif
2318# endif
2319 END DO
2320 END DO
2321 DO j=jstrv,jend
2322 DO i=istr,iend
2323 rvfrc(i,j)=rvfrc(i,j)-rhs_vbar(i,j)
2324 rhs_vbar(i,j)=rhs_vbar(i,j)+ &
2325 & 1.5_r8*rvfrc(i,j)-0.5_r8*rv(i,j,0,nnew)
2326 rv(i,j,0,nstp)=rvfrc(i,j)
2327# ifdef DIAGNOSTICS_UV
2328 DO idiag=1,m2pgrd
2329 diarvfrc(i,j,3,idiag)=diarvfrc(i,j,3,idiag)- &
2330 & diav2rhs(i,j,idiag)
2331 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)+ &
2332 & 1.5_r8*diarvfrc(i,j,3,idiag)- &
2333 & 0.5_r8*diarvfrc(i,j,nnew,idiag)
2334 diarvfrc(i,j,nstp,idiag)=diarvfrc(i,j,3,idiag)
2335 END DO
2336 diav2rhs(i,j,m2sstr)=1.5_r8*diarvfrc(i,j,3,m2sstr)- &
2337 & 0.5_r8*diarvfrc(i,j,nnew,m2sstr)
2338 diarvfrc(i,j,nstp,m2sstr)=diarvfrc(i,j,3,m2sstr)
2339 diav2rhs(i,j,m2bstr)=1.5_r8*diarvfrc(i,j,3,m2bstr)- &
2340 & 0.5_r8*diarvfrc(i,j,nnew,m2bstr)
2341 diarvfrc(i,j,nstp,m2bstr)=diarvfrc(i,j,3,m2bstr)
2342# ifdef WEC_VF
2343 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2zeta)+ &
2344 & 1.5_r8*diarvfrc(i,j,3,m2pgrd)- &
2345 & 0.5_r8*diarvfrc(i,j,nnew,m2pgrd)
2346# endif
2347# endif
2348 END DO
2349 END DO
2350 ELSE
2351 cff1=23.0_r8/12.0_r8
2352 cff2=16.0_r8/12.0_r8
2353 cff3= 5.0_r8/12.0_r8
2354 DO j=jstr,jend
2355 DO i=istru,iend
2356 rufrc(i,j)=rufrc(i,j)-rhs_ubar(i,j)
2357 rhs_ubar(i,j)=rhs_ubar(i,j)+ &
2358 & cff1*rufrc(i,j)- &
2359 & cff2*ru(i,j,0,nnew)+ &
2360 & cff3*ru(i,j,0,nstp)
2361 ru(i,j,0,nstp)=rufrc(i,j)
2362# ifdef DIAGNOSTICS_UV
2363 DO idiag=1,m2pgrd
2364 diarufrc(i,j,3,idiag)=diarufrc(i,j,3,idiag)- &
2365 & diau2rhs(i,j,idiag)
2366 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)+ &
2367 & cff1*diarufrc(i,j,3,idiag)- &
2368 & cff2*diarufrc(i,j,nnew,idiag)+ &
2369 & cff3*diarufrc(i,j,nstp,idiag)
2370 diarufrc(i,j,nstp,idiag)=diarufrc(i,j,3,idiag)
2371 END DO
2372 diau2rhs(i,j,m2sstr)=cff1*diarufrc(i,j,3,m2sstr)- &
2373 & cff2*diarufrc(i,j,nnew,m2sstr)+ &
2374 & cff3*diarufrc(i,j,nstp,m2sstr)
2375 diarufrc(i,j,nstp,m2sstr)=diarufrc(i,j,3,m2sstr)
2376 diau2rhs(i,j,m2bstr)=cff1*diarufrc(i,j,3,m2bstr)- &
2377 & cff2*diarufrc(i,j,nnew,m2bstr)+ &
2378 & cff3*diarufrc(i,j,nstp,m2bstr)
2379 diarufrc(i,j,nstp,m2bstr)=diarufrc(i,j,3,m2bstr)
2380# ifdef WEC_VF
2381 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2zeta)+ &
2382 & cff1*diarufrc(i,j,3,m2pgrd)- &
2383 & cff2*diarufrc(i,j,nnew,m2pgrd)+ &
2384 & cff3*diarufrc(i,j,nstp,m2pgrd)
2385# endif
2386# endif
2387 END DO
2388 END DO
2389 DO j=jstrv,jend
2390 DO i=istr,iend
2391 rvfrc(i,j)=rvfrc(i,j)-rhs_vbar(i,j)
2392 rhs_vbar(i,j)=rhs_vbar(i,j)+ &
2393 & cff1*rvfrc(i,j)- &
2394 & cff2*rv(i,j,0,nnew)+ &
2395 & cff3*rv(i,j,0,nstp)
2396 rv(i,j,0,nstp)=rvfrc(i,j)
2397# ifdef DIAGNOSTICS_UV
2398 DO idiag=1,m2pgrd
2399 diarvfrc(i,j,3,idiag)=diarvfrc(i,j,3,idiag)- &
2400 & diav2rhs(i,j,idiag)
2401 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)+ &
2402 & cff1*diarvfrc(i,j,3,idiag)- &
2403 & cff2*diarvfrc(i,j,nnew,idiag)+ &
2404 & cff3*diarvfrc(i,j,nstp,idiag)
2405 diarvfrc(i,j,nstp,idiag)=diarvfrc(i,j,3,idiag)
2406 END DO
2407 diav2rhs(i,j,m2sstr)=cff1*diarvfrc(i,j,3,m2sstr)- &
2408 & cff2*diarvfrc(i,j,nnew,m2sstr)+ &
2409 & cff3*diarvfrc(i,j,nstp,m2sstr)
2410 diarvfrc(i,j,nstp,m2sstr)=diarvfrc(i,j,3,m2sstr)
2411 diav2rhs(i,j,m2bstr)=cff1*diarvfrc(i,j,3,m2bstr)- &
2412 & cff2*diarvfrc(i,j,nnew,m2bstr)+ &
2413 & cff3*diarvfrc(i,j,nstp,m2bstr)
2414 diarvfrc(i,j,nstp,m2bstr)=diarvfrc(i,j,3,m2bstr)
2415# ifdef WEC_VF
2416 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2zeta)+ &
2417 & cff1*diarvfrc(i,j,3,m2pgrd)- &
2418 & cff2*diarvfrc(i,j,nnew,m2pgrd)+ &
2419 & cff3*diarvfrc(i,j,nstp,m2pgrd)
2420# endif
2421# endif
2422 END DO
2423 END DO
2424 END IF
2425 ELSE
2426 DO j=jstr,jend
2427 DO i=istru,iend
2428 rhs_ubar(i,j)=rhs_ubar(i,j)+rufrc(i,j)
2429# ifdef DIAGNOSTICS_UV
2430 DO idiag=1,m2pgrd
2431 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)+ &
2432 & diarufrc(i,j,3,idiag)
2433 END DO
2434 diau2rhs(i,j,m2sstr)=diarufrc(i,j,3,m2sstr)
2435 diau2rhs(i,j,m2bstr)=diarufrc(i,j,3,m2bstr)
2436# ifdef WEC_VF
2437 diau2rhs(i,j,m2zeta)=diau2rhs(i,j,m2zeta)+ &
2438 & diarufrc(i,j,3,m2pgrd)
2439# endif
2440# endif
2441 END DO
2442 END DO
2443 DO j=jstrv,jend
2444 DO i=istr,iend
2445 rhs_vbar(i,j)=rhs_vbar(i,j)+rvfrc(i,j)
2446# ifdef DIAGNOSTICS_UV
2447 DO idiag=1,m2pgrd
2448 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)+ &
2449 & diarvfrc(i,j,3,idiag)
2450 END DO
2451 diav2rhs(i,j,m2sstr)=diarvfrc(i,j,3,m2sstr)
2452 diav2rhs(i,j,m2bstr)=diarvfrc(i,j,3,m2bstr)
2453# ifdef WEC_VF
2454 diav2rhs(i,j,m2zeta)=diav2rhs(i,j,m2zeta)+ &
2455 & diarvfrc(i,j,3,m2pgrd)
2456# endif
2457# endif
2458 END DO
2459 END DO
2460 END IF
2461#else
2462!
2463!-----------------------------------------------------------------------
2464! Add in surface momentum stress.
2465!-----------------------------------------------------------------------
2466!
2467 DO j=jstr,jend
2468 DO i=istru,iend
2469 fac=sustr(i,j)*om_u(i,j)*on_u(i,j)
2470 rhs_ubar(i,j)=rhs_ubar(i,j)+fac
2471# ifdef DIAGNOSTICS_UV
2472 diau2rhs(i,j,m2sstr)=fac
2473# endif
2474 END DO
2475 END DO
2476 DO j=jstrv,jend
2477 DO i=istr,iend
2478 fac=svstr(i,j)*om_v(i,j)*on_v(i,j)
2479 rhs_vbar(i,j)=rhs_vbar(i,j)+fac
2480# ifdef DIAGNOSTICS_UV
2481 diav2rhs(i,j,m2sstr)=fac
2482# endif
2483 END DO
2484 END DO
2485#endif
2486!
2487!=======================================================================
2488! Time step 2D momentum equations.
2489!=======================================================================
2490!
2491! Compute total water column depth.
2492!
2493 DO j=jstrv-1,jend
2494 DO i=istru-1,iend
2495 dstp(i,j)=zeta(i,j,kstp)+h(i,j)
2496 END DO
2497 END DO
2498!
2499! During the first time-step, the predictor step is Forward-Euler
2500! and the corrector step is Backward-Euler. Otherwise, the predictor
2501! step is Leap-frog and the corrector step is Adams-Moulton.
2502!
2503 IF (first_2d_step) THEN
2504 cff1=0.5_r8*dtfast(ng)
2505#ifdef WET_DRY
2506 cff2=1.0_r8/cff1
2507#endif
2508 DO j=jstr,jend
2509 DO i=istru,iend
2510 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2511 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2512 ubar(i,j,knew)=(ubar(i,j,kstp)* &
2513 & (dstp(i,j)+dstp(i-1,j))+ &
2514 & cff*cff1*rhs_ubar(i,j))*fac
2515#ifdef MASKING
2516 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2517#endif
2518#ifdef WET_DRY
2519 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
2520 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2521 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2522 ubar(i,j,knew)=ubar(i,j,knew)*cff7
2523 rhs_ubar(i,j)=rhs_ubar(i,j)*cff7
2524# ifdef SOLVE3D
2525 IF (predictor_2d_step(ng)) THEN
2526 rufrc(i,j)=rufrc(i,j)*cff7
2527 ru(i,j,0,nstp)=rufrc(i,j)
2528 END IF
2529# endif
2530#endif
2531#if defined NESTING && !defined SOLVE3D
2532 du_flux(i,j)=ubar(i,j,knew)* &
2533 & 0.5_r8*(dnew(i,j)+dnew(i-1,j))*on_u(i,j)
2534#endif
2535 END DO
2536 END DO
2537 DO j=jstrv,jend
2538 DO i=istr,iend
2539 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2540 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2541 vbar(i,j,knew)=(vbar(i,j,kstp)* &
2542 & (dstp(i,j)+dstp(i,j-1))+ &
2543 & cff*cff1*rhs_vbar(i,j))*fac
2544#ifdef MASKING
2545 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2546#endif
2547#ifdef WET_DRY
2548 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
2549 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
2550 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2551 vbar(i,j,knew)=vbar(i,j,knew)*cff7
2552 rhs_vbar(i,j)=rhs_vbar(i,j)*cff7
2553# ifdef SOLVE3D
2554 IF (predictor_2d_step(ng)) THEN
2555 rvfrc(i,j)=rvfrc(i,j)*cff7
2556 rv(i,j,0,nstp)=rvfrc(i,j)
2557 END IF
2558# endif
2559#endif
2560#if defined NESTING && !defined SOLVE3D
2561 dv_flux(i,j)=vbar(i,j,knew)* &
2562 & 0.5_r8*(dnew(i,j)+dnew(i,j-1))*om_v(i,j)
2563#endif
2564 END DO
2565 END DO
2566 ELSE IF (predictor_2d_step(ng)) THEN
2567 cff1=dtfast(ng)
2568#ifdef WET_DRY
2569 cff2=1.0_r8/cff1
2570#endif
2571 DO j=jstr,jend
2572 DO i=istru,iend
2573 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2574 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2575 ubar(i,j,knew)=(ubar(i,j,kstp)* &
2576 & (dstp(i,j)+dstp(i-1,j))+ &
2577 & cff*cff1*rhs_ubar(i,j))*fac
2578#ifdef MASKING
2579 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2580#endif
2581#ifdef WET_DRY
2582 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
2583 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2584 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2585 ubar(i,j,knew)=ubar(i,j,knew)*cff7
2586 rhs_ubar(i,j)=rhs_ubar(i,j)*cff7
2587#endif
2588#if defined NESTING && !defined SOLVE3D
2589 du_flux(i,j)=ubar(i,j,knew)* &
2590 & 0.5_r8*(dnew(i,j)+dnew(i-1,j))*on_u(i,j)
2591#endif
2592 END DO
2593 END DO
2594 DO j=jstrv,jend
2595 DO i=istr,iend
2596 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2597 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2598 vbar(i,j,knew)=(vbar(i,j,kstp)* &
2599 & (dstp(i,j)+dstp(i,j-1))+ &
2600 & cff*cff1*rhs_vbar(i,j))*fac
2601#ifdef MASKING
2602 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2603#endif
2604#ifdef WET_DRY
2605 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
2606 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
2607 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2608 vbar(i,j,knew)=vbar(i,j,knew)*cff7
2609 rhs_vbar(i,j)=rhs_vbar(i,j)*cff7
2610#endif
2611#if defined NESTING && !defined SOLVE3D
2612 dv_flux(i,j)=vbar(i,j,knew)* &
2613 & 0.5_r8*(dnew(i,j)+dnew(i,j-1))*om_v(i,j)
2614#endif
2615 END DO
2616 END DO
2617 ELSE IF (corrector_2d_step) THEN
2618 cff1=0.5_r8*dtfast(ng)*5.0_r8/12.0_r8
2619 cff2=0.5_r8*dtfast(ng)*8.0_r8/12.0_r8
2620 cff3=0.5_r8*dtfast(ng)*1.0_r8/12.0_r8
2621#ifdef WET_DRY
2622 cff4=1.0_r8/cff1
2623#endif
2624 DO j=jstr,jend
2625 DO i=istru,iend
2626 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2627 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2628 ubar(i,j,knew)=(ubar(i,j,kstp)* &
2629 & (dstp(i,j)+dstp(i-1,j))+ &
2630 & cff*(cff1*rhs_ubar(i,j)+ &
2631 & cff2*rubar(i,j,kstp)- &
2632 & cff3*rubar(i,j,ptsk)))*fac
2633#ifdef MASKING
2634 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2635#endif
2636#ifdef WET_DRY
2637 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
2638 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2639 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2640 ubar(i,j,knew)=ubar(i,j,knew)*cff7
2641 rhs_ubar(i,j)=rhs_ubar(i,j)*cff7
2642#endif
2643#if defined NESTING && !defined SOLVE3D
2644 du_flux(i,j)=ubar(i,j,knew)* &
2645 & 0.5_r8*(dnew(i,j)+dnew(i-1,j))*on_u(i,j)
2646#endif
2647 END DO
2648 END DO
2649 DO j=jstrv,jend
2650 DO i=istr,iend
2651 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2652 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2653 vbar(i,j,knew)=(vbar(i,j,kstp)* &
2654 & (dstp(i,j)+dstp(i,j-1))+ &
2655 & cff*(cff1*rhs_vbar(i,j)+ &
2656 & cff2*rvbar(i,j,kstp)- &
2657 & cff3*rvbar(i,j,ptsk)))*fac
2658#ifdef MASKING
2659 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2660#endif
2661#ifdef WET_DRY
2662 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
2663 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
2664 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2665 vbar(i,j,knew)=vbar(i,j,knew)*cff7
2666 rhs_vbar(i,j)=rhs_vbar(i,j)*cff7
2667#endif
2668#if defined NESTING && !defined SOLVE3D
2669 dv_flux(i,j)=vbar(i,j,knew)* &
2670 & 0.5_r8*(dnew(i,j)+dnew(i,j-1))*om_v(i,j)
2671#endif
2672 END DO
2673 END DO
2674 END IF
2675
2676#ifdef DIAGNOSTICS_UV
2677!
2678!-----------------------------------------------------------------------
2679! Time step 2D momentum diagnostic terms.
2680!-----------------------------------------------------------------------
2681
2682# ifdef MASKING
2683!
2684! Apply land/sea mask.
2685!
2686 DO idiag=1,ndm2d-1
2687 DO j=jstr,jend
2688 DO i=istru,iend
2689 diau2rhs(i,j,idiag)=diau2rhs(i,j,idiag)*umask(i,j)
2690 END DO
2691 END DO
2692 DO j=jstrv,jend
2693 DO i=istr,iend
2694 diav2rhs(i,j,idiag)=diav2rhs(i,j,idiag)*vmask(i,j)
2695 END DO
2696 END DO
2697 END DO
2698# endif
2699# ifdef SOLVE3D
2700!
2701! The arrays "DiaU2rhs" and "DiaV2rhs" contain the contributions of
2702! each of the 2D right-hand-side terms for the momentum equations.
2703!
2704! These values are integrated, time-stepped and converted to mass flux
2705! units (m3 s-1) for coupling with the 3D diagnostic terms.
2706!
2707 fac=weight(1,iif(ng),ng)
2708 IF (first_2d_step.and.corrector_2d_step) THEN
2709 cff1=0.5_r8*dtfast(ng)
2710 DO idiag=1,ndm2d-1
2711 DO j=jstr,jend
2712 DO i=istru,iend
2713 diau2int(i,j,idiag)=cff1*diau2rhs(i,j,idiag)
2714 diau2wrk(i,j,idiag)=diau2int(i,j,idiag)* &
2715 & (pm(i-1,j)+pm(i,j))*fac
2716 END DO
2717 END DO
2718 DO j=jstrv,jend
2719 DO i=istr,iend
2720 diav2int(i,j,idiag)=cff1*diav2rhs(i,j,idiag)
2721 diav2wrk(i,j,idiag)=diav2int(i,j,idiag)* &
2722 & (pn(i,j)+pn(i,j-1))*fac
2723 END DO
2724 END DO
2725 END DO
2726 ELSE IF (corrector_2d_step) THEN
2727 cff1=0.5_r8*dtfast(ng)*5.0_r8/12.0_r8
2728 cff2=0.5_r8*dtfast(ng)*8.0_r8/12.0_r8
2729 cff3=0.5_r8*dtfast(ng)*1.0_r8/12.0_r8
2730 DO idiag=1,ndm2d-1
2731 DO j=jstr,jend
2732 DO i=istru,iend
2733 diau2int(i,j,idiag)=diau2int(i,j,idiag)+ &
2734 & (cff1*diau2rhs(i,j,idiag)+ &
2735 & cff2*diarubar(i,j,kstp,idiag)- &
2736 & cff3*diarubar(i,j,ptsk,idiag))
2737 diau2wrk(i,j,idiag)=diau2wrk(i,j,idiag)+ &
2738 & diau2int(i,j,idiag)* &
2739 & (pm(i-1,j)+pm(i,j))*fac
2740 END DO
2741 END DO
2742 DO j=jstrv,jend
2743 DO i=istr,iend
2744 diav2int(i,j,idiag)=diav2int(i,j,idiag)+ &
2745 & (cff1*diav2rhs(i,j,idiag)+ &
2746 & cff2*diarvbar(i,j,kstp,idiag)- &
2747 & cff3*diarvbar(i,j,ptsk,idiag))
2748 diav2wrk(i,j,idiag)=diav2wrk(i,j,idiag)+ &
2749 & diav2int(i,j,idiag)* &
2750 & (pn(i,j)+pn(i,j-1))*fac
2751 END DO
2752 END DO
2753 END DO
2754 END IF
2755# else
2756!
2757! Time-step the diagnostic terms.
2758!
2759 IF (first_2d_step.and.corrector_2d_step) THEN
2760 cff1=0.5_r8*dtfast(ng)
2761 DO idiag=1,ndm2d-1
2762 DO j=jstr,jend
2763 DO i=istru,iend
2764 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2765 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2766 diau2wrk(i,j,idiag)=cff*cff1*diau2rhs(i,j,idiag)*fac
2767 END DO
2768 END DO
2769 DO j=jstrv,jend
2770 DO i=istr,iend
2771 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2772 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2773 diav2wrk(i,j,idiag)=cff*cff1*diav2rhs(i,j,idiag)*fac
2774 END DO
2775 END DO
2776 END DO
2777 DO j=jstr,jend
2778 DO i=istru,iend
2779 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2780 diau2wrk(i,j,m2rate)=ubar(i,j,knew)-ubar(i,j,kstp)* &
2781 & (dstp(i,j)+dstp(i-1,j))*fac
2782 END DO
2783 END DO
2784 DO j=jstrv,jend
2785 DO i=istr,iend
2786 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2787 diav2wrk(i,j,m2rate)=vbar(i,j,knew)-vbar(i,j,kstp)* &
2788 & (dstp(i,j)+dstp(i,j-1))*fac
2789 END DO
2790 END DO
2791 ELSE IF (corrector_2d_step) THEN
2792 cff1=0.5_r8*dtfast(ng)*5.0_r8/12.0_r8
2793 cff2=0.5_r8*dtfast(ng)*8.0_r8/12.0_r8
2794 cff3=0.5_r8*dtfast(ng)*1.0_r8/12.0_r8
2795 DO idiag=1,ndm2d-1
2796 DO j=jstr,jend
2797 DO i=istru,iend
2798 cff=(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2799 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2800 diau2wrk(i,j,idiag)=cff*(cff1*diau2rhs(i,j,idiag)+ &
2801 & cff2*diarubar(i,j,kstp,idiag)- &
2802 & cff3*diarubar(i,j,ptsk,idiag))* &
2803 & fac
2804 END DO
2805 END DO
2806 DO j=jstrv,jend
2807 DO i=istr,iend
2808 cff=(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2809 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2810 diav2wrk(i,j,idiag)=cff*(cff1*diav2rhs(i,j,idiag)+ &
2811 & cff2*diarvbar(i,j,kstp,idiag)- &
2812 & cff3*diarvbar(i,j,ptsk,idiag))* &
2813 & fac
2814 END DO
2815 END DO
2816 END DO
2817 DO j=jstr,jend
2818 DO i=istru,iend
2819 fac=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2820 diau2wrk(i,j,m2rate)=ubar(i,j,knew)- &
2821 & ubar(i,j,kstp)* &
2822 & (dstp(i,j)+dstp(i-1,j))*fac
2823 END DO
2824 END DO
2825 DO j=jstrv,jend
2826 DO i=istr,iend
2827 fac=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2828 diav2wrk(i,j,m2rate)=vbar(i,j,knew)- &
2829 & vbar(i,j,kstp)* &
2830 & (dstp(i,j)+dstp(i,j-1))*fac
2831 END DO
2832 END DO
2833 END IF
2834# endif
2835#endif
2836!
2837! If predictor step, load right-side-term into shared arrays for
2838! future use during the subsequent corrector step.
2839!
2840 IF (predictor_2d_step(ng)) THEN
2841 DO j=jstr,jend
2842 DO i=istru,iend
2843 rubar(i,j,krhs)=rhs_ubar(i,j)
2844 END DO
2845 END DO
2846 DO j=jstrv,jend
2847 DO i=istr,iend
2848 rvbar(i,j,krhs)=rhs_vbar(i,j)
2849 END DO
2850 END DO
2851#ifdef DIAGNOSTICS_UV
2852 DO idiag=1,ndm2d-1
2853 DO j=jstr,jend
2854 DO i=istru,iend
2855 diarubar(i,j,krhs,idiag)=diau2rhs(i,j,idiag)
2856 END DO
2857 END DO
2858 DO j=jstrv,jend
2859 DO i=istr,iend
2860 diarvbar(i,j,krhs,idiag)=diav2rhs(i,j,idiag)
2861 END DO
2862 END DO
2863 END DO
2864#endif
2865 END IF
2866!
2867!-----------------------------------------------------------------------
2868! Apply lateral boundary conditions.
2869!-----------------------------------------------------------------------
2870!
2871 CALL u2dbc_tile (ng, tile, &
2872 & lbi, ubi, lbj, ubj, &
2873 & imins, imaxs, jmins, jmaxs, &
2874 & krhs, kstp, knew, &
2875 & ubar, vbar, zeta)
2876 CALL v2dbc_tile (ng, tile, &
2877 & lbi, ubi, lbj, ubj, &
2878 & imins, imaxs, jmins, jmaxs, &
2879 & krhs, kstp, knew, &
2880 & ubar, vbar, zeta)
2881!
2882! Compute integral mass flux across open boundaries and adjust
2883! for volume conservation.
2884!
2885 IF (any(volcons(:,ng))) THEN
2886 CALL obc_flux_tile (ng, tile, &
2887 & lbi, ubi, lbj, ubj, &
2888 & imins, imaxs, jmins, jmaxs, &
2889 & knew, &
2890#ifdef MASKING
2891 & umask, vmask, &
2892#endif
2893 & h, om_v, on_u, &
2894 & ubar, vbar, zeta)
2895 END IF
2896
2897#if defined NESTING && !defined SOLVE3D
2898!
2899!-----------------------------------------------------------------------
2900! Set barotropic fluxes along physical boundaries.
2901!-----------------------------------------------------------------------
2902!
2903 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
2904 IF (domain(ng)%Western_Edge(tile)) THEN
2905 DO j=jstr-1,jendr
2906 dnew(istr-1,j)=h(istr-1,j)+zeta_new(istr-1,j)
2907 END DO
2908 END IF
2909 END IF
2910 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
2911 IF (domain(ng)%Eastern_Edge(tile)) THEN
2912 DO j=jstr-1,jendr
2913 dnew(iend+1,j)=h(iend+1,j)+zeta_new(iend+1,j)
2914 END DO
2915 END IF
2916 END IF
2917 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
2918 IF (domain(ng)%Southern_Edge(tile)) THEN
2919 DO i=istr-1,iendr
2920 dnew(i,jstr-1)=h(i,jstr-1)+zeta_new(i,jstr-1)
2921 END DO
2922 END IF
2923 END IF
2924 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
2925 IF (domain(ng)%Northern_Edge(tile)) THEN
2926 DO i=istr-1,iendr
2927 dnew(i,jend+1)=h(i,jend+1)+zeta_new(i,jend+1)
2928 END DO
2929 END IF
2930 END IF
2931!
2932 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
2933 IF (domain(ng)%Western_Edge(tile)) THEN
2934 DO j=jstrr,jendr
2935 du_flux(istru-1,j)=ubar(istru-1,j,knew)*on_u(istru-1,j)* &
2936 & 0.5_r8*(dnew(istru-1,j)+dnew(istru-2,j))
2937 END DO
2938 DO j=jstrv,jend
2939 dv_flux(istr-1,j)=vbar(istr-1,j,knew)*om_v(istr-1,j)* &
2940 & 0.5_r8*(dnew(istr-1,j)+dnew(istr-1,j-1))
2941 END DO
2942 END IF
2943 END IF
2944 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
2945 IF (domain(ng)%Eastern_Edge(tile)) THEN
2946 DO j=jstrr,jendr
2947 du_flux(iend+1,j)=ubar(iend+1,j,knew)*on_u(iend+1,j)* &
2948 & 0.5_r8*(dnew(iend+1,j)+dnew(iend,j))
2949 END DO
2950 DO j=jstrv,jend
2951 dv_flux(iend+1,j)=vbar(iend+1,j,knew)*om_v(iend+1,j)* &
2952 & 0.5_r8*(dnew(iend+1,j)+dnew(iend+1,j-1))
2953 END DO
2954 END IF
2955 END IF
2956 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
2957 IF (domain(ng)%Southern_Edge(tile)) THEN
2958 DO i=istru,iend
2959 du_flux(i,jstr-1)=ubar(i,jstr-1,knew)*on_u(i,jstr-1)* &
2960 & 0.5_r8*(dnew(i,jstr-1)+dnew(i-1,jstr-1))
2961 END DO
2962 DO i=istrr,iendr
2963 dv_flux(i,jstrv-1)=vbar(i,jstrv-1,knew)*om_v(i,jstrv-1)* &
2964 & 0.5_r8*(dnew(i,jstrv-1)+dnew(i,jstrv-2))
2965 END DO
2966 END IF
2967 END IF
2968 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
2969 IF (domain(ng)%Northern_Edge(tile)) THEN
2970 DO i=istru,iend
2971 du_flux(i,jend+1)=ubar(i,jend+1,knew)*on_u(i,jend+1)* &
2972 & 0.5_r8*(dnew(i,jend+1)+dnew(i-1,jend+1))
2973 END DO
2974 DO i=istrr,iendr
2975 dv_flux(i,jend+1)=vbar(i,jend+1,knew)*om_v(i,jend+1)* &
2976 & 0.5_r8*(dnew(i,jend+1)+dnew(i,jend))
2977 END DO
2978 END IF
2979 END IF
2980#endif
2981!
2982!-----------------------------------------------------------------------
2983! Apply momentum transport point sources (like river runoff), if any.
2984!
2985! Dsrc(is) = 0, flow across grid cell u-face (positive or negative)
2986! Dsrc(is) = 1, flow across grid cell v-face (positive or negative)
2987!-----------------------------------------------------------------------
2988!
2989 IF (luvsrc(ng)) THEN
2990 DO is=1,nsrc(ng)
2991 i=sources(ng)%Isrc(is)
2992 j=sources(ng)%Jsrc(is)
2993 IF (((istrr.le.i).and.(i.le.iendr)).and. &
2994 & ((jstrr.le.j).and.(j.le.jendr))) THEN
2995 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
2996 cff=1.0_r8/(on_u(i,j)* &
2997 & 0.5_r8*(zeta(i-1,j,knew)+h(i-1,j)+ &
2998 & zeta(i ,j,knew)+h(i ,j)))
2999 ubar(i,j,knew)=sources(ng)%Qbar(is)*cff
3000#if defined NESTING && !defined SOLVE3D
3001 du_flux(i,j)=sources(ng)%Qbar(is)
3002#endif
3003 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
3004 cff=1.0_r8/(om_v(i,j)* &
3005 & 0.5_r8*(zeta(i,j-1,knew)+h(i,j-1)+ &
3006 & zeta(i,j ,knew)+h(i,j )))
3007 vbar(i,j,knew)=sources(ng)%Qbar(is)*cff
3008#if defined NESTING && !defined SOLVE3D
3009 dv_flux(i,j)=sources(ng)%Qbar(is)
3010#endif
3011 END IF
3012 END IF
3013 END DO
3014 END IF
3015!
3016!-----------------------------------------------------------------------
3017! Exchange boundary information.
3018!-----------------------------------------------------------------------
3019!
3020 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
3021 CALL exchange_u2d_tile (ng, tile, &
3022 & lbi, ubi, lbj, ubj, &
3023 & ubar(:,:,knew))
3024 CALL exchange_v2d_tile (ng, tile, &
3025 & lbi, ubi, lbj, ubj, &
3026 & vbar(:,:,knew))
3027
3028#if defined NESTING && !defined SOLVE3D
3029 CALL exchange_u2d_tile (ng, tile, &
3030 & lbi, ubi, lbj, ubj, &
3031 & du_flux)
3032 CALL exchange_v2d_tile (ng, tile, &
3033 & lbi, ubi, lbj, ubj, &
3034 & dv_flux)
3035#endif
3036 END IF
3037
3038#ifdef DISTRIBUTE
3039!
3040# if defined NESTING && !defined SOLVE3D
3041 CALL mp_exchange2d (ng, tile, inlm, 4, &
3042# else
3043 CALL mp_exchange2d (ng, tile, inlm, 2, &
3044# endif
3045 & lbi, ubi, lbj, ubj, &
3046 & nghostpoints, &
3047 & ewperiodic(ng), nsperiodic(ng), &
3048# if defined NESTING && !defined SOLVE3D
3049 & du_flux, dv_flux, &
3050# endif
3051 & ubar(:,:,knew), &
3052 & vbar(:,:,knew))
3053#endif
3054!
3055 RETURN
3056 END SUBROUTINE step2d_tile
3057!
3058 END MODULE step2d_mod
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_r2d_tile
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:253
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
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_ncparam
Definition mod_ncparam.F:2
mod_ncparam::isvbar
integer isvbar
Definition mod_ncparam.F:505
mod_ncparam::isubar
integer isubar
Definition mod_ncparam.F:504
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::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::lbc
type(t_lbc), dimension(:,:,:), allocatable lbc
Definition mod_param.F:375
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_scalars
Definition mod_scalars.F:2
mod_scalars::luvsrc
logical, dimension(:), allocatable luvsrc
Definition mod_scalars.F:542
mod_scalars::m2fcor
integer m2fcor
Definition mod_scalars.F:150
mod_scalars::lnudgem2clm
logical, dimension(:), allocatable lnudgem2clm
Definition mod_scalars.F:533
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::dcrit
real(r8), dimension(:), allocatable dcrit
Definition mod_scalars.F:782
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::predictor_2d_step
logical, dimension(:), allocatable predictor_2d_step
Definition mod_scalars.F:245
mod_scalars::weight
real(dp), dimension(:,:,:), allocatable weight
Definition mod_scalars.F:809
mod_scalars::volcons
logical, dimension(:,:), allocatable volcons
Definition mod_scalars.F:515
mod_scalars::nfast
integer, dimension(:), allocatable nfast
Definition mod_scalars.F:259
mod_scalars::ndtfast
integer, dimension(:), allocatable ndtfast
Definition mod_scalars.F:258
mod_scalars::lwsrc
logical, dimension(:), allocatable lwsrc
Definition mod_scalars.F:543
mod_scalars::gamma2
real(r8), dimension(:), allocatable gamma2
Definition mod_scalars.F:743
mod_scalars::m2hvis
integer m2hvis
Definition mod_scalars.F:168
mod_scalars::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::m2rate
integer m2rate
Definition mod_scalars.F:173
mod_scalars::m2yvis
integer m2yvis
Definition mod_scalars.F:170
mod_scalars::dtfast
real(dp), dimension(:), allocatable dtfast
Definition mod_scalars.F:270
mod_scalars::ntfirst
integer, dimension(:), allocatable ntfirst
Definition mod_scalars.F:377
mod_scalars::m2sstr
integer m2sstr
Definition mod_scalars.F:171
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::m2hadv
integer m2hadv
Definition mod_scalars.F:151
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::m2xvis
integer m2xvis
Definition mod_scalars.F:169
mod_scalars::m2bstr
integer m2bstr
Definition mod_scalars.F:172
mod_scalars::iif
integer, dimension(:), allocatable iif
Definition mod_scalars.F:250
mod_sedbed
Definition sedbed_mod.h:1
mod_sedbed::sedbed
type(t_sedbed), dimension(:), allocatable sedbed
Definition sedbed_mod.h:157
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::kstp
integer, dimension(:), allocatable kstp
Definition mod_stepping.F:66
mod_stepping::knew
integer, dimension(:), allocatable knew
Definition mod_stepping.F:64
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::krhs
integer, dimension(:), allocatable krhs
Definition mod_stepping.F:65
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
obc_volcons_mod
Definition obc_volcons.F:3
obc_volcons_mod::set_duv_bc_tile
subroutine, public set_duv_bc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, om_v, on_u, ubar, vbar, drhs, duon, dvom)
Definition obc_volcons.F:245
obc_volcons_mod::obc_flux_tile
subroutine, public obc_flux_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, h, om_v, on_u, ubar, vbar, zeta)
Definition obc_volcons.F:69
step2d_mod
Definition step2d_FB.h:3
step2d_mod::step2d_tile
subroutine step2d_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, krhs, kstp, knew, nstp, nnew, pmask, rmask, umask, vmask, pmask_wet, pmask_full, rmask_wet, rmask_full, umask_wet, umask_full, vmask_wet, vmask_full, rmask_wet_avg, fomn, h, om_u, om_v, on_u, on_v, pm, pn, dndx, dmde, rdrag, rdrag2, pmon_r, pnom_r, pmon_p, pnom_p, om_r, on_r, om_p, on_p, visc2_p, visc2_r, visc4_p, visc4_r, bed_thick, eq_tide, sustr, svstr,
Definition step2d_FB.h:219
step2d_mod::step2d
subroutine, public step2d(ng, tile)
Definition step2d_FB.h:75
u2dbc_mod
Definition u2dbc_im.F:2
u2dbc_mod::u2dbc_tile
subroutine, public u2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta)
Definition u2dbc_im.F:56
v2dbc_mod
Definition v2dbc_im.F:3
v2dbc_mod::v2dbc_tile
subroutine, public v2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta)
Definition v2dbc_im.F:57
wetdry_mod
Definition wetdry.F:2
wetdry_mod::wetdry_tile
subroutine wetdry_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, pmask, rmask, umask, vmask, h, zeta, du_avg1, dv_avg1, rmask_wet_avg, pmask_wet, pmask_full, rmask_wet, rmask_full, umask_wet, umask_full, vmask_wet, vmask_full)
Definition wetdry.F:108
zetabc_mod
Definition zetabc.F:2
zetabc_mod::zetabc_tile
subroutine, public zetabc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, zeta)
Definition zetabc.F:65
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