ROMS
Loading...
Searching...
No Matches
tl_step2d_FB.h
Go to the documentation of this file.
1#define DEBUG
2
3 MODULE tl_step2d_mod
4!
5!git $Id$
6!=======================================================================
7! !
8! Solves perturbation tangent linear model shallow-water primitive !
9! equations (barotropic mode) using the Generalized Forward-Backward !
10! 3rd-order Adams-Bashforth / 4th-order Adams-Moulton (FB AB3-AM4) !
11! time stepping algorithm (Shchepetkin and McWilliams, 2005); see !
12! section 2.3 starting with equation 2.49. In 3D applications, it !
13! perform fast-time averaging to interact with 3D momentum equations !
14! (baroclinic mode). !
15! !
16! Reference: !
17! !
18! Shchepetkin, A.F. and J.C. McWilliams, 2005: The regional oceanic !
19! modeling system (ROMS): a split-explicit, free-surface, !
20! topography-following-coordinate oceanic model, Ocean Modelling, !
21! 9, 347-404, doi:10.1016/j.ocemod.2004.08.002. !
22! !
23! Shchepetkin, A.F., and J.C. McWilliams, 2009: Computational kernel !
24! algorithms for fine-scale, multiprocess, longtime oceanic !
25! simulations, pp 121-183. In 'Handbook of Numerical Analysis: !
26! Computational Methods for the Atmosphere and Oceans', R.M. Teman !
27! and J.J. Tribbia, eds, Elsevier Science. !
28! !
29!=======================================================================
30!
31 USE mod_param
32 USE mod_parallel
33#ifdef SOLVE3D
34 USE mod_coupling
35#endif
36#ifdef DIAGNOSTICS_UV
37!! USE mod_diags
38#endif
39 USE mod_forces
40 USE mod_grid
41#if defined UV_VIS2 || defined UV_VIS4
42 USE mod_mixing
43#endif
44 USE mod_ocean
45 USE mod_scalars
46#if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET && \
47 defined solve3d
48 USE mod_sedbed
49#endif
50 USE mod_sources
51 USE mod_stepping
52!
53 USE exchange_2d_mod
54#ifdef DISTRIBUTE
55 USE mp_exchange_mod, ONLY : mp_exchange2d
56#endif
57 USE obc_volcons_mod, ONLY : obc_flux_tile, &
58 & set_duv_bc_tile
59 USE tl_obc_volcons_mod, ONLY : tl_obc_flux_tile, &
60 & tl_set_duv_bc_tile
61#ifdef SOLVE3D
62 USE tl_set_depth_mod, ONLY : tl_set_depth
63#endif
64 USE tl_u2dbc_mod, ONLY : tl_u2dbc_tile
65 USE tl_v2dbc_mod, ONLY : tl_v2dbc_tile
66 USE tl_zetabc_mod, ONLY : tl_zetabc_local
67#ifdef WET_DRY_NOT_YET
68 USE wetdry_mod, ONLY : wetdry_tile
69#endif
70!
71 implicit none
72!
73 PRIVATE
74 PUBLIC :: tl_step2d
75!
76 CONTAINS
77!
78!***********************************************************************
79 SUBROUTINE tl_step2d (ng, tile)
80!***********************************************************************
81!
82! Imported variable declarations.
83!
84 integer, intent(in) :: ng, tile
85!
86! Local variable declarations.
87!
88 character (len=*), parameter :: myfile = &
89 & __FILE__
90!
91#include "tile.h"
92!
93#ifdef PROFILE
94 CALL wclock_on (ng, itlm, 9, __line__, myfile)
95#endif
96 CALL tl_step2d_tile (ng, tile, &
97 & lbi, ubi, lbj, ubj, n(ng), &
98 & imins, imaxs, jmins, jmaxs, &
99 & krhs(ng), kstp(ng), knew(ng), &
100#ifdef SOLVE3D
101 & nstp(ng), nnew(ng), &
102#endif
103#ifdef MASKING
104 & grid(ng) % pmask, grid(ng) % rmask, &
105 & grid(ng) % umask, grid(ng) % vmask, &
106#endif
107#ifdef WET_DRY_NOT_YET
108 & grid(ng) % pmask_wet, grid(ng) % pmask_full, &
109 & grid(ng) % rmask_wet, grid(ng) % rmask_full, &
110 & grid(ng) % umask_wet, grid(ng) % umask_full, &
111 & grid(ng) % vmask_wet, grid(ng) % vmask_full, &
112# ifdef SOLVE3D
113 & grid(ng) % rmask_wet_avg, &
114# endif
115#endif
116#if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
117 & grid(ng) % fomn, &
118#endif
119 & grid(ng) % h, grid(ng) % tl_h, &
120 & grid(ng) % om_u, grid(ng) % om_v, &
121 & grid(ng) % on_u, grid(ng) % on_v, &
122 & grid(ng) % pm, grid(ng) % pn, &
123#if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
124 & grid(ng) % dndx, grid(ng) % dmde, &
125#endif
126 & grid(ng) % rdrag, &
127#if defined UV_QDRAG && !defined SOLVE3D
128 & grid(ng) % rdrag2, &
129#endif
130#if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
131 & grid(ng) % pmon_r, grid(ng) % pnom_r, &
132 & grid(ng) % pmon_p, grid(ng) % pnom_p, &
133 & grid(ng) % om_r, grid(ng) % on_r, &
134 & grid(ng) % om_p, grid(ng) % on_p, &
135# ifdef UV_VIS2
136 & mixing(ng) % visc2_p, &
137 & mixing(ng) % visc2_r, &
138# endif
139# ifdef UV_VIS4
140 & mixing(ng) % visc4_p, &
141 & mixing(ng) % visc4_r, &
142# endif
143#endif
144#if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
145 & sedbed(ng) % tl_bed_thick, &
146#endif
147#ifdef WEC_MELLOR
148 & mixing(ng) % tl_rustr2d, &
149 & mixing(ng) % tl_rvstr2d, &
150 & ocean(ng) % tl_rulag2d, &
151 & ocean(ng) % tl_rvlag2d, &
152 & ocean(ng) % ubar_stokes, &
153 & ocean(ng) % tl_ubar_stokes, &
154 & ocean(ng) % vbar_stokes, &
155 & ocean(ng) % tl_vbar_stokes, &
156#endif
157#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
158 & ocean(ng) % eq_tide, &
159 & ocean(ng) % tl_eq_tide, &
160#endif
161#ifndef SOLVE3D
162 & forces(ng) % sustr, forces(ng) % tl_sustr, &
163 & forces(ng) % svstr, forces(ng) % tl_svstr, &
164# ifdef ATM_PRESS
165 & forces(ng) % Pair, &
166# endif
167#else
168# ifdef VAR_RHO_2D
169 & coupling(ng) % rhoA, &
170 & coupling(ng) % tl_rhoA, &
171 & coupling(ng) % rhoS, &
172 & coupling(ng) % tl_rhoS, &
173# endif
174 & coupling(ng) % tl_DU_avg1, &
175 & coupling(ng) % tl_DU_avg2, &
176 & coupling(ng) % tl_DV_avg1, &
177 & coupling(ng) % tl_DV_avg2, &
178 & coupling(ng) % tl_Zt_avg1, &
179 & coupling(ng) % rufrc, &
180 & coupling(ng) % tl_rufrc, &
181 & coupling(ng) % rvfrc, &
182 & coupling(ng) % tl_rvfrc, &
183 & coupling(ng) % tl_rufrc_bak, &
184 & coupling(ng) % tl_rvfrc_bak, &
185#endif
186#if defined NESTING && !defined SOLVE3D
187 & ocean(ng) % tl_DU_flux, &
188 & ocean(ng) % tl_DV_flux, &
189#endif
190#ifdef DIAGNOSTICS_UV
191!! & DIAGS(ng) % DiaU2wrk, DIAGS(ng) % DiaV2wrk, &
192!! & DIAGS(ng) % DiaRUbar, DIAGS(ng) % DiaRVbar, &
193# ifdef SOLVE3D
194!! & DIAGS(ng) % DiaU2int, DIAGS(ng) % DiaV2int, &
195!! & DIAGS(ng) % DiaRUfrc, DIAGS(ng) % DiaRVfrc, &
196# endif
197#endif
198 & ocean(ng) % ubar, ocean(ng) % tl_ubar, &
199 & ocean(ng) % vbar, ocean(ng) % tl_vbar, &
200 & ocean(ng) % zeta, ocean(ng) % tl_zeta)
201#ifdef PROFILE
202 CALL wclock_off (ng, itlm, 9, __line__, myfile)
203#endif
204!
205 RETURN
206 END SUBROUTINE tl_step2d
207!
208!***********************************************************************
209 SUBROUTINE tl_step2d_tile (ng, tile, &
210 & LBi, UBi, LBj, UBj, UBk, &
211 & IminS, ImaxS, JminS, JmaxS, &
212 & krhs, kstp, knew, &
213#ifdef SOLVE3D
214 & nstp, nnew, &
215#endif
216#ifdef MASKING
217 & pmask, rmask, umask, vmask, &
218#endif
219#ifdef WET_DRY_NOT_YET
220 & pmask_wet, pmask_full, &
221 & rmask_wet, rmask_full, &
222 & umask_wet, umask_full, &
223 & vmask_wet, vmask_full, &
224# ifdef SOLVE3D
225 & rmask_wet_avg, &
226# endif
227#endif
228#if (defined UV_COR && !defined SOLVE3D) || defined step2d_coriolis
229 & fomn, &
230#endif
231 & h, tl_h, &
232 & om_u, om_v, on_u, on_v, pm, pn, &
233#if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
234 & dndx, dmde, &
235#endif
236 & rdrag, &
237#if defined UV_QDRAG && !defined SOLVE3D
238 & rdrag2, &
239#endif
240#if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
241 & pmon_r, pnom_r, pmon_p, pnom_p, &
242 & om_r, on_r, om_p, on_p, &
243# ifdef UV_VIS2
244 & visc2_p, visc2_r, &
245# endif
246# ifdef UV_VIS4
247 & visc4_p, visc4_r, &
248# endif
249#endif
250#if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
251 & tl_bed_thick, &
252#endif
253#ifdef WEC_MELLOR
254 & tl_rustr2d, tl_rvstr2d, &
255 & tl_rulag2d, tl_rvlag2d, &
256 & ubar_stokes, tl_ubar_stokes, &
257 & vbar_stokes, tl_vbar_stokes, &
258#endif
259#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
260 & eq_tide, tl_eq_tide, &
261#endif
262#ifndef SOLVE3D
263 & sustr, tl_sustr, &
264 & svstr, tl_svstr, &
265# ifdef ATM_PRESS
266 & pair, &
267# endif
268#else
269# ifdef VAR_RHO_2D
270 & rhoa, tl_rhoa, &
271 & rhos, tl_rhos, &
272# endif
273 & tl_du_avg1, tl_du_avg2, &
274 & tl_dv_avg1, tl_dv_avg2, &
275 & tl_zt_avg1, &
276 & rufrc, tl_rufrc, &
277 & rvfrc, tl_rvfrc, &
278 & tl_rufrc_bak, tl_rvfrc_bak, &
279#endif
280#if defined NESTING && !defined SOLVE3D
281 & tl_du_flux, tl_dv_flux, &
282#endif
283#ifdef DIAGNOSTICS_UV
284!! & DiaU2wrk, DiaV2wrk, &
285!! & DiaRUbar, DiaRVbar, &
286# ifdef SOLVE3D
287!! & DiaU2int, DiaV2int, &
288!! & DiaRUfrc, DiaRVfrc, &
289# endif
290#endif
291 & ubar, tl_ubar, &
292 & vbar, tl_vbar, &
293 & zeta, tl_zeta)
294!***********************************************************************
295!
296! Imported variable declarations.
297!
298 integer, intent(in ) :: ng, tile
299 integer, intent(in ) :: LBi, UBi, LBj, UBj, UBk
300 integer, intent(in ) :: IminS, ImaxS, JminS, JmaxS
301 integer, intent(in ) :: krhs, kstp, knew
302#ifdef SOLVE3D
303 integer, intent(in ) :: nstp, nnew
304#endif
305!
306#ifdef ASSUMED_SHAPE
307# ifdef MASKING
308 real(r8), intent(in ) :: pmask(LBi:,LBj:)
309 real(r8), intent(in ) :: rmask(LBi:,LBj:)
310 real(r8), intent(in ) :: umask(LBi:,LBj:)
311 real(r8), intent(in ) :: vmask(LBi:,LBj:)
312# endif
313# if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
314 real(r8), intent(in ) :: fomn(LBi:,LBj:)
315# endif
316 real(r8), intent(in ) :: h(LBi:,LBj:)
317 real(r8), intent(in ) :: om_u(LBi:,LBj:)
318 real(r8), intent(in ) :: om_v(LBi:,LBj:)
319 real(r8), intent(in ) :: on_u(LBi:,LBj:)
320 real(r8), intent(in ) :: on_v(LBi:,LBj:)
321 real(r8), intent(in ) :: pm(LBi:,LBj:)
322 real(r8), intent(in ) :: pn(LBi:,LBj:)
323# if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
324 real(r8), intent(in ) :: dndx(LBi:,LBj:)
325 real(r8), intent(in ) :: dmde(LBi:,LBj:)
326# endif
327 real(r8), intent(in ) :: rdrag(LBi:,LBj:)
328# if defined UV_QDRAG && !defined SOLVE3D
329 real(r8), intent(in ) :: rdrag2(LBi:,LBj:)
330# endif
331# if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
332 real(r8), intent(in ) :: pmon_r(LBi:,LBj:)
333 real(r8), intent(in ) :: pnom_r(LBi:,LBj:)
334 real(r8), intent(in ) :: pmon_p(LBi:,LBj:)
335 real(r8), intent(in ) :: pnom_p(LBi:,LBj:)
336 real(r8), intent(in ) :: om_r(LBi:,LBj:)
337 real(r8), intent(in ) :: on_r(LBi:,LBj:)
338 real(r8), intent(in ) :: om_p(LBi:,LBj:)
339 real(r8), intent(in ) :: on_p(LBi:,LBj:)
340# ifdef UV_VIS2
341 real(r8), intent(in ) :: visc2_p(LBi:,LBj:)
342 real(r8), intent(in ) :: visc2_r(LBi:,LBj:)
343# endif
344# ifdef UV_VIS4
345 real(r8), intent(in ) :: visc4_p(LBi:,LBj:)
346 real(r8), intent(in ) :: visc4_r(LBi:,LBj:)
347# endif
348# endif
349# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
350 real(r8), intent(in ) :: tl_bed_thick(LBi:,LBj:,:)
351# endif
352# ifdef WEC_MELLOR
353 real(r8), intent(in ) :: ubar_stokes(LBi:,LBj:)
354 real(r8), intent(in ) :: vbar_stokes(LBi:,LBj:)
355# endif
356# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
357 real(r8), intent(in ) :: eq_tide(LBi:,LBj:)
358 real(r8), intent(in ) :: tl_eq_tide(LBi:,LBj:)
359# endif
360 real(r8), intent(in ) :: ubar(LBi:,LBj:,:)
361 real(r8), intent(in ) :: vbar(LBi:,LBj:,:)
362 real(r8), intent(in ) :: zeta(LBi:,LBj:,:)
363# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
364 real(r8), intent(inout) :: tl_h(LBi:,LBj:)
365# else
366 real(r8), intent(in ) :: tl_h(LBi:,LBj:)
367# endif
368# ifndef SOLVE3D
369 real(r8), intent(in ) :: tl_sustr(LBi:,LBj:)
370 real(r8), intent(in ) :: tl_svstr(LBi:,LBj:)
371# ifdef ATM_PRESS
372 real(r8), intent(in ) :: Pair(LBi:,LBj:)
373# endif
374# else
375# ifdef VAR_RHO_2D
376 real(r8), intent(in ) :: rhoA(LBi:,LBj:)
377 real(r8), intent(in ) :: rhoS(LBi:,LBj:)
378 real(r8), intent(in ) :: tl_rhoA(LBi:,LBj:)
379 real(r8), intent(in ) :: tl_rhoS(LBi:,LBj:)
380# endif
381 real(r8), intent(in ) :: rufrc(LBi:,LBj:)
382 real(r8), intent(in ) :: rvfrc(LBi:,LBj:)
383
384 real(r8), intent(inout) :: tl_DU_avg1(LBi:,LBj:)
385 real(r8), intent(inout) :: tl_DU_avg2(LBi:,LBj:)
386 real(r8), intent(inout) :: tl_DV_avg1(LBi:,LBj:)
387 real(r8), intent(inout) :: tl_DV_avg2(LBi:,LBj:)
388 real(r8), intent(inout) :: tl_Zt_avg1(LBi:,LBj:)
389 real(r8), intent(inout) :: tl_rufrc(LBi:,LBj:)
390 real(r8), intent(inout) :: tl_rvfrc(LBi:,LBj:)
391 real(r8), intent(inout) :: tl_rufrc_bak(LBi:,LBj:,:)
392 real(r8), intent(inout) :: tl_rvfrc_bak(LBi:,LBj:,:)
393# endif
394# ifdef WEC_MELLOR
395 real(r8), intent(inout) :: tl_rustr2d(LBi:,LBj:)
396 real(r8), intent(inout) :: tl_rvstr2d(LBi:,LBj:)
397 real(r8), intent(inout) :: tl_rulag2d(LBi:,LBj:)
398 real(r8), intent(inout) :: tl_rvlag2d(LBi:,LBj:)
399 real(r8), intent(inout) :: tl_ubar_stokes(LBi:,LBj:)
400 real(r8), intent(inout) :: tl_vbar_stokes(LBi:,LBj:)
401# endif
402# ifdef WET_DRY_NOT_YET
403 real(r8), intent(inout) :: pmask_full(LBi:,LBj:)
404 real(r8), intent(inout) :: rmask_full(LBi:,LBj:)
405 real(r8), intent(inout) :: umask_full(LBi:,LBj:)
406 real(r8), intent(inout) :: vmask_full(LBi:,LBj:)
407
408 real(r8), intent(inout) :: pmask_wet(LBi:,LBj:)
409 real(r8), intent(inout) :: rmask_wet(LBi:,LBj:)
410 real(r8), intent(inout) :: umask_wet(LBi:,LBj:)
411 real(r8), intent(inout) :: vmask_wet(LBi:,LBj:)
412# ifdef SOLVE3D
413 real(r8), intent(inout) :: rmask_wet_avg(LBi:,LBj:)
414# endif
415# endif
416# ifdef DIAGNOSTICS_UV
417!! real(r8), intent(inout) :: DiaU2wrk(LBi:,LBj:,:)
418!! real(r8), intent(inout) :: DiaV2wrk(LBi:,LBj:,:)
419!! real(r8), intent(inout) :: DiaRUbar(LBi:,LBj:,:,:)
420!! real(r8), intent(inout) :: DiaRVbar(LBi:,LBj:,:,:)
421# ifdef SOLVE3D
422!! real(r8), intent(inout) :: DiaU2int(LBi:,LBj:,:)
423!! real(r8), intent(inout) :: DiaV2int(LBi:,LBj:,:)
424!! real(r8), intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
425!! real(r8), intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
426# endif
427# endif
428 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
429 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
430 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
431# if defined NESTING && !defined SOLVE3D
432 real(r8), intent(out ) :: tl_DU_flux(LBi:,LBj:)
433 real(r8), intent(out ) :: tl_DV_flux(LBi:,LBj:)
434# endif
435
436#else
437
438# ifdef MASKING
439 real(r8), intent(in ) :: pmask(LBi:UBi,LBj:UBj)
440 real(r8), intent(in ) :: rmask(LBi:UBi,LBj:UBj)
441 real(r8), intent(in ) :: umask(LBi:UBi,LBj:UBj)
442 real(r8), intent(in ) :: vmask(LBi:UBi,LBj:UBj)
443# endif
444# if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
445 real(r8), intent(in ) :: fomn(LBi:UBi,LBj:UBj)
446# endif
447 real(r8), intent(in ) :: h(LBi:UBi,LBj:UBj)
448 real(r8), intent(in ) :: om_u(LBi:UBi,LBj:UBj)
449 real(r8), intent(in ) :: om_v(LBi:UBi,LBj:UBj)
450 real(r8), intent(in ) :: on_u(LBi:UBi,LBj:UBj)
451 real(r8), intent(in ) :: on_v(LBi:UBi,LBj:UBj)
452 real(r8), intent(in ) :: pm(LBi:UBi,LBj:UBj)
453 real(r8), intent(in ) :: pn(LBi:UBi,LBj:UBj)
454# if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
455 real(r8), intent(in ) :: dndx(LBi:UBi,LBj:UBj)
456 real(r8), intent(in ) :: dmde(LBi:UBi,LBj:UBj)
457# endif
458 real(r8), intent(in ) :: rdrag(LBi:UBi,LBj:UBj)
459# if defined UV_QDRAG && !defined SOLVE3D
460 real(r8), intent(in ) :: rdrag2(LBi:UBi,LBj:UBj)
461# endif
462# if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
463 real(r8), intent(in ) :: pmon_r(LBi:UBi,LBj:UBj)
464 real(r8), intent(in ) :: pnom_r(LBi:UBi,LBj:UBj)
465 real(r8), intent(in ) :: pmon_p(LBi:UBi,LBj:UBj)
466 real(r8), intent(in ) :: pnom_p(LBi:UBi,LBj:UBj)
467 real(r8), intent(in ) :: om_r(LBi:UBi,LBj:UBj)
468 real(r8), intent(in ) :: on_r(LBi:UBi,LBj:UBj)
469 real(r8), intent(in ) :: om_p(LBi:UBi,LBj:UBj)
470 real(r8), intent(in ) :: on_p(LBi:UBi,LBj:UBj)
471# ifdef UV_VIS2
472 real(r8), intent(in ) :: visc2_p(LBi:UBi,LBj:UBj)
473 real(r8), intent(in ) :: visc2_r(LBi:UBi,LBj:UBj)
474# endif
475# ifdef UV_VIS4
476 real(r8), intent(in ) :: visc4_p(LBi:UBi,LBj:UBj)
477 real(r8), intent(in ) :: visc4_r(LBi:UBi,LBj:UBj)
478# endif
479# endif
480# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
481 real(r8), intent(in ) :: tl_bed_thick(LBi:UBi,LBj:UBj,3)
482# endif
483# ifdef WEC_MELLOR
484 real(r8), intent(in ) :: ubar_stokes(LBi:UBi,LBj:UBj)
485 real(r8), intent(in ) :: vbar_stokes(LBi:UBi,LBj:UBj)
486# endif
487# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
488 real(r8), intent(in ) :: eq_tide(LBi:UBi,LBj:UBj)
489 real(r8), intent(in ) :: tl_eq_tide(LBi:UBi,LBj:UBj)
490# endif
491 real(r8), intent(in ) :: ubar(LBi:UBi,LBj:UBj,:)
492 real(r8), intent(in ) :: vbar(LBi:UBi,LBj:UBj,:)
493 real(r8), intent(in ) :: zeta(LBi:UBi,LBj:UBj,:)
494# if defined SEDIMENT_NOT_YET && defined SED_MORPH_NOT_YET
495 real(r8), intent(inout) :: tl_h(LBi:UBi,LBj:UBj)
496# else
497 real(r8), intent(in ) :: tl_h(LBi:UBi,LBj:UBj)
498# endif
499# ifndef SOLVE3D
500 real(r8), intent(in ) :: tl_sustr(LBi:UBi,LBj:UBj)
501 real(r8), intent(in ) :: tl_svstr(LBi:UBi,LBj:UBj)
502# ifdef ATM_PRESS
503 real(r8), intent(in ) :: Pair(LBi:UBi,LBj:UBj)
504# endif
505# else
506# ifdef VAR_RHO_2D
507 real(r8), intent(in ) :: rhoA(LBi:UBi,LBj:UBj)
508 real(r8), intent(in ) :: rhoS(LBi:UBi,LBj:UBj)
509 real(r8), intent(in ) :: tl_rhoA(LBi:UBi,LBj:UBj)
510 real(r8), intent(in ) :: tl_rhoS(LBi:UBi,LBj:UBj)
511# endif
512 real(r8), intent(in ) :: rufrc(LBi:UBi,LBj:UBj)
513 real(r8), intent(in ) :: rvfrc(LBi:UBi,LBj:UBj)
514
515 real(r8), intent(inout) :: tl_DU_avg1(LBi:UBi,LBj:UBj)
516 real(r8), intent(inout) :: tl_DU_avg2(LBi:UBi,LBj:UBj)
517 real(r8), intent(inout) :: tl_DV_avg1(LBi:UBi,LBj:UBj)
518 real(r8), intent(inout) :: tl_DV_avg2(LBi:UBi,LBj:UBj)
519 real(r8), intent(inout) :: tl_Zt_avg1(LBi:UBi,LBj:UBj)
520 real(r8), intent(inout) :: tl_rufrc(LBi:UBi,LBj:UBj)
521 real(r8), intent(inout) :: tl_rvfrc(LBi:UBi,LBj:UBj)
522 real(r8), intent(inout) :: tl_rufrc_bak(LBi:UBi,LBj:UBj,2)
523 real(r8), intent(inout) :: tl_rvfrc_bak(LBi:UBi,LBj:UBj,2)
524# endif
525# ifdef WEC_MELLOR
526 real(r8), intent(inout) :: tl_rustr2d(LBi:UBi,LBj:UBj)
527 real(r8), intent(inout) :: tl_rvstr2d(LBi:UBi,LBj:UBj)
528 real(r8), intent(inout) :: tl_rulag2d(LBi:UBi,LBj:UBj)
529 real(r8), intent(inout) :: tl_rvlag2d(LBi:UBi,LBj:UBj)
530 real(r8), intent(inout) :: tl_ubar_stokes(LBi:UBi,LBj:UBj)
531 real(r8), intent(inout) :: tl_vbar_stokes(LBi:UBi,LBj:UBj)
532# endif
533# ifdef WET_DRY_NOT_YET
534 real(r8), intent(inout) :: pmask_full(LBi:UBi,LBj:UBj)
535 real(r8), intent(inout) :: rmask_full(LBi:UBi,LBj:UBj)
536 real(r8), intent(inout) :: umask_full(LBi:UBi,LBj:UBj)
537 real(r8), intent(inout) :: vmask_full(LBi:UBi,LBj:UBj)
538
539 real(r8), intent(inout) :: pmask_wet(LBi:UBi,LBj:UBj)
540 real(r8), intent(inout) :: rmask_wet(LBi:UBi,LBj:UBj)
541 real(r8), intent(inout) :: umask_wet(LBi:UBi,LBj:UBj)
542 real(r8), intent(inout) :: vmask_wet(LBi:UBi,LBj:UBj)
543# ifdef SOLVE3D
544 real(r8), intent(inout) :: rmask_wet_avg(LBi:UBi,LBj:UBj)
545# endif
546# endif
547# ifdef DIAGNOSTICS_UV
548!! real(r8), intent(inout) :: DiaU2wrk(LBi:UBi,LBj:UBj,NDM2d)
549!! real(r8), intent(inout) :: DiaV2wrk(LBi:UBi,LBj:UBj,NDM2d)
550!! real(r8), intent(inout) :: DiaRUbar(LBi:UBi,LBj:UBj,2,NDM2d-1)
551!! real(r8), intent(inout) :: DiaRVbar(LBi:UBi,LBj:UBj,2,NDM2d-1)
552# ifdef SOLVE3D
553!! real(r8), intent(inout) :: DiaU2int(LBi:UBi,LBj:UBj,NDM2d)
554!! real(r8), intent(inout) :: DiaV2int(LBi:UBi,LBj:UBj,NDM2d)
555!! real(r8), intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
556!! real(r8), intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
557# endif
558# endif
559 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
560 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
561 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
562# if defined NESTING && !defined SOLVE3D
563 real(r8), intent(out ) :: tl_DU_flux(LBi:UBi,LBj:UBj)
564 real(r8), intent(out ) :: tl_DV_flux(LBi:UBi,LBj:UBj)
565# endif
566#endif
567!
568! Local variable declarations.
569!
570 integer :: i, is, j
571 integer :: kbak, kold
572#ifdef DIAGNOSTICS_UV
573 integer :: idiag
574#endif
575!
576 real(r8) :: bkw0, bkw1, bkw2, bkw_new
577 real(r8) :: fwd0, fwd1, fwd2
578#ifdef SOLVE3D
579 real(r8) :: cfwd0, cfwd1, cfwd2
580#endif
581 real(r8) :: cff, cff1, cff2, cff3, cff4
582#ifdef WET_DRY_NOT_YET
583 real(r8) :: cff5, cff6, cff7
584#endif
585 real(r8) :: fac, fac1, fac2
586 real(r8) :: tl_cff, tl_cff1, tl_cff2, tl_cff3, tl_cff4
587#ifdef WET_DRY_NOT_YET
588 real(r8) :: tl_cff5, tl_cff6, tl_cff7
589#endif
590 real(r8) :: tl_fac, tl_fac1, tl_fac2
591#ifdef DEBUG
592 real(r8), parameter :: IniVal = 0.0_r8
593#endif
594!
595#if defined UV_C4ADVECTION && !defined SOLVE3D
596 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dgrad
597#endif
598 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dnew
599 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dnew_rd
600 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs
601#if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
602 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs_p
603#endif
604 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dstp
605 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUon
606 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVom
607#ifdef WEC_MELLOR
608 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUSon
609 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVSom
610#endif
611#if defined STEP2D_CORIOLIS || !defined SOLVE3D
612 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
613 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
614#endif
615#if !defined SOLVE3D
616 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
617 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
618#endif
619#if defined UV_C4ADVECTION && !defined SOLVE3D
620 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
621#endif
622 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzeta2
623#if defined VAR_RHO_2D && defined SOLVE3D
624 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzetaSA
625#endif
626 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rubar
627 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rvbar
628 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzeta
629 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: urhs
630 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: vrhs
631 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zeta_new
632 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zwrk
633#ifdef WET_DRY_NOT_YET
634 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: wetdry
635#endif
636#ifdef DIAGNOSTICS_UV
637!! real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Uwrk
638!! real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vwrk
639!! real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaU2rhs
640!! real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaV2rhs
641#endif
642!
643#if defined UV_C4ADVECTION && !defined SOLVE3D
644 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Dgrad
645#endif
646 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Dnew
647 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Dnew_rd
648 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Drhs
649#if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
650 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Drhs_p
651#endif
652 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_Dstp
653 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_DUon
654 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_DVom
655#ifdef WEC_MELLOR
656 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_DUSon
657 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_DVSom
658#endif
659#if defined STEP2D_CORIOLIS || !defined SOLVE3D
660 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_UFx
661 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_VFe
662#endif
663#if !defined SOLVE3D
664 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_UFe
665 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_VFx
666#endif
667#if defined UV_C4ADVECTION && !defined SOLVE3D
668 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_grad
669#endif
670 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_rzeta2
671#if defined VAR_RHO_2D && defined SOLVE3D
672 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_rzetaSA
673#endif
674 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_rzeta
675 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_rubar
676 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_rvbar
677 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_urhs
678 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_vrhs
679 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_zwrk
680!
681 real(r8), allocatable :: tl_zeta_new(:,:)
682
683#include "set_bounds.h"
684
685#ifdef DEBUG
686!
687!-----------------------------------------------------------------------
688! Initialize private arrays for debugging.
689!-----------------------------------------------------------------------
690!
691# if defined UV_C4ADVECTION && !defined SOLVE3D
692 dgrad=inival
693# endif
694 dnew=inival
695 dnew_rd=inival
696 drhs=inival
697# if (defined UV_VIS2 || defined UV_VIS4) && !defined SOLVE3D
698 drhs_p=inival
699# endif
700 dstp=inival
701 duon=inival
702 dvom=inival
703# if defined STEP2D_CORIOLIS || !defined SOLVE3D
704 ufx=inival
705 vfe=inival
706# endif
707# if !defined SOLVE3D
708 ufe=inival
709 vfx=inival
710# endif
711# if defined UV_C4ADVECTION && !defined SOLVE3D
712 grad=inival
713# endif
714 rzeta2=inival
715# if defined VAR_RHO_2D && defined SOLVE3D
716 rzetasa=inival
717# endif
718 rzeta=inival
719 rubar=inival
720 rvbar=inival
721 urhs=inival
722 vrhs=inival
723 zwrk=inival
724# ifdef WET_DRY_NOT_YET
725 wetdry=inival
726# endif
727# ifdef DIAGNOSTICS_UV
728 uwrk=inival
729 vwrk=inival
730 diau2rhs=inival
731 diav2rhs=inival
732# endif
733#endif
734!
735!-----------------------------------------------------------------------
736! Set coefficients for AB3-AM4 forward-backward algorithm.
737!-----------------------------------------------------------------------
738!
739! Because the Forward Euler step is used to update "zeta" during the
740! first barotropic step, the pressure-gradient term in the momentum
741! equation must be computed via the Backward step to keep it
742! numerically stable. However, this interferes with the computation
743! of forcing terms "rufrc" and "rvfrc" because the free surface in
744! pressure gradient computation in 3D is exactly at the time
745! corresponding to baroclinic step "nstp" (rather than ahead by one
746! barotropic step after it updated by a normal forward-backward step).
747! To resolve this conflict, the pressure gradient term is computed in
748! two stages during the first barotropic step. It uses zeta(:,:,kstp)
749! at first to ensure exact consistency with 3D model. Then, after
750! vertical integrals of 3D right-hand-side "rufrc" and "rvfrc" are
751! converted into forcing terms, add correction based on the difference
752! zeta_new(:,:)-zeta(:,:,kstp) to "rubar" and "rvbar" to make them
753! consistent with the Backward step for pressure gradient.
754! For pressure gradient terms, search for the label PGF_FB_CORRECTION
755! below.
756!
757 IF (first_2d_step) THEN ! Meaning of time indices
758 kbak=kstp !------------------------
759 kold=kstp ! m-2 m-1 m m+1
760 fwd0=1.0_r8 ! kold kbak kstp knew
761 fwd1=0.0_r8 ! fwd2 fwd1 fwd0
762 fwd2=0.0_r8 ! bkw2 bkw1 bkw0 bkw_new
763#ifdef SOLVE3D
764 bkw_new=0.0_r8
765 bkw0=1.0_r8
766#else
767 bkw_new=1.0_r8
768 bkw0=0.0_r8
769#endif
770 bkw1=0.0_r8
771 bkw2=0.0_r8
772 ELSE IF (first_2d_step+1) THEN
773 kbak=kstp-1
774 IF (kbak.lt.1) kbak=4
775 kold=kbak
776 fwd0=1.0_r8 ! Logically AB2-AM3 forward-
777 fwd1=0.0_r8 ! backward scheme with maximum
778 fwd2=0.0_r8 ! stability coefficients while
779 bkw_new=1.0833333333333_r8 ! maintaining third-order
780 bkw0=-0.1666666666666_r8 ! accuracy, alpha_max=1.73
781 bkw1= 0.0833333333333_r8
782 bkw2= 0.0_r8
783 ELSE
784 kbak=kstp-1
785 IF (kbak.lt.1) kbak=4
786 kold=kbak-1
787 IF (kold.lt.1) kold=4
788 fwd0=1.781105_r8
789 fwd1=-1.06221_r8
790 fwd2=0.281105_r8
791 bkw_new=0.614_r8
792 bkw0=0.285_r8
793 bkw1=0.0880_r8
794 bkw2=0.013_r8
795 END IF
796
797#ifdef DEBUG
798!
799 IF (master) THEN
800 WRITE (20,10) iic(ng), iif(ng), kold, kbak, kstp, knew
801 10 FORMAT (' iic = ',i5.5,' iif = ',i3.3, &
802 & ' kold = ',i1,' kbak = ',i1,' kstp = ',i1,' knew = ',i1)
803 END IF
804#endif
805!
806!-----------------------------------------------------------------------
807! Preliminary steps.
808!-----------------------------------------------------------------------
809!
810! Compute total depth of water column and vertically integrated fluxes
811! needed for computing horizontal divergence to advance free surface
812! and for advection terms for the barotropic momentum equations.
813! Set BASIC STATE for 'Drhs', 'urhs', and 'vrhs'.
814!
815#if defined DISTRIBUTE && !defined NESTING
816# define IR_RANGE IstrUm2-1,Iendp2
817# define JR_RANGE JstrVm2-1,Jendp2
818# define IU_RANGE IstrUm1-1,Iendp2
819# define JU_RANGE Jstrm1-1,Jendp2
820# define IV_RANGE Istrm1-1,Iendp2
821# define JV_RANGE JstrVm1-1,Jendp2
822#else
823# define IR_RANGE IstrUm2-1,Iendp2
824# define JR_RANGE JstrVm2-1,Jendp2
825# define IU_RANGE IstrUm2,Iendp2
826# define JU_RANGE JstrVm2-1,Jendp2
827# define IV_RANGE IstrUm2-1,Iendp2
828# define JV_RANGE JstrVm2,Jendp2
829#endif
830
831 DO j=jr_range
832 DO i=ir_range
833!^ Drhs(i,j)=h(i,j)+fwd0*zeta(i,j,kstp)+ &
834!^ & fwd1*zeta(i,j,kbak)+ &
835!^ & fwd2*zeta(i,j,kold)
836!^
837 drhs(i,j)=h(i,j)+zeta(i,j,kstp) ! using background value
838 tl_drhs(i,j)=tl_h(i,j)+fwd0*tl_zeta(i,j,kstp)+ &
839 & fwd1*tl_zeta(i,j,kbak)+ &
840 & fwd2*tl_zeta(i,j,kold)
841 END DO
842 END DO
843!
844 DO j=ju_range
845 DO i=iu_range
846 cff=0.5_r8*on_u(i,j)
847 cff1=cff*(drhs(i,j)+drhs(i-1,j))
848 tl_cff1=cff*(tl_drhs(i,j)+tl_drhs(i-1,j))
849!^ urhs(i,j)=fwd0*ubar(i,j,kstp)+ &
850!^ & fwd1*ubar(i,j,kbak)+ &
851!^ & fwd2*ubar(i,j,kold)
852!^
853 urhs(i,j)=ubar(i,j,kstp) ! using background value
854 tl_urhs(i,j)=fwd0*tl_ubar(i,j,kstp)+ &
855 & fwd1*tl_ubar(i,j,kbak)+ &
856 & fwd2*tl_ubar(i,j,kold)
857 duon(i,j)=urhs(i,j)*cff1
858 tl_duon(i,j)=tl_urhs(i,j)*cff1+ &
859 & urhs(i,j)*tl_cff1
860 END DO
861 END DO
862!
863 DO j=jv_range
864 DO i=iv_range
865 cff=0.5_r8*om_v(i,j)
866 cff1=cff*(drhs(i,j)+drhs(i,j-1))
867 tl_cff1=cff*(tl_drhs(i,j)+tl_drhs(i,j-1))
868!^ vrhs(i,j)=fwd0*vbar(i,j,kstp)+ &
869!^ & fwd1*vbar(i,j,kbak)+ &
870!^ & fwd2*vbar(i,j,kold)
871!^
872 vrhs(i,j)=vbar(i,j,kstp) ! using background value
873 tl_vrhs(i,j)=fwd0*tl_vbar(i,j,kstp)+ &
874 & fwd1*tl_vbar(i,j,kbak)+ &
875 & fwd2*tl_vbar(i,j,kold)
876 dvom(i,j)=vrhs(i,j)*cff1
877 tl_dvom(i,j)=tl_vrhs(i,j)*cff1+ &
878 & vrhs(i,j)*tl_cff1
879 END DO
880 END DO
881
882#undef IU_RANGE
883#undef IV_RANGE
884#undef JU_RANGE
885#undef JV_RANGE
886
887#if defined DISTRIBUTE && \
888 defined uv_adv && defined uv_c4advection && !defined SOLVE3D
889!
890! In distributed-memory, the I- and J-ranges are different and a
891! special exchange is done here to avoid having three ghost points
892! for high-order numerical stencils. Notice that a private array is
893! passed below to the exchange routine. It also applies periodic
894! boundary conditions, if appropriate and no partitions in I- or
895! J-directions.
896!
897 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
898 CALL exchange_u2d_tile (ng, tile, &
899 & imins, imaxs, jmins, jmaxs, &
900 & duon)
901 CALL exchange_u2d_tile (ng, tile, &
902 & imins, imaxs, jmins, jmaxs, &
903 & tl_duon)
904 CALL exchange_v2d_tile (ng, tile, &
905 & imins, imaxs, jmins, jmaxs, &
906 & dvom)
907 CALL exchange_v2d_tile (ng, tile, &
908 & imins, imaxs, jmins, jmaxs, &
909 & tl_dvom)
910 END IF
911 CALL mp_exchange2d (ng, tile, itlm, 4, &
912 & imins, imaxs, jmins, jmaxs, &
913 & nghostpoints, &
914 & ewperiodic(ng), nsperiodic(ng), &
915 & duon, dvom, &
916 & tl_duon, tl_dvom)
917#endif
918!
919! Compute integral mass flux across open boundaries and adjust
920! for volume conservation. Compute BASIC STATE value.
921! HGA: Need to resolve 'krhs' index here.
922!
923 IF (any(volcons(:,ng))) THEN
924 CALL obc_flux_tile (ng, tile, &
925 & lbi, ubi, lbj, ubj, &
926 & imins, imaxs, jmins, jmaxs, &
927 & knew, &
928#ifdef MASKING
929 & umask, vmask, &
930#endif
931 & h, om_v, on_u, &
932 & ubar, vbar, zeta)
933!
934 CALL tl_obc_flux_tile (ng, tile, &
935 & lbi, ubi, lbj, ubj, &
936 & imins, imaxs, jmins, jmaxs, &
937 & knew, &
938# ifdef MASKING
939 & umask, vmask, &
940# endif
941 & h, tl_h, om_v, on_u, &
942 & ubar, vbar, zeta, &
943 & tl_ubar, tl_vbar, tl_zeta)
944!
945! Set vertically integrated mass fluxes DUon and DVom along the open
946! boundaries in such a way that the integral volume is conserved.
947!
948 CALL set_duv_bc_tile (ng, tile, &
949 & lbi, ubi, lbj, ubj, &
950 & imins, imaxs, jmins, jmaxs, &
951 & krhs, &
952#ifdef MASKING
953 & umask, vmask, &
954#endif
955 & om_v, on_u, &
956 & ubar, vbar, &
957 & drhs, duon, dvom)
958 CALL tl_set_duv_bc_tile (ng, tile, &
959 & lbi, ubi, lbj, ubj, &
960 & imins, imaxs, jmins, jmaxs, &
961 & krhs, &
962#ifdef MASKING
963 & umask, vmask, &
964#endif
965 & om_v, on_u, &
966 & ubar, vbar, &
967 & tl_ubar, tl_vbar, &
968 & drhs, duon, dvom, &
969 & tl_drhs, tl_duon, tl_dvom)
970 END IF
971!
972!-----------------------------------------------------------------------
973! Advance free-surface.
974!-----------------------------------------------------------------------
975!
976! Notice that the new local free-surface is allocated so it can be
977! passed as an argumment to "zetabc_local". An automatic array cannot
978! be used here because of weird memory problems.
979!
980 allocate ( tl_zeta_new(imins:imaxs,jmins:jmaxs) )
981 tl_zeta_new = 0.0_r8
982!
983! Get background zeta_new from BASIC state. Notice the I- and J-range
984! used to avoid calling nonlinear 'zetabc_local' routine.
985!
986 DO j=jr_range
987 DO i=ir_range
988 zeta_new(i,j)=zeta(i,j,knew)
989#ifdef MASKING
990 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
991# ifdef WET_DRY_NOT_YET
992!! zeta_new(i,j)=zeta_new(i,j)+ &
993!! & (Dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
994# endif
995#endif
996 dnew(i,j)=h(i,j)+zeta_new(i,j)
997 dnew_rd(i,j)=dnew(i,j)
998 dstp(i,j)=h(i,j)+zeta(i,j,kstp)
999 END DO
1000 END DO
1001
1002#undef IR_RANGE
1003#undef JR_RANGE
1004!
1005! Compute "zeta_new" at new time step and interpolate it half-step
1006! backward, "zwrk" for the subsequent computation of the tangent
1007! linear barotropic pressure gradient. Here, we use the BASIC STATE
1008! values. Thus, the nonlinear correction to the pressure-gradient
1009! term from "kstp" to "knew" is not needed for Forward-Euler to
1010! Forward-Backward steps (PGF_FB_CORRECTION method).
1011!
1012 DO j=jstrv-1,jend
1013 DO i=istru-1,iend
1014 fac=dtfast(ng)*pm(i,j)*pn(i,j)
1015!^ zeta_new(i,j)=zeta(i,j,kstp)+ &
1016!^ & fac*(DUon(i,j)-DUon(i+1,j)+ &
1017!^ & DVom(i,j)-DVom(i,j+1))
1018!^
1019 tl_zeta_new(i,j)=tl_zeta(i,j,kstp)+ &
1020 & fac*(tl_duon(i,j)-tl_duon(i+1,j)+ &
1021 & tl_dvom(i,j)-tl_dvom(i,j+1))
1022#ifdef MASKING
1023!^ zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
1024!^
1025 tl_zeta_new(i,j)=tl_zeta_new(i,j)*rmask(i,j)
1026# ifdef WET_DRY_NOT_YET
1027!! zeta_new(i,j)=zeta_new(i,j)+ &
1028!! & (Dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
1029# endif
1030#endif
1031 zwrk(i,j)=bkw_new*zeta_new(i,j)+ &
1032 & bkw0*zeta(i,j,kstp)+ &
1033 & bkw1*zeta(i,j,kbak)+ &
1034 & bkw2*zeta(i,j,kold)
1035 tl_zwrk(i,j)=bkw_new*tl_zeta_new(i,j)+ &
1036 & bkw0*tl_zeta(i,j,kstp)+ &
1037 & bkw1*tl_zeta(i,j,kbak)+ &
1038 & bkw2*tl_zeta(i,j,kold)
1039
1040#if defined VAR_RHO_2D && defined SOLVE3D
1041 rzeta(i,j)=(1.0_r8+rhos(i,j))*zwrk(i,j)
1042 tl_rzeta(i,j)=(1.0_r8+rhos(i,j))*tl_zwrk(i,j)+ &
1043 & tl_rhos(i,j)*zwrk(i,j)
1044 rzeta2(i,j)=rzeta(i,j)*zwrk(i,j)
1045 tl_rzeta2(i,j)=tl_rzeta(i,j)*zwrk(i,j)+ &
1046 & rzeta(i,j)*tl_zwrk(i,j)
1047 rzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
1048 tl_rzetasa(i,j)=tl_zwrk(i,j)*(rhos(i,j)-rhoa(i,j))+ &
1049 & zwrk(i,j)*(tl_rhos(i,j)-tl_rhoa(i,j))
1050#else
1051 rzeta(i,j)=zwrk(i,j)
1052 tl_rzeta(i,j)=tl_zwrk(i,j)
1053 rzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
1054 tl_rzeta2(i,j)=2.0_r8*tl_zwrk(i,j)*zwrk(i,j)
1055#endif
1056 END DO
1057 END DO
1058!
1059! Apply mass point sources (volume vertical influx), if any.
1060!
1061! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
1062!
1063 IF (lwsrc(ng)) THEN
1064 DO is=1,nsrc(ng)
1065 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
1066 i=sources(ng)%Isrc(is)
1067 j=sources(ng)%Jsrc(is)
1068 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1069 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1070!^ zeta_new(i,j)=zeta_new(i,j)+ &
1071!^ & SOURCES(ng)%Qbar(is)* &
1072!^ & pm(i,j)*pn(i,j)*dtfast(ng)
1073!^
1074! tl_zeta_new(i,j)=tl_zeta_new(i,j)+0.0_r8
1075 END IF
1076 END IF
1077 END DO
1078 END IF
1079!
1080! Apply boundary conditions to newly computed free-surface "zeta_new".
1081! Notice that we are using the local routine, which passes the private
1082! "zeta_new" array as argument.
1083!
1084! Here, we use the local "zetabc" since the private array "zeta_new"
1085! is passed as an argument to allow computing the lateral boundary
1086! conditions on the range IstrU-1:Iend and JstrV-1:Jend, so parallel
1087! tile exchanges are avoided.
1088!
1089!^ CALL zetabc_local (ng, tile, &
1090!^ & LBi, UBi, LBj, UBj, &
1091!^ & IminS, ImaxS, JminS, JmaxS, &
1092!^ & kstp, &
1093!^ & zeta, &
1094!^ & zeta_new)
1095!^
1096 CALL tl_zetabc_local (ng, tile, &
1097 & lbi, ubi, lbj, ubj, &
1098 & imins, imaxs, jmins, jmaxs, &
1099 & kstp, &
1100 & zeta, tl_zeta, &
1101 & zeta_new, tl_zeta_new)
1102!
1103! Load new computed free-surface into global state array.
1104!
1105 DO j=jstrr,jendr
1106 DO i=istrr,iendr
1107!^ zeta(i,j,knew)=zeta_new(i,j)
1108!^
1109 tl_zeta(i,j,knew)=tl_zeta_new(i,j)
1110 END DO
1111 END DO
1112
1113#ifdef SOLVE3D
1114!
1115!-----------------------------------------------------------------------
1116! Compute fast-time-averaged fields over all barotropic time steps.
1117!-----------------------------------------------------------------------
1118!
1119! Reset/initialize arrays for averaged fields during the first
1120! barotropic time step. Then, accumulate it time average. Include
1121! physical boundary points, but not periodic ghost points or
1122! computation distributed-memory computational margins.
1123!
1124 cff1=weight(1,iif(ng),ng)
1125 cff2=weight(2,iif(ng),ng)
1126!
1127 IF (first_2d_step) THEN
1128 DO j=jstrr,jendr
1129 DO i=istrr,iendr
1130!^ Zt_avg1(i,j)=cff1*zeta(i,j,knew)
1131!^
1132 tl_zt_avg1(i,j)=cff1*tl_zeta(i,j,knew)
1133 IF (i.ge.istr) THEN
1134!^ DU_avg1(i,j)=0.0_r8
1135!^
1136 tl_du_avg1(i,j)=0.0_r8
1137!^ DU_avg2(i,j)=cff2*DUon(i,j)
1138!^
1139 tl_du_avg2(i,j)=cff2*tl_duon(i,j)
1140 END IF
1141 IF (j.ge.jstr) THEN
1142!^ DV_avg1(i,j)=0.0_r8
1143!^
1144 tl_dv_avg1(i,j)=0.0_r8
1145!^ DV_avg2(i,j)=cff2*DVom(i,j)
1146!^
1147 tl_dv_avg2(i,j)=cff2*tl_dvom(i,j)
1148 END IF
1149 END DO
1150 END DO
1151 ELSE
1152 DO j=jstrr,jendr
1153 DO i=istrr,iendr
1154!^ Zt_avg1(i,j)=Zt_avg1(i,j)+cff1*zeta(i,j,knew)
1155!^
1156 tl_zt_avg1(i,j)=tl_zt_avg1(i,j)+cff1*tl_zeta(i,j,knew)
1157 IF (i.ge.istr) THEN
1158!^ DU_avg2(i,j)=DU_avg2(i,j)+cff2*DUon(i,j)
1159!^
1160 tl_du_avg2(i,j)=tl_du_avg2(i,j)+cff2*tl_duon(i,j)
1161 END IF
1162 IF (j.ge.jstr) THEN
1163!^ DV_avg2(i,j)=DV_avg2(i,j)+cff2*DVom(i,j)
1164!^
1165 tl_dv_avg2(i,j)=tl_dv_avg2(i,j)+cff2*tl_dvom(i,j)
1166 END IF
1167 END DO
1168 END DO
1169 END IF
1170#endif
1171!
1172!=======================================================================
1173! Compute right-hand-side for the 2D momentum equations.
1174!=======================================================================
1175#ifdef SOLVE3D
1176!
1177! Notice that we are suppressing the computation of momentum advection,
1178! Coriolis, and lateral viscosity terms in 3D Applications because
1179! these terms are already included in the baroclinic-to-barotropic
1180! forcing arrays "rufrc" and "rvfrc". It does not mean we are entirely
1181! omitting them, but it is a choice between recomputing them at every
1182! barotropic step or keeping them "frozen" during the fast-time
1183! stepping.
1184# ifdef STEP2D_CORIOLIS
1185! However, in some coarse grid applications with larger baroclinic
1186! timestep (say, DT around 20 minutes or larger), adding the Coriolis
1187! term in the barotropic equations is useful since f*DT is no longer
1188! small.
1189# endif
1190#endif
1191!
1192!-----------------------------------------------------------------------
1193! Compute pressure-gradient terms.
1194!-----------------------------------------------------------------------
1195!
1196 cff1=0.5_r8*g
1197#if defined VAR_RHO_2D && defined SOLVE3D
1198 cff2=0.333333333333_r8
1199#endif
1200#if defined ATM_PRESS && !defined SOLVE3D
1201 cff3=0.5_r8*100.0_r8/rho0
1202#endif
1203 DO j=jstr,jend
1204 DO i=istr,iend
1205 IF (i.ge.istru) THEN
1206!^ rubar(i,j)=cff1*on_u(i,j)* &
1207!^ & ((h(i-1,j)+ &
1208!^ & h(i ,j))* &
1209!^ & (rzeta(i-1,j)- &
1210!^ & rzeta(i ,j))+ &
1211#if defined VAR_RHO_2D && defined SOLVE3D
1212!^ & (h(i-1,j)- &
1213!^ & h(i ,j))* &
1214!^ & (rzetaSA(i-1,j)+ &
1215!^ & rzetaSA(i ,j)+ &
1216!^ & cff2*(rhoA(i-1,j)- &
1217!^ & rhoA(i ,j))* &
1218!^ & (zwrk(i-1,j)- &
1219!^ & zwrk(i,j)))+ &
1220#endif
1221!^ & (rzeta2(i-1,j)- &
1222!^ & rzeta2(i ,j)))
1223!^
1224 tl_rubar(i,j)=cff1*on_u(i,j)* &
1225 & ((tl_h(i-1,j)+ &
1226 & tl_h(i ,j))* &
1227 & (rzeta(i-1,j)- &
1228 & rzeta(i ,j))+ &
1229 & (h(i-1,j)+ &
1230 & h(i ,j))* &
1231 & (tl_rzeta(i-1,j)- &
1232 & tl_rzeta(i ,j))+ &
1233#if defined VAR_RHO_2D && defined SOLVE3D
1234 & (tl_h(i-1,j)- &
1235 & tl_h(i ,j))* &
1236 & (rzetasa(i-1,j)+ &
1237 & rzetasa(i ,j)+ &
1238 & cff2*(rhoa(i-1,j)- &
1239 & rhoa(i ,j))* &
1240 & (zwrk(i-1,j)- &
1241 & zwrk(i ,j)))+ &
1242 & (h(i-1,j)- &
1243 & h(i ,j))* &
1244 & (tl_rzetasa(i-1,j)+ &
1245 & tl_rzetasa(i ,j)+ &
1246 & cff2*((tl_rhoa(i-1,j)- &
1247 & tl_rhoa(i ,j))* &
1248 & (zwrk(i-1,j)- &
1249 & zwrk(i ,j))+ &
1250 & (rhoa(i-1,j)- &
1251 & rhoa(i ,j))* &
1252 & (tl_zwrk(i-1,j)- &
1253 & tl_zwrk(i ,j))))+ &
1254#endif
1255 & (tl_rzeta2(i-1,j)- &
1256 & tl_rzeta2(i ,j)))
1257#if defined ATM_PRESS && !defined SOLVE3D
1258!^ rubar(i,j)=rubar(i,j)- &
1259!^ & cff3*on_u(i,j)* &
1260!^ & (h(i-1,j)+h(i,j)+ &
1261!^ & rzeta(i-1,j)+rzeta(i,j))* &
1262!^ & (Pair(i,j)-Pair(i-1,j))
1263!^
1264 tl_rubar(i,j)=tl_rubar(i,j)- &
1265 & cff3*on_u(i,j)* &
1266 & (tl_h(i-1,j)+tl_h(i,j)+ &
1267 & tl_rzeta(i-1,j)+tl_rzeta(i,j))* &
1268 & (pair(i,j)-pair(i-1,j))
1269#endif
1270#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
1271!^ rubar(i,j)=rubar(i,j)- &
1272!^ & cff1*on_u(i,j)* &
1273!^ & (h(i-1,j)+h(i,j)+ &
1274!^ & rzeta(i-1,j)+rzeta(i,j))* &
1275!^ & (eq_tide(i,j)-eq_tide(i-1,j))
1276!^
1277 tl_rubar(i,j)=tl_rubar(i,j)- &
1278 & cff1*on_u(i,j)* &
1279 & ((tl_h(i-1,j)+tl_h(i,j)+ &
1280 & tl_rzeta(i-1,j)+tl_rzeta(i,j))* &
1281 & (eq_tide(i,j)-eq_tide(i-1,j))+ &
1282 & (h(i-1,j)+h(i,j)+ &
1283 & rzeta(i-1,j)+rzeta(i,j))* &
1284 & (tl_eq_tide(i,j)-tl_eq_tide(i-1,j)))
1285#endif
1286#ifdef DIAGNOSTICS_UV
1287!! DiaU2rhs(i,j,M2pgrd)=rubar(i,j)
1288#endif
1289 END IF
1290!
1291 IF (j.ge.jstrv) THEN
1292!^ rvbar(i,j)=cff1*om_v(i,j)* &
1293!^ & ((h(i,j-1)+ &
1294!^ & h(i,j ))* &
1295!^ & (rzeta(i,j-1)- &
1296!^ & rzeta(i,j ))+ &
1297#if defined VAR_RHO_2D && defined SOLVE3D
1298!^ & (h(i,j-1)- &
1299!^ & h(i,j ))* &
1300!^ & (rzetaSA(i,j-1)+ &
1301!^ & rzetaSA(i,j )+ &
1302!^ & cff2*(rhoA(i,j-1)- &
1303!^ & rhoA(i,j ))* &
1304!^ & (zwrk(i,j-1)- &
1305!^ & zwrk(i,j )))+ &
1306#endif
1307!^ & (rzeta2(i,j-1)- &
1308!^ & rzeta2(i,j )))
1309!^
1310 tl_rvbar(i,j)=cff1*om_v(i,j)* &
1311 & ((tl_h(i,j-1)+ &
1312 & tl_h(i,j ))* &
1313 & (rzeta(i,j-1)- &
1314 & rzeta(i,j ))+ &
1315 & (h(i,j-1)+ &
1316 & h(i,j ))* &
1317 & (tl_rzeta(i,j-1)- &
1318 & tl_rzeta(i,j ))+ &
1319#if defined VAR_RHO_2D && defined SOLVE3D
1320 & (tl_h(i,j-1)- &
1321 & tl_h(i,j ))* &
1322 & (rzetasa(i,j-1)+ &
1323 & rzetasa(i,j )+ &
1324 & cff2*(rhoa(i,j-1)- &
1325 & rhoa(i,j ))* &
1326 & (zwrk(i,j-1)- &
1327 & zwrk(i,j )))+ &
1328 & (h(i,j-1)- &
1329 & h(i,j ))* &
1330 & (tl_rzetasa(i,j-1)+ &
1331 & tl_rzetasa(i,j )+ &
1332 & cff2*((tl_rhoa(i,j-1)- &
1333 & tl_rhoa(i,j ))* &
1334 & (zwrk(i,j-1)- &
1335 & zwrk(i,j ))+ &
1336 & (rhoa(i,j-1)- &
1337 & rhoa(i,j ))* &
1338 & (tl_zwrk(i,j-1)- &
1339 & tl_zwrk(i,j ))))+ &
1340#endif
1341 & (tl_rzeta2(i,j-1)- &
1342 & tl_rzeta2(i,j )))
1343#if defined ATM_PRESS && !defined SOLVE3D
1344!^ rvbar(i,j)=rvbar(i,j)- &
1345!^ & cff3*om_v(i,j)* &
1346!^ & (h(i,j-1)+h(i,j)+ &
1347!^ & rzeta(i,j-1)+rzeta(i,j))* &
1348!^ & (Pair(i,j)-Pair(i,j-1))
1349!^
1350 tl_rvbar(i,j)=tl_rvbar(i,j)- &
1351 & cff3*om_v(i,j)* &
1352 & (tl_h(i,j-1)+tl_h(i,j)+ &
1353 & tl_rzeta(i,j-1)+tl_rzeta(i,j))* &
1354 & (pair(i,j)-pair(i,j-1))
1355#endif
1356#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
1357!^ rvbar(i,j)=rvbar(i,j)- &
1358!^ & cff1*om_v(i,j)* &
1359!^ & (h(i,j-1)+h(i,j)+ &
1360!^ & rzeta(i,j-1)+rzeta(i,j))* &
1361!^ & (eq_tide(i,j)-eq_tide(i,j-1))
1362!^
1363 tl_rvbar(i,j)=tl_rvbar(i,j)- &
1364 & cff1*om_v(i,j)* &
1365 & ((tl_h(i,j-1)+tl_h(i,j)+ &
1366 & tl_rzeta(i,j-1)+tl_rzeta(i,j))* &
1367 & (eq_tide(i,j)-eq_tide(i,j-1))+ &
1368 & (h(i,j-1)+h(i,j)+ &
1369 & rzeta(i,j-1)+rzeta(i,j))* &
1370 & (tl_eq_tide(i,j)-tl_eq_tide(i,j-1)))
1371#endif
1372#ifdef DIAGNOSTICS_UV
1373!! DiaV2rhs(i,j,M2pgrd)=rvbar(i,j)
1374#endif
1375 END IF
1376 END DO
1377 END DO
1378
1379#if defined UV_ADV && !defined SOLVE3D
1380!
1381!-----------------------------------------------------------------------
1382! Add in horizontal advection of momentum.
1383!-----------------------------------------------------------------------
1384
1385# ifdef UV_C2ADVECTION
1386!
1387! Second-order, centered differences advection fluxes.
1388!
1389 DO j=jstr,jend
1390 DO i=istr-1,iend
1391 IF (i.ge.istru-1) THEN
1392!^ UFx(i,j)=0.25_r8* &
1393!^ & (DUon(i,j)+DUon(i+1,j))* &
1394!^ & (urhs(i ,j)+ &
1395# ifdef WEC_MELLOR
1396!^ & ubar_stokes(i ,j)+ &
1397!^ & ubar_stokes(i+1,j)+ &
1398# endif
1399!^ & urhs(i+1,j))
1400!^
1401 tl_ufx(i,j)=0.25_r8* &
1402 & ((tl_duon(i,j)+tl_duon(i+1,j))* &
1403 & (urhs(i ,j)+ &
1404# ifdef WEC_MELLOR
1405 & ubar_stokes(i ,j)+ &
1406 & ubar_stokes(i+1,j)+ &
1407# endif
1408 & urhs(i+1,j))+ &
1409 & (duon(i,j)+duon(i+1,j))* &
1410 & (tl_urhs(i ,j)+ &
1411# ifdef WEC_MELLOR
1412 & tl_ubar_stokes(i ,j)+ &
1413 & tl_ubar_stokes(i+1,j)+ &
1414# endif
1415 & tl_urhs(i+1,j)))
1416 END IF
1417!
1418!^ VFx(i+1,j)=0.25_r8* &
1419# ifdef MASKING
1420!^ pmask(i+1,j)* &
1421# endif
1422!^ & (DUon(i+1,j)+DUon(i+1,j-1))* &
1423!^ & (vrhs(i+1,j)+ &
1424# ifdef WEC_MELLOR
1425!^ & vbar_stokes(i ,j)+ &
1426!^ & vbar_stokes(i-1,j)+ &
1427# endif
1428!^ & vrhs(i ,j))
1429!^
1430 tl_vfx(i+1,j)=0.25_r8* &
1431# ifdef MASKING
1432 & pmask(i+1,j)* &
1433# endif
1434 & ((tl_duon(i+1,j)+tl_duon(i+1,j-1))* &
1435 & (vrhs(i+1,j)+ &
1436# ifdef WEC_MELLOR
1437 & vbar_stokes(i ,j)+ &
1438 & vbar_stokes(i-1,j)+ &
1439# endif
1440 & vrhs(i ,j))+ &
1441 & (duon(i+1,j)+duon(i+1,j-1))* &
1442 & (tl_vrhs(i+1,j)+ &
1443# ifdef WEC_MELLOR
1444 & tl_vbar_stokes(i ,j)+ &
1445 & tl_vbar_stokes(i-1,j)+ &
1446# endif
1447 & tl_vrhs(i ,j)))
1448 END DO
1449 END DO
1450!
1451 DO j=jstr-1,jend
1452 DO i=istr,iend
1453 IF (j.ge.jstrv-1) THEN
1454!^ VFe(i,j)=0.25_r8* &
1455!^ & (DVom(i,j)+DVom(i,j+1))* &
1456!^ & (vrhs(i,j )+ &
1457# ifdef WEC_MELLOR
1458!^ & vbar_stokes(i,j )+ &
1459!^ & vbar_stokes(i,j+1)+ &
1460# endif
1461!^ & vrhs(i,j+1))
1462!^
1463 tl_vfe(i,j)=0.25_r8* &
1464 & ((tl_dvom(i,j)+tl_dvom(i,j+1))* &
1465 & (vrhs(i,j )+ &
1466# ifdef WEC_MELLOR
1467 & vbar_stokes(i,j )+ &
1468 & vbar_stokes(i,j+1)+ &
1469# endif
1470 & vrhs(i,j+1))+ &
1471 & (dvom(i,j)+dvom(i,j+1))* &
1472 & (tl_vrhs(i,j )+ &
1473# ifdef WEC_MELLOR
1474 & tl_vbar_stokes(i,j )+ &
1475 & tl_vbar_stokes(i,j+1)+ &
1476# endif
1477 & tl_vrhs(i,j+1)))
1478 END IF
1479!
1480!^ UFe(i,j+1)=0.25_r8* &
1481# ifdef MASKING
1482!^ & pmask(i,j+1)* &
1483# endif
1484!^ & (DVom(i,j+1)+DVom(i-1,j+1))* &
1485!^ & (urhs(i,j+1)+ &
1486# ifdef WEC_MELLOR
1487!^ & ubar_stokes(i,j+1)+ &
1488!^ & ubar_stokes(i,j )+ &
1489# endif
1490!^ & urhs(i,j ))
1491!^
1492 tl_ufe(i,j+1)=0.25_r8* &
1493# ifdef MASKING
1494 & pmask(i,j+1)* &
1495# endif
1496 & ((tl_dvom(i,j+1)+tl_dvom(i-1,j+1))* &
1497 & (urhs(i,j+1)+ &
1498# ifdef WEC_MELLOR
1499 & ubar_stokes(i,j+1)+ &
1500 & ubar_stokes(i,j )+ &
1501# endif
1502 & urhs(i,j ))+ &
1503 & (dvom(i,j+1)+dvom(i-1,j+1))* &
1504 & (tl_urhs(i,j+1)+ &
1505# ifdef WEC_MELLOR
1506 & tl_ubar_stokes(i,j+1)+ &
1507 & tl_ubar_stokes(i,j )+ &
1508# endif
1509 & tl_urhs(i,j )))
1510 END DO
1511 END DO
1512
1513# elif defined UV_C4ADVECTION
1514!
1515! Fourth-order, centered differences u-momentum advection fluxes.
1516!
1517 DO j=jstr,jend
1518 DO i=istrum1,iendp1
1519 grad(i,j)=urhs(i-1,j)-2.0_r8*urhs(i,j)+ &
1520# ifdef WEC_MELLOR
1521 & ubar_stokes(i-1,j)-2.0_r8*ubar_stokes(i,j)+ &
1522 & ubar_stokes(i+1,j)+ &
1523# endif
1524 & urhs(i+1,j)
1525 tl_grad(i,j)=tl_urhs(i-1,j)-2.0_r8*tl_urhs(i,j)+ &
1526# ifdef WEC_MELLOR
1527 & tl_ubar_stokes(i-1,j)-2.0_r8*tl_ubar_stokes(i,j)+&
1528 & tl_ubar_stokes(i+1,j)+ &
1529# endif
1530 & tl_urhs(i+1,j)
1531 dgrad(i,j)=duon(i-1,j)-2.0_r8*duon(i,j)+duon(i+1,j)
1532 tl_dgrad(i,j)=tl_duon(i-1,j)-2.0_r8*tl_duon(i,j)+ &
1533 & tl_duon(i+1,j)
1534 END DO
1535 END DO
1536 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1537 IF (domain(ng)%Western_Edge(tile)) THEN
1538 DO j=jstr,jend
1539 grad(istr,j)=grad(istr+1,j)
1540 tl_grad(istr,j)=tl_grad(istr+1,j)
1541 dgrad(istr,j)=dgrad(istr+1,j)
1542 tl_dgrad(istr,j)=tl_dgrad(istr+1,j)
1543 END DO
1544 END IF
1545 END IF
1546 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1547 IF (domain(ng)%Eastern_Edge(tile)) THEN
1548 DO j=jstr,jend
1549 grad(iend+1,j)=grad(iend,j)
1550 tl_grad(iend+1,j)=tl_grad(iend,j)
1551 dgrad(iend+1,j)=dgrad(iend,j)
1552 tl_dgrad(iend+1,j)=tl_dgrad(iend,j)
1553 END DO
1554 END IF
1555 END IF
1556! d/dx(Duu/n)
1557 cff=1.0_r8/6.0_r8
1558 DO j=jstr,jend
1559 DO i=istru-1,iend
1560!^ UFx(i,j)=0.25_r8*(urhs(i ,j)+ &
1561# ifdef WEC_MELLOR
1562!^ & ubar_stokes(i ,j)+ &
1563!^ & ubar_stokes(i+1,j)+ &
1564# endif
1565!^ & urhs(i+1,j)- &
1566!^ & cff*(grad (i,j)+grad (i+1,j)))* &
1567!^ & (DUon(i,j)+DUon(i+1,j)- &
1568!^ & cff*(Dgrad(i,j)+Dgrad(i+1,j)))
1569!^
1570 tl_ufx(i,j)=0.25_r8* &
1571 & ((urhs(i ,j)+ &
1572# ifdef WEC_MELLOR
1573 & ubar_stokes(i ,j)+ &
1574 & ubar_stokes(i+1,j)+ &
1575# endif
1576 & urhs(i+1,j)- &
1577 & cff*(grad(i,j)+grad(i+1,j)))* &
1578 & (tl_duon(i,j)+tl_duon(i+1,j)- &
1579 & cff*(tl_dgrad(i,j)+tl_dgrad(i+1,j)))+ &
1580 & (tl_urhs(i ,j)+ &
1581# ifdef WEC_MELLOR
1582 & tl_ubar_stokes(i ,j)+ &
1583 & tl_ubar_stokes(i+1,j)+ &
1584# endif
1585 & tl_urhs(i+1,j)- &
1586 & cff*(tl_grad(i,j)+tl_grad(i+1,j)))* &
1587 & (duon(i,j)+duon(i+1,j)- &
1588 & cff*(dgrad(i,j)+dgrad(i+1,j))))
1589 END DO
1590 END DO
1591!
1592 DO j=jstrm1,jendp1
1593 DO i=istru,iend
1594 grad(i,j)=urhs(i,j-1)-2.0_r8*urhs(i,j)+ &
1595# ifdef WEC_MELLOR
1596 & ubar_stokes(i,j-1)-2.0_r8*ubar_stokes(i,j)+ &
1597 & ubar_stokes(i,j+1)+ &
1598# endif
1599 & urhs(i,j+1)
1600 tl_grad(i,j)=tl_urhs(i,j-1)-2.0_r8*tl_urhs(i,j)+ &
1601# ifdef WEC_MELLOR
1602 & tl_ubar_stokes(i,j-1)-2.0_r8*tl_ubar_stokes(i,j)+&
1603 & tl_ubar_stokes(i,j+1)+ &
1604# endif
1605 & tl_urhs(i,j+1)
1606 END DO
1607 END DO
1608 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1609 IF (domain(ng)%Southern_Edge(tile)) THEN
1610 DO i=istru,iend
1611 grad(i,jstr-1)=grad(i,jstr)
1612 tl_grad(i,jstr-1)=tl_grad(i,jstr)
1613 END DO
1614 END IF
1615 END IF
1616 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1617 IF (domain(ng)%Northern_Edge(tile)) THEN
1618 DO i=istru,iend
1619 grad(i,jend+1)=grad(i,jend)
1620 tl_grad(i,jend+1)=tl_grad(i,jend)
1621 END DO
1622 END IF
1623 END IF
1624 DO j=jstr,jend+1
1625 DO i=istru-1,iend
1626 dgrad(i,j)=dvom(i-1,j)-2.0_r8*dvom(i,j)+dvom(i+1,j)
1627 tl_dgrad(i,j)=tl_dvom(i-1,j)-2.0_r8*tl_dvom(i,j)+ &
1628 & tl_dvom(i+1,j)
1629 END DO
1630 END DO
1631! d/dy(Duv/m)
1632 cff=1.0_r8/6.0_r8
1633 DO j=jstr,jend+1
1634 DO i=istru,iend
1635!^ UFe(i,j)=0.25_r8*(urhs(i,j )+ &
1636# ifdef WEC_MELLOR
1637!^ & ubar_stokes(i,j )+ &
1638!^ & ubar_stokes(i,j-1)+ &
1639# endif
1640!^ & urhs(i,j-1)- &
1641!^ & cff*(grad (i,j)+grad (i,j-1)))* &
1642!^ & (DVom(i,j)+DVom(i-1,j)- &
1643!^ & cff*(Dgrad(i,j)+Dgrad(i-1,j)))
1644!^
1645 tl_ufe(i,j)=0.25_r8* &
1646 & ((tl_urhs(i,j )+ &
1647# ifdef WEC_MELLOR
1648 & tl_ubar_stokes(i,j )+ &
1649 & tl_ubar_stokes(i,j-1)+ &
1650# endif
1651 & tl_urhs(i,j-1)- &
1652 & cff*(tl_grad(i,j)+tl_grad(i,j-1)))* &
1653 & (dvom(i,j)+dvom(i-1,j)- &
1654 & cff*(dgrad(i,j)+dgrad(i-1,j)))+ &
1655 & (urhs(i,j )+ &
1656# ifdef WEC_MELLOR
1657 & ubar_stokes(i,j )+ &
1658 & ubar_stokes(i,j-1)+ &
1659# endif
1660 & urhs(i,j-1)- &
1661 & cff*(grad(i,j)+grad(i,j-1)))* &
1662 & (tl_dvom(i,j)+tl_dvom(i-1,j)- &
1663 & cff*(tl_dgrad(i,j)+tl_dgrad(i-1,j))))
1664 END DO
1665 END DO
1666!
1667! Fourth-order, centered differences v-momentum advection fluxes.
1668!
1669 DO j=jstrv,jend
1670 DO i=istrm1,iendp1
1671 grad(i,j)=vrhs(i-1,j)-2.0_r8*vrhs(i,j)+ &
1672# ifdef WEC_MELLOR
1673 & vbar_stokes(i-1,j)-2.0_r8*vbar_stokes(i,j)+ &
1674 & vbar_stokes(i+1,j)+ &
1675# endif
1676 & vrhs(i+1,j)
1677 tl_grad(i,j)=tl_vrhs(i-1,j)-2.0_r8*tl_vrhs(i,j)+ &
1678# ifdef WEC_MELLOR
1679 & tl_vbar_stokes(i-1,j)-2.0_r8*tl_vbar_stokes(i,j)+&
1680 & tl_vbar_stokes(i+1,j)+ &
1681# endif
1682 & tl_vrhs(i+1,j)
1683 END DO
1684 END DO
1685 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1686 IF (domain(ng)%Western_Edge(tile)) THEN
1687 DO j=jstrv,jend
1688 grad(istr-1,j)=grad(istr,j)
1689 tl_grad(istr-1,j)=tl_grad(istr,j)
1690 END DO
1691 END IF
1692 END IF
1693 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1694 IF (domain(ng)%Eastern_Edge(tile)) THEN
1695 DO j=jstrv,jend
1696 grad(iend+1,j)=grad(iend,j)
1697 tl_grad(iend+1,j)=tl_grad(iend,j)
1698 END DO
1699 END IF
1700 END IF
1701 DO j=jstrv-1,jend
1702 DO i=istr,iend+1
1703 dgrad(i,j)=duon(i,j-1)-2.0_r8*duon(i,j)+duon(i,j+1)
1704 tl_dgrad(i,j)=tl_duon(i,j-1)-2.0_r8*tl_duon(i,j)+ &
1705 & tl_duon(i,j+1)
1706 END DO
1707 END DO
1708! d/dx(Duv/n)
1709 cff=1.0_r8/6.0_r8
1710 DO j=jstrv,jend
1711 DO i=istr,iend+1
1712!^ VFx(i,j)=0.25_r8*(vrhs(i ,j)+ &
1713# ifdef WEC_MELLOR
1714!^ & vbar_stokes(i ,j)+ &
1715!^ & vbar_stokes(i-1,j)+ &
1716# endif
1717!^ & vrhs(i-1,j)- &
1718!^ & cff*(grad (i,j)+grad (i-1,j)))* &
1719!^ & (DUon(i,j)+DUon(i,j-1)- &
1720!^ & cff*(Dgrad(i,j)+Dgrad(i,j-1)))
1721!^
1722 tl_vfx(i,j)=0.25_r8* &
1723 & ((tl_vrhs(i ,j)+ &
1724# ifdef WEC_MELLOR
1725 & tl_vbar_stokes(i ,j)+ &
1726 & tl_vbar_stokes(i-1,j)+ &
1727# endif
1728 & tl_vrhs(i-1,j)- &
1729 & cff*(tl_grad(i,j)+tl_grad(i-1,j)))* &
1730 & (duon(i,j)+duon(i,j-1)- &
1731 & cff*(dgrad(i,j)+dgrad(i,j-1)))+ &
1732 & (vrhs(i ,j)+ &
1733# ifdef WEC_MELLOR
1734 & vbar_stokes(i ,j)+ &
1735 & vbar_stokes(i-1,j)+ &
1736# endif
1737 & vrhs(i-1,j)- &
1738 & cff*(grad(i,j)+grad(i-1,j)))* &
1739 & (tl_duon(i,j)+tl_duon(i,j-1)- &
1740 & cff*(tl_dgrad(i,j)+tl_dgrad(i,j-1))))
1741 END DO
1742 END DO
1743!
1744 DO j=jstrvm1,jendp1
1745 DO i=istr,iend
1746 grad(i,j)=vrhs(i,j-1)-2.0_r8*vrhs(i,j)+ &
1747# ifdef WEC_MELLOR
1748 & vbar_stokes(i,j-1)-2.0_r8*vbar_stokes(i,j)+ &
1749 & vbar_stokes(i,j+1)+ &
1750# endif
1751 & vrhs(i,j+1)
1752 tl_grad(i,j)=tl_vrhs(i,j-1)-2.0_r8*tl_vrhs(i,j)+ &
1753# ifdef WEC_MELLOR
1754 & tl_vbar_stokes(i,j-1)-2.0_r8*tl_vbar_stokes(i,j)+&
1755 & tl_vbar_stokes(i,j+1)+ &
1756# endif
1757 & tl_vrhs(i,j+1)
1758 dgrad(i,j)=dvom(i,j-1)-2.0_r8*dvom(i,j)+dvom(i,j+1)
1759 tl_dgrad(i,j)=tl_dvom(i,j-1)-2.0_r8*tl_dvom(i,j)+ &
1760 & tl_dvom(i,j+1)
1761 END DO
1762 END DO
1763 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1764 IF (domain(ng)%Southern_Edge(tile)) THEN
1765 DO i=istr,iend
1766 grad(i,jstr)=grad(i,jstr+1)
1767 tl_grad(i,jstr)=tl_grad(i,jstr+1)
1768 dgrad(i,jstr)=dgrad(i,jstr+1)
1769 tl_dgrad(i,jstr)=tl_dgrad(i,jstr+1)
1770 END DO
1771 END IF
1772 END IF
1773 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1774 IF (domain(ng)%Northern_Edge(tile)) THEN
1775 DO i=istr,iend
1776 grad(i,jend+1)=grad(i,jend)
1777 tl_grad(i,jend+1)=tl_grad(i,jend)
1778 dgrad(i,jend+1)=dgrad(i,jend)
1779 tl_dgrad(i,jend+1)=tl_dgrad(i,jend)
1780 END DO
1781 END IF
1782 END IF
1783! d/dy(Dvv/m)
1784 cff=1.0_r8/6.0_r8
1785 DO j=jstrv-1,jend
1786 DO i=istr,iend
1787!^ VFe(i,j)=0.25_r8*(vrhs(i,j )+ &
1788# ifdef WEC_MELLOR
1789!^ & vbar_stokes(i,j )+ &
1790!^ & vbar_stokes(i,j+1)+ &
1791# endif
1792!^ & vrhs(i,j+1)- &
1793!^ & cff*(grad (i,j)+grad (i,j+1)))* &
1794!^ & (DVom(i,j)+DVom(i,j+1)- &
1795!^ & cff*(Dgrad(i,j)+Dgrad(i,j+1)))
1796!^
1797 tl_vfe(i,j)=0.25_r8* &
1798 & ((tl_vrhs(i,j )+ &
1799# ifdef WEC_MELLOR
1800 & tl_vbar_stokes(i,j )+ &
1801 & tl_vbar_stokes(i,j+1)+ &
1802# endif
1803 & tl_vrhs(i,j+1)- &
1804 & cff*(tl_grad(i,j)+tl_grad(i,j+1)))* &
1805 & (dvom(i,j)+dvom(i,j+1)- &
1806 & cff*(dgrad(i,j)+dgrad(i,j+1)))+ &
1807 & (vrhs(i,j )+ &
1808# ifdef WEC_MELLOR
1809 & vbar_stokes(i,j )+ &
1810 & vbar_stokes(i,j+1)+ &
1811# endif
1812 & vrhs(i,j+1)- &
1813 & cff*(grad(i,j)+grad(i,j+1)))* &
1814 & (tl_dvom(i,j)+tl_dvom(i,j+1)- &
1815 & cff*(tl_dgrad(i,j)+tl_dgrad(i,j+1))))
1816 END DO
1817 END DO
1818# endif
1819!
1820! Add advection to RHS terms.
1821!
1822 DO j=jstr,jend
1823 DO i=istr,iend
1824 IF (i.ge.istru) THEN
1825!^ cff1=UFx(i,j)-UFx(i-1,j)
1826!^
1827 tl_cff1=tl_ufx(i,j)-tl_ufx(i-1,j)
1828!^ cff2=UFe(i,j+1)-UFe(i,j)
1829!^
1830 tl_cff2=tl_ufe(i,j+1)-tl_ufe(i,j)
1831!^ fac1=cff1+cff2
1832!^
1833 tl_fac1=tl_cff1+tl_cff2
1834!^ rubar(i,j)=rubar(i,j)-fac1
1835!^
1836 tl_rubar(i,j)=tl_rubar(i,j)-tl_fac1
1837# if defined DIAGNOSTICS_UV
1838!! DiaU2rhs(i,j,M2xadv)=-cff1
1839!! DiaU2rhs(i,j,M2yadv)=-cff2
1840!! DiaU2rhs(i,j,M2hadv)=-fac1
1841# endif
1842 END IF
1843!
1844 IF (j.ge.jstrv) THEN
1845!^ cff3=VFx(i+1,j)-VFx(i,j)
1846!^
1847 tl_cff3=tl_vfx(i+1,j)-tl_vfx(i,j)
1848!^ cff4=VFe(i,j)-VFe(i,j-1)
1849!^
1850 tl_cff4=tl_vfe(i,j)-tl_vfe(i,j-1)
1851!^ fac2=cff3+cff4
1852!^
1853 tl_fac2=tl_cff3+tl_cff4
1854!^ rvbar(i,j)=rvbar(i,j)-fac2
1855!^
1856 tl_rvbar(i,j)=tl_rvbar(i,j)-tl_fac2
1857# if defined DIAGNOSTICS_UV
1858!! DiaV2rhs(i,j,M2xadv)=-cff3
1859!! DiaV2rhs(i,j,M2yadv)=-cff4
1860!! DiaV2rhs(i,j,M2hadv)=-fac2
1861# endif
1862 END IF
1863 END DO
1864 END DO
1865#endif
1866
1867#if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
1868!
1869!-----------------------------------------------------------------------
1870! Add in Coriolis term.
1871!-----------------------------------------------------------------------
1872!
1873 DO j=jstrv-1,jend
1874 DO i=istru-1,iend
1875 cff=0.5_r8*drhs(i,j)*fomn(i,j)
1876 tl_cff=0.5_r8*tl_drhs(i,j)*fomn(i,j)
1877!^ UFx(i,j)=cff*(vrhs(i,j )+ &
1878# ifdef WEC_MELLOR
1879!^ & vbar_stokes(i,j )+ &
1880!^ & vbar_stokes(i,j+1)+ &
1881# endif
1882!^ & vrhs(i,j+1))
1883!^
1884 tl_ufx(i,j)=tl_cff*(vrhs(i,j )+ &
1885# ifdef WEC_MELLOR
1886 & vbar_stokes(i,j )+ &
1887 & vbar_stokes(i,j+1)+ &
1888# endif
1889 & vrhs(i,j+1))+ &
1890 & cff*(tl_vrhs(i,j )+ &
1891# ifdef WEC_MELLOR
1892 & tl_vbar_stokes(i,j )+ &
1893 & tl_vbar_stokes(i,j+1)+ &
1894# endif
1895 & tl_vrhs(i,j+1))
1896!
1897!^ VFe(i,j)=cff*(urhs(i ,j)+ &
1898# ifdef WEC_MELLOR
1899!^ & ubar_stokes(i ,j)+ &
1900!^ & ubar_stokes(i+1,j)+ &
1901# endif
1902!^ & urhs(i+1,j))
1903!^
1904 tl_vfe(i,j)=tl_cff*(urhs(i ,j)+ &
1905# ifdef WEC_MELLOR
1906 & ubar_stokes(i ,j)+ &
1907 & ubar_stokes(i+1,j)+ &
1908# endif
1909 & urhs(i+1,j))+ &
1910 & cff*(tl_urhs(i ,j)+ &
1911# ifdef WEC_MELLOR
1912 & tl_ubar_stokes(i ,j)+ &
1913 & tl_ubar_stokes(i+1,j)+ &
1914# endif
1915 & tl_urhs(i+1,j))
1916 END DO
1917 END DO
1918!
1919 DO j=jstr,jend
1920 DO i=istr,iend
1921 IF (i.ge.istru) THEN
1922!^ fac1=0.5_r8*(UFx(i,j)+UFx(i-1,j))
1923!^
1924 tl_fac1=0.5_r8*(tl_ufx(i,j)+tl_ufx(i-1,j))
1925!^ rubar(i,j)=rubar(i,j)+fac1
1926!^
1927 tl_rubar(i,j)=tl_rubar(i,j)+tl_fac1
1928# if defined DIAGNOSTICS_UV
1929!! DiaU2rhs(i,j,M2fcor)=fac1
1930# endif
1931 END IF
1932!
1933 IF (j.ge.jstrv) THEN
1934!^ fac2=0.5_r8*(VFe(i,j)+VFe(i,j-1))
1935!^
1936 tl_fac2=0.5_r8*(tl_vfe(i,j)+tl_vfe(i,j-1))
1937!^ rvbar(i,j)=rvbar(i,j)-fac2
1938!^
1939 tl_rvbar(i,j)=tl_rvbar(i,j)-tl_fac2
1940# if defined DIAGNOSTICS_UV
1941!! DiaV2rhs(i,j,M2fcor)=-fac2
1942# endif
1943 END IF
1944 END DO
1945 END DO
1946#endif
1947
1948#if (defined CURVGRID && defined UV_ADV) && !defined SOLVE3D
1949!
1950!-----------------------------------------------------------------------
1951! Add in curvilinear transformation terms.
1952!-----------------------------------------------------------------------
1953!
1954 DO j=jstrv-1,jend
1955 DO i=istru-1,iend
1956 cff1=0.5_r8*(vrhs(i,j )+ &
1957# ifdef WEC_MELLOR
1958 & vbar_stokes(i,j )+ &
1959 & vbar_stokes(i,j+1)+ &
1960# endif
1961 & vrhs(i,j+1))
1962 tl_cff1=0.5_r8*(tl_vrhs(i,j )+ &
1963# ifdef WEC_MELLOR
1964 & tl_vbar_stokes(i,j )+ &
1965 & tl_vbar_stokes(i,j+1)+ &
1966# endif
1967 & tl_vrhs(i,j+1))
1968!
1969 cff2=0.5_r8*(urhs(i ,j)+
1970# ifdef WEC_MELLOR
1971 & ubar_stokes(i ,j)+ &
1972 & ubar_stokes(i+1,j)+ &
1973# endif
1974 & urhs(i+1,j))
1975 tl_cff2=0.5_r8*(tl_urhs(i ,j)+ &
1976# ifdef WEC_MELLOR
1977 & tl_ubar_stokes(i ,j)+ &
1978 & tl_ubar_stokes(i+1,j)+ &
1979# endif
1980 & tl_urhs(i+1,j))
1981!
1982 cff3=cff1*dndx(i,j)
1983 tl_cff3=tl_cff1*dndx(i,j)
1984 cff4=cff2*dmde(i,j)
1985 tl_cff4=tl_cff2*dmde(i,j)
1986 cff=drhs(i,j)*(cff3-cff4)
1987 tl_cff=tl_drhs(i,j)*(cff3-cff4)+ &
1988 & drhs(i,j)*(tl_cff3-tl_cff4)
1989!^ UFx(i,j)=cff*cff1
1990!^
1991 tl_ufx(i,j)=tl_cff*cff1+cff*tl_cff1
1992!^ VFe(i,j)=cff*cff2
1993!^
1994 tl_vfe(i,j)=tl_cff*cff2+cff*tl_cff2
1995# if defined DIAGNOSTICS_UV
1996!! cff=Drhs(i,j)*cff4
1997!! Uwrk(i,j)=-cff*cff1 ! ubar equation, ETA-term
1998!! Vwrk(i,j)=-cff*cff2 ! vbar equation, ETA-term
1999# endif
2000 END DO
2001 END DO
2002!
2003 DO j=jstr,jend
2004 DO i=istr,iend
2005 IF (i.ge.istru) THEN
2006!^ fac1=0.5_r8*(UFx(i,j)+UFx(i-1,j))
2007!^
2008 tl_fac1=0.5_r8*(tl_ufx(i,j)+tl_ufx(i-1,j))
2009!^ rubar(i,j)=rubar(i,j)+fac1
2010!^
2011 tl_rubar(i,j)=tl_rubar(i,j)+tl_fac1
2012# if defined DIAGNOSTICS_UV
2013!! fac2=0.5_r8*(Uwrk(i,j)+Uwrk(i-1,j))
2014!! DiaU2rhs(i,j,M2xadv)=DiaU2rhs(i,j,M2xadv)+fac1-fac2
2015!! DiaU2rhs(i,j,M2yadv)=DiaU2rhs(i,j,M2yadv)+fac2
2016!! DiaU2rhs(i,j,M2hadv)=DiaU2rhs(i,j,M2hadv)+fac1
2017# endif
2018 END IF
2019!
2020 IF (j.ge.jstrv) THEN
2021!^ fac1=0.5_r8*(VFe(i,j)+VFe(i,j-1))
2022!^
2023 tl_fac1=0.5_r8*(tl_vfe(i,j)+tl_vfe(i,j-1))
2024!^ rvbar(i,j)=rvbar(i,j)-fac1
2025!^
2026 tl_rvbar(i,j)=tl_rvbar(i,j)-tl_fac1
2027# if defined DIAGNOSTICS_UV
2028!! fac2=0.5_r8*(Vwrk(i,j)+Vwrk(i,j-1))
2029!! DiaV2rhs(i,j,M2xadv)=DiaV2rhs(i,j,M2xadv)-fac1+fac2
2030!! DiaV2rhs(i,j,M2yadv)=DiaV2rhs(i,j,M2yadv)-fac2
2031!! DiaV2rhs(i,j,M2hadv)=DiaV2rhs(i,j,M2hadv)-fac1
2032# endif
2033 END IF
2034 END DO
2035 END DO
2036#endif
2037
2038#if defined UV_VIS2 && !defined SOLVE3D
2039!
2040!-----------------------------------------------------------------------
2041! Add in horizontal harmonic viscosity.
2042!-----------------------------------------------------------------------
2043!
2044! Compute total depth at PSI-points.
2045!
2046 DO j=jstr,jend+1
2047 DO i=istr,iend+1
2048 drhs_p(i,j)=0.25_r8*(drhs(i,j )+drhs(i-1,j )+ &
2049 & drhs(i,j-1)+drhs(i-1,j-1))
2050 tl_drhs_p(i,j)=0.25_r8*(tl_drhs(i,j )+tl_drhs(i-1,j )+ &
2051 & tl_drhs(i,j-1)+tl_drhs(i-1,j-1))
2052 END DO
2053 END DO
2054!
2055! Compute flux-components of the horizontal divergence of the stress
2056! tensor (m5/s2) in XI- and ETA-directions.
2057!
2058 DO j=jstrv-1,jend
2059 DO i=istru-1,iend
2060!^ cff=visc2_r(i,j)*Drhs(i,j)*0.5_r8* &
2061!^ & (pmon_r(i,j)* &
2062!^ & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,kstp)- &
2063!^ & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,kstp))- &
2064!^ & pnom_r(i,j)* &
2065!^ & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,kstp)- &
2066!^ & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,kstp)))
2067!^
2068 tl_cff=visc2_r(i,j)*0.5_r8* &
2069 & (tl_drhs(i,j)* &
2070 & (pmon_r(i,j)* &
2071 & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,kstp)- &
2072 & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,kstp))- &
2073 & pnom_r(i,j)* &
2074 & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,kstp)- &
2075 & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,kstp)))+ &
2076 & drhs(i,j)* &
2077 & (pmon_r(i,j)* &
2078 & ((pn(i ,j)+pn(i+1,j))*tl_ubar(i+1,j,kstp)- &
2079 & (pn(i-1,j)+pn(i ,j))*tl_ubar(i ,j,kstp))- &
2080 & pnom_r(i,j)* &
2081 & ((pm(i,j )+pm(i,j+1))*tl_vbar(i,j+1,kstp)- &
2082 & (pm(i,j-1)+pm(i,j ))*tl_vbar(i,j ,kstp))))
2083!^ UFx(i,j)=on_r(i,j)*on_r(i,j)*cff
2084!^
2085 tl_ufx(i,j)=on_r(i,j)*on_r(i,j)*tl_cff
2086!^ VFe(i,j)=om_r(i,j)*om_r(i,j)*cff
2087!^
2088 tl_vfe(i,j)=om_r(i,j)*om_r(i,j)*tl_cff
2089 END DO
2090 END DO
2091!
2092 DO j=jstr,jend+1
2093 DO i=istr,iend+1
2094!^ cff=visc2_p(i,j)*Drhs_p(i,j)*0.5_r8* &
2095!^ & (pmon_p(i,j)* &
2096!^ & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,kstp)- &
2097!^ & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,kstp))+ &
2098!^ & pnom_p(i,j)* &
2099!^ & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,kstp)- &
2100!^ & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,kstp)))
2101!^
2102 tl_cff=visc2_p(i,j)*0.5_r8* &
2103 & (tl_drhs_p(i,j)* &
2104 & (pmon_p(i,j)* &
2105 & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,kstp)- &
2106 & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,kstp))+ &
2107 & pnom_p(i,j)* &
2108 & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,kstp)- &
2109 & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,kstp)))+ &
2110 & drhs_p(i,j)* &
2111 & (pmon_p(i,j)* &
2112 & ((pn(i ,j-1)+pn(i ,j))*tl_vbar(i ,j,kstp)- &
2113 & (pn(i-1,j-1)+pn(i-1,j))*tl_vbar(i-1,j,kstp))+ &
2114 & pnom_p(i,j)* &
2115 & ((pm(i-1,j )+pm(i,j ))*tl_ubar(i,j ,kstp)- &
2116 & (pm(i-1,j-1)+pm(i,j-1))*tl_ubar(i,j-1,kstp))))
2117# ifdef MASKING
2118!^ cff=cff*pmask(i,j)
2119!^
2120 tl_cff=tl_cff*pmask(i,j)
2121# endif
2122# ifdef WET_DRY_NOT_YET
2123!^ cff=cff*pmask_wet(i,j)
2124!^
2125 tl_cff=tl_cff*pmask_wet(i,j)
2126# endif
2127!^ UFe(i,j)=om_p(i,j)*om_p(i,j)*cff
2128!^
2129 tl_ufe(i,j)=om_p(i,j)*om_p(i,j)*tl_cff
2130!^ VFx(i,j)=on_p(i,j)*on_p(i,j)*cff
2131!^
2132 tl_vfx(i,j)=on_p(i,j)*on_p(i,j)*tl_cff
2133 END DO
2134 END DO
2135!
2136! Add in harmonic viscosity.
2137!
2138 DO j=jstr,jend
2139 DO i=istr,iend
2140 IF (i.ge.istru) THEN
2141!^ cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(UFx(i,j )-UFx(i-1,j))
2142!^
2143 tl_cff1=0.5_r8*(pn(i-1,j)+pn(i,j))* &
2144 & (tl_ufx(i,j )-tl_ufx(i-1,j))
2145!^ cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(UFe(i,j+1)-UFe(i ,j))
2146!^
2147 tl_cff2=0.5_r8*(pm(i-1,j)+pm(i,j))* &
2148 & (tl_ufe(i,j+1)-tl_ufe(i ,j))
2149!^ fac1=cff1+cff2
2150!^
2151 tl_fac1=tl_cff1+tl_cff2
2152!^ rubar(i,j)=rubar(i,j)+fac1
2153!^
2154 tl_rubar(i,j)=tl_rubar(i,j)+tl_fac1
2155# if defined DIAGNOSTICS_UV
2156!! DiaU2rhs(i,j,M2hvis)=fac1
2157!! DiaU2rhs(i,j,M2xvis)=cff1
2158!! DiaU2rhs(i,j,M2yvis)=cff2
2159# endif
2160 END IF
2161!
2162 IF (j.ge.jstrv) THEN
2163!^ cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(VFx(i+1,j)-VFx(i,j ))
2164!^
2165 tl_cff1=0.5_r8*(pn(i,j-1)+pn(i,j))* &
2166 & (tl_vfx(i+1,j)-tl_vfx(i,j ))
2167!^ cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(VFe(i ,j)-VFe(i,j-1))
2168!^
2169 tl_cff2=0.5_r8*(pm(i,j-1)+pm(i,j))* &
2170 & (tl_vfe(i ,j)-tl_vfe(i,j-1))
2171!^ fac1=cff1-cff2
2172!^
2173 tl_fac1=tl_cff1-tl_cff2
2174!^ rvbar(i,j)=rvbar(i,j)+fac1
2175!^
2176 tl_rvbar(i,j)=tl_rvbar(i,j)+tl_fac1
2177# if defined DIAGNOSTICS_UV
2178!! DiaV2rhs(i,j,M2hvis)=fac1
2179!! DiaV2rhs(i,j,M2xvis)= cff1
2180!! DiaV2rhs(i,j,M2yvis)=-cff2
2181# endif
2182 END IF
2183 END DO
2184 END DO
2185#endif
2186
2187#ifdef SOLVE3D
2188!
2189!-----------------------------------------------------------------------
2190! Coupling between 2D and 3D equations.
2191!-----------------------------------------------------------------------
2192!
2193! Before the first barotropic time step, arrays "rufrc" and "rvfrc"
2194! contain vertical integrals of the 3D right-hand-side terms for the
2195! momentum equations (including surface and bottom stresses). During
2196! the first barotropic time step, convert them into forcing terms by
2197! subtracting the fast-time "rubar" and "rvbar" from them.
2198!
2199! In the predictor-coupled mode, the resultant forcing terms "rufrc"
2200! and "rvfrc" are extrapolated forward in time, so they become
2201! centered effectively at time n+1/2. This is done using optimized
2202! Adams-Bashforth weights. In the code below, rufrc_bak(:,:,nstp) is
2203! at (n-1)time step, while rufrc_bak(:,:,3-nstp) is at (n-2). After
2204! its use as input, the latter is overwritten by the value at time
2205! step "nstp" (mathematically "n") during the next step.
2206!
2207! From now on, the computed forcing terms "rufrc" and "rvfrc" will
2208! remain constant during the fast-time stepping and will be added
2209! to "rubar" and "rvbar" during all subsequent barotropic steps.
2210!
2211 coupled_step : IF (first_2d_step) THEN
2212!
2213! Predictor coupled barotropic mode: Set coefficients for AB3-like
2214! forward-in-time extrapolation of 3D to 2D forcing terms "rufrc" and
2215! "rvfrc".
2216!
2217 IF (iic(ng).eq.ntstart(ng)) THEN
2218 cfwd0=1.0_r8
2219 cfwd1=0.0_r8
2220 cfwd2=0.0_r8
2221 ELSE IF (iic(ng).eq.ntstart(ng)+1) THEN
2222 cfwd0=1.5_r8
2223 cfwd1=-0.5_r8
2224 cfwd2=0.0_r8
2225 ELSE
2226 cfwd2=0.281105_r8
2227 cfwd1=-0.5_r8-2.0_r8*cfwd2
2228 cfwd0=1.5_r8+cfwd2
2229 END IF
2230!
2231 DO j=jstr,jend
2232 DO i=istr,iend
2233!
2234! Compensate for (cancel) bottom drag terms: at input into step2d
2235! "rufrc" and "rvfrc" contain bottom drag terms computed by 3D mode.
2236! However, there are no 2D counterparts in "rubar" and "rvbar" because
2237! 2D bottom drag will be computed implicitly during the final stage of
2238! updating ubar(:,:,knew) and vbar(:,:,knew) below. Note that unlike
2239! the other terms, the bottom drag should not be extrapolated forward,
2240! if "rufrc" and "rvfrc" are, so this cancelation needs to be done
2241! right now rather than at the bottom of this loop.
2242!
2243 IF (i.ge.istru) THEN
2244!^ rufrc(i,j)=rufrc(i,j)+ &
2245!^ & 0.5_r8*(rdrag(i,j)+rdrag(i-1,j))* &
2246!^ & om_u(i,j)*on_u(i,j)*ubar(i,j,kstp)
2247!^
2248 tl_rufrc(i,j)=tl_rufrc(i,j)+ &
2249 & 0.5_r8*(rdrag(i,j)+rdrag(i-1,j))* &
2250 & om_u(i,j)*on_u(i,j)*tl_ubar(i,j,kstp)
2251 END IF
2252!
2253 IF (j.ge.jstrv) THEN
2254!^ rvfrc(i,j)=rvfrc(i,j)+ &
2255!^ & 0.5_r8*(rdrag(i,j)+rdrag(i,j-1))* &
2256!^ & om_v(i,j)*on_v(i,j)*vbar(i,j,kstp)
2257!^
2258 tl_rvfrc(i,j)=tl_rvfrc(i,j)+ &
2259 & 0.5_r8*(rdrag(i,j)+rdrag(i,j-1))* &
2260 & om_v(i,j)*on_v(i,j)*tl_vbar(i,j,kstp)
2261 END IF
2262!
2263! Barotropic mode running predictor stage: forward extrapolation.
2264!
2265 IF (i.ge.istru) THEN
2266!^ cff1=rufrc(i,j)-rubar(i,j)
2267!^
2268 tl_cff1=tl_rufrc(i,j)-tl_rubar(i,j)
2269!^ rufrc(i,j)=cfwd0*cff1+ &
2270!^ & cfwd1*rufrc_bak(i,j, nstp)+ &
2271!^ & cfwd2*rufrc_bak(i,j,3-nstp)
2272!^
2273 tl_rufrc(i,j)=cfwd0*tl_cff1+ &
2274 & cfwd1*tl_rufrc_bak(i,j, nstp)+ &
2275 & cfwd2*tl_rufrc_bak(i,j,3-nstp)
2276!^ rufrc_bak(i,j,3-nstp)=cff1
2277!^
2278 tl_rufrc_bak(i,j,3-nstp)=tl_cff1
2279 END IF
2280!
2281 IF (j.ge.jstrv) THEN
2282!^ cff2=rvfrc(i,j)-rvbar(i,j)
2283!^
2284 tl_cff2=tl_rvfrc(i,j)-tl_rvbar(i,j)
2285!^ rvfrc(i,j)=cfwd0*cff2+ &
2286!^ & cfwd1*rvfrc_bak(i,j, nstp)+ &
2287!^ & cfwd2*rvfrc_bak(i,j,3-nstp)
2288!^
2289 tl_rvfrc(i,j)=cfwd0*tl_cff2+ &
2290 & cfwd1*tl_rvfrc_bak(i,j, nstp)+ &
2291 & cfwd2*tl_rvfrc_bak(i,j,3-nstp)
2292!^ rvfrc_bak(i,j,3-nstp)=cff2
2293!^
2294 tl_rvfrc_bak(i,j,3-nstp)=tl_cff2
2295 END IF
2296 END DO
2297 END DO
2298!
2299! Add correction term to shift pressure-gradient terms from "kstp" to
2300! "knew". That is, it converts the first 2D step from Forward-Euler
2301! to Forward-Backward (this is PGF_FB_CORRECTION mentioned above).
2302!
2303 DO j=jstrv-1,jend
2304 DO i=istru-1,iend
2305!^ zwrk(i,j)=zeta_new(i,j)-zeta(i,j,kstp)
2306!^
2307 tl_zwrk(i,j)=tl_zeta_new(i,j)-tl_zeta(i,j,kstp)
2308# if defined VAR_RHO_2D && defined SOLVE3D
2309!^ rzeta(i,j)=(1.0_r8+rhoS(i,j))*zwrk(i,j)
2310!^
2311 tl_rzeta(i,j)=(1.0_r8+rhos(i,j))*tl_zwrk(i,j)+ &
2312 & tl_rhos(i,j)*zwrk(i,j)
2313!^ rzeta2(i,j)=rzeta(i,j)*(zeta_new(i,j)+zeta(i,j,kstp))
2314!^
2315 tl_rzeta2(i,j)=tl_rzeta(i,j)* &
2316 & (zeta_new(i,j)+zeta(i,j,kstp))+ &
2317 & rzeta(i,j)* &
2318 & (tl_zeta_new(i,j)+tl_zeta(i,j,kstp))
2319!^ rzetaSA(i,j)=zwrk(i,j)*(rhoS(i,j)-rhoA(i,j))
2320!^
2321 tl_rzetasa(i,j)=tl_zwrk(i,j)* &
2322 & (rhos(i,j)-rhoa(i,j))+ &
2323 & zwrk(i,j)* &
2324 & (tl_rhos(i,j)-tl_rhoa(i,j))
2325# else
2326!^ rzeta(i,j)=zwrk(i,j)
2327!^
2328 tl_rzeta(i,j)=tl_zwrk(i,j)
2329!^ rzeta2(i,j)=zwrk(i,j)*(zeta_new(i,j)+zeta(i,j,kstp))
2330!^
2331 tl_rzeta2(i,j)=tl_zwrk(i,j)* &
2332 & (zeta_new(i,j)+zeta(i,j,kstp))+ &
2333 & zwrk(i,j)* &
2334 & (tl_zeta_new(i,j)+tl_zeta(i,j,kstp))
2335# endif
2336 END DO
2337 END DO
2338!
2339 cff1=0.5*g
2340# if defined VAR_RHO_2D && defined SOLVE3D
2341 cff2=0.333333333333_r8
2342# endif
2343 DO j=jstr,jend
2344 DO i=istr,iend
2345 IF (i.ge.istru) THEN
2346!^ rubar(i,j)=rubar(i,j)+ &
2347!^ & cff1*on_u(i,j)* &
2348!^ & ((h(i-1,j)+ &
2349!^ & h(i ,j))* &
2350!^ & (rzeta(i-1,j)- &
2351!^ & rzeta(i ,j))+ &
2352# if defined VAR_RHO_2D && defined SOLVE3D
2353!^ & (h(i-1,j)- &
2354!^ & h(i ,j))* &
2355!^ & (rzetaSA(i-1,j)+ &
2356!^ & rzetaSA(i ,j)+ &
2357!^ & cff2*(rhoA(i-1,j)- &
2358!^ & rhoA(i ,j))* &
2359!^ & (zwrk(i-1,j)- &
2360!^ & zwrk(i ,j)))+ &
2361# endif
2362!^ & (rzeta2(i-1,j)- &
2363!^ & rzeta2(i ,j)))
2364!^
2365 tl_rubar(i,j)=tl_rubar(i,j)+ &
2366 & cff1*on_u(i,j)* &
2367 & ((tl_h(i-1,j)+ &
2368 & tl_h(i ,j))* &
2369 & (rzeta(i-1,j)- &
2370 & rzeta(i ,j))+ &
2371 & (h(i-1,j)+ &
2372 & h(i ,j))* &
2373 & (tl_rzeta(i-1,j)- &
2374 & tl_rzeta(i ,j))+ &
2375# if defined VAR_RHO_2D && defined SOLVE3D
2376 & (tl_h(i-1,j)- &
2377 & tl_h(i ,j))* &
2378 & (rzetasa(i-1,j)+ &
2379 & rzetasa(i ,j)+ &
2380 & cff2*(rhoa(i-1,j)- &
2381 & rhoa(i ,j))* &
2382 & (zwrk(i-1,j)- &
2383 & zwrk(i ,j)))+ &
2384 & (h(i-1,j)- &
2385 & h(i ,j))* &
2386 & (tl_rzetasa(i-1,j)+ &
2387 & tl_rzetasa(i ,j)+ &
2388 & cff2*((tl_rhoa(i-1,j)- &
2389 & tl_rhoa(i ,j))* &
2390 & (zwrk(i-1,j)- &
2391 & zwrk(i ,j))+ &
2392 & (rhoa(i-1,j)- &
2393 & rhoa(i ,j))* &
2394 & (tl_zwrk(i-1,j)- &
2395 & tl_zwrk(i ,j))))+ &
2396# endif
2397 & (tl_rzeta2(i-1,j)- &
2398 & tl_rzeta2(i ,j)))
2399# ifdef DIAGNOSTICS_UV
2400!! DiaU2rhs(i,j,M2pgrd)=DiaU2rhs(i,j,M2pgrd)+ &
2401!! & rubar(i,j)
2402# endif
2403 END IF
2404!
2405 IF (j.ge.jstrv) THEN
2406!^ rvbar(i,j)=rvbar(i,j)+ &
2407!^ & cff1*om_v(i,j)* &
2408!^ & ((h(i,j-1)+ &
2409!^ & h(i,j ))* &
2410!^ & (rzeta(i,j-1)- &
2411!^ & rzeta(i,j ))+ &
2412# if defined VAR_RHO_2D && defined SOLVE3D
2413!^ & (h(i,j-1)- &
2414!^ & h(i,j ))* &
2415!^ & (rzetaSA(i,j-1)+ &
2416!^ & rzetaSA(i,j )+ &
2417!^ & cff2*(rhoA(i,j-1)- &
2418!^ & rhoA(i,j ))* &
2419!^ & (zwrk(i,j-1)- &
2420!^ & zwrk(i,j )))+ &
2421# endif
2422!^ & (rzeta2(i,j-1)- &
2423!^ & rzeta2(i,j )))
2424!^
2425 tl_rvbar(i,j)=tl_rvbar(i,j)+ &
2426 & cff1*om_v(i,j)* &
2427 & ((tl_h(i,j-1)+ &
2428 & tl_h(i,j ))* &
2429 & (rzeta(i,j-1)- &
2430 & rzeta(i,j ))+ &
2431 & (h(i,j-1)+ &
2432 & h(i,j ))* &
2433 & (tl_rzeta(i,j-1)- &
2434 & tl_rzeta(i,j ))+ &
2435# if defined VAR_RHO_2D && defined SOLVE3D
2436 & (tl_h(i,j-1)- &
2437 & tl_h(i,j ))* &
2438 & (rzetasa(i,j-1)+ &
2439 & rzetasa(i,j )+ &
2440 & cff2*(rhoa(i,j-1)- &
2441 & rhoa(i,j ))* &
2442 & (zwrk(i,j-1)- &
2443 & zwrk(i,j )))+ &
2444 & (h(i,j-1)- &
2445 & h(i,j ))* &
2446 & (tl_rzetasa(i,j-1)+ &
2447 & tl_rzetasa(i,j )+ &
2448 & cff2*((tl_rhoa(i,j-1)- &
2449 & tl_rhoa(i,j ))* &
2450 & (zwrk(i,j-1)- &
2451 & zwrk(i,j ))+ &
2452 & (rhoa(i,j-1)- &
2453 & rhoa(i,j ))* &
2454 & (tl_zwrk(i,j-1)- &
2455 & tl_zwrk(i,j ))))+ &
2456# endif
2457 & (tl_rzeta2(i,j-1)- &
2458 & tl_rzeta2(i,j )))
2459# ifdef DIAGNOSTICS_UV
2460!! DiaV2rhs(i,j,M2pgrd)=DiaV2rhs(i,j,M2pgrd)+ &
2461!! & rvbar(i,j)
2462# endif
2463 END IF
2464 END DO
2465 END DO
2466 END IF coupled_step
2467#endif
2468!
2469!-----------------------------------------------------------------------
2470! Time step 2D momentum equations.
2471!-----------------------------------------------------------------------
2472!
2473! Advance 2D momentum components while simultaneously adding them to
2474! accumulate fast-time averages to compute barotropic fluxes. Doing so
2475! straight away yields a more computationally dense code. However, the
2476! fast-time averaged fluxes (DU_avg1 and DV_avg1) are needed both at
2477! the interior and physical boundary points. Thus, we need separate
2478! loops along the domain boundaries after setting "ubar" and "vbar"
2479! lateral boundary conditions. Also, note that bottom drag is treated
2480! implicitly:
2481!
2482! Dnew*ubar(:,:,m+1) = Dold*ubar(:,:,m) +
2483! dtfast(ng)*rhs2D(:,:) -
2484! dtfast(ng)*rdrag(:,:)*ubar(:,:,m+1)
2485! hence
2486!
2487! ubar(:,:,m+1)=[Dold * ubar(..,m) + dtfast(ng) * rhs2D(:,:)] /
2488! [Dnew + dtfast(ng) * rdrag(:,:)]
2489!
2490! DU_avg1 = DU_avg1 +
2491! weight(m+1) * Dnew * ubar(:,:,m+1) * on_u(:,:)
2492!
2493! where it should be noted that Dnew .ne. Dnew + dtfast * rdrag
2494!
2495 DO j=jstrv-1,jend
2496 DO i=istru-1,iend
2497!^ Dnew(i,j)=h(i,j)+zeta_new(i,j)
2498!^
2499 tl_dnew(i,j)=tl_h(i,j)+tl_zeta_new(i,j)
2500!^ Dnew_rd(i,j)=Dnew(i,j)+dtfast(ng)*rdrag(i,j)
2501!^
2502 tl_dnew_rd(i,j)=tl_dnew(i,j)
2503!^ Dstp(i,j)=h(i,j)+zeta(i,j,kstp)
2504!^
2505 tl_dstp(i,j)=tl_h(i,j)+tl_zeta(i,j,kstp)
2506 END DO
2507 END DO
2508
2509#if defined UV_QDRAG && !defined SOLVE3D
2510!
2511! Add quadratic drag term associated in shallow-water applications.
2512!
2513! Here, the SQRT(3) is due to a linear interpolation with second order
2514! accuaracy that ensures positive and negative values of the velocity
2515! components:
2516!
2517! u^2(i+1/2) = (1/3)*[u(i)*u(i) + u(i)*u(i+1) + u(i+1)*u(i+1)]
2518!
2519! If u(i)=1 and u(i+1)=-1, then u^2(i+1/2)=1/3 as it should be.
2520!
2521 cff=dtfast(ng)/sqrt(3.0_r8)
2522 DO j=jstrv-1,jend
2523 DO i=istru-1,iend
2524 cff1=ubar(i ,j,kstp)**2+ &
2525 & ubar(i+1,j,kstp)**2+ &
2526 & ubar(i ,j,kstp)*ubar(i+1,j,kstp)+ &
2527 & vbar(i,j ,kstp)**2+ &
2528 & vbar(i,j+1,kstp)**2+ &
2529 & vbar(i,j ,kstp)*vbar(i,j+1,kstp)
2530 tl_cff1=2.0_r8*ubar(i ,j,kstp)*tl_ubar(i ,j,kstp)+ &
2531 & 2.0_r8*ubar(i+1,j,kstp)*tl_ubar(i+1,j,kstp)+ &
2532 & tl_ubar(i ,j,kstp)*ubar(i+1,j,kstp)+ &
2533 & tl_ubar(i+1,j,kstp)*ubar(i ,j,kstp)+ &
2534 & 2.0_r8*vbar(i,j ,kstp)*tl_vbar(i,j ,kstp)+ &
2535 & 2.0_r8*vbar(i,j+1,kstp)*tl_ vbar(i,j+1,kstp)+ &
2536 & tl_vbar(i,j ,kstp)*vbar(i,j+1,kstp)+ &
2537 & tl_vbar(i,j+1,kstp)*vbar(i,j ,kstp)
2538 cff2=sqrt(cff1)
2539 tl_cff2=0.5_r8*tl_cff1/cff2
2540!^ Dnew_rd(i,j)=Dnew_rd(i,j)+ &
2541!^ & cff*rdrag2(i,j)*cff2
2542!^
2543 tl_dnew_rd(i,j)=tl_dnew_rd(i,j)+ &
2544 & cff*rdrag2(i,j)*tl_cff2
2545 END DO
2546 END DO
2547#endif
2548!
2549! Step 2D momentum equations.
2550!
2551 cff=0.5_r8*dtfast(ng)
2552#ifdef SOLVE3D
2553 cff1=0.5_r8*weight(1,iif(ng),ng)
2554#else
2555 cff2=2.0_r8*dtfast(ng)
2556#endif
2557 DO j=jstr,jend
2558 DO i=istru,iend
2559 cff3=cff*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2560 fac1=1.0_r8/(dnew_rd(i,j)+dnew_rd(i-1,j))
2561 tl_fac1=-fac1*fac1*(tl_dnew_rd(i,j)+tl_dnew_rd(i-1,j))
2562!^ ubar(i,j,knew)=fac1* &
2563!^ & ((Dstp(i,j)+Dstp(i-1,j))*ubar(i,j,kstp)+ &
2564#ifdef SOLVE3D
2565!^ & cff3*(rubar(i,j)+rufrc(i,j)))
2566#else
2567!^ & cff3*rubar(i,j)+cff2*sustr(i,j))
2568#endif
2569!^
2570 tl_ubar(i,j,knew)=tl_fac1* &
2571 & ((dstp(i,j)+dstp(i-1,j))*ubar(i,j,kstp)+ &
2572#ifdef SOLVE3D
2573 & cff3*(rubar(i,j)+rufrc(i,j)))+ &
2574#else
2575 & cff3*rubar(i,j)+cff2*sustr(i,j))+ &
2576#endif
2577 & fac1* &
2578 & ((dstp(i,j)+dstp(i-1,j))* &
2579 & tl_ubar(i,j,kstp)+ &
2580 & (tl_dstp(i,j)+tl_dstp(i-1,j))* &
2581 & ubar(i,j,kstp)+ &
2582#ifdef SOLVE3D
2583 & cff3*(tl_rubar(i,j)+tl_rufrc(i,j)))
2584#else
2585 & cff3*tl_rubar(i,j)+cff2*tl_sustr(i,j))
2586#endif
2587#ifdef MASKING
2588!^ ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2589!^
2590 tl_ubar(i,j,knew)=tl_ubar(i,j,knew)*umask(i,j)
2591#endif
2592#ifdef WET_DRY_NOT_YET
2593!^ cff5=ABS(ABS(umask_wet(i,j))-1.0_r8)
2594!^ cff6=0.5_r8+DSIGN(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2595!^ cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2596!^ ubar(i,j,knew)=ubar(i,j,knew)*cff7
2597!^
2598!^ HGA: TLM code needed here.
2599!^
2600#endif
2601#ifdef SOLVE3D
2602!^ DU_avg1(i,j)=DU_avg1(i,j)+ &
2603!^ & cff1*on_u(i,j)* &
2604!^ & (Dnew(i,j)+Dnew(i-1,j))*ubar(i,j,knew)
2605!^
2606 tl_du_avg1(i,j)=tl_du_avg1(i,j)+ &
2607 & cff1*on_u(i,j)* &
2608 & ((dnew(i,j)+dnew(i-1,j))* &
2609 & tl_ubar(i,j,knew)+ &
2610 & (tl_dnew(i,j)+tl_dnew(i-1,j))* &
2611 & ubar(i,j,knew))
2612#endif
2613#if defined NESTING && !defined SOLVE3D
2614!^ DU_flux(i,j)=0.5_r8*on_u(i,j)* &
2615!^ & (Dnew(i,j)+Dnew(i-1,j))*ubar(i,j,knew)
2616!^
2617 tl_du_flux(i,j)=0.5_r8*on_u(i,j)* &
2618 & ((dnew(i,j)+dnew(i-1,j))* &
2619 & tl_ubar(i,j,knew)+ &
2620 & (tl_dnew(i,j)+tl_dnew(i-1,j))* &
2621 & ubar(i,j,knew))
2622#endif
2623 END DO
2624 END DO
2625!
2626 DO j=jstrv,jend
2627 DO i=istr,iend
2628 cff3=cff*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2629 fac2=1.0_r8/(dnew_rd(i,j)+dnew_rd(i,j-1))
2630 tl_fac2=-fac2*fac2*(tl_dnew_rd(i,j)+tl_dnew_rd(i,j-1))
2631!^ vbar(i,j,knew)=fac2* &
2632!^ & ((Dstp(i,j)+Dstp(i,j-1))*vbar(i,j,kstp)+ &
2633#ifdef SOLVE3D
2634!^ & cff3*(rvbar(i,j)+rvfrc(i,j)))
2635#else
2636!^ & cff3*rvbar(i,j)+cff2*svstr(i,j))
2637#endif
2638!^
2639 tl_vbar(i,j,knew)=tl_fac2* &
2640 & ((dstp(i,j)+dstp(i,j-1))*vbar(i,j,kstp)+ &
2641#ifdef SOLVE3D
2642 & cff3*(rvbar(i,j)+rvfrc(i,j)))+ &
2643#else
2644 & cff3*rvbar(i,j)+cff2*svstr(i,j))+ &
2645#endif
2646 & fac2* &
2647 & ((dstp(i,j)+dstp(i,j-1))* &
2648 & tl_vbar(i,j,kstp)+ &
2649 & (tl_dstp(i,j)+tl_dstp(i,j-1))* &
2650 & vbar(i,j,kstp)+ &
2651#ifdef SOLVE3D
2652 & cff3*(tl_rvbar(i,j)+tl_rvfrc(i,j)))
2653#else
2654 & cff3*tl_rvbar(i,j)+cff2*tl_svstr(i,j))
2655#endif
2656#ifdef MASKING
2657!^ vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2658!^
2659 tl_vbar(i,j,knew)=tl_vbar(i,j,knew)*vmask(i,j)
2660#endif
2661#ifdef WET_DRY_NOT_YET
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!^
2667!^ HGA: TLM code needed here.
2668!^
2669#endif
2670#ifdef SOLVE3D
2671!^ DV_avg1(i,j)=DV_avg1(i,j)+ &
2672!^ & cff1*om_v(i,j)* &
2673!^ & (Dnew(i,j)+Dnew(i,j-1))*vbar(i,j,knew)
2674!^
2675 tl_dv_avg1(i,j)=tl_dv_avg1(i,j)+ &
2676 & cff1*om_v(i,j)* &
2677 & ((dnew(i,j)+dnew(i,j-1))* &
2678 & tl_vbar(i,j,knew)+ &
2679 & (tl_dnew(i,j)+tl_dnew(i,j-1))* &
2680 & vbar(i,j,knew))
2681#endif
2682#if defined NESTING && !defined SOLVE3D
2683!^ DV_flux(i,j)=0.5_r8*om_v(i,j)* &
2684!^ & (Dnew(i,j)+Dnew(i,j-1))*vbar(i,j,knew)
2685!^
2686 tl_dv_flux(i,j)=0.5_r8*om_v(i,j)* &
2687 & ((dnew(i,j)+dnew(i,j-1))* &
2688 & tl_vbar(i,j,knew)+ &
2689 & (tl_dnew(i,j)+tl_dnew(i,j-1))* &
2690 & vbar(i,j,knew))
2691#endif
2692 END DO
2693 END DO
2694!
2695! Apply lateral boundary conditions.
2696!
2697!^ CALL u2dbc_tile (ng, tile, &
2698!^ & LBi, UBi, LBj, UBj, &
2699!^ & IminS, ImaxS, JminS, JmaxS, &
2700!^ & krhs, kstp, knew, &
2701!^ & ubar, vbar, zeta)
2702!^
2703 CALL tl_u2dbc_tile (ng, tile, &
2704 & lbi, ubi, lbj, ubj, &
2705 & imins, imaxs, jmins, jmaxs, &
2706 & krhs, kstp, knew, &
2707 & ubar, vbar, zeta, &
2708 & tl_ubar, tl_vbar, tl_zeta)
2709!^ CALL v2dbc_tile (ng, tile, &
2710!^ & LBi, UBi, LBj, UBj, &
2711!^ & IminS, ImaxS, JminS, JmaxS, &
2712!^ & krhs, kstp, knew, &
2713!^ & ubar, vbar, zeta)
2714!^
2715 CALL tl_v2dbc_tile (ng, tile, &
2716 & lbi, ubi, lbj, ubj, &
2717 & imins, imaxs, jmins, jmaxs, &
2718 & krhs, kstp, knew, &
2719 & ubar, vbar, zeta, &
2720 & tl_ubar, tl_vbar, tl_zeta)
2721!
2722! Compute integral mass flux across open boundaries and adjust
2723! for volume conservation.
2724!
2725 IF (any(volcons(:,ng))) THEN
2726!^ CALL obc_flux_tile (ng, tile, &
2727!^ & LBi, UBi, LBj, UBj, &
2728!^ & IminS, ImaxS, JminS, JmaxS, &
2729!^ & knew, &
2730#ifdef MASKING
2731!^ & umask, vmask, &
2732#endif
2733!^ & h, om_v, on_u, &
2734!^ & ubar, vbar, zeta)
2735!^
2736 CALL tl_obc_flux_tile (ng, tile, &
2737 & lbi, ubi, lbj, ubj, &
2738 & imins, imaxs, jmins, jmaxs, &
2739 & knew, &
2740#ifdef MASKING
2741 & umask, vmask, &
2742#endif
2743 & h, tl_h, om_v, on_u, &
2744 & ubar, vbar, zeta, &
2745 & tl_ubar, tl_vbar, tl_zeta)
2746 END IF
2747
2748#if defined SOLVE3D || (defined NESTING && !defined SOLVE3D)
2749!
2750! Set barotropic fluxes along physical boundaries.
2751!
2752 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
2753 IF (domain(ng)%Western_Edge(tile)) THEN
2754 DO j=jstr-1,jendr
2755!^ Dnew(Istr-1,j)=h(Istr-1,j)+zeta_new(Istr-1,j)
2756!^
2757 tl_dnew(istr-1,j)=tl_h(istr-1,j)+tl_zeta_new(istr-1,j)
2758 END DO
2759 END IF
2760 END IF
2761 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
2762 IF (domain(ng)%Eastern_Edge(tile)) THEN
2763 DO j=jstr-1,jendr
2764!^ Dnew(Iend+1,j)=h(Iend+1,j)+zeta_new(Iend+1,j)
2765!^
2766 tl_dnew(iend+1,j)=tl_h(iend+1,j)+tl_zeta_new(iend+1,j)
2767 END DO
2768 END IF
2769 END IF
2770 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
2771 IF (domain(ng)%Southern_Edge(tile)) THEN
2772 DO i=istr-1,iendr
2773!^ Dnew(i,Jstr-1)=h(i,Jstr-1)+zeta_new(i,Jstr-1)
2774!^
2775 tl_dnew(i,jstr-1)=tl_h(i,jstr-1)+tl_zeta_new(i,jstr-1)
2776 END DO
2777 END IF
2778 END IF
2779 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
2780 IF (domain(ng)%Northern_Edge(tile)) THEN
2781 DO i=istr-1,iendr
2782!^ Dnew(i,Jend+1)=h(i,Jend+1)+zeta_new(i,Jend+1)
2783!^
2784 tl_dnew(i,jend+1)=tl_h(i,jend+1)+tl_zeta_new(i,jend+1)
2785 END DO
2786 END IF
2787 END IF
2788
2789# ifdef SOLVE3D
2790!
2791 cff1=0.5*weight(1,iif(ng),ng)
2792# endif
2793!
2794 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
2795 IF (domain(ng)%Western_Edge(tile)) THEN
2796 DO j=jstrr,jendr
2797# if defined NESTING && !defined SOLVE3D
2798!^ DU_flux(IstrU-1,j)=0.5_r8*on_u(IstrU-1,j)* &
2799!^ & (Dnew(IstrU-1,j)+Dnew(IstrU-2,j))* &
2800!^ & ubar(IstrU-1,j,knew)
2801!^
2802 tl_du_flux(istru-1,j)=0.5_r8*on_u(istru-1,j)* &
2803 & ((dnew(istru-1,j)+ &
2804 & dnew(istru-2,j))* &
2805 & tl_ubar(istru-1,j,knew)+ &
2806 & (tl_dnew(istru-1,j)+ &
2807 & tl_dnew(istru-2,j))* &
2808 & ubar(istru-1,j,knew))
2809# else
2810!^ DU_avg1(IstrU-1,j)=DU_avg1(IstrU-1,j)+ &
2811!^ & cff1*on_u(IstrU-1,j)* &
2812!^ & (Dnew(IstrU-1,j)+Dnew(IstrU-2,j))* &
2813!^ & ubar(IstrU-1,j,knew)
2814!^
2815 tl_du_avg1(istru-1,j)=tl_du_avg1(istru-1,j)+ &
2816 & cff1*on_u(istru-1,j)* &
2817 & ((dnew(istru-1,j)+ &
2818 & dnew(istru-2,j))* &
2819 & tl_ubar(istru-1,j,knew)+ &
2820 & (tl_dnew(istru-1,j)+ &
2821 & tl_dnew(istru-2,j))* &
2822 & ubar(istru-1,j,knew))
2823# endif
2824 END DO
2825 DO j=jstrv,jend
2826# if defined NESTING && !defined SOLVE3D
2827!^ DV_flux(Istr-1,j)=0.5_r8*om_v(Istr-1,j)* &
2828!^ & (Dnew(Istr-1,j)+Dnew(Istr-1,j-1))* &
2829!^ & vbar(Istr-1,j,knew)
2830!^
2831 tl_dv_flux(istr-1,j)=0.5_r8*om_v(istr-1,j)* &
2832 & ((dnew(istr-1,j )+ &
2833 & dnew(istr-1,j-1))* &
2834 & tl_vbar(istr-1,j,knew)+ &
2835 & (tl_dnew(istr-1,j )+ &
2836 & tl_dnew(istr-1,j-1))* &
2837 & vbar(istr-1,j,knew))
2838# else
2839!^ DV_avg1(Istr-1,j)=DV_avg1(Istr-1,j)+ &
2840!^ & cff1*om_v(Istr-1,j)* &
2841!^ & (Dnew(Istr-1,j)+Dnew(Istr-1,j-1))* &
2842!^ & vbar(Istr-1,j,knew)
2843!^
2844 tl_dv_avg1(istr-1,j)=tl_dv_avg1(istr-1,j)+ &
2845 & cff1*om_v(istr-1,j)* &
2846 & ((dnew(istr-1,j )+ &
2847 & dnew(istr-1,j-1))* &
2848 & tl_vbar(istr-1,j,knew)+ &
2849 & (tl_dnew(istr-1,j )+ &
2850 & tl_dnew(istr-1,j-1))* &
2851 & vbar(istr-1,j,knew))
2852# endif
2853 END DO
2854 END IF
2855 END IF
2856 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
2857 IF (domain(ng)%Eastern_Edge(tile)) THEN
2858 DO j=jstrr,jendr
2859# if defined NESTING && !defined SOLVE3D
2860!^ DU_flux(Iend+1,j)=0.5_r8*on_u(Iend+1,j)* &
2861!^ & (Dnew(Iend+1,j)+Dnew(Iend,j))* &
2862!^ & ubar(Iend+1,j,knew)
2863!^
2864 tl_du_flux(iend+1,j)=0.5_r8*on_u(iend+1,j)* &
2865 & ((dnew(iend+1,j)+ &
2866 & dnew(iend ,j))* &
2867 & tl_ubar(iend+1,j,knew)+ &
2868 & (tl_dnew(iend+1,j)+ &
2869 & tl_dnew(iend ,j))* &
2870 & ubar(iend+1,j,knew))
2871# else
2872!^ DU_avg1(Iend+1,j)=DU_avg1(Iend+1,j)+ &
2873!^ & cff1*on_u(Iend+1,j)* &
2874!^ & (Dnew(Iend+1,j)+Dnew(Iend,j))* &
2875!^ & ubar(Iend+1,j,knew)
2876!^
2877 tl_du_avg1(iend+1,j)=tl_du_avg1(iend+1,j)+ &
2878 & cff1*on_u(iend+1,j)* &
2879 & ((dnew(iend+1,j)+ &
2880 & dnew(iend ,j))* &
2881 & tl_ubar(iend+1,j,knew)+ &
2882 & (tl_dnew(iend+1,j)+ &
2883 & tl_dnew(iend ,j))* &
2884 & ubar(iend+1,j,knew))
2885# endif
2886 END DO
2887 DO j=jstrv,jend
2888# if defined NESTING && !defined SOLVE3D
2889!^ DV_flux(Iend+1,j)=0.5_r8*om_v(Iend+1,j)* &
2890!^ & (Dnew(Iend+1,j)+Dnew(Iend+1,j-1))* &
2891!^ & vbar(Iend+1,j,knew)
2892!^
2893 tl_dv_flux(iend+1,j)=0.5_r8*om_v(iend+1,j)* &
2894 & ((dnew(iend+1,j )+ &
2895 & dnew(iend+1,j-1))* &
2896 & tl_vbar(iend+1,j,knew)+ &
2897 & (tl_dnew(iend+1,j )+ &
2898 & tl_dnew(iend+1,j-1))* &
2899 & vbar(iend+1,j,knew))
2900# else
2901!^ DV_avg1(Iend+1,j)=DV_avg1(Iend+1,j)+ &
2902!^ & cff1*om_v(Iend+1,j)* &
2903!^ & (Dnew(Iend+1,j)+Dnew(Iend+1,j-1))* &
2904!^ & vbar(Iend+1,j,knew)
2905!^
2906 tl_dv_avg1(iend+1,j)=tl_dv_avg1(iend+1,j)+ &
2907 & cff1*om_v(iend+1,j)* &
2908 & ((dnew(iend+1,j )+ &
2909 & dnew(iend+1,j-1))* &
2910 & tl_vbar(iend+1,j,knew)+ &
2911 & (tl_dnew(iend+1,j )+ &
2912 & tl_dnew(iend+1,j-1))* &
2913 & vbar(iend+1,j,knew))
2914# endif
2915 END DO
2916 END IF
2917 END IF
2918 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
2919 IF (domain(ng)%Southern_Edge(tile)) THEN
2920 DO i=istru,iend
2921# if defined NESTING && !defined SOLVE3D
2922!^ DU_flux(i,Jstr-1)=0.5_r8*on_u(i,Jstr-1)* &
2923!^ & (Dnew(i,Jstr-1)+Dnew(i-1,Jstr-1))* &
2924!^ & ubar(i,Jstr-1,knew)
2925!^
2926 tl_du_flux(i,jstr-1)=0.5_r8*on_u(i,jstr-1)* &
2927 & ((dnew(i ,jstr-1)+ &
2928 & dnew(i-1,jstr-1))* &
2929 & tl_ubar(i,jstr-1,knew)+ &
2930 & (tl_dnew(i ,jstr-1)+ &
2931 & tl_dnew(i-1,jstr-1))* &
2932 & ubar(i,jstr-1,knew))
2933# else
2934!^ DU_avg1(i,Jstr-1)=DU_avg1(i,Jstr-1)+ &
2935!^ & cff1*on_u(i,Jstr-1)* &
2936!^ & (Dnew(i,Jstr-1)+Dnew(i-1,Jstr-1))* &
2937!^ & ubar(i,Jstr-1,knew)
2938!^
2939 tl_du_avg1(i,jstr-1)=tl_du_avg1(i,jstr-1)+ &
2940 & cff1*on_u(i,jstr-1)* &
2941 & ((dnew(i ,jstr-1)+ &
2942 & dnew(i-1,jstr-1))* &
2943 & tl_ubar(i,jstr-1,knew)+ &
2944 & (tl_dnew(i ,jstr-1)+ &
2945 & tl_dnew(i-1,jstr-1))* &
2946 & ubar(i,jstr-1,knew))
2947# endif
2948 END DO
2949 DO i=istrr,iendr
2950# if defined NESTING && !defined SOLVE3D
2951!^ DV_flux(i,JstrV-1)=0.5_r8*om_v(i,JstrV-1)* &
2952!^ & (Dnew(i,JstrV-1)+Dnew(i,JstrV-2))* &
2953!^ & vbar(i,JstrV-1,knew)
2954!^
2955 tl_dv_flux(i,jstrv-1)=0.5_r8*om_v(i,jstrv-1)* &
2956 & ((dnew(i,jstrv-1)+ &
2957 & dnew(i,jstrv-2))* &
2958 & tl_vbar(i,jstrv-1,knew)+ &
2959 & (tl_dnew(i,jstrv-1)+ &
2960 & tl_dnew(i,jstrv-2))* &
2961 & vbar(i,jstrv-1,knew))
2962# else
2963!^ DV_avg1(i,JstrV-1)=DV_avg1(i,JstrV-1)+ &
2964!^ & cff1*om_v(i,JstrV-1)* &
2965!^ & (Dnew(i,JstrV-1)+Dnew(i,JstrV-2))* &
2966!^ & vbar(i,JstrV-1,knew)
2967!^
2968 tl_dv_avg1(i,jstrv-1)=tl_dv_avg1(i,jstrv-1)+ &
2969 & cff1*om_v(i,jstrv-1)* &
2970 & ((dnew(i,jstrv-1)+ &
2971 & dnew(i,jstrv-2))* &
2972 & tl_vbar(i,jstrv-1,knew)+ &
2973 & (tl_dnew(i,jstrv-1)+ &
2974 & tl_dnew(i,jstrv-2))* &
2975 & vbar(i,jstrv-1,knew))
2976# endif
2977 END DO
2978 END IF
2979 END IF
2980 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
2981 IF (domain(ng)%Northern_Edge(tile)) THEN
2982 DO i=istru,iend
2983# if defined NESTING && !defined SOLVE3D
2984!^ DU_flux(i,Jend+1)=0.5_r8*on_u(i,Jend+1)* &
2985!^ & (Dnew(i,Jend+1)+Dnew(i-1,Jend+1))* &
2986!^ & ubar(i,Jend+1,knew)
2987!^
2988 tl_du_flux(i,jend+1)=0.5_r8*on_u(i,jend+1)* &
2989 & ((dnew(i ,jend+1)+ &
2990 & dnew(i-1,jend+1))* &
2991 & tl_ubar(i,jend+1,knew)+ &
2992 & (tl_dnew(i ,jend+1)+ &
2993 & tl_dnew(i-1,jend+1))* &
2994 & ubar(i,jend+1,knew))
2995# else
2996!^ DU_avg1(i,Jend+1)=DU_avg1(i,Jend+1)+ &
2997!^ & cff1*on_u(i,Jend+1)* &
2998!^ & (Dnew(i,Jend+1)+Dnew(i-1,Jend+1))* &
2999!^ & ubar(i,Jend+1,knew)
3000!^
3001 tl_du_avg1(i,jend+1)=tl_du_avg1(i,jend+1)+ &
3002 & cff1*on_u(i,jend+1)* &
3003 & ((dnew(i ,jend+1)+ &
3004 & dnew(i-1,jend+1))* &
3005 & tl_ubar(i,jend+1,knew)+ &
3006 & (tl_dnew(i ,jend+1)+ &
3007 & tl_dnew(i-1,jend+1))* &
3008 & ubar(i,jend+1,knew))
3009# endif
3010 END DO
3011 DO i=istrr,iendr
3012# if defined NESTING && !defined SOLVE3D
3013!^ DV_flux(i,Jend+1)=0.5_r8*om_v(i,Jend+1)* &
3014!^ & (Dnew(i,Jend+1)+Dnew(i,Jend))* &
3015!^ & vbar(i,Jend+1,knew)
3016!^
3017 tl_dv_flux(i,jend+1)=0.5_r8*om_v(i,jend+1)* &
3018 & ((dnew(i,jend+1)+ &
3019 & dnew(i,jend ))* &
3020 & tl_vbar(i,jend+1,knew)+ &
3021 & (tl_dnew(i,jend+1)+ &
3022 & tl_dnew(i,jend ))* &
3023 & vbar(i,jend+1,knew))
3024# else
3025!^ DV_avg1(i,Jend+1)=DV_avg1(i,Jend+1)+ &
3026!^ & cff1*om_v(i,Jend+1)* &
3027!^ & (Dnew(i,Jend+1)+Dnew(i,Jend))* &
3028!^ & vbar(i,Jend+1,knew)
3029!^
3030 tl_dv_avg1(i,jend+1)=tl_dv_avg1(i,jend+1)+ &
3031 & cff1*om_v(i,jend+1)* &
3032 & ((dnew(i,jend+1)+ &
3033 & dnew(i,jend ))* &
3034 & tl_vbar(i,jend+1,knew)+ &
3035 & (tl_dnew(i,jend+1)+ &
3036 & tl_dnew(i,jend ))* &
3037 & vbar(i,jend+1,knew))
3038# endif
3039 END DO
3040 END IF
3041 END IF
3042#endif
3043!
3044!-----------------------------------------------------------------------
3045! Apply momentum transport point sources (like river runoff), if any.
3046!
3047! Dsrc(is) = 0, flow across grid cell u-face (positive or negative)
3048! Dsrc(is) = 1, flow across grid cell v-face (positive or negative)
3049!-----------------------------------------------------------------------
3050!
3051 IF (luvsrc(ng)) THEN
3052 DO is=1,nsrc(ng)
3053 i=sources(ng)%Isrc(is)
3054 j=sources(ng)%Jsrc(is)
3055 IF (((istrr.le.i).and.(i.le.iendr)).and. &
3056 & ((jstrr.le.j).and.(j.le.jendr))) THEN
3057 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
3058 cff=1.0_r8/(on_u(i,j)* &
3059 & 0.5_r8*(dnew(i-1,j)+dnew(i,j)))
3060 tl_cff=-cff*cff*on_u(i,j)* &
3061 & 0.5_r8*(tl_dnew(i-1,j)+tl_dnew(i ,j))
3062!^ ubar(i,j,knew)=SOURCES(ng)%Qbar(is)*cff
3063!^
3064 tl_ubar(i,j,knew)=sources(ng)%tl_Qbar(is)*cff+ &
3065 & sources(ng)%Qbar(is)*tl_cff
3066#ifdef SOLVE3D
3067!^ DU_avg1(i,j)=SOURCES(ng)%Qbar(is)
3068!^
3069 tl_du_avg1(i,j)=sources(ng)%tl_Qbar(is)
3070#endif
3071#if defined NESTING && !defined SOLVE3D
3072!^ DU_flux(i,j)=SOURCES(ng)%Qbar(is)
3073!^
3074 tl_du_flux(i,j)=sources(ng)%tl_Qbar(is)
3075#endif
3076 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
3077 cff=1.0_r8/(om_v(i,j)* &
3078 & 0.5_r8*(dnew(i,j-1)+dnew(i,j)))
3079 tl_cff=-cff*cff*om_v(i,j)* &
3080 & 0.5_r8*(tl_dnew(i,j-1)+tl_dnew(i,j))
3081!^ vbar(i,j,knew)=SOURCES(ng)%Qbar(is)*cff
3082!^
3083 tl_vbar(i,j,knew)=sources(ng)%tl_Qbar(is)*cff+ &
3084 & sources(ng)%Qbar(is)*tl_cff
3085#ifdef SOLVE3D
3086!^ DV_avg1(i,j)=SOURCES(ng)%Qbar(is)
3087!^
3088 tl_dv_avg1(i,j)=sources(ng)%tl_Qbar(is)
3089#endif
3090#if defined NESTING && !defined SOLVE3D
3091!^ DV_flux(i,j)=SOURCES(ng)%Qbar(is)
3092!^
3093 tl_dv_flux(i,j)=sources(ng)%tl_Qbar(is)
3094#endif
3095 END IF
3096 END IF
3097 END DO
3098 END IF
3099!
3100! Deallocate local new free-surface.
3101!
3102 deallocate ( tl_zeta_new )
3103
3104#ifdef WET_DRY_NOT_YET
3105!
3106!-----------------------------------------------------------------------
3107! Compute new wet/dry masks.
3108!-----------------------------------------------------------------------
3109!
3110!^ CALL wetdry_tile (ng, tile, &
3111!^ & LBi, UBi, LBj, UBj, &
3112!^ & IminS, ImaxS, JminS, JmaxS, &
3113# ifdef MASKING
3114!^ & pmask, rmask, umask, vmask, &
3115# endif
3116!^ & h, zeta(:,:,knew), &
3117# ifdef SOLVE3D
3118!^ & DU_avg1, DV_avg1, &
3119!^ & rmask_wet_avg, &
3120# endif
3121!^ & pmask_wet, pmask_full, &
3122!^ & rmask_wet, rmask_full, &
3123!^ & umask_wet, umask_full, &
3124!^ & vmask_wet, vmask_full)
3125!^
3126!^ HGA: Need the TLM code here.
3127!^
3128#endif
3129
3130#ifdef SOLVE3D
3131!
3132!-----------------------------------------------------------------------
3133! At the end of the last 2D time step replace the new free-surface
3134! zeta(:,:,knew) with it fast time-averaged value, Zt_avg1. Recall
3135! this is state variable is the one that communicates with the 3D
3136! kernel. Compute time-dependent depths.
3137!-----------------------------------------------------------------------
3138!
3139 IF (iif(ng).eq.nfast(ng)) THEN
3140 DO j=jstrr,jendr
3141 DO i=istrr,iendr
3142!^ zeta(i,j,knew)=Zt_avg1(i,j)
3143!^
3144 tl_zeta(i,j,knew)=tl_zt_avg1(i,j)
3145 END DO
3146 END DO
3147!^ CALL set_depth (ng, tile, iNLM)
3148!^
3149 CALL tl_set_depth (ng, tile, itlm)
3150 END IF
3151#endif
3152
3153#ifdef NESTING
3154# ifdef SOLVE3D
3155!
3156!-----------------------------------------------------------------------
3157! If nesting and after all fast time steps are completed, exchange
3158! halo information to time averaged fields.
3159!-----------------------------------------------------------------------
3160!
3161 IF (iif(ng).eq.nfast(ng)) THEN
3162!
3163! In nesting applications with refinement grids, we need to exchange
3164! the DU_avg2 and DV_avg2 fluxes boundary information for the case
3165! that a contact point is at a tile partition. Notice that in such
3166! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
3167!
3168 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
3169!^ CALL exchange_r2d_tile (ng, tile, &
3170!^ & LBi, UBi, LBj, UBj, &
3171!^ & Zt_avg1)
3172!^
3173 CALL exchange_r2d_tile (ng, tile, &
3174 & lbi, ubi, lbj, ubj, &
3175 & tl_zt_avg1)
3176!^ CALL exchange_u2d_tile (ng, tile, &
3177!^ & LBi, UBi, LBj, UBj, &
3178!^ & DU_avg1)
3179!^
3180 CALL exchange_u2d_tile (ng, tile, &
3181 & lbi, ubi, lbj, ubj, &
3182 & tl_du_avg1)
3183!^ CALL exchange_v2d_tile (ng, tile, &
3184!^ & LBi, UBi, LBj, UBj, &
3185!^ & DV_avg1)
3186!^
3187 CALL exchange_v2d_tile (ng, tile, &
3188 & lbi, ubi, lbj, ubj, &
3189 & tl_dv_avg1)
3190!^ CALL exchange_u2d_tile (ng, tile, &
3191!^ & LBi, UBi, LBj, UBj, &
3192!^ & DU_avg2)
3193!^
3194 CALL exchange_u2d_tile (ng, tile, &
3195 & lbi, ubi, lbj, ubj, &
3196 & tl_du_avg2)
3197!^ CALL exchange_v2d_tile (ng, tile, &
3198!^ & LBi, UBi, LBj, UBj, &
3199!^ & DV_avg2)
3200!^
3201 CALL exchange_v2d_tile (ng, tile, &
3202 & lbi, ubi, lbj, ubj, &
3203 & tl_dv_avg2)
3204 END IF
3205
3206# ifdef DISTRIBUTE
3207!
3208!^ CALL mp_exchange2d (ng, tile, iNLM, 3, &
3209!^ & LBi, UBi, LBj, UBj, &
3210!^ & NghostPoints, &
3211!^ & EWperiodic(ng), NSperiodic(ng), &
3212!^ & Zt_avg1, DU_avg1, DV_avg1)
3213!^
3214 CALL mp_exchange2d (ng, tile, itlm, 3, &
3215 & lbi, ubi, lbj, ubj, &
3216 & nghostpoints, &
3217 & ewperiodic(ng), nsperiodic(ng), &
3218 & tl_zt_avg1, tl_du_avg1, tl_dv_avg1)
3219!^ CALL mp_exchange2d (ng, tile, iNLM, 2, &
3220!^ & LBi, UBi, LBj, UBj, &
3221!^ & NghostPoints, &
3222!^ & EWperiodic(ng), NSperiodic(ng), &
3223!^ & DU_avg2, DV_avg2)
3224!^
3225 CALL mp_exchange2d (ng, tile, itlm, 2, &
3226 & lbi, ubi, lbj, ubj, &
3227 & nghostpoints, &
3228 & ewperiodic(ng), nsperiodic(ng), &
3229 & tl_du_avg2, tl_dv_avg2)
3230# endif
3231 END IF
3232# else
3233!
3234! In nesting applications with refinement grids, we need to exchange
3235! the DU_flux and DV_flux fluxes boundary information for the case
3236! that a contact point is at a tile partition. Notice that in such
3237! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
3238!
3239 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
3240!^ CALL exchange_u2d_tile (ng, tile, &
3241!^ & LBi, UBi, LBj, UBj, &
3242!^ & DU_flux)
3243!^
3244 CALL exchange_u2d_tile (ng, tile, &
3245 & lbi, ubi, lbj, ubj, &
3246 & tl_du_flux)
3247!^ CALL exchange_v2d_tile (ng, tile, &
3248!^ & LBi, UBi, LBj, UBj, &
3249!^ & DV_flux)
3250!^
3251 CALL exchange_v2d_tile (ng, tile, &
3252 & lbi, ubi, lbj, ubj, &
3253 & tl_dv_flux)
3254 END IF
3255
3256# ifdef DISTRIBUTE
3257!
3258!^ CALL mp_exchange2d (ng, tile, iNLM, 2, &
3259!^ & LBi, UBi, LBj, UBj, &
3260!^ & NghostPoints, &
3261!^ & EWperiodic(ng), NSperiodic(ng), &
3262!^ & DU_flux, DV_flux)
3263!^
3264 CALL mp_exchange2d (ng, tile, itlm, 2, &
3265 & lbi, ubi, lbj, ubj, &
3266 & nghostpoints, &
3267 & ewperiodic(ng), nsperiodic(ng), &
3268 & tl_du_flux, tl_dv_flux)
3269# endif
3270# endif
3271#endif
3272!
3273!-----------------------------------------------------------------------
3274! Exchange halo tile information.
3275!-----------------------------------------------------------------------
3276!
3277 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
3278!^ CALL exchange_r2d_tile (ng, tile, &
3279!^ & LBi, UBi, LBj, UBj, &
3280!^ & zeta(:,:,knew))
3281!^
3282 CALL exchange_r2d_tile (ng, tile, &
3283 & lbi, ubi, lbj, ubj, &
3284 & tl_zeta(:,:,knew))
3285!^ CALL exchange_u2d_tile (ng, tile, &
3286!^ & LBi, UBi, LBj, UBj, &
3287!^ & ubar(:,:,knew))
3288!^
3289 CALL exchange_u2d_tile (ng, tile, &
3290 & lbi, ubi, lbj, ubj, &
3291 & tl_ubar(:,:,knew))
3292!^ CALL exchange_v2d_tile (ng, tile, &
3293!^ & LBi, UBi, LBj, UBj, &
3294!^ & vbar(:,:,knew))
3295!^
3296 CALL exchange_v2d_tile (ng, tile, &
3297 & lbi, ubi, lbj, ubj, &
3298 & tl_vbar(:,:,knew))
3299 END IF
3300
3301#ifdef DISTRIBUTE
3302!
3303!^ CALL mp_exchange2d (ng, tile, iNLM, 3, &
3304!^ & LBi, UBi, LBj, UBj, &
3305!^ & NghostPoints, &
3306!^ & EWperiodic(ng), NSperiodic(ng), &
3307!^ & zeta(:,:,knew), &
3308!^ & ubar(:,:,knew), &
3309!^ & vbar(:,:,knew))
3310!^
3311 CALL mp_exchange2d (ng, tile, itlm, 3, &
3312 & lbi, ubi, lbj, ubj, &
3313 & nghostpoints, &
3314 & ewperiodic(ng), nsperiodic(ng), &
3315 & tl_zeta(:,:,knew), &
3316 & tl_ubar(:,:,knew), &
3317 & tl_vbar(:,:,knew))
3318#endif
3319!
3320 RETURN
3321 END SUBROUTINE tl_step2d_tile
3322
3323 END MODULE tl_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_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::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_param::itlm
integer, parameter itlm
Definition mod_param.F:663
mod_scalars
Definition mod_scalars.F:2
mod_scalars::luvsrc
logical, dimension(:), allocatable luvsrc
Definition mod_scalars.F:542
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
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::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::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
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::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::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
tl_obc_volcons_mod
Definition tl_obc_volcons.F:2
tl_obc_volcons_mod::tl_obc_flux_tile
subroutine, public tl_obc_flux_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, h, tl_h, om_v, on_u, ubar, vbar, zeta, tl_ubar, tl_vbar, tl_zeta)
Definition tl_obc_volcons.F:76
tl_obc_volcons_mod::tl_set_duv_bc_tile
subroutine, public tl_set_duv_bc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, om_v, on_u, ubar, vbar, tl_ubar, tl_vbar, drhs, duon, dvom, tl_drhs, tl_duon, tl_dvom)
Definition tl_obc_volcons.F:322
tl_set_depth_mod
Definition tl_set_depth.F:2
tl_set_depth_mod::tl_set_depth
subroutine, public tl_set_depth(ng, tile, model)
Definition tl_set_depth.F:41
tl_step2d_mod
Definition tl_step2d_FB.h:3
tl_step2d_mod::tl_step2d
subroutine, public tl_step2d(ng, tile)
Definition tl_step2d_FB.h:80
tl_step2d_mod::tl_step2d_tile
subroutine tl_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, tl_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, tl_bed_thick, tl_rustr2d, tl_rvstr2d, tl_rulag2d, tl_rvlag2d, ubar_stokes, tl_ubar_stokes, vbar_stokes, tl_vbar_stokes, eq_tide, tl_eq_tide, sustr, tl_sustr, svstr, tl_svstr, pair, rhoa, tl_rhoa, rhos, tl_rhos, tl_du_avg1, tl_du_avg2, tl_dv_avg1, tl_dv_avg2, tl_zt_avg1, rufrc, tl_rufrc, rvfrc, tl_rvfrc, tl_rufrc_bak, tl_rvfrc_bak, tl_du_flux, tl_dv_flux, ubar, tl_ubar, vbar, tl_vbar, zeta, tl_zeta)
Definition tl_step2d_FB.h:294
tl_u2dbc_mod
Definition tl_u2dbc_im.F:2
tl_u2dbc_mod::tl_u2dbc_tile
subroutine, public tl_u2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta, tl_ubar, tl_vbar, tl_zeta)
Definition tl_u2dbc_im.F:64
tl_v2dbc_mod
Definition tl_v2dbc_im.F:2
tl_v2dbc_mod::tl_v2dbc_tile
subroutine, public tl_v2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta, tl_ubar, tl_vbar, tl_zeta)
Definition tl_v2dbc_im.F:64
tl_zetabc_mod
Definition tl_zetabc.F:2
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
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