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