ROMS
Loading...
Searching...
No Matches
step2d_FB.h
Go to the documentation of this file.
1#define DEBUG
2
3 MODULE step2d_mod
4!
5!git $Id$
6!=======================================================================
7! !
8! Solves nonlinear shallow-water primitive equations (barotropic mode)!
9! using the Generalized Forward-Backward 3rd-order Adams-Bashforth / !
10! 4th-order Adams-Moulton (FB AB3-AM4) time stepping algorithm !
11! (Shchepetkin and McWilliams, 2005); see section 2.3 starting with !
12! equation 2.49. In 3D applications, it perform fast-time averaging !
13! to interact with 3D momentum equations (baroclinic mode). !
14! !
15! Reference: !
16! !
17! Shchepetkin, A.F. and J.C. McWilliams, 2005: The regional oceanic !
18! modeling system (ROMS): a split-explicit, free-surface, !
19! topography-following-coordinate oceanic model, Ocean Modelling, !
20! 9, 347-404, doi:10.1016/j.ocemod.2004.08.002. !
21! !
22! Shchepetkin, A.F., and J.C. McWilliams, 2009: Computational kernel !
23! algorithms for fine-scale, multiprocess, longtime oceanic !
24! simulations, pp 121-183. In 'Handbook of Numerical Analysis: !
25! Computational Methods for the Atmosphere and Oceans', R.M. Teman !
26! and J.J. Tribbia, eds, Elsevier Science. !
27! !
28! Adapted from A.F. Shchepetkin routine "step2d_FB.F" (07-26-2022) !
29! !
30!=======================================================================
31!
32 USE mod_param
33 USE mod_parallel
34#ifdef SOLVE3D
35 USE mod_coupling
36#endif
37#ifdef DIAGNOSTICS_UV
38 USE mod_diags
39#endif
40 USE mod_forces
41 USE mod_grid
42 USE mod_mixing
43 USE mod_ocean
44 USE mod_scalars
45#if defined SEDIMENT && defined SED_MORPH && defined SOLVE3D
46 USE mod_sedbed
47#endif
48 USE mod_sources
49 USE mod_stepping
50!
51 USE exchange_2d_mod
52#ifdef DISTRIBUTE
53 USE mp_exchange_mod, ONLY : mp_exchange2d
54#endif
55 USE obc_volcons_mod, ONLY : obc_flux_tile, set_duv_bc_tile
56#ifdef SOLVE3D
57 USE set_depth_mod, ONLY : set_depth
58#endif
59 USE u2dbc_mod, ONLY : u2dbc_tile
60 USE v2dbc_mod, ONLY : v2dbc_tile
61#ifdef WET_DRY
62 USE wetdry_mod, ONLY : wetdry_tile
63#endif
64 USE zetabc_mod, ONLY : zetabc_local
65!
66 implicit none
67!
68 PRIVATE
69 PUBLIC :: step2d
70!
71 CONTAINS
72!
73!***********************************************************************
74 SUBROUTINE step2d (ng, tile)
75!***********************************************************************
76!
77! Imported variable declarations.
78!
79 integer, intent(in) :: ng, tile
80!
81! Local variable declarations.
82!
83 character (len=*), parameter :: myfile = &
84 & __FILE__
85!
86#include "tile.h"
87!
88#ifdef PROFILE
89 CALL wclock_on (ng, inlm, 9, __line__, myfile)
90#endif
91 CALL step2d_tile (ng, tile, &
92 & lbi, ubi, lbj, ubj, n(ng), &
93 & imins, imaxs, jmins, jmaxs, &
94 & krhs(ng), kstp(ng), knew(ng), &
95#ifdef SOLVE3D
96 & nstp(ng), nnew(ng), &
97#endif
98#ifdef MASKING
99 & grid(ng) % pmask, grid(ng) % rmask, &
100 & grid(ng) % umask, grid(ng) % vmask, &
101#endif
102#ifdef WET_DRY
103 & grid(ng) % pmask_wet, grid(ng) % pmask_full, &
104 & grid(ng) % rmask_wet, grid(ng) % rmask_full, &
105 & grid(ng) % umask_wet, grid(ng) % umask_full, &
106 & grid(ng) % vmask_wet, grid(ng) % vmask_full, &
107# ifdef SOLVE3D
108 & grid(ng) % rmask_wet_avg, &
109# endif
110#endif
111#if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
112 & grid(ng) % fomn, &
113#endif
114 & grid(ng) % h, &
115 & grid(ng) % om_u, grid(ng) % om_v, &
116 & grid(ng) % on_u, grid(ng) % on_v, &
117 & grid(ng) % pm, grid(ng) % pn, &
118#if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
119 & grid(ng) % dndx, grid(ng) % dmde, &
120#endif
121 & grid(ng) % rdrag, &
122#if defined UV_QDRAG && !defined SOLVE3D
123 & grid(ng) % rdrag2, &
124#endif
125#if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
126 & grid(ng) % pmon_r, grid(ng) % pnom_r, &
127 & grid(ng) % pmon_p, grid(ng) % pnom_p, &
128 & grid(ng) % om_r, grid(ng) % on_r, &
129 & grid(ng) % om_p, grid(ng) % on_p, &
130# ifdef UV_VIS2
131 & mixing(ng) % visc2_p, mixing(ng) % visc2_r, &
132# endif
133# ifdef UV_VIS4
134 & mixing(ng) % visc4_p, mixing(ng) % visc4_r, &
135# endif
136#endif
137#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
138 & ocean(ng) % eq_tide, &
139#endif
140#ifndef SOLVE3D
141 & forces(ng) % sustr, forces(ng) % svstr, &
142# ifdef ATM_PRESS
143 & forces(ng) % Pair, &
144# endif
145#else
146# ifdef VAR_RHO_2D
147 & coupling(ng) % rhoA, coupling(ng) % rhoS, &
148# endif
149 & coupling(ng) % DU_avg1, coupling(ng) % DU_avg2, &
150 & coupling(ng) % DV_avg1, coupling(ng) % DV_avg2, &
151 & coupling(ng) % Zt_avg1, &
152 & coupling(ng) % rufrc, &
153 & coupling(ng) % rvfrc, &
154 & coupling(ng) % rufrc_bak, &
155 & coupling(ng) % rvfrc_bak, &
156#endif
157#if defined NESTING && !defined SOLVE3D
158 & ocean(ng) % DU_flux, ocean(ng) % DV_flux, &
159#endif
160 & ocean(ng) % ubar, ocean(ng) % vbar, &
161 & ocean(ng) % zeta)
162#ifdef PROFILE
163 CALL wclock_off (ng, inlm, 9, __line__, myfile)
164#endif
165!
166 RETURN
167 END SUBROUTINE step2d
168!
169!***********************************************************************
170 SUBROUTINE step2d_tile (ng, tile, &
171 & LBi, UBi, LBj, UBj, UBk, &
172 & IminS, ImaxS, JminS, JmaxS, &
173 & krhs, kstp, knew, &
174#ifdef SOLVE3D
175 & nstp, nnew, &
176#endif
177#ifdef MASKING
178 & pmask, rmask, umask, vmask, &
179#endif
180#ifdef WET_DRY
181 & pmask_wet, pmask_full, &
182 & rmask_wet, rmask_full, &
183 & umask_wet, umask_full, &
184 & vmask_wet, vmask_full, &
185# ifdef SOLVE3D
186 & rmask_wet_avg, &
187# endif
188#endif
189#if (defined UV_COR && !defined SOLVE3D) || defined step2d_coriolis
190 & fomn, &
191#endif
192 & h, &
193 & om_u, om_v, on_u, on_v, pm, pn, &
194#if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
195 & dndx, dmde, &
196#endif
197 & rdrag, &
198#if defined UV_QDRAG && !defined SOLVE3D
199 & rdrag2, &
200#endif
201#if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
202 & pmon_r, pnom_r, pmon_p, pnom_p, &
203 & om_r, on_r, om_p, on_p, &
204# ifdef UV_VIS2
205 & visc2_p, visc2_r, &
206# endif
207# ifdef UV_VIS4
208 & visc4_p, visc4_r, &
209# endif
210#endif
211#if defined SEDIMENT && defined SED_MORPH
212 & bed_thick, &
213#endif
214#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
215 & eq_tide, &
216#endif
217#ifndef SOLVE3D
218 & sustr, svstr,
219# ifdef ATM_PRESS
220 & pair, &
221# endif
222#else
223# ifdef VAR_RHO_2D
224 & rhoa, rhos, &
225# endif
226 & du_avg1, du_avg2, &
227 & dv_avg1, dv_avg2, &
228 & zt_avg1, &
229 & rufrc, rvfrc, &
230 & rufrc_bak, rvfrc_bak, &
231#endif
232#ifdef DIAGNOSTICS_UV
233 & diau2wrk, diav2wrk, &
234 & diarubar, diarvbar, &
235# ifdef SOLVE3D
236 & diau2int, diav2int, &
237 & diarufrc, diarvfrc, &
238# endif
239#endif
240#if defined NESTING && !defined SOLVE3D
241 & du_flux, dv_flux, &
242#endif
243 & ubar, vbar, zeta)
244!***********************************************************************
245!
246! Imported variable declarations.
247!
248 integer, intent(in ) :: ng, tile
249 integer, intent(in ) :: LBi, UBi, LBj, UBj, UBk
250 integer, intent(in ) :: IminS, ImaxS, JminS, JmaxS
251 integer, intent(in ) :: krhs, kstp, knew
252#ifdef SOLVE3D
253 integer, intent(in ) :: nstp, nnew
254#endif
255!
256#ifdef ASSUMED_SHAPE
257# ifdef MASKING
258 real(r8), intent(in ) :: pmask(LBi:,LBj:)
259 real(r8), intent(in ) :: rmask(LBi:,LBj:)
260 real(r8), intent(in ) :: umask(LBi:,LBj:)
261 real(r8), intent(in ) :: vmask(LBi:,LBj:)
262# endif
263# if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
264 real(r8), intent(in ) :: fomn(LBi:,LBj:)
265# endif
266# if defined SEDIMENT && defined SED_MORPH
267 real(r8), intent(inout) :: h(LBi:,LBj:)
268# else
269 real(r8), intent(in ) :: h(LBi:,LBj:)
270# endif
271 real(r8), intent(in ) :: om_u(LBi:,LBj:)
272 real(r8), intent(in ) :: om_v(LBi:,LBj:)
273 real(r8), intent(in ) :: on_u(LBi:,LBj:)
274 real(r8), intent(in ) :: on_v(LBi:,LBj:)
275 real(r8), intent(in ) :: pm(LBi:,LBj:)
276 real(r8), intent(in ) :: pn(LBi:,LBj:)
277# if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
278 real(r8), intent(in ) :: dndx(LBi:,LBj:)
279 real(r8), intent(in ) :: dmde(LBi:,LBj:)
280# endif
281 real(r8), intent(in ) :: rdrag(LBi:,LBj:)
282# if defined UV_QDRAG && !defined SOLVE3D
283 real(r8), intent(in ) :: rdrag2(LBi:,LBj:)
284# endif
285# if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
286 real(r8), intent(in ) :: pmon_r(LBi:,LBj:)
287 real(r8), intent(in ) :: pnom_r(LBi:,LBj:)
288 real(r8), intent(in ) :: pmon_p(LBi:,LBj:)
289 real(r8), intent(in ) :: pnom_p(LBi:,LBj:)
290 real(r8), intent(in ) :: om_r(LBi:,LBj:)
291 real(r8), intent(in ) :: on_r(LBi:,LBj:)
292 real(r8), intent(in ) :: om_p(LBi:,LBj:)
293 real(r8), intent(in ) :: on_p(LBi:,LBj:)
294# ifdef UV_VIS2
295 real(r8), intent(in ) :: visc2_p(LBi:,LBj:)
296 real(r8), intent(in ) :: visc2_r(LBi:,LBj:)
297# endif
298# ifdef UV_VIS4
299 real(r8), intent(in ) :: visc4_p(LBi:,LBj:)
300 real(r8), intent(in ) :: visc4_r(LBi:,LBj:)
301# endif
302# endif
303# if defined SEDIMENT && defined SED_MORPH
304 real(r8), intent(in ) :: bed_thick(LBi:,LBj:,:)
305# endif
306# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
307 real(r8), intent(in ) :: eq_tide(LBi:,LBj:)
308# endif
309# ifndef SOLVE3D
310 real(r8), intent(in ) :: sustr(LBi:,LBj:)
311 real(r8), intent(in ) :: svstr(LBi:,LBj:)
312# ifdef ATM_PRESS
313 real(r8), intent(in ) :: Pair(LBi:,LBj:)
314# endif
315# else
316# ifdef VAR_RHO_2D
317 real(r8), intent(in ) :: rhoA(LBi:,LBj:)
318 real(r8), intent(in ) :: rhoS(LBi:,LBj:)
319# endif
320 real(r8), intent(inout) :: DU_avg1(LBi:,LBj:)
321 real(r8), intent(inout) :: DU_avg2(LBi:,LBj:)
322 real(r8), intent(inout) :: DV_avg1(LBi:,LBj:)
323 real(r8), intent(inout) :: DV_avg2(LBi:,LBj:)
324 real(r8), intent(inout) :: Zt_avg1(LBi:,LBj:)
325 real(r8), intent(inout) :: rufrc(LBi:,LBj:)
326 real(r8), intent(inout) :: rvfrc(LBi:,LBj:)
327 real(r8), intent(inout) :: rufrc_bak(LBi:,LBj:,:)
328 real(r8), intent(inout) :: rvfrc_bak(LBi:,LBj:,:)
329# endif
330# ifdef WET_DRY
331 real(r8), intent(inout) :: pmask_full(LBi:,LBj:)
332 real(r8), intent(inout) :: rmask_full(LBi:,LBj:)
333 real(r8), intent(inout) :: umask_full(LBi:,LBj:)
334 real(r8), intent(inout) :: vmask_full(LBi:,LBj:)
335
336 real(r8), intent(inout) :: pmask_wet(LBi:,LBj:)
337 real(r8), intent(inout) :: rmask_wet(LBi:,LBj:)
338 real(r8), intent(inout) :: umask_wet(LBi:,LBj:)
339 real(r8), intent(inout) :: vmask_wet(LBi:,LBj:)
340# ifdef SOLVE3D
341 real(r8), intent(inout) :: rmask_wet_avg(LBi:,LBj:)
342# endif
343# endif
344 real(r8), intent(inout) :: ubar(LBi:,LBj:,:)
345 real(r8), intent(inout) :: vbar(LBi:,LBj:,:)
346 real(r8), intent(inout) :: zeta(LBi:,LBj:,:)
347# if defined NESTING && !defined SOLVE3D
348 real(r8), intent(out ) :: DU_flux(LBi:,LBj:)
349 real(r8), intent(out ) :: DV_flux(LBi:,LBj:)
350# endif
351
352#else
353
354# ifdef MASKING
355 real(r8), intent(in ) :: pmask(LBi:UBi,LBj:UBj)
356 real(r8), intent(in ) :: rmask(LBi:UBi,LBj:UBj)
357 real(r8), intent(in ) :: umask(LBi:UBi,LBj:UBj)
358 real(r8), intent(in ) :: vmask(LBi:UBi,LBj:UBj)
359# endif
360# if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
361 real(r8), intent(in ) :: fomn(LBi:UBi,LBj:UBj)
362# endif
363# if defined SEDIMENT && defined SED_MORPH
364 real(r8), intent(inout) :: h(LBi:UBi,LBj:UBj)
365# else
366 real(r8), intent(in ) :: h(LBi:UBi,LBj:UBj)
367# endif
368 real(r8), intent(in ) :: om_u(LBi:UBi,LBj:UBj)
369 real(r8), intent(in ) :: om_v(LBi:UBi,LBj:UBj)
370 real(r8), intent(in ) :: on_u(LBi:UBi,LBj:UBj)
371 real(r8), intent(in ) :: on_v(LBi:UBi,LBj:UBj)
372 real(r8), intent(in ) :: pm(LBi:UBi,LBj:UBj)
373 real(r8), intent(in ) :: pn(LBi:UBi,LBj:UBj)
374# if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
375 real(r8), intent(in ) :: dndx(LBi:UBi,LBj:UBj)
376 real(r8), intent(in ) :: dmde(LBi:UBi,LBj:UBj)
377# endif
378 real(r8), intent(in ) :: rdrag(LBi:UBi,LBj:UBj)
379# if defined UV_QDRAG && !defined SOLVE3D
380 real(r8), intent(in ) :: rdrag2(LBi:UBi,LBj:UBj)
381# endif
382# if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
383 real(r8), intent(in ) :: pmon_r(LBi:UBi,LBj:UBj)
384 real(r8), intent(in ) :: pnom_r(LBi:UBi,LBj:UBj)
385 real(r8), intent(in ) :: pmon_p(LBi:UBi,LBj:UBj)
386 real(r8), intent(in ) :: pnom_p(LBi:UBi,LBj:UBj)
387 real(r8), intent(in ) :: om_r(LBi:UBi,LBj:UBj)
388 real(r8), intent(in ) :: on_r(LBi:UBi,LBj:UBj)
389 real(r8), intent(in ) :: om_p(LBi:UBi,LBj:UBj)
390 real(r8), intent(in ) :: on_p(LBi:UBi,LBj:UBj)
391# ifdef UV_VIS2
392 real(r8), intent(in ) :: visc2_p(LBi:UBi,LBj:UBj)
393 real(r8), intent(in ) :: visc2_r(LBi:UBi,LBj:UBj)
394# endif
395# ifdef UV_VIS4
396 real(r8), intent(in ) :: visc4_p(LBi:UBi,LBj:UBj)
397 real(r8), intent(in ) :: visc4_r(LBi:UBi,LBj:UBj)
398# endif
399# endif
400# if defined SEDIMENT && defined SED_MORPH
401 real(r8), intent(in ) :: bed_thick(LBi:UBi,LBj:UBj,1:3)
402# endif
403# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
404 real(r8), intent(in ) :: eq_tide(LBi:UBi,LBj:UBj)
405# endif
406# ifndef SOLVE3D
407 real(r8), intent(in ) :: sustr(LBi:UBi,LBj:UBj)
408 real(r8), intent(in ) :: svstr(LBi:UBi,LBj:UBj)
409# ifdef ATM_PRESS
410 real(r8), intent(in ) :: Pair(LBi:UBi,LBj:UBj)
411# endif
412# else
413# ifdef VAR_RHO_2D
414 real(r8), intent(in ) :: rhoA(LBi:UBi,LBj:UBj)
415 real(r8), intent(in ) :: rhoS(LBi:UBi,LBj:UBj)
416# endif
417 real(r8), intent(inout) :: DU_avg1(LBi:UBi,LBj:UBj)
418 real(r8), intent(inout) :: DU_avg2(LBi:UBi,LBj:UBj)
419 real(r8), intent(inout) :: DV_avg1(LBi:UBi,LBj:UBj)
420 real(r8), intent(inout) :: DV_avg2(LBi:UBi,LBj:UBj)
421 real(r8), intent(inout) :: Zt_avg1(LBi:UBi,LBj:UBj)
422 real(r8), intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
423 real(r8), intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
424 real(r8), intent(inout) :: rufrc_bak(LBi:UBi,LBj:UBj,2)
425 real(r8), intent(inout) :: rvfrc_bak(LBi:UBi,LBj:UBj,2)
426# endif
427# ifdef WET_DRY
428 real(r8), intent(inout) :: pmask_full(LBi:UBi,LBj:UBj)
429 real(r8), intent(inout) :: rmask_full(LBi:UBi,LBj:UBj)
430 real(r8), intent(inout) :: umask_full(LBi:UBi,LBj:UBj)
431 real(r8), intent(inout) :: vmask_full(LBi:UBi,LBj:UBj)
432
433 real(r8), intent(inout) :: pmask_wet(LBi:UBi,LBj:UBj)
434 real(r8), intent(inout) :: rmask_wet(LBi:UBi,LBj:UBj)
435 real(r8), intent(inout) :: umask_wet(LBi:UBi,LBj:UBj)
436 real(r8), intent(inout) :: vmask_wet(LBi:UBi,LBj:UBj)
437# ifdef SOLVE3D
438 real(r8), intent(inout) :: rmask_wet_avg(LBi:UBi,LBj:UBj)
439# endif
440# endif
441 real(r8), intent(inout) :: ubar(LBi:UBi,LBj:UBj,:)
442 real(r8), intent(inout) :: vbar(LBi:UBi,LBj:UBj,:)
443 real(r8), intent(inout) :: zeta(LBi:UBi,LBj:UBj,:)
444# if defined NESTING && !defined SOLVE3D
445 real(r8), intent(out ) :: DU_flux(LBi:UBi,LBj:UBj)
446 real(r8), intent(out ) :: DV_flux(LBi:UBi,LBj:UBj)
447# endif
448#endif
449!
450! Local variable declarations.
451!
452 integer :: i, is, j
453 integer :: kbak, kold
454#ifdef DIAGNOSTICS_UV
455 integer :: idiag
456#endif
457!
458 real(r8) :: bkw0, bkw1, bkw2, bkw_new
459 real(r8) :: fwd0, fwd1, fwd2
460#ifdef SOLVE3D
461 real(r8) :: cfwd0, cfwd1, cfwd2
462#endif
463 real(r8) :: cff, cff1, cff2, cff3, cff4
464#ifdef WET_DRY
465 real(r8) :: cff5, cff6, cff7
466#endif
467 real(r8) :: fac, fac1, fac2
468#ifdef DEBUG
469 real(r8), parameter :: IniVal = 0.0_r8
470#endif
471!
472#if defined UV_C4ADVECTION && !defined SOLVE3D
473 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dgrad
474#endif
475 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dnew
476 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dnew_rd
477 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs
478#if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
479 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs_p
480#endif
481 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dstp
482 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUon
483 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVom
484#if defined STEP2D_CORIOLIS || !defined SOLVE3D
485 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
486 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
487#endif
488#if !defined SOLVE3D
489 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
490 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
491#endif
492#if defined UV_C4ADVECTION && !defined SOLVE3D
493 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
494#endif
495 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzeta2
496#if defined VAR_RHO_2D && defined SOLVE3D
497 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzetaSA
498#endif
499 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rubar
500 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rvbar
501 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzeta
502 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: urhs
503 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: vrhs
504 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zwrk
505#ifdef WET_DRY
506 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: wetdry
507#endif
508#ifdef DIAGNOSTICS_UV
509 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Uwrk
510 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vwrk
511 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaU2rhs
512 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaV2rhs
513#endif
514!
515 real(r8), allocatable :: zeta_new(:,:)
516
517#include "set_bounds.h"
518
519#ifdef DEBUG
520!
521!-----------------------------------------------------------------------
522! Initialize private arrays for debugging.
523!-----------------------------------------------------------------------
524!
525# if defined UV_C4ADVECTION && !defined SOLVE3D
526 dgrad=inival
527# endif
528 dnew=inival
529 dnew_rd=inival
530 drhs=inival
531# if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
532 drhs_p=inival
533# endif
534 dstp=inival
535 duon=inival
536 dvom=inival
537# if defined STEP2D_CORIOLIS || !defined SOLVE3D
538 ufx=inival
539 vfe=inival
540# endif
541# if !defined SOLVE3D
542 ufe=inival
543 vfx=inival
544# endif
545# if defined UV_C4ADVECTION && !defined SOLVE3D
546 grad=inival
547# endif
548 rzeta2=inival
549# if defined VAR_RHO_2D && defined SOLVE3D
550 rzetasa=inival
551# endif
552 rzeta=inival
553 rubar=inival
554 rvbar=inival
555 urhs=inival
556 vrhs=inival
557 zwrk=inival
558# ifdef WET_DRY
559 wetdry=inival
560# endif
561# ifdef DIAGNOSTICS_UV
562 uwrk=inival
563 vwrk=inival
564 diau2rhs=inival
565 diav2rhs=inival
566# endif
567#endif
568!
569!-----------------------------------------------------------------------
570! Set coefficients for AB3-AM4 forward-backward algorithm.
571!-----------------------------------------------------------------------
572!
573! Because the Forward Euler step is used to update "zeta" during the
574! first barotropic step, the pressure-gradient term in the momentum
575! equation must be computed via the Backward step to keep it
576! numerically stable. However, this interferes with the computation
577! of forcing terms "rufrc" and "rvfrc" because the free surface in
578! pressure gradient computation in 3D is exactly at the time
579! corresponding to baroclinic step "nstp" (rather than ahead by one
580! barotropic step after it updated by a normal forward-backward step).
581! To resolve this conflict, the pressure gradient term is computed in
582! two stages during the first barotropic step. It uses zeta(:,:,kstp)
583! at first to ensure exact consistency with 3D model. Then, after
584! vertical integrals of 3D right-hand-side "rufrc" and "rvfrc" are
585! converted into forcing terms, add correction based on the difference
586! zeta_new(:,:)-zeta(:,:,kstp) to "rubar" and "rvbar" to make them
587! consistent with the Backward step for pressure gradient.
588! For pressure gradient terms, search for the label PGF_FB_CORRECTION
589! below.
590!
591 IF (first_2d_step) THEN ! Meaning of time indices
592 kbak=kstp !------------------------
593 kold=kstp ! m-2 m-1 m m+1
594 fwd0=1.0_r8 ! kold kbak kstp knew
595 fwd1=0.0_r8 ! fwd2 fwd1 fwd0
596 fwd2=0.0_r8 ! bkw2 bkw1 bkw0 bkw_new
597#ifdef SOLVE3D
598 bkw_new=0.0_r8
599 bkw0=1.0_r8
600#else
601 bkw_new=1.0_r8
602 bkw0=0.0_r8
603#endif
604 bkw1=0.0_r8
605 bkw2=0.0_r8
606 ELSE IF (first_2d_step+1) THEN
607 kbak=kstp-1
608 IF (kbak.lt.1) kbak=4
609 kold=kbak
610 fwd0=1.0_r8 ! Logically AB2-AM3 forward-
611 fwd1=0.0_r8 ! backward scheme with maximum
612 fwd2=0.0_r8 ! stability coefficients while
613 bkw_new=1.0833333333333_r8 ! maintaining third-order
614 bkw0=-0.1666666666666_r8 ! accuracy, alpha_max=1.73
615 bkw1= 0.0833333333333_r8
616 bkw2= 0.0_r8
617 ELSE
618 kbak=kstp-1
619 IF (kbak.lt.1) kbak=4
620 kold=kbak-1
621 IF (kold.lt.1) kold=4
622 fwd0=1.781105_r8
623 fwd1=-1.06221_r8
624 fwd2=0.281105_r8
625 bkw_new=0.614_r8
626 bkw0=0.285_r8
627 bkw1=0.0880_r8
628 bkw2=0.013_r8
629 END IF
630
631#ifdef DEBUG
632!
633 IF (master) THEN
634 WRITE (20,10) iic(ng), iif(ng), kold, kbak, kstp, knew
635 10 FORMAT (' iic = ',i5.5,' iif = ',i3.3, &
636 & ' kold = ',i1,' kbak = ',i1,' kstp = ',i1,' knew = ',i1)
637 END IF
638#endif
639!
640!-----------------------------------------------------------------------
641! Preliminary steps.
642!-----------------------------------------------------------------------
643!
644! Compute total depth of water column and vertically integrated fluxes
645! needed for computing horizontal divergence to advance free surface
646! and for nonlinear advection terms for the barotropic momentum
647! equations.
648!
649#if defined DISTRIBUTE && !defined NESTING
650# define IR_RANGE IstrUm2-1,Iendp2
651# define JR_RANGE JstrVm2-1,Jendp2
652# define IU_RANGE IstrUm1-1,Iendp2
653# define JU_RANGE Jstrm1-1,Jendp2
654# define IV_RANGE Istrm1-1,Iendp2
655# define JV_RANGE JstrVm1-1,Jendp2
656#else
657# define IR_RANGE IstrUm2-1,Iendp2
658# define JR_RANGE JstrVm2-1,Jendp2
659# define IU_RANGE IstrUm2,Iendp2
660# define JU_RANGE JstrVm2-1,Jendp2
661# define IV_RANGE IstrUm2-1,Iendp2
662# define JV_RANGE JstrVm2,Jendp2
663#endif
664
665 DO j=jr_range
666 DO i=ir_range
667 drhs(i,j)=h(i,j)+fwd0*zeta(i,j,kstp)+ &
668 & fwd1*zeta(i,j,kbak)+ &
669 & fwd2*zeta(i,j,kold)
670 END DO
671 END DO
672!
673 DO j=ju_range
674 DO i=iu_range
675 cff=0.5_r8*on_u(i,j)
676 cff1=cff*(drhs(i,j)+drhs(i-1,j))
677 urhs(i,j)=fwd0*ubar(i,j,kstp)+ &
678 & fwd1*ubar(i,j,kbak)+ &
679 & fwd2*ubar(i,j,kold)
680 duon(i,j)=urhs(i,j)*cff1
681 END DO
682 END DO
683!
684 DO j=jv_range
685 DO i=iv_range
686 cff=0.5_r8*om_v(i,j)
687 cff1=cff*(drhs(i,j)+drhs(i,j-1))
688 vrhs(i,j)=fwd0*vbar(i,j,kstp)+ &
689 & fwd1*vbar(i,j,kbak)+ &
690 & fwd2*vbar(i,j,kold)
691 dvom(i,j)=vrhs(i,j)*cff1
692 END DO
693 END DO
694
695#undef IR_RANGE
696#undef IU_RANGE
697#undef IV_RANGE
698#undef JR_RANGE
699#undef JU_RANGE
700#undef JV_RANGE
701
702#if defined DISTRIBUTE && \
703 defined uv_adv && defined uv_c4advection && !defined SOLVE3D
704!
705! In distributed-memory, the I- and J-ranges are different and a
706! special exchange is done here to avoid having three ghost points
707! for high-order numerical stencils. Notice that a private array is
708! passed below to the exchange routine. It also applies periodic
709! boundary conditions, if appropriate and no partitions in I- or
710! J-directions.
711!
712 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
713 CALL exchange_u2d_tile (ng, tile, &
714 & imins, imaxs, jmins, jmaxs, &
715 & duon)
716 CALL exchange_v2d_tile (ng, tile, &
717 & imins, imaxs, jmins, jmaxs, &
718 & dvom)
719 END IF
720 CALL mp_exchange2d (ng, tile, inlm, 2, &
721 & imins, imaxs, jmins, jmaxs, &
722 & nghostpoints, &
723 & ewperiodic(ng), nsperiodic(ng), &
724 & duon, dvom)
725#endif
726!
727! Set vertically integrated mass fluxes DUon and DVom along the open
728! boundaries in such a way that the integral volume is conserved.
729! HGA: Need to resolve 'krhs' index here.
730!
731 IF (any(volcons(:,ng))) THEN
732 CALL set_duv_bc_tile (ng, tile, &
733 & lbi, ubi, lbj, ubj, &
734 & imins, imaxs, jmins, jmaxs, &
735 & krhs, &
736#ifdef MASKING
737 & umask, vmask, &
738#endif
739 & om_v, on_u, &
740 & ubar, vbar, &
741 & drhs, duon, dvom)
742 END IF
743!
744!-----------------------------------------------------------------------
745! Advance free-surface.
746!-----------------------------------------------------------------------
747!
748! Notice that the new local free-surface is allocated so it can be
749! passed as an argumment to "zetabc_local". An automatic array cannot
750! be used here because of weird memory problems.
751!
752 allocate ( zeta_new(imins:imaxs,jmins:jmaxs) )
753 zeta_new = 0.0_r8
754!
755! Compute "zeta_new" at new time step and interpolate it half-step
756! backward, "zwrk" for the subsequent computation of the tangent
757! linear barotropic pressure gradient. Here, we use the BASIC STATE
758! values. Thus, the nonlinear correction to the pressure-gradient
759! term from "kstp" to "knew" is not needed for Forward-Euler to
760! Forward-Backward steps (PGF_FB_CORRECTION method).
761!
762 DO j=jstrv-1,jend
763 DO i=istru-1,iend
764 fac=dtfast(ng)*pm(i,j)*pn(i,j)
765 zeta_new(i,j)=zeta(i,j,kstp)+ &
766 & fac*(duon(i,j)-duon(i+1,j)+ &
767 & dvom(i,j)-dvom(i,j+1))
768#ifdef MASKING
769 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
770# ifdef WET_DRY
771 zeta_new(i,j)=zeta_new(i,j)+ &
772 & (dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
773# endif
774#endif
775 zwrk(i,j)=bkw_new*zeta_new(i,j)+ &
776 & bkw0*zeta(i,j,kstp)+ &
777 & bkw1*zeta(i,j,kbak)+ &
778 & bkw2*zeta(i,j,kold)
779
780#if defined VAR_RHO_2D && defined SOLVE3D
781 rzeta(i,j)=(1.0_r8+rhos(i,j))*zwrk(i,j)
782 rzeta2(i,j)=rzeta(i,j)*zwrk(i,j)
783 rzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
784#else
785 rzeta(i,j)=zwrk(i,j)
786 rzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
787#endif
788 END DO
789 END DO
790!
791! Apply mass point sources (volume vertical influx), if any.
792!
793! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
794!
795 IF (lwsrc(ng)) THEN
796 DO is=1,nsrc(ng)
797 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
798 i=sources(ng)%Isrc(is)
799 j=sources(ng)%Jsrc(is)
800 IF (((istrr.le.i).and.(i.le.iendr)).and. &
801 & ((jstrr.le.j).and.(j.le.jendr))) THEN
802 zeta_new(i,j)=zeta_new(i,j)+ &
803 & sources(ng)%Qbar(is)* &
804 & pm(i,j)*pn(i,j)*dtfast(ng)
805 END IF
806 END IF
807 END DO
808 END IF
809!
810! Apply boundary conditions to newly computed free-surface "zeta_new".
811! Notice that we are using the local routine, which passes the private
812! "zeta_new" array as argument.
813!
814! Here, we use the local "zetabc" since the private array "zeta_new"
815! is passed as an argument to allow computing the lateral boundary
816! conditions on the range IstrU-1:Iend and JstrV-1:Jend, so parallel
817! tile exchanges are avoided.
818!
819 CALL zetabc_local (ng, tile, &
820 & lbi, ubi, lbj, ubj, &
821 & imins, imaxs, jmins, jmaxs, &
822 & kstp, &
823 & zeta, &
824 & zeta_new)
825!
826! Load new computed free-surface into global state array.
827!
828 DO j=jstrr,jendr
829 DO i=istrr,iendr
830 zeta(i,j,knew)=zeta_new(i,j)
831 END DO
832 END DO
833
834#ifdef SOLVE3D
835!
836!-----------------------------------------------------------------------
837! Compute fast-time-averaged fields over all barotropic time steps.
838!-----------------------------------------------------------------------
839!
840! Reset/initialize arrays for averaged fields during the first
841! barotropic time step. Then, accumulate it time average. Include
842! physical boundary points, but not periodic ghost points or
843! computation distributed-memory computational margins.
844!
845 cff1=weight(1,iif(ng),ng)
846 cff2=weight(2,iif(ng),ng)
847!
848 IF (first_2d_step) THEN
849 DO j=jstrr,jendr
850 DO i=istrr,iendr
851 zt_avg1(i,j)=cff1*zeta(i,j,knew)
852 IF (i.ge.istr) THEN
853 du_avg1(i,j)=0.0_r8
854 du_avg2(i,j)=cff2*duon(i,j)
855 END IF
856 IF (j.ge.jstr) THEN
857 dv_avg1(i,j)=0.0_r8
858 dv_avg2(i,j)=cff2*dvom(i,j)
859 END IF
860 END DO
861 END DO
862 ELSE
863 DO j=jstrr,jendr
864 DO i=istrr,iendr
865 zt_avg1(i,j)=zt_avg1(i,j)+cff1*zeta(i,j,knew)
866 IF (i.ge.istr) THEN
867 du_avg2(i,j)=du_avg2(i,j)+cff2*duon(i,j)
868 END IF
869 IF (j.ge.jstr) THEN
870 dv_avg2(i,j)=dv_avg2(i,j)+cff2*dvom(i,j)
871 END IF
872 END DO
873 END DO
874 END IF
875#endif
876!
877!=======================================================================
878! Compute right-hand-side for the 2D momentum equations.
879!=======================================================================
880#ifdef SOLVE3D
881!
882! Notice that we are suppressing the computation of momentum advection,
883! Coriolis, and lateral viscosity terms in 3D Applications because
884! these terms are already included in the baroclinic-to-barotropic
885! forcing arrays "rufrc" and "rvfrc". It does not mean we are entirely
886! omitting them, but it is a choice between recomputing them at every
887! barotropic step or keeping them "frozen" during the fast-time
888! stepping.
889# ifdef STEP2D_CORIOLIS
890! However, in some coarse grid applications with larger baroclinic
891! timestep (say, DT around 20 minutes or larger), adding the Coriolis
892! term in the barotropic equations is useful since f*DT is no longer
893! small.
894# endif
895#endif
896!
897!-----------------------------------------------------------------------
898! Compute pressure-gradient terms.
899!-----------------------------------------------------------------------
900!
901 cff1=0.5_r8*g
902#if defined VAR_RHO_2D && defined SOLVE3D
903 cff2=0.333333333333_r8
904#endif
905#if defined ATM_PRESS && !defined SOLVE3D
906 cff3=0.5_r8*100.0_r8/rho0
907#endif
908 DO j=jstr,jend
909 DO i=istr,iend
910 IF (i.ge.istru) THEN
911 rubar(i,j)=cff1*on_u(i,j)* &
912 & ((h(i-1,j)+ &
913 & h(i ,j))* &
914 & (rzeta(i-1,j)- &
915 & rzeta(i ,j))+ &
916#if defined VAR_RHO_2D && defined SOLVE3D
917 & (h(i-1,j)- &
918 & h(i ,j))* &
919 & (rzetasa(i-1,j)+ &
920 & rzetasa(i ,j)+ &
921 & cff2*(rhoa(i-1,j)- &
922 & rhoa(i ,j))* &
923 & (zwrk(i-1,j)- &
924 & zwrk(i,j)))+ &
925#endif
926 & (rzeta2(i-1,j)- &
927 & rzeta2(i ,j)))
928#if defined ATM_PRESS && !defined SOLVE3D
929 rubar(i,j)=rubar(i,j)- &
930 & cff3*on_u(i,j)* &
931 & (h(i-1,j)+h(i,j)+ &
932 & rzeta(i-1,j)+rzeta(i,j))* &
933 & (pair(i,j)-pair(i-1,j))
934#endif
935#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
936 rubar(i,j)=rubar(i,j)- &
937 & cff1*on_u(i,j)* &
938 & (h(i-1,j)+h(i,j)+ &
939 & rzeta(i-1,j)+rzeta(i,j))* &
940 & (eq_tide(i,j)-eq_tide(i-1,j))
941#endif
942#ifdef DIAGNOSTICS_UV
943 diau2rhs(i,j,m2pgrd)=rubar(i,j)
944#endif
945 END IF
946!
947 IF (j.ge.jstrv) THEN
948 rvbar(i,j)=cff1*om_v(i,j)* &
949 & ((h(i,j-1)+ &
950 & h(i,j ))* &
951 & (rzeta(i,j-1)- &
952 & rzeta(i,j ))+ &
953#if defined VAR_RHO_2D && defined SOLVE3D
954 & (h(i,j-1)- &
955 & h(i,j ))* &
956 & (rzetasa(i,j-1)+ &
957 & rzetasa(i,j )+ &
958 & cff2*(rhoa(i,j-1)- &
959 & rhoa(i,j ))* &
960 & (zwrk(i,j-1)- &
961 & zwrk(i,j )))+ &
962#endif
963 & (rzeta2(i,j-1)- &
964 & rzeta2(i,j )))
965#if defined ATM_PRESS && !defined SOLVE3D
966 rvbar(i,j)=rvbar(i,j)- &
967 & cff3*om_v(i,j)* &
968 & (h(i,j-1)+h(i,j)+ &
969 & rzeta(i,j-1)+rzeta(i,j))* &
970 & (pair(i,j)-pair(i,j-1))
971#endif
972#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
973 rvbar(i,j)=rvbar(i,j)- &
974 & cff1*om_v(i,j)* &
975 & (h(i,j-1)+h(i,j)+ &
976 & rzeta(i,j-1)+rzeta(i,j))* &
977 & (eq_tide(i,j)-eq_tide(i,j-1))
978#endif
979#ifdef DIAGNOSTICS_UV
980 diav2rhs(i,j,m2pgrd)=rvbar(i,j)
981#endif
982 END IF
983 END DO
984 END DO
985
986#if defined UV_ADV && !defined SOLVE3D
987!
988!-----------------------------------------------------------------------
989! Add in horizontal advection of momentum.
990!-----------------------------------------------------------------------
991
992# ifdef UV_C2ADVECTION
993!
994! Second-order, centered differences advection fluxes.
995!
996 DO j=jstr,jend
997 DO i=istr-1,iend
998 IF (i.ge.istru-1) THEN
999 ufx(i,j)=0.25_r8* &
1000 & (duon(i,j)+duon(i+1,j))* &
1001 & (urhs(i ,j)+ &
1002 & urhs(i+1,j))
1003 END IF
1004!
1005 vfx(i+1,j)=0.25_r8* &
1006# ifdef MASKING
1007 & pmask(i+1,j)* &
1008# endif
1009 & (duon(i+1,j)+duon(i+1,j-1))* &
1010 & (vrhs(i+1,j)+ &
1011 & vrhs(i ,j))
1012 END DO
1013 END DO
1014!
1015 DO j=jstr-1,jend
1016 DO i=istr,iend
1017 IF (j.ge.jstrv-1) THEN
1018 vfe(i,j)=0.25_r8* &
1019 & (dvom(i,j)+dvom(i,j+1))* &
1020 & (vrhs(i,j )+ &
1021 & vrhs(i,j+1))
1022 END IF
1023!
1024 ufe(i,j+1)=0.25_r8* &
1025# ifdef MASKING
1026 & pmask(i,j+1)* &
1027# endif
1028 & (dvom(i,j+1)+dvom(i-1,j+1))* &
1029 & (urhs(i,j+1)+ &
1030 & urhs(i,j ))
1031 END DO
1032 END DO
1033
1034# elif defined UV_C4ADVECTION
1035!
1036! Fourth-order, centered differences u-momentum advection fluxes.
1037!
1038 DO j=jstr,jend
1039 DO i=istrum1,iendp1
1040 grad(i,j)=urhs(i-1,j)-2.0_r8*urhs(i,j)+ &
1041 & urhs(i+1,j)
1042 dgrad(i,j)=duon(i-1,j)-2.0_r8*duon(i,j)+duon(i+1,j)
1043 END DO
1044 END DO
1045 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1046 IF (domain(ng)%Western_Edge(tile)) THEN
1047 DO j=jstr,jend
1048 grad(istr,j)=grad(istr+1,j)
1049 dgrad(istr,j)=dgrad(istr+1,j)
1050 END DO
1051 END IF
1052 END IF
1053 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1054 IF (domain(ng)%Eastern_Edge(tile)) THEN
1055 DO j=jstr,jend
1056 grad(iend+1,j)=grad(iend,j)
1057 dgrad(iend+1,j)=dgrad(iend,j)
1058 END DO
1059 END IF
1060 END IF
1061! d/dx(Duu/n)
1062 cff=1.0_r8/6.0_r8
1063 DO j=jstr,jend
1064 DO i=istru-1,iend
1065 ufx(i,j)=0.25_r8*(urhs(i ,j)+ &
1066 & urhs(i+1,j)- &
1067 & cff*(grad(i,j)+grad(i+1,j)))* &
1068 & (duon(i,j)+duon(i+1,j)- &
1069 & cff*(dgrad(i,j)+dgrad(i+1,j)))
1070 END DO
1071 END DO
1072!
1073 DO j=jstrm1,jendp1
1074 DO i=istru,iend
1075 grad(i,j)=urhs(i,j-1)-2.0_r8*urhs(i,j)+ &
1076 & urhs(i,j+1)
1077 END DO
1078 END DO
1079 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1080 IF (domain(ng)%Southern_Edge(tile)) THEN
1081 DO i=istru,iend
1082 grad(i,jstr-1)=grad(i,jstr)
1083 END DO
1084 END IF
1085 END IF
1086 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1087 IF (domain(ng)%Northern_Edge(tile)) THEN
1088 DO i=istru,iend
1089 grad(i,jend+1)=grad(i,jend)
1090 END DO
1091 END IF
1092 END IF
1093 DO j=jstr,jend+1
1094 DO i=istru-1,iend
1095 dgrad(i,j)=dvom(i-1,j)-2.0_r8*dvom(i,j)+dvom(i+1,j)
1096 END DO
1097 END DO
1098! d/dy(Duv/m)
1099 cff=1.0_r8/6.0_r8
1100 DO j=jstr,jend+1
1101 DO i=istru,iend
1102 ufe(i,j)=0.25_r8*(urhs(i,j )+ &
1103 & urhs(i,j-1)- &
1104 & cff*(grad(i,j)+grad(i,j-1)))* &
1105 & (dvom(i,j)+dvom(i-1,j)- &
1106 & cff*(dgrad(i,j)+dgrad(i-1,j)))
1107 END DO
1108 END DO
1109!
1110! Fourth-order, centered differences v-momentum advection fluxes.
1111!
1112 DO j=jstrv,jend
1113 DO i=istrm1,iendp1
1114 grad(i,j)=vrhs(i-1,j)-2.0_r8*vrhs(i,j)+ &
1115 & vrhs(i+1,j)
1116 END DO
1117 END DO
1118 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1119 IF (domain(ng)%Western_Edge(tile)) THEN
1120 DO j=jstrv,jend
1121 grad(istr-1,j)=grad(istr,j)
1122 END DO
1123 END IF
1124 END IF
1125 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1126 IF (domain(ng)%Eastern_Edge(tile)) THEN
1127 DO j=jstrv,jend
1128 grad(iend+1,j)=grad(iend,j)
1129 END DO
1130 END IF
1131 END IF
1132 DO j=jstrv-1,jend
1133 DO i=istr,iend+1
1134 dgrad(i,j)=duon(i,j-1)-2.0_r8*duon(i,j)+duon(i,j+1)
1135 END DO
1136 END DO
1137! d/dx(Duv/n)
1138 cff=1.0_r8/6.0_r8
1139 DO j=jstrv,jend
1140 DO i=istr,iend+1
1141 vfx(i,j)=0.25_r8*(vrhs(i ,j)+ &
1142 & vrhs(i-1,j)- &
1143 & cff*(grad(i,j)+grad(i-1,j)))* &
1144 & (duon(i,j)+duon(i,j-1)- &
1145 & cff*(dgrad(i,j)+dgrad(i,j-1)))
1146 END DO
1147 END DO
1148!
1149 DO j=jstrvm1,jendp1
1150 DO i=istr,iend
1151 grad(i,j)=vrhs(i,j-1)-2.0_r8*vrhs(i,j)+ &
1152 & vrhs(i,j+1)
1153 dgrad(i,j)=dvom(i,j-1)-2.0_r8*dvom(i,j)+dvom(i,j+1)
1154 END DO
1155 END DO
1156 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1157 IF (domain(ng)%Southern_Edge(tile)) THEN
1158 DO i=istr,iend
1159 grad(i,jstr)=grad(i,jstr+1)
1160 dgrad(i,jstr)=dgrad(i,jstr+1)
1161 END DO
1162 END IF
1163 END IF
1164 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1165 IF (domain(ng)%Northern_Edge(tile)) THEN
1166 DO i=istr,iend
1167 grad(i,jend+1)=grad(i,jend)
1168 dgrad(i,jend+1)=dgrad(i,jend)
1169 END DO
1170 END IF
1171 END IF
1172! d/dy(Dvv/m)
1173 cff=1.0_r8/6.0_r8
1174 DO j=jstrv-1,jend
1175 DO i=istr,iend
1176 vfe(i,j)=0.25_r8*(vrhs(i,j )+ &
1177 & vrhs(i,j+1)- &
1178 & cff*(grad(i,j)+grad(i,j+1)))* &
1179 & (dvom(i,j)+dvom(i,j+1)- &
1180 & cff*(dgrad(i,j)+dgrad(i,j+1)))
1181 END DO
1182 END DO
1183# endif
1184!
1185! Add advection to RHS terms.
1186!
1187 DO j=jstr,jend
1188 DO i=istr,iend
1189 IF (i.ge.istru) THEN
1190 cff1=ufx(i,j)-ufx(i-1,j)
1191 cff2=ufe(i,j+1)-ufe(i,j)
1192 fac1=cff1+cff2
1193 rubar(i,j)=rubar(i,j)-fac1
1194# if defined DIAGNOSTICS_UV
1195 diau2rhs(i,j,m2xadv)=-cff1
1196 diau2rhs(i,j,m2yadv)=-cff2
1197 diau2rhs(i,j,m2hadv)=-fac1
1198# endif
1199 END IF
1200!
1201 IF (j.ge.jstrv) THEN
1202 cff3=vfx(i+1,j)-vfx(i,j)
1203 cff4=vfe(i,j)-vfe(i,j-1)
1204 fac2=cff3+cff4
1205 rvbar(i,j)=rvbar(i,j)-fac2
1206# if defined DIAGNOSTICS_UV
1207 diav2rhs(i,j,m2xadv)=-cff3
1208 diav2rhs(i,j,m2yadv)=-cff4
1209 diav2rhs(i,j,m2hadv)=-fac2
1210# endif
1211 END IF
1212 END DO
1213 END DO
1214#endif
1215
1216#if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
1217!
1218!-----------------------------------------------------------------------
1219! Add in Coriolis term.
1220!-----------------------------------------------------------------------
1221!
1222 DO j=jstrv-1,jend
1223 DO i=istru-1,iend
1224 cff=0.5_r8*drhs(i,j)*fomn(i,j)
1225 ufx(i,j)=cff*(vrhs(i,j )+ &
1226 & vrhs(i,j+1))
1227 vfe(i,j)=cff*(urhs(i ,j)+ &
1228 & urhs(i+1,j))
1229 END DO
1230 END DO
1231!
1232 DO j=jstr,jend
1233 DO i=istr,iend
1234 IF (i.ge.istru) THEN
1235 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1236 rubar(i,j)=rubar(i,j)+fac1
1237# if defined DIAGNOSTICS_UV
1238 diau2rhs(i,j,m2fcor)=fac1
1239# endif
1240
1241 END IF
1242!
1243 IF (j.ge.jstrv) THEN
1244 fac2=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1245 rvbar(i,j)=rvbar(i,j)-fac2
1246# if defined DIAGNOSTICS_UV
1247 diav2rhs(i,j,m2fcor)=-fac2
1248# endif
1249 END IF
1250 END DO
1251 END DO
1252#endif
1253
1254#if (defined CURVGRID && defined UV_ADV) && !defined SOLVE3D
1255!
1256!-----------------------------------------------------------------------
1257! Add in curvilinear transformation terms.
1258!-----------------------------------------------------------------------
1259!
1260 DO j=jstrv-1,jend
1261 DO i=istru-1,iend
1262 cff1=0.5_r8*(vrhs(i,j )+ &
1263 & vrhs(i,j+1))
1264 cff2=0.5_r8*(urhs(i ,j)+
1265 & urhs(i+1,j))
1266 cff3=cff1*dndx(i,j)
1267 cff4=cff2*dmde(i,j)
1268 cff=drhs(i,j)*(cff3-cff4)
1269 ufx(i,j)=cff*cff1
1270 vfe(i,j)=cff*cff2
1271# if defined DIAGNOSTICS_UV
1272 cff=drhs(i,j)*cff4
1273 uwrk(i,j)=-cff*cff1 ! ubar equation, ETA-term
1274 vwrk(i,j)=-cff*cff2 ! vbar equation, ETA-term
1275# endif
1276 END DO
1277 END DO
1278!
1279 DO j=jstr,jend
1280 DO i=istr,iend
1281 IF (i.ge.istru) THEN
1282 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1283 rubar(i,j)=rubar(i,j)+fac1
1284# if defined DIAGNOSTICS_UV
1285 fac2=0.5_r8*(uwrk(i,j)+uwrk(i-1,j))
1286 diau2rhs(i,j,m2xadv)=diau2rhs(i,j,m2xadv)+fac1-fac2
1287 diau2rhs(i,j,m2yadv)=diau2rhs(i,j,m2yadv)+fac2
1288 diau2rhs(i,j,m2hadv)=diau2rhs(i,j,m2hadv)+fac1
1289# endif
1290 END IF
1291!
1292 IF (j.ge.jstrv) THEN
1293 fac1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1294 rvbar(i,j)=rvbar(i,j)-fac1
1295# if defined DIAGNOSTICS_UV
1296 fac2=0.5_r8*(vwrk(i,j)+vwrk(i,j-1))
1297 diav2rhs(i,j,m2xadv)=diav2rhs(i,j,m2xadv)-fac1+fac2
1298 diav2rhs(i,j,m2yadv)=diav2rhs(i,j,m2yadv)-fac2
1299 diav2rhs(i,j,m2hadv)=diav2rhs(i,j,m2hadv)-fac1
1300# endif
1301 END IF
1302 END DO
1303 END DO
1304#endif
1305
1306#if defined UV_VIS2 && !defined SOLVE3D
1307!
1308!-----------------------------------------------------------------------
1309! Add in horizontal harmonic viscosity.
1310!-----------------------------------------------------------------------
1311!
1312! Compute total depth at PSI-points.
1313!
1314 DO j=jstr,jend+1
1315 DO i=istr,iend+1
1316 drhs_p(i,j)=0.25_r8*(drhs(i,j )+drhs(i-1,j )+ &
1317 & drhs(i,j-1)+drhs(i-1,j-1))
1318 END DO
1319 END DO
1320!
1321! Compute flux-components of the horizontal divergence of the stress
1322! tensor (m5/s2) in XI- and ETA-directions.
1323!
1324 DO j=jstrv-1,jend
1325 DO i=istru-1,iend
1326 cff=visc2_r(i,j)*drhs(i,j)*0.5_r8* &
1327 & (pmon_r(i,j)* &
1328 & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,kstp)- &
1329 & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,kstp))- &
1330 & pnom_r(i,j)* &
1331 & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,kstp)- &
1332 & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,kstp)))
1333 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
1334 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
1335 END DO
1336 END DO
1337!
1338 DO j=jstr,jend+1
1339 DO i=istr,iend+1
1340 cff=visc2_p(i,j)*drhs_p(i,j)*0.5_r8* &
1341 & (pmon_p(i,j)* &
1342 & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,kstp)- &
1343 & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,kstp))+ &
1344 & pnom_p(i,j)* &
1345 & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,kstp)- &
1346 & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,kstp)))
1347# ifdef MASKING
1348 cff=cff*pmask(i,j)
1349# endif
1350# ifdef WET_DRY
1351 cff=cff*pmask_wet(i,j)
1352# endif
1353 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
1354 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
1355 END DO
1356 END DO
1357!
1358! Add in harmonic viscosity.
1359!
1360 DO j=jstr,jend
1361 DO i=istr,iend
1362 IF (i.ge.istru) THEN
1363 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
1364 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
1365 fac1=cff1+cff2
1366 rubar(i,j)=rubar(i,j)+fac1
1367# if defined DIAGNOSTICS_UV
1368 diau2rhs(i,j,m2hvis)=fac1
1369 diau2rhs(i,j,m2xvis)=cff1
1370 diau2rhs(i,j,m2yvis)=cff2
1371# endif
1372 END IF
1373!
1374 IF (j.ge.jstrv) THEN
1375 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
1376 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
1377 fac1=cff1-cff2
1378 rvbar(i,j)=rvbar(i,j)+fac1
1379# if defined DIAGNOSTICS_UV
1380 diav2rhs(i,j,m2hvis)=fac1
1381 diav2rhs(i,j,m2xvis)= cff1
1382 diav2rhs(i,j,m2yvis)=-cff2
1383# endif
1384 END IF
1385 END DO
1386 END DO
1387#endif
1388
1389#ifdef SOLVE3D
1390!
1391!-----------------------------------------------------------------------
1392! Coupling between 2D and 3D equations.
1393!-----------------------------------------------------------------------
1394!
1395! Before the first barotropic time step, arrays "rufrc" and "rvfrc"
1396! contain vertical integrals of the 3D right-hand-side terms for the
1397! momentum equations (including surface and bottom stresses). During
1398! the first barotropic time step, convert them into forcing terms by
1399! subtracting the fast-time "rubar" and "rvbar" from them.
1400!
1401! In the predictor-coupled mode, the resultant forcing terms "rufrc"
1402! and "rvfrc" are extrapolated forward in time, so they become
1403! centered effectively at time n+1/2. This is done using optimized
1404! Adams-Bashforth weights. In the code below, rufrc_bak(:,:,nstp) is
1405! at (n-1)time step, while rufrc_bak(:,:,3-nstp) is at (n-2). After
1406! its use as input, the latter is overwritten by the value at time
1407! step "nstp" (mathematically "n") during the next step.
1408!
1409! From now on, the computed forcing terms "rufrc" and "rvfrc" will
1410! remain constant during the fast-time stepping and will be added
1411! to "rubar" and "rvbar" during all subsequent barotropic steps.
1412!
1413 coupled_step : IF (first_2d_step) THEN
1414!
1415! Predictor coupled barotropic mode: Set coefficients for AB3-like
1416! forward-in-time extrapolation of 3D to 2D forcing terms "rufrc" and
1417! "rvfrc".
1418!
1419 IF (iic(ng).eq.ntstart(ng)) THEN
1420 cfwd0=1.0_r8
1421 cfwd1=0.0_r8
1422 cfwd2=0.0_r8
1423 ELSE IF (iic(ng).eq.ntstart(ng)+1) THEN
1424 cfwd0=1.5_r8
1425 cfwd1=-0.5_r8
1426 cfwd2=0.0_r8
1427 ELSE
1428 cfwd2=0.281105_r8
1429 cfwd1=-0.5_r8-2.0_r8*cfwd2
1430 cfwd0=1.5_r8+cfwd2
1431 END IF
1432!
1433 DO j=jstr,jend
1434 DO i=istr,iend
1435!
1436! Compensate for (cancel) bottom drag terms: at input into step2d
1437! "rufrc" and "rvfrc" contain bottom drag terms computed by 3D mode.
1438! However, there are no 2D counterparts in "rubar" and "rvbar" because
1439! 2D bottom drag will be computed implicitly during the final stage of
1440! updating ubar(:,:,knew) and vbar(:,:,knew) below. Note that unlike
1441! the other terms, the bottom drag should not be extrapolated forward,
1442! if "rufrc" and "rvfrc" are, so this cancelation needs to be done
1443! right now rather than at the bottom of this loop.
1444!
1445 IF (i.ge.istru) THEN
1446 rufrc(i,j)=rufrc(i,j)+ &
1447 & 0.5_r8*(rdrag(i,j)+rdrag(i-1,j))* &
1448 & om_u(i,j)*on_u(i,j)*ubar(i,j,kstp)
1449 END IF
1450!
1451 IF (j.ge.jstrv) THEN
1452 rvfrc(i,j)=rvfrc(i,j)+ &
1453 & 0.5_r8*(rdrag(i,j)+rdrag(i,j-1))* &
1454 & om_v(i,j)*on_v(i,j)*vbar(i,j,kstp)
1455 END IF
1456!
1457! Barotropic mode running predictor stage: forward extrapolation.
1458!
1459 IF (i.ge.istru) THEN
1460 cff1=rufrc(i,j)-rubar(i,j)
1461 rufrc(i,j)=cfwd0*cff1+ &
1462 & cfwd1*rufrc_bak(i,j, nstp)+ &
1463 & cfwd2*rufrc_bak(i,j,3-nstp)
1464 rufrc_bak(i,j,3-nstp)=cff1
1465 END IF
1466!
1467 IF (j.ge.jstrv) THEN
1468 cff2=rvfrc(i,j)-rvbar(i,j)
1469 rvfrc(i,j)=cfwd0*cff2+ &
1470 & cfwd1*rvfrc_bak(i,j, nstp)+ &
1471 & cfwd2*rvfrc_bak(i,j,3-nstp)
1472 rvfrc_bak(i,j,3-nstp)=cff2
1473 END IF
1474 END DO
1475 END DO
1476!
1477! Add correction term to shift pressure-gradient terms from "kstp" to
1478! "knew". That is, it converts the first 2D step from Forward-Euler
1479! to Forward-Backward (this is PGF_FB_CORRECTION mentioned above).
1480!
1481 DO j=jstrv-1,jend
1482 DO i=istru-1,iend
1483 zwrk(i,j)=zeta_new(i,j)-zeta(i,j,kstp)
1484# if defined VAR_RHO_2D && defined SOLVE3D
1485 rzeta(i,j)=(1.0_r8+rhos(i,j))*zwrk(i,j)
1486 rzeta2(i,j)=rzeta(i,j)*(zeta_new(i,j)+zeta(i,j,kstp))
1487 rzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
1488# else
1489 rzeta(i,j)=zwrk(i,j)
1490 rzeta2(i,j)=zwrk(i,j)*(zeta_new(i,j)+zeta(i,j,kstp))
1491# endif
1492 END DO
1493 END DO
1494!
1495 cff1=0.5*g
1496# if defined VAR_RHO_2D && defined SOLVE3D
1497 cff2=0.333333333333_r8
1498# endif
1499 DO j=jstr,jend
1500 DO i=istr,iend
1501 IF (i.ge.istru) THEN
1502 rubar(i,j)=rubar(i,j)+ &
1503 & cff1*on_u(i,j)* &
1504 & ((h(i-1,j)+ &
1505 & h(i ,j))* &
1506 & (rzeta(i-1,j)- &
1507 & rzeta(i ,j))+ &
1508# if defined VAR_RHO_2D && defined SOLVE3D
1509 & (h(i-1,j)- &
1510 & h(i ,j))* &
1511 & (rzetasa(i-1,j)+ &
1512 & rzetasa(i ,j)+ &
1513 & cff2*(rhoa(i-1,j)- &
1514 & rhoa(i ,j))* &
1515 & (zwrk(i-1,j)- &
1516 & zwrk(i ,j)))+ &
1517# endif
1518 & (rzeta2(i-1,j)- &
1519 & rzeta2(i ,j)))
1520# ifdef DIAGNOSTICS_UV
1521 diau2rhs(i,j,m2pgrd)=diau2rhs(i,j,m2pgrd)+ &
1522 & rubar(i,j)
1523# endif
1524 END IF
1525!
1526 IF (j.ge.jstrv) THEN
1527 rvbar(i,j)=rvbar(i,j)+ &
1528 & cff1*om_v(i,j)* &
1529 & ((h(i,j-1)+ &
1530 & h(i,j ))* &
1531 & (rzeta(i,j-1)- &
1532 & rzeta(i,j ))+ &
1533# if defined VAR_RHO_2D && defined SOLVE3D
1534 & (h(i,j-1)- &
1535 & h(i,j ))* &
1536 & (rzetasa(i,j-1)+ &
1537 & rzetasa(i,j )+ &
1538 & cff2*(rhoa(i,j-1)- &
1539 & rhoa(i,j ))* &
1540 & (zwrk(i,j-1)- &
1541 & zwrk(i,j )))+ &
1542# endif
1543 & (rzeta2(i,j-1)- &
1544 & rzeta2(i,j )))
1545# ifdef DIAGNOSTICS_UV
1546 diav2rhs(i,j,m2pgrd)=diav2rhs(i,j,m2pgrd)+ &
1547 & rvbar(i,j)
1548# endif
1549 END IF
1550 END DO
1551 END DO
1552 END IF coupled_step
1553#endif
1554!
1555!-----------------------------------------------------------------------
1556! Time step 2D momentum equations.
1557!-----------------------------------------------------------------------
1558!
1559! Advance 2D momentum components while simultaneously adding them to
1560! accumulate fast-time averages to compute barotropic fluxes. Doing so
1561! straight away yields a more computationally dense code. However, the
1562! fast-time averaged fluxes (DU_avg1 and DV_avg1) are needed both at
1563! the interior and physical boundary points. Thus, we need separate
1564! loops along the domain boundaries after setting "ubar" and "vbar"
1565! lateral boundary conditions. Also, note that bottom drag is treated
1566! implicitly:
1567!
1568! Dnew*ubar(:,:,m+1) = Dold*ubar(:,:,m) +
1569! dtfast(ng)*rhs2D(:,:) -
1570! dtfast(ng)*rdrag(:,:)*ubar(:,:,m+1)
1571! hence
1572!
1573! ubar(:,:,m+1)=[Dold * ubar(..,m) + dtfast(ng) * rhs2D(:,:)] /
1574! [Dnew + dtfast(ng) * rdrag(:,:)]
1575!
1576! DU_avg1 = DU_avg1 +
1577! weight(m+1) * Dnew * ubar(:,:,m+1) * on_u(:,:)
1578!
1579! where it should be noted that Dnew .ne. Dnew + dtfast * rdrag
1580!
1581 DO j=jstrv-1,jend
1582 DO i=istru-1,iend
1583 dnew(i,j)=h(i,j)+zeta_new(i,j)
1584 dnew_rd(i,j)=dnew(i,j)+dtfast(ng)*rdrag(i,j)
1585 dstp(i,j)=h(i,j)+zeta(i,j,kstp)
1586 END DO
1587 END DO
1588
1589#if defined UV_QDRAG && !defined SOLVE3D
1590!
1591! Add quadratic drag term associated in shallow-water applications.
1592!
1593! Here, the SQRT(3) is due to a linear interpolation with second order
1594! accuaracy that ensures positive and negative values of the velocity
1595! components:
1596!
1597! u^2(i+1/2) = (1/3)*[u(i)*u(i) + u(i)*u(i+1) + u(i+1)*u(i+1)]
1598!
1599! If u(i)=1 and u(i+1)=-1, then u^2(i+1/2)=1/3 as it should be.
1600!
1601 cff=dtfast(ng)/sqrt(3.0_r8)
1602 DO j=jstrv-1,jend
1603 DO i=istru-1,iend
1604 cff1=ubar(i ,j,kstp)**2+ &
1605 & ubar(i+1,j,kstp)**2+ &
1606 & ubar(i ,j,kstp)*ubar(i+1,j,kstp)+ &
1607 & vbar(i,j ,kstp)**2+ &
1608 & vbar(i,j+1,kstp)**2+ &
1609 & vbar(i,j ,kstp)*vbar(i,j+1,kstp)
1610 cff2=sqrt(cff1)
1611 dnew_rd(i,j)=dnew_rd(i,j)+ &
1612 & cff*rdrag2(i,j)*cff2
1613 END DO
1614 END DO
1615#endif
1616!
1617! Step 2D momentum equations.
1618!
1619 cff=0.5_r8*dtfast(ng)
1620#ifdef SOLVE3D
1621 cff1=0.5_r8*weight(1,iif(ng),ng)
1622#else
1623 cff2=2.0_r8*dtfast(ng)
1624#endif
1625 DO j=jstr,jend
1626 DO i=istru,iend
1627 cff3=cff*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
1628 fac1=1.0_r8/(dnew_rd(i,j)+dnew_rd(i-1,j))
1629 ubar(i,j,knew)=fac1*((dstp(i,j)+dstp(i-1,j))*ubar(i,j,kstp)+ &
1630#ifdef SOLVE3D
1631 & cff3*(rubar(i,j)+rufrc(i,j)))
1632#else
1633 & cff3*rubar(i,j)+cff2*sustr(i,j))
1634#endif
1635#ifdef MASKING
1636 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
1637#endif
1638#ifdef WET_DRY
1639 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
1640 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
1641 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
1642 ubar(i,j,knew)=ubar(i,j,knew)*cff7
1643#endif
1644#ifdef SOLVE3D
1645 du_avg1(i,j)=du_avg1(i,j)+ &
1646 & cff1*on_u(i,j)* &
1647 & (dnew(i,j)+dnew(i-1,j))*ubar(i,j,knew)
1648#endif
1649#if defined NESTING && !defined SOLVE3D
1650 du_flux(i,j)=0.5_r8*on_u(i,j)* &
1651 & (dnew(i,j)+dnew(i-1,j))*ubar(i,j,knew)
1652#endif
1653 END DO
1654 END DO
1655!
1656 DO j=jstrv,jend
1657 DO i=istr,iend
1658 cff3=cff*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
1659 fac2=1.0_r8/(dnew_rd(i,j)+dnew_rd(i,j-1))
1660 vbar(i,j,knew)=fac2*((dstp(i,j)+dstp(i,j-1))*vbar(i,j,kstp)+ &
1661#ifdef SOLVE3D
1662 & cff3*(rvbar(i,j)+rvfrc(i,j)))
1663#else
1664 & cff3*rvbar(i,j)+cff2*svstr(i,j))
1665#endif
1666#ifdef MASKING
1667 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
1668#endif
1669#ifdef WET_DRY
1670 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
1671 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
1672 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
1673 vbar(i,j,knew)=vbar(i,j,knew)*cff7
1674#endif
1675#ifdef SOLVE3D
1676 dv_avg1(i,j)=dv_avg1(i,j)+ &
1677 & cff1*om_v(i,j)* &
1678 & (dnew(i,j)+dnew(i,j-1))*vbar(i,j,knew)
1679#endif
1680#if defined NESTING && !defined SOLVE3D
1681 dv_flux(i,j)=0.5_r8*om_v(i,j)* &
1682 & (dnew(i,j)+dnew(i,j-1))*vbar(i,j,knew)
1683#endif
1684 END DO
1685 END DO
1686!
1687! Apply lateral boundary conditions.
1688!
1689 CALL u2dbc_tile (ng, tile, &
1690 & lbi, ubi, lbj, ubj, &
1691 & imins, imaxs, jmins, jmaxs, &
1692 & krhs, kstp, knew, &
1693 & ubar, vbar, zeta)
1694 CALL v2dbc_tile (ng, tile, &
1695 & lbi, ubi, lbj, ubj, &
1696 & imins, imaxs, jmins, jmaxs, &
1697 & krhs, kstp, knew, &
1698 & ubar, vbar, zeta)
1699!
1700! Compute integral mass flux across open boundaries and adjust
1701! for volume conservation.
1702!
1703 IF (any(volcons(:,ng))) THEN
1704 CALL obc_flux_tile (ng, tile, &
1705 & lbi, ubi, lbj, ubj, &
1706 & imins, imaxs, jmins, jmaxs, &
1707 & knew, &
1708#ifdef MASKING
1709 & umask, vmask, &
1710#endif
1711 & h, om_v, on_u, &
1712 & ubar, vbar, zeta)
1713 END IF
1714
1715#if defined SOLVE3D || (defined NESTING && !defined SOLVE3D)
1716!
1717! Set barotropic fluxes along physical boundaries.
1718!
1719 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1720 IF (domain(ng)%Western_Edge(tile)) THEN
1721 DO j=jstr-1,jendr
1722 dnew(istr-1,j)=h(istr-1,j)+zeta_new(istr-1,j)
1723 END DO
1724 END IF
1725 END IF
1726 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1727 IF (domain(ng)%Eastern_Edge(tile)) THEN
1728 DO j=jstr-1,jendr
1729 dnew(iend+1,j)=h(iend+1,j)+zeta_new(iend+1,j)
1730 END DO
1731 END IF
1732 END IF
1733 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1734 IF (domain(ng)%Southern_Edge(tile)) THEN
1735 DO i=istr-1,iendr
1736 dnew(i,jstr-1)=h(i,jstr-1)+zeta_new(i,jstr-1)
1737 END DO
1738 END IF
1739 END IF
1740 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1741 IF (domain(ng)%Northern_Edge(tile)) THEN
1742 DO i=istr-1,iendr
1743 dnew(i,jend+1)=h(i,jend+1)+zeta_new(i,jend+1)
1744 END DO
1745 END IF
1746 END IF
1747
1748# ifdef SOLVE3D
1749!
1750 cff1=0.5*weight(1,iif(ng),ng)
1751# endif
1752!
1753 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1754 IF (domain(ng)%Western_Edge(tile)) THEN
1755 DO j=jstrr,jendr
1756# if defined NESTING && !defined SOLVE3D
1757 du_flux(istru-1,j)=0.5_r8*on_u(istru-1,j)* &
1758 & (dnew(istru-1,j)+dnew(istru-2,j))* &
1759 & ubar(istru-1,j,knew)
1760# else
1761 du_avg1(istru-1,j)=du_avg1(istru-1,j)+ &
1762 & cff1*on_u(istru-1,j)* &
1763 & (dnew(istru-1,j)+dnew(istru-2,j))* &
1764 & ubar(istru-1,j,knew)
1765# endif
1766 END DO
1767 DO j=jstrv,jend
1768# if defined NESTING && !defined SOLVE3D
1769 dv_flux(istr-1,j)=0.5_r8*om_v(istr-1,j)* &
1770 & (dnew(istr-1,j)+dnew(istr-1,j-1))* &
1771 & vbar(istr-1,j,knew)
1772# else
1773 dv_avg1(istr-1,j)=dv_avg1(istr-1,j)+ &
1774 & cff1*om_v(istr-1,j)* &
1775 & (dnew(istr-1,j)+dnew(istr-1,j-1))* &
1776 & vbar(istr-1,j,knew)
1777# endif
1778 END DO
1779 END IF
1780 END IF
1781 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1782 IF (domain(ng)%Eastern_Edge(tile)) THEN
1783 DO j=jstrr,jendr
1784# if defined NESTING && !defined SOLVE3D
1785 du_flux(iend+1,j)=0.5_r8*on_u(iend+1,j)* &
1786 & (dnew(iend+1,j)+dnew(iend,j))* &
1787 & ubar(iend+1,j,knew)
1788# else
1789 du_avg1(iend+1,j)=du_avg1(iend+1,j)+ &
1790 & cff1*on_u(iend+1,j)* &
1791 & (dnew(iend+1,j)+dnew(iend,j))* &
1792 & ubar(iend+1,j,knew)
1793# endif
1794 END DO
1795 DO j=jstrv,jend
1796# if defined NESTING && !defined SOLVE3D
1797 dv_flux(iend+1,j)=0.5_r8*om_v(iend+1,j)* &
1798 & (dnew(iend+1,j)+dnew(iend+1,j-1))* &
1799 & vbar(iend+1,j,knew)
1800# else
1801 dv_avg1(iend+1,j)=dv_avg1(iend+1,j)+ &
1802 & cff1*om_v(iend+1,j)* &
1803 & (dnew(iend+1,j)+dnew(iend+1,j-1))* &
1804 & vbar(iend+1,j,knew)
1805# endif
1806 END DO
1807 END IF
1808 END IF
1809 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1810 IF (domain(ng)%Southern_Edge(tile)) THEN
1811 DO i=istru,iend
1812# if defined NESTING && !defined SOLVE3D
1813 du_flux(i,jstr-1)=0.5_r8*on_u(i,jstr-1)* &
1814 & (dnew(i,jstr-1)+dnew(i-1,jstr-1))* &
1815 & ubar(i,jstr-1,knew)
1816# else
1817 du_avg1(i,jstr-1)=du_avg1(i,jstr-1)+ &
1818 & cff1*on_u(i,jstr-1)* &
1819 & (dnew(i,jstr-1)+dnew(i-1,jstr-1))* &
1820 & ubar(i,jstr-1,knew)
1821# endif
1822 END DO
1823 DO i=istrr,iendr
1824# if defined NESTING && !defined SOLVE3D
1825 dv_flux(i,jstrv-1)=0.5_r8*om_v(i,jstrv-1)* &
1826 & (dnew(i,jstrv-1)+dnew(i,jstrv-2))* &
1827 & vbar(i,jstrv-1,knew)
1828# else
1829 dv_avg1(i,jstrv-1)=dv_avg1(i,jstrv-1)+ &
1830 & cff1*om_v(i,jstrv-1)* &
1831 & (dnew(i,jstrv-1)+dnew(i,jstrv-2))* &
1832 & vbar(i,jstrv-1,knew)
1833# endif
1834 END DO
1835 END IF
1836 END IF
1837 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1838 IF (domain(ng)%Northern_Edge(tile)) THEN
1839 DO i=istru,iend
1840# if defined NESTING && !defined SOLVE3D
1841 du_flux(i,jend+1)=0.5_r8*on_u(i,jend+1)* &
1842 & (dnew(i,jend+1)+dnew(i-1,jend+1))* &
1843 & ubar(i,jend+1,knew)
1844# else
1845 du_avg1(i,jend+1)=du_avg1(i,jend+1)+ &
1846 & cff1*on_u(i,jend+1)* &
1847 & (dnew(i,jend+1)+dnew(i-1,jend+1))* &
1848 & ubar(i,jend+1,knew)
1849# endif
1850 END DO
1851 DO i=istrr,iendr
1852# if defined NESTING && !defined SOLVE3D
1853 dv_flux(i,jend+1)=0.5_r8*om_v(i,jend+1)* &
1854 & (dnew(i,jend+1)+dnew(i,jend))* &
1855 & vbar(i,jend+1,knew)
1856# else
1857 dv_avg1(i,jend+1)=dv_avg1(i,jend+1)+ &
1858 & cff1*om_v(i,jend+1)* &
1859 & (dnew(i,jend+1)+dnew(i,jend))* &
1860 & vbar(i,jend+1,knew)
1861# endif
1862 END DO
1863 END IF
1864 END IF
1865#endif
1866!
1867!-----------------------------------------------------------------------
1868! Apply momentum transport point sources (like river runoff), if any.
1869!
1870! Dsrc(is) = 0, flow across grid cell u-face (positive or negative)
1871! Dsrc(is) = 1, flow across grid cell v-face (positive or negative)
1872!-----------------------------------------------------------------------
1873!
1874 IF (luvsrc(ng)) THEN
1875 DO is=1,nsrc(ng)
1876 i=sources(ng)%Isrc(is)
1877 j=sources(ng)%Jsrc(is)
1878 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1879 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1880 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
1881 cff=1.0_r8/(on_u(i,j)* &
1882 & 0.5_r8*(dnew(i-1,j)+dnew(i,j)))
1883 ubar(i,j,knew)=sources(ng)%Qbar(is)*cff
1884#ifdef SOLVE3D
1885 du_avg1(i,j)=sources(ng)%Qbar(is)
1886#endif
1887#if defined NESTING && !defined SOLVE3D
1888 du_flux(i,j)=sources(ng)%Qbar(is)
1889#endif
1890 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
1891 cff=1.0_r8/(om_v(i,j)* &
1892 & 0.5_r8*(dnew(i,j-1)+dnew(i,j)))
1893 vbar(i,j,knew)=sources(ng)%Qbar(is)*cff
1894#ifdef SOLVE3D
1895 dv_avg1(i,j)=sources(ng)%Qbar(is)
1896#endif
1897#if defined NESTING && !defined SOLVE3D
1898 dv_flux(i,j)=sources(ng)%Qbar(is)
1899#endif
1900 END IF
1901 END IF
1902 END DO
1903 END IF
1904!
1905! Deallocate local new free-surface.
1906!
1907 deallocate ( zeta_new )
1908
1909#ifdef WET_DRY
1910!
1911!-----------------------------------------------------------------------
1912! Compute new wet/dry masks.
1913!-----------------------------------------------------------------------
1914!
1915 CALL wetdry_tile (ng, tile, &
1916 & lbi, ubi, lbj, ubj, &
1917 & imins, imaxs, jmins, jmaxs, &
1918# ifdef MASKING
1919 & pmask, rmask, umask, vmask, &
1920# endif
1921 & h, zeta(:,:,knew), &
1922# ifdef SOLVE3D
1923 & du_avg1, dv_avg1, &
1924 & rmask_wet_avg, &
1925# endif
1926 & pmask_wet, pmask_full, &
1927 & rmask_wet, rmask_full, &
1928 & umask_wet, umask_full, &
1929 & vmask_wet, vmask_full)
1930#endif
1931
1932#ifdef SOLVE3D
1933!
1934!-----------------------------------------------------------------------
1935! At the end of the last 2D time step replace the new free-surface
1936! zeta(:,:,knew) with it fast time-averaged value, Zt_avg1. Recall
1937! this is state variable is the one that communicates with the 3D
1938! kernel. Then, compute time-dependent depths.
1939!-----------------------------------------------------------------------
1940!
1941 IF (iif(ng).eq.nfast(ng)) THEN
1942 DO j=jstrr,jendr
1943 DO i=istrr,iendr
1944 zeta(i,j,knew)=zt_avg1(i,j)
1945 END DO
1946 END DO
1947 CALL set_depth (ng, tile, inlm)
1948 END IF
1949#endif
1950
1951#ifdef NESTING
1952# ifdef SOLVE3D
1953!
1954!-----------------------------------------------------------------------
1955! If nesting and after all fast time steps are completed, exchange
1956! halo information to time averaged fields.
1957!-----------------------------------------------------------------------
1958!
1959 IF (iif(ng).eq.nfast(ng)) THEN
1960!
1961! In nesting applications with refinement grids, we need to exchange
1962! the DU_avg2 and DV_avg2 fluxes boundary information for the case
1963! that a contact point is at a tile partition. Notice that in such
1964! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
1965!
1966 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1967 CALL exchange_r2d_tile (ng, tile, &
1968 & lbi, ubi, lbj, ubj, &
1969 & zt_avg1)
1970 CALL exchange_u2d_tile (ng, tile, &
1971 & lbi, ubi, lbj, ubj, &
1972 & du_avg1)
1973 CALL exchange_v2d_tile (ng, tile, &
1974 & lbi, ubi, lbj, ubj, &
1975 & dv_avg1)
1976 CALL exchange_u2d_tile (ng, tile, &
1977 & lbi, ubi, lbj, ubj, &
1978 & du_avg2)
1979 CALL exchange_v2d_tile (ng, tile, &
1980 & lbi, ubi, lbj, ubj, &
1981 & dv_avg2)
1982 END IF
1983
1984# ifdef DISTRIBUTE
1985!
1986 CALL mp_exchange2d (ng, tile, inlm, 3, &
1987 & lbi, ubi, lbj, ubj, &
1988 & nghostpoints, &
1989 & ewperiodic(ng), nsperiodic(ng), &
1990 & zt_avg1, du_avg1, dv_avg1)
1991 CALL mp_exchange2d (ng, tile, inlm, 2, &
1992 & lbi, ubi, lbj, ubj, &
1993 & nghostpoints, &
1994 & ewperiodic(ng), nsperiodic(ng), &
1995 & du_avg2, dv_avg2)
1996# endif
1997 END IF
1998# else
1999!
2000! In nesting applications with refinement grids, we need to exchange
2001! the DU_flux and DV_flux fluxes boundary information for the case
2002! that a contact point is at a tile partition. Notice that in such
2003! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
2004!
2005 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
2006 CALL exchange_u2d_tile (ng, tile, &
2007 & lbi, ubi, lbj, ubj, &
2008 & du_flux)
2009 CALL exchange_v2d_tile (ng, tile, &
2010 & lbi, ubi, lbj, ubj, &
2011 & dv_flux)
2012 END IF
2013
2014# ifdef DISTRIBUTE
2015!
2016 CALL mp_exchange2d (ng, tile, inlm, 2, &
2017 & lbi, ubi, lbj, ubj, &
2018 & nghostpoints, &
2019 & ewperiodic(ng), nsperiodic(ng), &
2020 & du_flux, dv_flux)
2021# endif
2022# endif
2023#endif
2024!
2025!-----------------------------------------------------------------------
2026! Exchange halo tile information.
2027!-----------------------------------------------------------------------
2028!
2029 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
2030 CALL exchange_r2d_tile (ng, tile, &
2031 & lbi, ubi, lbj, ubj, &
2032 & zeta(:,:,knew))
2033 CALL exchange_u2d_tile (ng, tile, &
2034 & lbi, ubi, lbj, ubj, &
2035 & ubar(:,:,knew))
2036 CALL exchange_v2d_tile (ng, tile, &
2037 & lbi, ubi, lbj, ubj, &
2038 & vbar(:,:,knew))
2039 END IF
2040
2041#ifdef DISTRIBUTE
2042!
2043 CALL mp_exchange2d (ng, tile, inlm, 3, &
2044 & lbi, ubi, lbj, ubj, &
2045 & nghostpoints, &
2046 & ewperiodic(ng), nsperiodic(ng), &
2047 & zeta(:,:,knew), &
2048 & ubar(:,:,knew), &
2049 & vbar(:,:,knew))
2050#endif
2051!
2052 RETURN
2053 END SUBROUTINE step2d_tile
2054
2055 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_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_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_parallel
Definition mod_parallel.F:2
mod_parallel::master
logical master
Definition mod_parallel.F:40
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::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::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
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::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::lwsrc
logical, dimension(:), allocatable lwsrc
Definition mod_scalars.F:543
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::m2yvis
integer m2yvis
Definition mod_scalars.F:170
mod_scalars::dtfast
real(dp), dimension(:), allocatable dtfast
Definition mod_scalars.F:270
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::ntstart
integer, dimension(:), allocatable ntstart
Definition mod_scalars.F:378
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::m2xvis
integer m2xvis
Definition mod_scalars.F:169
mod_scalars::iif
integer, dimension(:), allocatable iif
Definition mod_scalars.F:250
mod_sedbed
Definition sedbed_mod.h:1
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
set_depth_mod
Definition set_depth.F:2
set_depth_mod::set_depth
subroutine, public set_depth(ng, tile, model)
Definition set_depth.F:34
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
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