ROMS
Loading...
Searching...
No Matches
tl_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 tl_step2d_mod
5!
6!git $Id$
7!=======================================================================
8! !
9! It solves perturbation 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 tl_obc_volcons_mod, ONLY : tl_obc_flux_tile, &
55 & tl_set_duv_bc_tile
56#ifdef SOLVE3D
57 USE tl_set_depth_mod, ONLY : tl_set_depth
58#endif
59 USE tl_u2dbc_mod, ONLY : tl_u2dbc_tile
60 USE tl_v2dbc_mod, ONLY : tl_v2dbc_tile
61 USE tl_zetabc_mod, ONLY : tl_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 :: tl_step2d
70!
71 CONTAINS
72!
73!************************************************************************
74 SUBROUTINE tl_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, itlm, 9, __line__, myfile)
90#endif
91 CALL tl_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, itlm, 9, __line__, myfile)
191#endif
192!
193 RETURN
194 END SUBROUTINE tl_step2d
195!
196!***********************************************************************
197 SUBROUTINE tl_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 END DO
716 END DO
717 DO j=jv_range
718 DO i=iv_range
719 cff=0.5_r8*om_v(i,j)
720 cff1=cff*(drhs(i,j)+drhs(i,j-1))
721 tl_cff1=cff*(tl_drhs(i,j)+tl_drhs(i,j-1))
722 dvom(i,j)=vbar(i,j,krhs)*cff1
723 tl_dvom(i,j)=tl_vbar(i,j,krhs)*cff1+ &
724 & vbar(i,j,krhs)*tl_cff1
725 END DO
726 END DO
727
728#undef IR_RANGE
729#undef IU_RANGE
730#undef IV_RANGE
731#undef JR_RANGE
732#undef JU_RANGE
733#undef JV_RANGE
734
735#if defined DISTRIBUTE && \
736 defined uv_adv && defined uv_c4advection && !defined SOLVE3D
737!
738! In distributed-memory, the I- and J-ranges are different and a
739! special exchange is done here to avoid having three ghost points
740! for high-order numerical stencils. Notice that a private array is
741! passed below to the exchange routine. It also applies periodic
742! boundary conditions, if appropriate and no partitions in I- or
743! J-directions.
744!
745 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
746 CALL exchange_u2d_tile (ng, tile, &
747 & imins, imaxs, jmins, jmaxs, &
748 & duon)
749 CALL exchange_u2d_tile (ng, tile, &
750 & imins, imaxs, jmins, jmaxs, &
751 & tl_duon)
752 CALL exchange_v2d_tile (ng, tile, &
753 & imins, imaxs, jmins, jmaxs, &
754 & dvom)
755 CALL exchange_v2d_tile (ng, tile, &
756 & imins, imaxs, jmins, jmaxs, &
757 & tl_dvom)
758 END IF
759 CALL mp_exchange2d (ng, tile, itlm, 4, &
760 & imins, imaxs, jmins, jmaxs, &
761 & nghostpoints, &
762 & ewperiodic(ng), nsperiodic(ng), &
763 & duon, dvom, &
764 & tl_duon, tl_dvom)
765#endif
766!
767! Compute integral mass flux across open boundaries and adjust
768! for volume conservation. Compute BASIC STATE value.
769!
770 IF (any(volcons(:,ng))) THEN
771 CALL obc_flux_tile (ng, tile, &
772 & lbi, ubi, lbj, ubj, &
773 & imins, imaxs, jmins, jmaxs, &
774 & knew, &
775#ifdef MASKING
776 & umask, vmask, &
777#endif
778 & h, om_v, on_u, &
779 & ubar, vbar, zeta)
780!
781! Set vertically integrated mass fluxes DUon and DVom along the open
782! boundaries in such a way that the integral volume is conserved.
783!
784 CALL set_duv_bc_tile (ng, tile, &
785 & lbi, ubi, lbj, ubj, &
786 & imins, imaxs, jmins, jmaxs, &
787 & krhs, &
788#ifdef MASKING
789 & umask, vmask, &
790#endif
791 & om_v, on_u, &
792 & ubar, vbar, &
793 & drhs, duon, dvom)
794 CALL tl_set_duv_bc_tile (ng, tile, &
795 & lbi, ubi, lbj, ubj, &
796 & imins, imaxs, jmins, jmaxs, &
797 & krhs, &
798#ifdef MASKING
799 & umask, vmask, &
800#endif
801 & om_v, on_u, &
802 & ubar, vbar, &
803 & tl_ubar, tl_vbar, &
804 & drhs, duon, dvom, &
805 & tl_drhs, tl_duon, tl_dvom)
806 END IF
807
808#ifdef SOLVE3D
809!
810!-----------------------------------------------------------------------
811! Fields averaged over all barotropic time steps.
812!-----------------------------------------------------------------------
813!
814! Notice that the index ranges here are designed to include physical
815! boundaries only. Periodic ghost points and internal mpi computational
816! margins are NOT included.
817!
818! Reset all barotropic mode time-averaged arrays during the first
819! predictor step. At all subsequent time steps, accumulate averages
820! of the first kind using the DELAYED way. For example, "Zt_avg1" is
821! not summed immediately after the corrector step when computed but
822! during the subsequent predictor substep. It allows saving operations
823! because "DUon" and "DVom" are calculated anyway. The last time step
824! has a special code to add all three barotropic variables after the
825! last corrector substep.
826!
827 IF (predictor_2d_step) THEN ! PREDICTOR STEP
828 IF (first_2d_step) THEN
829 DO j=jstrr,jendr
830 DO i=istrr,iendr
831!^ Zt_avg1(i,j)=0.0_r8
832!^
833 tl_zt_avg1(i,j)=0.0_r8
834!^ DU_avg1(i,j)=0.0_r8
835!^
836 tl_du_avg1(i,j)=0.0_r8
837!^ DV_avg1(i,j)=0.0_r8
838!^
839 tl_dv_avg1(i,j)=0.0_r8
840!^ DU_avg2(i,j)=0.0_r8
841!^
842 tl_du_avg2(i,j)=0.0_r8
843!^ DV_avg2(i,j)=0.0_r8
844!^
845 tl_dv_avg2(i,j)=0.0_r8
846 END DO
847 END DO
848 ELSE
849 cff=weight(1,iif(ng)-1,ng)
850 DO j=jstrr,jendr
851 DO i=istrr,iendr
852!^ Zt_avg1(i,j)=Zt_avg1(i,j)+cff*zeta(i,j,krhs)
853!^
854 tl_zt_avg1(i,j)=tl_zt_avg1(i,j)+cff*tl_zeta(i,j,krhs)
855 IF (i.ge.istr) THEN
856!^ DU_avg1(i,j)=DU_avg1(i,j)+cff*DUon(i,j)
857!^
858 tl_du_avg1(i,j)=tl_du_avg1(i,j)+cff*tl_duon(i,j)
859 END IF
860 IF (j.ge.jstr) THEN
861!^ DV_avg1(i,j)=DV_avg1(i,j)+cff*DVom(i,j)
862!^
863 tl_dv_avg1(i,j)=tl_dv_avg1(i,j)+cff*tl_dvom(i,j)
864 END IF
865 END DO
866 END DO
867 END IF
868 ELSE ! CORRECTOR STEP
869 cff=weight(2,iif(ng),ng)
870 DO j=jstrr,jendr
871 DO i=istrr,iendr
872 IF (i.ge.istr) THEN
873!^ DU_avg2(i,j)=DU_avg2(i,j)+cff*DUon(i,j)
874!^
875 tl_du_avg2(i,j)=tl_du_avg2(i,j)+cff*tl_duon(i,j)
876 END IF
877 IF (j.ge.jstr) THEN
878!^ DV_avg2(i,j)=DV_avg2(i,j)+cff*DVom(i,j)
879!^
880 tl_dv_avg2(i,j)=tl_dv_avg2(i,j)+cff*tl_dvom(i,j)
881 END IF
882 END DO
883 END DO
884 END IF
885#endif
886!
887!-----------------------------------------------------------------------
888! Tangent linear of advance free-surface.
889!-----------------------------------------------------------------------
890!
891! Notice that the new local free-surface is allocated so it can be
892! passed as an argumment to "zetabc" to avoid memory issues.
893!
894 allocate ( tl_zeta_new(imins:imaxs,jmins:jmaxs) )
895 tl_zeta_new = 0.0_r8
896!
897! Get background "zeta_new" from BASIC state. Notice the I- and J-range
898! used to avoid calling nonlinear 'zetabc_local' routine.
899!
900 DO j=lbj,ubj
901 DO i=lbi,ubi
902 zeta_new(i,j)=zeta(i,j,knew)
903#ifdef MASKING
904 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
905# ifdef WET_DRY_NOT_YET
906!! zeta_new(i,j)=zeta_new(i,j)+ &
907!! & (Dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
908# endif
909#endif
910 dnew(i,j)=h(i,j)+zeta_new(i,j)
911 dstp(i,j)=h(i,j)+zeta(i,j,kstp)
912 END DO
913 END DO
914!
915! Compute "zeta_new" at the new time step and interpolate backward for
916! the subsequent computation of barotropic pressure-gradient terms.
917! Notice that during the predictor of the first 2D step in 3D mode,
918! the pressure gradient terms are computed using just zeta(:,:,kstp),
919! i.e., like in the Forward Euler step, rather than the more accurate
920! predictor of generalized RK2. This is to keep it consistent with the
921! computation of pressure gradient in 3D mode, which uses precisely
922! the initial value of "zeta" rather than the value changed by the
923! first barotropic predictor step. Later in this code, just after
924! "rufrc, rvfrc" are finalized, a correction term based on the
925! difference zeta_new(:,:)-zeta(:,:,kstp) to "rubar, rvbar" to make
926! them consistent with generalized RK2 stepping for pressure gradient
927! terms.
928!
929 IF (predictor_2d_step) THEN
930 IF (first_2d_step) THEN ! Modified RK2 time step (with
931 cff=dtfast(ng) ! Forward-Backward feedback with
932#ifdef SOLVE3D
933 cff1=0.0_r8 !==> Forward Euler
934 cff2=1.0_r8
935#else
936 cff1=0.333333333333_r8 ! optimally chosen beta=1/3 and
937 cff2=0.666666666667_r8 ! epsilon=2/3, see below) is used
938#endif
939 cff3=0.0_r8 ! here for the start up.
940 ELSE
941 cff=2.0_r8*dtfast(ng) ! In the code below "zwrk" is
942 cff1=beta ! time-centered at time step "n"
943 cff2=1.0_r8-2.0_r8*beta ! in the case of LF (for all but
944 cff3=beta ! the first time step)
945 END IF
946!
947 DO j=jstrv-1,jend
948 DO i=istru-1,iend
949 fac=cff*pm(i,j)*pn(i,j)
950!^ zeta_new(i,j)=zeta(i,j,kbak)+ &
951!^ & fac*(DUon(i,j)-DUon(i+1,j)+ &
952!^ & DVom(i,j)-DVom(i,j+1))
953!^
954 tl_zeta_new(i,j)=tl_zeta(i,j,kbak)+ &
955 & fac*(tl_duon(i,j)-tl_duon(i+1,j)+ &
956 & tl_dvom(i,j)-tl_dvom(i,j+1))
957#ifdef MASKING
958!^ zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
959!^
960 tl_zeta_new(i,j)=tl_zeta_new(i,j)*rmask(i,j)
961# ifdef WET_DRY_NOT_YET
962!! zeta_new(i,j)=zeta_new(i,j)+ &
963!! & (Dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
964# endif
965#endif
966!^ Dnew(i,j)=zeta_new(i,j)+h(i,j)
967!^
968 tl_dnew(i,j)=tl_zeta_new(i,j)+tl_h(i,j)
969
970 zwrk(i,j)=cff1*zeta_new(i,j)+ &
971 & cff2*zeta(i,j,kstp)+ &
972 & cff3*zeta(i,j,kbak)
973 tl_zwrk(i,j)=cff1*tl_zeta_new(i,j)+ &
974 & cff2*tl_zeta(i,j,kstp)+ &
975 & cff3*tl_zeta(i,j,kbak)
976#if defined VAR_RHO_2D && defined SOLVE3D
977 rzeta(i,j)=(1.0_r8+rhos(i,j))*zwrk(i,j)
978 tl_rzeta(i,j)=(1.0_r8+rhos(i,j))*tl_zwrk(i,j)+ &
979 & tl_rhos(i,j)*zwrk(i,j)
980 rzeta2(i,j)=rzeta(i,j)*zwrk(i,j)
981 tl_rzeta2(i,j)=tl_rzeta(i,j)*zwrk(i,j)+ &
982 & rzeta(i,j)*tl_zwrk(i,j)
983 rzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
984 tl_rzetasa(i,j)=tl_zwrk(i,j)*(rhos(i,j)-rhoa(i,j))+ &
985 & zwrk(i,j)*(tl_rhos(i,j)-tl_rhoa(i,j))
986#else
987 rzeta(i,j)=zwrk(i,j)
988 tl_rzeta(i,j)=tl_zwrk(i,j)
989 rzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
990 tl_rzeta2(i,j)=2.0_r8*tl_zwrk(i,j)*zwrk(i,j)
991#endif
992 END DO
993 END DO
994 ELSE !--> CORRECTOR STEP
995 IF (first_2d_step) THEN
996 cff =0.333333333333_r8 ! Modified RK2 weighting:
997 cff1=0.333333333333_r8 ! here "zwrk" is time-
998 cff2=0.333333333333_r8 ! centered at "n+1/2".
999 cff3=0.0_r8
1000 ELSE
1001 cff =1.0_r8-epsil ! zwrk is always time-
1002 cff1=(0.5_r8-gamma)*epsil ! centered at n+1/2
1003 cff2=(0.5_r8+2.0_r8*gamma)*epsil ! during corrector sub-
1004 cff3=-gamma *epsil ! step.
1005 END IF
1006!
1007 DO j=jstrv-1,jend
1008 DO i=istru-1,iend
1009 fac=dtfast(ng)*pm(i,j)*pn(i,j)
1010!^ zeta_new(i,j)=zeta(i,j,kstp)+ &
1011!^ & fac*(DUon(i,j)-DUon(i+1,j)+ &
1012!^ & DVom(i,j)-DVom(i,j+1))
1013!^
1014 tl_zeta_new(i,j)=tl_zeta(i,j,kstp)+ &
1015 & fac*(tl_duon(i,j)-tl_duon(i+1,j)+ &
1016 & tl_dvom(i,j)-tl_dvom(i,j+1))
1017#ifdef MASKING
1018!^ zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
1019!^
1020 tl_zeta_new(i,j)=tl_zeta_new(i,j)*rmask(i,j)
1021# ifdef WET_DRY_NOT_YET
1022!! zeta_new(i,j)=zeta_new(i,j)+ &
1023!! & (Dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
1024# endif
1025#endif
1026!^ Dnew(i,j)=zeta_new(i,j)+h(i,j)
1027!^
1028 tl_dnew(i,j)=tl_zeta_new(i,j)+tl_h(i,j)
1029
1030 zwrk(i,j)=cff *zeta(i,j,krhs)+ &
1031 & cff1*zeta_new(i,j)+ &
1032 & cff2*zeta(i,j,kstp)+ &
1033 & cff3*zeta(i,j,kbak)
1034 tl_zwrk(i,j)=cff *tl_zeta(i,j,krhs)+ &
1035 & cff1*tl_zeta_new(i,j)+ &
1036 & cff2*tl_zeta(i,j,kstp)+ &
1037 & cff3*tl_zeta(i,j,kbak)
1038#if defined VAR_RHO_2D && defined SOLVE3D
1039 rzeta(i,j)=(1.0_r8+rhos(i,j))*zwrk(i,j)
1040 tl_rzeta(i,j)=(1.0_r8+rhos(i,j))*tl_zwrk(i,j)+ &
1041 & tl_rhos(i,j)*zwrk(i,j)
1042 rzeta2(i,j)=rzeta(i,j)*zwrk(i,j)
1043 tl_rzeta2(i,j)=tl_rzeta(i,j)*zwrk(i,j)+ &
1044 & rzeta(i,j)*tl_zwrk(i,j)
1045 rzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
1046 tl_rzetasa(i,j)=tl_zwrk(i,j)*(rhos(i,j)-rhoa(i,j))+ &
1047 & zwrk(i,j)*(tl_rhos(i,j)-tl_rhoa(i,j))
1048#else
1049 rzeta(i,j)=zwrk(i,j)
1050 tl_rzeta(i,j)=tl_zwrk(i,j)
1051 rzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
1052 tl_rzeta2(i,j)=2.0_r8*tl_zwrk(i,j)*zwrk(i,j)
1053#endif
1054 END DO
1055 END DO
1056 END IF
1057!
1058! Apply mass point sources (volume vertical influx), if any.
1059!
1060! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
1061!
1062 IF (lwsrc(ng)) THEN
1063 DO is=1,nsrc(ng)
1064 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
1065 i=sources(ng)%Isrc(is)
1066 j=sources(ng)%Jsrc(is)
1067 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1068 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1069!^ zeta_new(i,j)=zeta_new(i,j)+ &
1070!^ & SOURCES(ng)%Qbar(is)* &
1071!^ & pm(i,j)*pn(i,j)*dtfast(ng)
1072!^
1073! tl_zeta_new(i,j)=tl_zeta_new(i,j)+0.0_r8
1074 END IF
1075 END IF
1076 END DO
1077 END IF
1078!
1079! Apply boundary conditions to newly computed free-surface "zeta_new"
1080! and load into global state array. Notice that "zeta_new" is always
1081! centered at time step "m+1", while zeta(:,:,knew) should be centered
1082! either at "m+1/2" after predictor step and at "m+1" after corrector.
1083! Chosing it to be this way makes it possible avoid storing RHS for
1084! zeta, ubar, and vbar between predictor and corrector sub-steps.
1085!
1086! Here, we use the local "zetabc" since the private array "zeta_new"
1087! is passed as an argument to allow computing the lateral boundary
1088! conditions on the range IstrU-1:Iend and JstrV-1:Jend, so parallel
1089! tile exchanges are avoided.
1090!
1091!^ CALL zetabc_local (ng, tile, &
1092!^ & LBi, UBi, LBj, UBj, &
1093!^ & IminS, ImaxS, JminS, JmaxS, &
1094!^ & kstp, &
1095!^ & zeta, &
1096!^ & zeta_new)
1097!^
1098 CALL tl_zetabc_local (ng, tile, &
1099 & lbi, ubi, lbj, ubj, &
1100 & imins, imaxs, jmins, jmaxs, &
1101 & kstp, &
1102 & zeta, tl_zeta, &
1103 & zeta_new, tl_zeta_new)
1104!
1105 IF (predictor_2d_step) THEN
1106 IF (first_2d_step) THEN
1107 cff1=0.5_r8
1108 cff2=0.5_r8
1109 cff3=0.0_r8
1110 ELSE
1111 cff1=0.5_r8-gamma
1112 cff2=0.5_r8+2.0_r8*gamma
1113 cff3=-gamma
1114 END IF
1115 DO j=jstrr,jendr
1116 DO i=istrr,iendr
1117!^ zeta(i,j,knew)=cff1*zeta_new(i,j)+ &
1118!^ & cff2*zeta(i,j,kstp)+ &
1119!^ & cff3*zeta(i,j,kbak)
1120!^
1121 tl_zeta(i,j,knew)=cff1*tl_zeta_new(i,j)+ &
1122 & cff2*tl_zeta(i,j,kstp)+ &
1123 & cff3*tl_zeta(i,j,kbak)
1124 END DO
1125 END DO
1126 ELSE
1127 DO j=jstrr,jendr
1128 DO i=istrr,iendr
1129!^ zeta(i,j,knew)=zeta_new(i,j)
1130!^
1131 tl_zeta(i,j,knew)=tl_zeta_new(i,j)
1132 END DO
1133 END DO
1134 END IF
1135!
1136!=======================================================================
1137! Compute right-hand-side for the 2D momentum equations.
1138!=======================================================================
1139#ifdef SOLVE3D
1140!
1141! Notice that we are suppressing the computation of momentum advection,
1142! Coriolis, and lateral viscosity terms in 3D Applications because
1143! these terms are already included in the baroclinic-to-barotropic
1144! forcing arrays "rufrc" and "rvfrc". It does not mean we are entirely
1145! omitting them, but it is a choice between recomputing them at every
1146! barotropic step or keeping them "frozen" during the fast-time
1147! stepping.
1148# ifdef STEP2D_CORIOLIS
1149! However, in some coarse grid applications with larger baroclinic
1150! timestep (say, DT around 20 minutes or larger), adding the Coriolis
1151! term in the barotropic equations is useful since f*DT is no longer
1152! small.
1153# endif
1154#endif
1155!
1156!-----------------------------------------------------------------------
1157! Compute pressure-gradient terms.
1158!-----------------------------------------------------------------------
1159!
1160! Notice that "rubar" and "rvbar" are computed within the same to allow
1161! shared references to array elements (i,j), which increases the
1162! computational density by almost a factor of 1.5 resulting in overall
1163! more efficient code.
1164!
1165 cff1=0.5*g
1166 cff2=0.333333333333_r8
1167#if !defined SOLVE3D && defined ATM_PRESS
1168 fac=0.5_r8*100.0_r8/rho0
1169#endif
1170 DO j=jstr,jend
1171 DO i=istr,iend
1172 IF (i.ge.istru) THEN
1173!^ rubar(i,j)=cff1*on_u(i,j)* &
1174!^ & ((h(i-1,j)+ &
1175!^ & h(i ,j))* &
1176!^ & (rzeta(i-1,j)- &
1177!^ & rzeta(i ,j))+ &
1178#if defined VAR_RHO_2D && defined SOLVE3D
1179!^ & (h(i-1,j)- &
1180!^ & h(i ,j))* &
1181!^ & (rzetaSA(i-1,j)+ &
1182!^ & rzetaSA(i ,j)+ &
1183!^ & cff2*(rhoA(i-1,j)- &
1184!^ & rhoA(i ,j))* &
1185!^ & (zwrk(i-1,j)- &
1186!^ & zwrk(i ,j)))+ &
1187#endif
1188!^ & (rzeta2(i-1,j)- &
1189!^ & rzeta2(i ,j)))
1190!^
1191 tl_rubar(i,j)=cff1*on_u(i,j)* &
1192 & ((tl_h(i-1,j)+ &
1193 & tl_h(i ,j))* &
1194 & (rzeta(i-1,j)- &
1195 & rzeta(i ,j))+ &
1196 & (h(i-1,j)+ &
1197 & h(i ,j))* &
1198 & (tl_rzeta(i-1,j)- &
1199 & tl_rzeta(i ,j))+ &
1200#if defined VAR_RHO_2D && defined SOLVE3D
1201 & (tl_h(i-1,j)- &
1202 & tl_h(i ,j))* &
1203 & (rzetasa(i-1,j)+ &
1204 & rzetasa(i ,j)+ &
1205 & cff2*(rhoa(i-1,j)- &
1206 & rhoa(i ,j))* &
1207 & (zwrk(i-1,j)- &
1208 & zwrk(i ,j)))+ &
1209 & (h(i-1,j)- &
1210 & h(i ,j))* &
1211 & (tl_rzetasa(i-1,j)+ &
1212 & tl_rzetasa(i ,j)+ &
1213 & cff2*((tl_rhoa(i-1,j)- &
1214 & tl_rhoa(i ,j))* &
1215 & (zwrk(i-1,j)- &
1216 & zwrk(i ,j))+ &
1217 & (rhoa(i-1,j)- &
1218 & rhoa(i ,j))* &
1219 & (tl_zwrk(i-1,j)- &
1220 & tl_zwrk(i ,j))))+ &
1221#endif
1222 & (tl_rzeta2(i-1,j)- &
1223 & tl_rzeta2(i ,j)))
1224#if defined ATM_PRESS && !defined SOLVE3D
1225!^ rubar(i,j)=rubar(i,j)- &
1226!^ & fac*on_u(i,j)* &
1227!^ & (h(i-1,j)+h(i,j)+ &
1228!^ & rzeta(i-1,j)+rzeta(i,j))* &
1229!^ & (Pair(i,j)-Pair(i-1,j))
1230!^
1231 tl_rubar(i,j)=tl_rubar(i,j)- &
1232 & fac*on_u(i,j)* &
1233 & (tl_h(i-1,j)+tl_h(i,j)+ &
1234 & tl_rzeta(i-1,j)+tl_rzeta(i,j))* &
1235 & (pair(i,j)-pair(i-1,j))
1236#endif
1237#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
1238!^ rubar(i,j)=rubar(i,j)- &
1239!^ & cff1*on_u(i,j)* &
1240!^ & (h(i-1,j)+h(i,j)+ &
1241!^ & rzeta(i-1,j)+rzeta(i,j))* &
1242!^ & (eq_tide(i,j)-eq_tide(i-1,j))
1243!^
1244 tl_rubar(i,j)=tl_rubar(i,j)- &
1245 & cff1*on_u(i,j)* &
1246 & ((tl_h(i-1,j)+tl_h(i,j)+ &
1247 & tl_rzeta(i-1,j)+tl_rzeta(i,j))* &
1248 & (eq_tide(i,j)-eq_tide(i-1,j))+ &
1249 & (h(i-1,j)+h(i,j)+ &
1250 & rzeta(i-1,j)+rzeta(i,j))* &
1251 & (tl_eq_tide(i,j)-tl_eq_tide(i-1,j)))
1252#endif
1253#ifdef DIAGNOSTICS_UV
1254!! DiaU2rhs(i,j,M2pgrd)=rubar(i,j)
1255#endif
1256 END IF
1257!
1258 IF (j.ge.jstrv) THEN
1259!^ rvbar(i,j)=cff1*om_v(i,j)* &
1260!^ & ((h(i,j-1)+ &
1261!^ & h(i,j ))* &
1262!^ & (rzeta(i,j-1)- &
1263!^ & rzeta(i,j ))+ &
1264#if defined VAR_RHO_2D && defined SOLVE3D
1265!^ & (h(i,j-1)- &
1266!^ & h(i,j ))* &
1267!^ & (rzetaSA(i,j-1)+ &
1268!^ & rzetaSA(i,j )+ &
1269!^ & cff2*(rhoA(i,j-1)- &
1270!^ & rhoA(i,j ))* &
1271!^ & (zwrk(i,j-1)- &
1272!^ & zwrk(i,j )))+ &
1273#endif
1274!^ & (rzeta2(i,j-1)- &
1275!^ & rzeta2(i,j )))
1276!^
1277 tl_rvbar(i,j)=cff1*om_v(i,j)* &
1278 & ((tl_h(i,j-1)+ &
1279 & tl_h(i,j ))* &
1280 & (rzeta(i,j-1)- &
1281 & rzeta(i,j ))+ &
1282 & (h(i,j-1)+ &
1283 & h(i,j ))* &
1284 & (tl_rzeta(i,j-1)- &
1285 & tl_rzeta(i,j ))+ &
1286#if defined VAR_RHO_2D && defined SOLVE3D
1287 & (tl_h(i,j-1)- &
1288 & tl_h(i,j ))* &
1289 & (rzetasa(i,j-1)+ &
1290 & rzetasa(i,j )+ &
1291 & cff2*(rhoa(i,j-1)- &
1292 & rhoa(i,j ))* &
1293 & (zwrk(i,j-1)- &
1294 & zwrk(i,j )))+ &
1295 & (h(i,j-1)- &
1296 & h(i,j ))* &
1297 & (tl_rzetasa(i,j-1)+ &
1298 & tl_rzetasa(i,j )+ &
1299 & cff2*((tl_rhoa(i,j-1)- &
1300 & tl_rhoa(i,j ))* &
1301 & (zwrk(i,j-1)- &
1302 & zwrk(i,j ))+ &
1303 & (rhoa(i,j-1)- &
1304 & rhoa(i,j ))* &
1305 & (tl_zwrk(i,j-1)- &
1306 & tl_zwrk(i,j ))))+ &
1307#endif
1308 & (tl_rzeta2(i,j-1)- &
1309 & tl_rzeta2(i,j )))
1310#if defined ATM_PRESS && !defined SOLVE3D
1311!^ rvbar(i,j)=rvbar(i,j)- &
1312!^ & fac*om_v(i,j)* &
1313!^ & (h(i,j-1)+h(i,j)+ &
1314!^ & rzeta(i,j-1)+rzeta(i,j))* &
1315!^ & (Pair(i,j)-Pair(i,j-1))
1316!^
1317 tl_rvbar(i,j)=tl_rvbar(i,j)- &
1318 & fac*om_v(i,j)* &
1319 & (tl_h(i,j-1)+tl_h(i,j)+ &
1320 & tl_rzeta(i,j-1)+tl_rzeta(i,j))* &
1321 & (pair(i,j)-pair(i,j-1))
1322#endif
1323#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
1324!^ rvbar(i,j)=rvbar(i,j)- &
1325!^ & cff1*om_v(i,j)* &
1326!^ & (h(i,j-1)+h(i,j)+ &
1327!^ & rzeta(i,j-1)+rzeta(i,j))* &
1328!^ & (eq_tide(i,j)-eq_tide(i,j-1))
1329!^
1330 tl_rvbar(i,j)=tl_rvbar(i,j)- &
1331 & cff1*om_v(i,j)* &
1332 & ((tl_h(i,j-1)+tl_h(i,j)+ &
1333 & tl_rzeta(i,j-1)+tl_rzeta(i,j))* &
1334 & (eq_tide(i,j)-eq_tide(i,j-1))+ &
1335 & (h(i,j-1)+h(i,j)+ &
1336 & rzeta(i,j-1)+rzeta(i,j))* &
1337 & (tl_eq_tide(i,j)-tl_eq_tide(i,j-1)))
1338#endif
1339#ifdef DIAGNOSTICS_UV
1340!! DiaV2rhs(i,j,M2pgrd)=rvbar(i,j)
1341#endif
1342 END IF
1343 END DO
1344 END DO
1345
1346#if defined UV_ADV && !defined SOLVE3D
1347!
1348!-----------------------------------------------------------------------
1349! Add in horizontal advection of momentum.
1350!-----------------------------------------------------------------------
1351
1352# ifdef UV_C2ADVECTION
1353!
1354! Second-order, centered differences advection fluxes.
1355!
1356 DO j=jstr,jend
1357 DO i=istru-1,iend
1358!^ UFx(i,j)=0.25_r8* &
1359!^ & (DUon(i,j)+DUon(i+1,j))* &
1360!^ & (ubar(i ,j,krhs)+ &
1361# ifdef WEC_MELLOR
1362!^ & ubar_stokes(i ,j)+ &
1363!^ & ubar_stokes(i+1,j)+ &
1364# endif
1365!^ & ubar(i+1,j,krhs))
1366!^
1367 tl_ufx(i,j)=0.25_r8* &
1368 & ((tl_duon(i,j)+tl_duon(i+1,j))* &
1369 & (ubar(i ,j,krhs)+ &
1370# ifdef WEC_MELLOR
1371 & ubar_stokes(i ,j)+ &
1372 & ubar_stokes(i+1,j)+ &
1373# endif
1374 & ubar(i+1,j,krhs))+ &
1375 & (duon(i,j)+duon(i+1,j))* &
1376 & (tl_ubar(i ,j,krhs)+ &
1377# ifdef WEC_MELLOR
1378 & tl_ubar_stokes(i ,j)+ &
1379 & tl_ubar_stokes(i+1,j)+ &
1380# endif
1381 & tl_ubar(i+1,j,krhs)))
1382 END DO
1383 END DO
1384!
1385 DO j=jstr,jend+1
1386 DO i=istru,iend
1387!^ UFe(i,j)=0.25_r8* &
1388!^ & (DVom(i,j)+DVom(i-1,j))* &
1389!^ & (ubar(i,j ,krhs)+ &
1390# ifdef WEC_MELLOR
1391!^ & ubar_stokes(i,j )+ &
1392!^ & ubar_stokes(i,j-1)+ &
1393# endif
1394!^ & ubar(i,j-1,krhs))
1395!^
1396 tl_ufe(i,j)=0.25_r8* &
1397 & ((tl_dvom(i,j)+tl_dvom(i-1,j))* &
1398 & (ubar(i,j ,krhs)+ &
1399# ifdef WEC_MELLOR
1400 & ubar_stokes(i,j )+ &
1401 & ubar_stokes(i,j-1)+ &
1402# endif
1403 & ubar(i,j-1,krhs))+ &
1404 & (dvom(i,j)+dvom(i-1,j))* &
1405 & (tl_ubar(i,j ,krhs)+ &
1406# ifdef WEC_MELLOR
1407 & tl_ubar_stokes(i,j )+ &
1408 & tl_ubar_stokes(i,j-1)+ &
1409# endif
1410 & tl_ubar(i,j-1,krhs)))
1411 END DO
1412 END DO
1413!
1414 DO j=jstrv,jend
1415 DO i=istr,iend+1
1416!^ VFx(i,j)=0.25_r8* &
1417!^ & (DUon(i,j)+DUon(i,j-1))* &
1418!^ & (vbar(i ,j,krhs)+ &
1419# ifdef WEC_MELLOR
1420!^ & vbar_stokes(i ,j)+ &
1421!^ & vbar_stokes(i-1,j)+ &
1422# endif
1423!^ & vbar(i-1,j,krhs))
1424!^
1425 tl_vfx(i,j)=0.25_r8* &
1426 & ((tl_duon(i,j)+tl_duon(i,j-1))* &
1427 & (vbar(i ,j,krhs)+ &
1428# ifdef WEC_MELLOR
1429 & vbar_stokes(i ,j)+ &
1430 & vbar_stokes(i-1,j)+ &
1431# endif
1432 & vbar(i-1,j,krhs))+ &
1433 & (duon(i,j)+duon(i,j-1))* &
1434 & (tl_vbar(i ,j,krhs)+ &
1435# ifdef WEC_MELLOR
1436 & tl_vbar_stokes(i ,j)+ &
1437 & tl_vbar_stokes(i-1,j)+ &
1438# endif
1439 & tl_vbar(i-1,j,krhs)))
1440 END DO
1441 END DO
1442!
1443 DO j=jstrv-1,jend
1444 DO i=istr,iend
1445!^ VFe(i,j)=0.25_r8* &
1446!^ & (DVom(i,j)+DVom(i,j+1))* &
1447!^ & (vbar(i,j ,krhs)+ &
1448# ifdef WEC_MELLOR
1449!^ & vbar_stokes(i,j )+ &
1450!^ & vbar_stokes(i,j+1)+ &
1451# endif
1452!^ & vbar(i,j+1,krhs))
1453!^
1454 tl_vfe(i,j)=0.25_r8* &
1455 & ((tl_dvom(i,j)+tl_dvom(i,j+1))* &
1456 & (vbar(i,j ,krhs)+ &
1457# ifdef WEC_MELLOR
1458 & vbar_stokes(i,j )+ &
1459 & vbar_stokes(i,j+1)+ &
1460# endif
1461 & vbar(i,j+1,krhs))+ &
1462 & (dvom(i,j)+dvom(i,j+1))* &
1463 & (tl_vbar(i,j ,krhs)+ &
1464# ifdef WEC_MELLOR
1465 & tl_vbar_stokes(i,j )+ &
1466 & tl_vbar_stokes(i,j+1)+ &
1467# endif
1468 & tl_vbar(i,j+1,krhs)))
1469 END DO
1470 END DO
1471
1472# elif defined UV_C4ADVECTION
1473!
1474! Fourth-order, centered differences u-momentum advection fluxes.
1475!
1476 DO j=jstr,jend
1477 DO i=istrum1,iendp1
1478 grad(i,j)=ubar(i-1,j,krhs)-2.0_r8*ubar(i,j,krhs)+ &
1479# ifdef WEC_MELLOR
1480 & ubar_stokes(i-1,j)-2.0_r8*ubar_stokes(i,j)+ &
1481 & ubar_stokes(i+1,j)+ &
1482# endif
1483 & ubar(i+1,j,krhs)
1484 tl_grad(i,j)=tl_ubar(i-1,j,krhs)-2.0_r8*tl_ubar(i,j,krhs)+ &
1485# ifdef WEC_MELLOR
1486 & tl_ubar_stokes(i-1,j)-2.0_r8*tl_ubar_stokes(i,j)+&
1487 & tl_ubar_stokes(i+1,j)+ &
1488# endif
1489 & tl_ubar(i+1,j,krhs)
1490 dgrad(i,j)=duon(i-1,j)-2.0_r8*duon(i,j)+duon(i+1,j)
1491 tl_dgrad(i,j)=tl_duon(i-1,j)-2.0_r8*tl_duon(i,j)+ &
1492 & tl_duon(i+1,j)
1493 END DO
1494 END DO
1495 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1496 IF (domain(ng)%Western_Edge(tile)) THEN
1497 DO j=jstr,jend
1498 grad(istr,j)=grad(istr+1,j)
1499 tl_grad(istr,j)=tl_grad(istr+1,j)
1500 dgrad(istr,j)=dgrad(istr+1,j)
1501 tl_dgrad(istr,j)=tl_dgrad(istr+1,j)
1502 END DO
1503 END IF
1504 END IF
1505 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1506 IF (domain(ng)%Eastern_Edge(tile)) THEN
1507 DO j=jstr,jend
1508 grad(iend+1,j)=grad(iend,j)
1509 tl_grad(iend+1,j)=tl_grad(iend,j)
1510 dgrad(iend+1,j)=dgrad(iend,j)
1511 tl_dgrad(iend+1,j)=tl_dgrad(iend,j)
1512 END DO
1513 END IF
1514 END IF
1515! d/dx(Duu/n)
1516 cff=1.0_r8/6.0_r8
1517 DO j=jstr,jend
1518 DO i=istru-1,iend
1519!^ UFx(i,j)=0.25_r8*(ubar(i ,j,krhs)+ &
1520# ifdef WEC_MELLOR
1521!^ & ubar_stokes(i ,j)+ &
1522!^ & ubar_stokes(i+1,j)+ &
1523# endif
1524!^ & ubar(i+1,j,krhs)- &
1525!^ & cff*(grad (i,j)+grad (i+1,j)))* &
1526!^ & (DUon(i,j)+DUon(i+1,j)- &
1527!^ & cff*(Dgrad(i,j)+Dgrad(i+1,j)))
1528!^
1529 tl_ufx(i,j)=0.25_r8* &
1530 & ((tl_ubar(i ,j,krhs)+ &
1531# ifdef WEC_MELLOR
1532 & tl_ubar_stokes(i ,j)+ &
1533 & tl_ubar_stokes(i+1,j)+ &
1534# endif
1535 & tl_ubar(i+1,j,krhs)- &
1536 & cff*(tl_grad(i,j)+tl_grad(i+1,j)))* &
1537 & (duon(i,j)+duon(i+1,j)- &
1538 & cff*(dgrad(i,j)+dgrad(i+1,j)))+ &
1539 & (ubar(i ,j,krhs)+ &
1540# ifdef WEC_MELLOR
1541 & ubar_stokes(i ,j)+ &
1542 & ubar_stokes(i+1,j)+ &
1543# endif
1544 & ubar(i+1,j,krhs)- &
1545 & cff*(grad(i,j)+grad(i+1,j)))* &
1546 & (tl_duon(i,j)+tl_duon(i+1,j)- &
1547 & cff*(tl_dgrad(i,j)+tl_dgrad(i+1,j))))
1548!
1549 DO j=jstrm1,jendp1
1550 DO i=istru,iend
1551 grad(i,j)=ubar(i,j-1,krhs)-2.0_r8*ubar(i,j,krhs)+ &
1552# ifdef WEC_MELLOR
1553 & ubar_stokes(i,j-1)-2.0_r8*ubar_stokes(i,j)+ &
1554 & ubar_stokes(i,j+1)+ &
1555# endif
1556 & ubar(i,j+1,krhs)
1557 tl_grad(i,j)=tl_ubar(i,j-1,krhs)-2.0_r8*tl_ubar(i,j,krhs)+ &
1558# ifdef WEC_MELLOR
1559 & tl_ubar_stokes(i,j-1)-2.0_r8*tl_ubar_stokes(i,j)+&
1560 & tl_ubar_stokes(i,j+1)+ &
1561# endif
1562 & tl_ubar(i,j+1,krhs)
1563 END DO
1564 END DO
1565!
1566 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1567 IF (domain(ng)%Southern_Edge(tile)) THEN
1568 DO i=istru,iend
1569 grad(i,jstr-1)=grad(i,jstr)
1570 tl_grad(i,jstr-1)=tl_grad(i,jstr)
1571 END DO
1572 END IF
1573 END IF
1574 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1575 IF (domain(ng)%Northern_Edge(tile)) THEN
1576 DO i=istru,iend
1577 grad(i,jend+1)=grad(i,jend)
1578 tl_grad(i,jend+1)=tl_grad(i,jend)
1579 END DO
1580 END IF
1581 END IF
1582 DO j=jstr,jend+1
1583 DO i=istru-1,iend
1584 dgrad(i,j)=dvom(i-1,j)-2.0_r8*dvom(i,j)+dvom(i+1,j)
1585 tl_dgrad(i,j)=tl_dvom(i-1,j)-2.0_r8*tl_dvom(i,j)+ &
1586 & tl_dvom(i+1,j)
1587 END DO
1588 END DO
1589! d/dy(Duv/m)
1590 cff=1.0_r8/6.0_r8
1591 DO j=jstr,jend+1
1592 DO i=istru,iend
1593!^ UFe(i,j)=0.25_r8*(ubar(i,j ,krhs)+ &
1594# ifdef WEC_MELLOR
1595!^ & ubar_stokes(i,j )+ &
1596!^ & ubar_stokes(i,j-1)+ &
1597# endif
1598!^ & ubar(i,j-1,krhs)- &
1599!^ & cff*(grad (i,j)+grad (i,j-1)))* &
1600!^ & (DVom(i,j)+DVom(i-1,j)- &
1601!^ & cff*(Dgrad(i,j)+Dgrad(i-1,j)))
1602!^
1603 tl_ufe(i,j)=0.25_r8* &
1604 & ((tl_ubar(i,j ,krhs)+ &
1605# ifdef WEC_MELLOR
1606 & tl_ubar_stokes(i,j )+ &
1607 & tl_ubar_stokes(i,j-1)+ &
1608# endif
1609 & tl_ubar(i,j-1,krhs)- &
1610 & cff*(tl_grad(i,j)+tl_grad(i,j-1)))* &
1611 & (dvom(i,j)+dvom(i-1,j)- &
1612 & cff*(dgrad(i,j)+dgrad(i-1,j)))+ &
1613 & (ubar(i,j ,krhs)+ &
1614# ifdef WEC_MELLOR
1615 & ubar_stokes(i,j )+ &
1616 & ubar_stokes(i,j-1)+ &
1617# endif
1618 & ubar(i,j-1,krhs)- &
1619 & cff*(grad(i,j)+grad(i,j-1)))* &
1620 & (tl_dvom(i,j)+tl_dvom(i-1,j)- &
1621 & cff*(tl_dgrad(i,j)+tl_dgrad(i-1,j))))
1622 END DO
1623 END DO
1624!
1625! Fourth-order, centered differences v-momentum advection fluxes.
1626!
1627 DO j=jstrv,jend
1628 DO i=istrm1,iendp1
1629 grad(i,j)=vbar(i-1,j,krhs)-2.0_r8*vbar(i,j,krhs)+ &
1630# ifdef WEC_MELLOR
1631 & vbar_stokes(i-1,j)-2.0_r8*vbar_stokes(i,j)+ &
1632 & vbar_stokes(i+1,j)+ &
1633# endif
1634 & vbar(i+1,j,krhs)
1635 tl_grad(i,j)=tl_vbar(i-1,j,krhs)-2.0_r8*tl_vbar(i,j,krhs)+ &
1636# ifdef WEC_MELLOR
1637 & tl_vbar_stokes(i-1,j)-2.0_r8*tl_vbar_stokes(i,j)+&
1638 & tl_vbar_stokes(i+1,j)+ &
1639# endif
1640 & tl_vbar(i+1,j,krhs)
1641 END DO
1642 END DO
1643!
1644 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1645 IF (domain(ng)%Western_Edge(tile)) THEN
1646 DO j=jstrv,jend
1647 grad(istr-1,j)=grad(istr,j)
1648 tl_grad(istr-1,j)=tl_grad(istr,j)
1649 END DO
1650 END IF
1651 END IF
1652 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1653 IF (domain(ng)%Eastern_Edge(tile)) THEN
1654 DO j=jstrv,jend
1655 grad(iend+1,j)=grad(iend,j)
1656 tl_grad(iend+1,j)=tl_grad(iend,j)
1657 END DO
1658 END IF
1659 END IF
1660 DO j=jstrv-1,jend
1661 DO i=istr,iend+1
1662 dgrad(i,j)=duon(i,j-1)-2.0_r8*duon(i,j)+duon(i,j+1)
1663 tl_dgrad(i,j)=tl_duon(i,j-1)-2.0_r8*tl_duon(i,j)+ &
1664 & tl_duon(i,j+1)
1665 END DO
1666 END DO
1667! d/dx(Duv/n)
1668 cff=1.0_r8/6.0_r8
1669 DO j=jstrv,jend
1670 DO i=istr,iend+1
1671!^ VFx(i,j)=0.25_r8*(vbar(i ,j,krhs)+ &
1672# ifdef WEC_MELLOR
1673!^ & vbar_stokes(i ,j)+ &
1674!^ & vbar_stokes(i-1,j)+ &
1675# endif
1676!^ & vbar(i-1,j,krhs)- &
1677!^ & cff*(grad (i,j)+grad (i-1,j)))* &
1678!^ & (DUon(i,j)+DUon(i,j-1)- &
1679!^ & cff*(Dgrad(i,j)+Dgrad(i,j-1)))
1680!^
1681 tl_vfx(i,j)=0.25_r8* &
1682 & ((tl_vbar(i ,j,krhs)+ &
1683# ifdef WEC_MELLOR
1684 & tl_vbar_stokes(i ,j)+ &
1685 & tl_vbar_stokes(i-1,j)+ &
1686# endif
1687 & tl_vbar(i-1,j,krhs)- &
1688 & cff*(tl_grad(i,j)+tl_grad(i-1,j)))* &
1689 & (duon(i,j)+duon(i,j-1)- &
1690 & cff*(dgrad(i,j)+dgrad(i,j-1)))+ &
1691 & (vbar(i ,j,krhs)+ &
1692# ifdef WEC_MELLOR
1693 & vbar_stokes(i ,j)+ &
1694 & vbar_stokes(i-1,j)+ &
1695# endif
1696 & vbar(i-1,j,krhs)- &
1697 & cff*(grad(i,j)+grad(i-1,j)))* &
1698 & (tl_duon(i,j)+tl_duon(i,j-1)- &
1699 & cff*(tl_dgrad(i,j)+tl_dgrad(i,j-1))))
1700 END DO
1701 END DO
1702!
1703 DO j=jstrvm1,jendp1
1704 DO i=istr,iend
1705 grad(i,j)=vbar(i,j-1,krhs)-2.0_r8*vbar(i,j,krhs)+ &
1706# ifdef WEC_MELLOR
1707 & vbar_stokes(i,j-1)-2.0_r8*vbar_stokes(i,j)+ &
1708 & vbar_stokes(i,j+1)+ &
1709# endif
1710 & vbar(i,j+1,krhs)
1711 tl_grad(i,j)=tl_vbar(i,j-1,krhs)-2.0_r8*tl_vbar(i,j,krhs)+ &
1712# ifdef WEC_MELLOR
1713 & tl_vbar_stokes(i,j-1)-2.0_r8*tl_vbar_stokes(i,j)+&
1714 & tl_vbar_stokes(i,j+1)+ &
1715# endif
1716 & tl_vbar(i,j+1,krhs)
1717 dgrad(i,j)=dvom(i,j-1)-2.0_r8*dvom(i,j)+dvom(i,j+1)
1718 tl_dgrad(i,j)=tl_dvom(i,j-1)-2.0_r8*tl_dvom(i,j)+ &
1719 & tl_dvom(i,j+1)
1720 END DO
1721 END DO
1722 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1723 IF (domain(ng)%Southern_Edge(tile)) THEN
1724 DO i=istr,iend
1725 grad(i,jstr)=grad(i,jstr+1)
1726 tl_grad(i,jstr)=tl_grad(i,jstr+1)
1727 dgrad(i,jstr)=dgrad(i,jstr+1)
1728 tl_dgrad(i,jstr)=tl_dgrad(i,jstr+1)
1729 END DO
1730 END IF
1731 END IF
1732 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1733 IF (domain(ng)%Northern_Edge(tile)) THEN
1734 DO i=istr,iend
1735 grad(i,jend+1)=grad(i,jend)
1736 tl_grad(i,jend+1)=tl_grad(i,jend)
1737 dgrad(i,jend+1)=dgrad(i,jend)
1738 tl_dgrad(i,jend+1)=tl_dgrad(i,jend)
1739 END DO
1740 END IF
1741 END IF
1742! d/dy(Dvv/m)
1743 cff=1.0_r8/6.0_r8
1744 DO j=jstrv-1,jend
1745 DO i=istr,iend
1746!^ VFe(i,j)=0.25_r8*(vbar(i,j ,krhs)+ &
1747# ifdef WEC_MELLOR
1748!^ & vbar_stokes(i,j )+ &
1749!^ & vbar_stokes(i,j+1)+ &
1750# endif
1751!^ & vbar(i,j+1,krhs)- &
1752!^ & cff*(grad (i,j)+grad (i,j+1)))* &
1753!^ & (DVom(i,j)+DVom(i,j+1)- &
1754!^ & cff*(Dgrad(i,j)+Dgrad(i,j+1)))
1755!^
1756 tl_vfe(i,j)=0.25_r8* &
1757 & ((tl_vbar(i,j ,krhs)+ &
1758# ifdef WEC_MELLOR
1759 & tl_vbar_stokes(i,j )+ &
1760 & tl_vbar_stokes(i,j+1)+ &
1761# endif
1762 & tl_vbar(i,j+1,krhs)- &
1763 & cff*(tl_grad(i,j)+tl_grad(i,j+1)))* &
1764 & (dvom(i,j)+dvom(i,j+1)- &
1765 & cff*(dgrad(i,j)+dgrad(i,j+1)))+ &
1766 & (vbar(i,j ,krhs)+ &
1767# ifdef WEC_MELLOR
1768 & vbar_stokes(i,j )+ &
1769 & vbar_stokes(i,j+1)+ &
1770# endif
1771 & vbar(i,j+1,krhs)- &
1772 & cff*(grad(i,j)+grad(i,j+1)))* &
1773 & (tl_dvom(i,j)+tl_dvom(i,j+1)- &
1774 & cff*(tl_dgrad(i,j)+tl_dgrad(i,j+1))))
1775 END DO
1776 END DO
1777# endif
1778!
1779! Add advection to RHS terms.
1780!
1781 DO j=jstr,jend
1782 DO i=istr,iend
1783 IF (i.ge.istru) THEN
1784!^ cff1=UFx(i,j)-UFx(i-1,j)
1785!^
1786 tl_cff1=tl_ufx(i,j)-tl_ufx(i-1,j)
1787!^ cff2=UFe(i,j+1)-UFe(i,j)
1788!^
1789 tl_cff2=tl_ufe(i,j+1)-tl_ufe(i,j)
1790!^ fac=cff1+cff2
1791!^
1792 tl_fac=tl_cff1+tl_cff2
1793!^ rubar(i,j)=rubar(i,j)-fac
1794!^
1795 tl_rubar(i,j)=tl_rubar(i,j)-tl_fac
1796# if defined DIAGNOSTICS_UV
1797!! DiaU2rhs(i,j,M2xadv)=-cff1
1798!! DiaU2rhs(i,j,M2yadv)=-cff2
1799!! DiaU2rhs(i,j,M2hadv)=-fac
1800# endif
1801 END IF
1802!
1803 IF (j.ge.jstrv) THEN
1804!^ cff1=VFx(i+1,j)-VFx(i,j)
1805!^
1806 tl_cff1=tl_vfx(i+1,j)-tl_vfx(i,j)
1807!^ cff2=VFe(i,j)-VFe(i,j-1)
1808!^
1809 tl_cff2=tl_vfe(i,j)-tl_vfe(i,j-1)
1810!^ fac=cff1+cff2
1811!^
1812 tl_fac=tl_cff1+tl_cff2
1813!^ rvbar(i,j)=rvbar(i,j)-fac
1814!^
1815 tl_rvbar(i,j)=tl_rvbar(i,j)-tl_fac
1816# if defined DIAGNOSTICS_UV
1817!! DiaV2rhs(i,j,M2xadv)=-cff1
1818!! DiaV2rhs(i,j,M2yadv)=-cff2
1819!! DiaV2rhs(i,j,M2hadv)=-fac
1820# endif
1821 END IF
1822 END DO
1823 END DO
1824#endif
1825
1826#if (defined UV_COR & !defined SOLVE3D) || defined STEP2D_CORIOLIS
1827!
1828!-----------------------------------------------------------------------
1829! Add in Coriolis term.
1830!-----------------------------------------------------------------------
1831!
1832 DO j=jstrv-1,jend
1833 DO i=istru-1,iend
1834 cff=0.5_r8*drhs(i,j)*fomn(i,j)
1835 tl_cff=0.5_r8*tl_drhs(i,j)*fomn(i,j)
1836!^ UFx(i,j)=cff*(vbar(i,j ,krhs)+ &
1837# ifdef WEC_MELLOR
1838!^ & vbar_stokes(i,j )+ &
1839!^ & vbar_stokes(i,j+1)+ &
1840# endif
1841!^ & vbar(i,j+1,krhs))
1842!^
1843 tl_ufx(i,j)=tl_cff*(vbar(i,j ,krhs)+ &
1844# ifdef WEC_MELLOR
1845 & vbar_stokes(i,j )+ &
1846 & vbar_stokes(i,j+1)+ &
1847# endif
1848 & vbar(i,j+1,krhs))+ &
1849 & cff*(tl_vbar(i,j ,krhs)+ &
1850# ifdef WEC_MELLOR
1851 & tl_vbar_stokes(i,j )+ &
1852 & tl_vbar_stokes(i,j+1)+ &
1853# endif
1854 & tl_vbar(i,j+1,krhs))
1855!^ VFe(i,j)=cff*(ubar(i ,j,krhs)+ &
1856# ifdef WEC_MELLOR
1857!^ & ubar_stokes(i ,j)+ &
1858!^ & ubar_stokes(i+1,j)+ &
1859# endif
1860!^ & ubar(i+1,j,krhs))
1861!^
1862 tl_vfe(i,j)=tl_cff*(ubar(i ,j,krhs)+ &
1863# ifdef WEC_MELLOR
1864 & ubar_stokes(i ,j)+ &
1865 & ubar_stokes(i+1,j)+ &
1866# endif
1867 & ubar(i+1,j,krhs))+ &
1868 & cff*(tl_ubar(i ,j,krhs)+ &
1869# ifdef WEC_MELLOR
1870 & tl_ubar_stokes(i ,j)+ &
1871 & tl_ubar_stokes(i+1,j)+ &
1872# endif
1873 & tl_ubar(i+1,j,krhs))
1874 END DO
1875 END DO
1876!
1877 DO j=jstr,jend
1878 DO i=istr,iend
1879 IF (i.ge.istru) THEN
1880!^ fac1=0.5_r8*(UFx(i,j)+UFx(i-1,j))
1881!^
1882 tl_fac1=0.5_r8*(tl_ufx(i,j)+tl_ufx(i-1,j))
1883!^ rubar(i,j)=rubar(i,j)+fac1
1884!^
1885 tl_rubar(i,j)=tl_rubar(i,j)+tl_fac1
1886# if defined DIAGNOSTICS_UV
1887!! DiaU2rhs(i,j,M2fcor)=fac1
1888# endif
1889 END IF
1890!
1891 IF (j.ge.jstrv) THEN
1892!^ fac2=0.5_r8*(VFe(i,j)+VFe(i,j-1))
1893!^
1894 tl_fac2=0.5_r8*(tl_vfe(i,j)+tl_vfe(i,j-1))
1895!^ rvbar(i,j)=rvbar(i,j)-fac2
1896!^
1897 tl_rvbar(i,j)=tl_rvbar(i,j)-tl_fac2
1898# if defined DIAGNOSTICS_UV
1899!! DiaV2rhs(i,j,M2fcor)=-fac2
1900# endif
1901 END IF
1902 END DO
1903 END DO
1904#endif
1905
1906#if (defined CURVGRID && defined UV_ADV) && !defined SOLVE3D
1907!
1908!-----------------------------------------------------------------------
1909! Add in curvilinear transformation terms.
1910!-----------------------------------------------------------------------
1911!
1912 DO j=jstrv-1,jend
1913 DO i=istru-1,iend
1914 cff1=0.5_r8*(vbar(i,j ,krhs)+ &
1915# ifdef WEC_MELLOR
1916 & vbar_stokes(i,j )+ &
1917 & vbar_stokes(i,j+1)+ &
1918# endif
1919 & vbar(i,j+1,krhs))
1920 tl_cff1=0.5_r8*(tl_vbar(i,j ,krhs)+ &
1921# ifdef WEC_MELLOR
1922 & tl_vbar_stokes(i,j )+ &
1923 & tl_vbar_stokes(i,j+1)+ &
1924# endif
1925 & tl_vbar(i,j+1,krhs))
1926 cff2=0.5_r8*(ubar(i ,j,krhs)+ &
1927# ifdef WEC_MELLOR
1928 & ubar_stokes(i ,j)+ &
1929 & ubar_stokes(i+1,j)+ &
1930# endif
1931 & ubar(i+1,j,krhs))
1932 tl_cff2=0.5_r8*(tl_ubar(i ,j,krhs)+ &
1933# ifdef WEC_MELLOR
1934 & tl_ubar_stokes(i ,j)+ &
1935 & tl_ubar_stokes(i+1,j)+ &
1936# endif
1937 & tl_ubar(i+1,j,krhs))
1938 cff3=cff1*dndx(i,j)
1939 tl_cff3=tl_cff1*dndx(i,j)
1940 cff4=cff2*dmde(i,j)
1941 tl_cff4=tl_cff2*dmde(i,j)
1942 cff=drhs(i,j)*(cff3-cff4)
1943 tl_cff=tl_drhs(i,j)*(cff3-cff4)+ &
1944 & drhs(i,j)*(tl_cff3-tl_cff4)
1945!^ UFx(i,j)=cff*cff1
1946!^
1947 tl_ufx(i,j)=tl_cff*cff1+cff*tl_cff1
1948!^ VFe(i,j)=cff*cff2
1949!^
1950 tl_vfe(i,j)=tl_cff*cff2+cff*tl_cff2
1951# if defined DIAGNOSTICS_UV
1952!! cff=Drhs(i,j)*cff4
1953!! Uwrk(i,j)=-cff*cff1 ! ubar equation, ETA-term
1954!! Vwrk(i,j)=-cff*cff2 ! vbar equation, ETA-term
1955# endif
1956 END DO
1957 END DO
1958!
1959 DO j=jstr,jend
1960 DO i=istr,iend
1961 IF (i.ge.istru) THEN
1962!^ fac1=0.5_r8*(UFx(i,j)+UFx(i-1,j))
1963!^
1964 tl_fac1=0.5_r8*(tl_ufx(i,j)+tl_ufx(i-1,j))
1965!^ rubar(i,j)=rubar(i,j)+fac1
1966!^
1967 tl_rubar(i,j)=tl_rubar(i,j)+tl_fac1
1968# if defined DIAGNOSTICS_UV
1969!! fac2=0.5_r8*(Uwrk(i,j)+Uwrk(i-1,j))
1970!! DiaU2rhs(i,j,M2xadv)=DiaU2rhs(i,j,M2xadv)+fac1-fac2
1971!! DiaU2rhs(i,j,M2yadv)=DiaU2rhs(i,j,M2yadv)+fac2
1972!! DiaU2rhs(i,j,M2hadv)=DiaU2rhs(i,j,M2hadv)+fac1
1973# endif
1974 END IF
1975!
1976 IF (j.ge.jstrv) THEN
1977!^ fac1=0.5_r8*(VFe(i,j)+VFe(i,j-1))
1978!^
1979 tl_fac1=0.5_r8*(tl_vfe(i,j)+tl_vfe(i,j-1))
1980!^ rvbar(i,j)=rvbar(i,j)-fac1
1981!^
1982 tl_rvbar(i,j)=tl_rvbar(i,j)-tl_fac1
1983# if defined DIAGNOSTICS_UV
1984!! fac2=0.5_r8*(Vwrk(i,j)+Vwrk(i,j-1))
1985!! DiaV2rhs(i,j,M2xadv)=DiaV2rhs(i,j,M2xadv)-fac1+fac2
1986!! DiaV2rhs(i,j,M2yadv)=DiaV2rhs(i,j,M2yadv)-fac2
1987!! DiaV2rhs(i,j,M2hadv)=DiaV2rhs(i,j,M2hadv)-fac1
1988# endif
1989 END IF
1990 END DO
1991 END DO
1992#endif
1993
1994#if defined UV_VIS2 && !defined SOLVE3D
1995!
1996!-----------------------------------------------------------------------
1997! Add in horizontal harmonic viscosity.
1998!-----------------------------------------------------------------------
1999!
2000! Compute total depth at PSI-points.
2001!
2002 DO j=jstr,jend+1
2003 DO i=istr,iend+1
2004 drhs_p(i,j)=0.25_r8*(drhs(i,j )+drhs(i-1,j )+ &
2005 & drhs(i,j-1)+drhs(i-1,j-1))
2006 tl_drhs_p(i,j)=0.25_r8*(tl_drhs(i,j )+tl_drhs(i-1,j )+ &
2007 & tl_drhs(i,j-1)+tl_drhs(i-1,j-1))
2008 END DO
2009 END DO
2010!
2011! Compute flux-components of the horizontal divergence of the stress
2012! tensor (m5/s2) in XI- and ETA-directions.
2013!
2014 DO j=jstrv-1,jend
2015 DO i=istru-1,iend
2016!^ cff=visc2_r(i,j)*Drhs(i,j)*0.5_r8* &
2017!^ & (pmon_r(i,j)* &
2018!^ & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,krhs)- &
2019!^ & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,krhs))- &
2020!^ & pnom_r(i,j)* &
2021!^ & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,krhs)- &
2022!^ & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,krhs)))
2023!^
2024 tl_cff=visc2_r(i,j)*0.5_r8* &
2025 & (tl_drhs(i,j)* &
2026 & (pmon_r(i,j)* &
2027 & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,krhs)- &
2028 & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,krhs))- &
2029 & pnom_r(i,j)* &
2030 & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,krhs)- &
2031 & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,krhs)))+ &
2032 & drhs(i,j)* &
2033 & (pmon_r(i,j)* &
2034 & ((pn(i ,j)+pn(i+1,j))*tl_ubar(i+1,j,krhs)- &
2035 & (pn(i-1,j)+pn(i ,j))*tl_ubar(i ,j,krhs))- &
2036 & pnom_r(i,j)* &
2037 & ((pm(i,j )+pm(i,j+1))*tl_vbar(i,j+1,krhs)- &
2038 & (pm(i,j-1)+pm(i,j ))*tl_vbar(i,j ,krhs))))
2039!^ UFx(i,j)=on_r(i,j)*on_r(i,j)*cff
2040!^
2041 tl_ufx(i,j)=on_r(i,j)*on_r(i,j)*tl_cff
2042!^ VFe(i,j)=om_r(i,j)*om_r(i,j)*cff
2043!^
2044 tl_vfe(i,j)=om_r(i,j)*om_r(i,j)*tl_cff
2045 END DO
2046 END DO
2047!
2048 DO j=jstr,jend+1
2049 DO i=istr,iend+1
2050!^ cff=visc2_p(i,j)*Drhs_p(i,j)*0.5_r8* &
2051!^ & (pmon_p(i,j)* &
2052!^ & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,krhs)- &
2053!^ & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,krhs))+ &
2054!^ & pnom_p(i,j)* &
2055!^ & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,krhs)- &
2056!^ & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,krhs)))
2057!^
2058 tl_cff=visc2_p(i,j)*0.5_r8* &
2059 & (tl_drhs_p(i,j)* &
2060 & (pmon_p(i,j)* &
2061 & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,krhs)- &
2062 & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,krhs))+ &
2063 & pnom_p(i,j)* &
2064 & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,krhs)- &
2065 & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,krhs)))+ &
2066 & drhs_p(i,j)* &
2067 & (pmon_p(i,j)* &
2068 & ((pn(i ,j-1)+pn(i ,j))*tl_vbar(i ,j,krhs)- &
2069 & (pn(i-1,j-1)+pn(i-1,j))*tl_vbar(i-1,j,krhs))+ &
2070 & pnom_p(i,j)* &
2071 & ((pm(i-1,j )+pm(i,j ))*tl_ubar(i,j ,krhs)- &
2072 & (pm(i-1,j-1)+pm(i,j-1))*tl_ubar(i,j-1,krhs))))
2073# ifdef MASKING
2074!^ cff=cff*pmask(i,j)
2075!^
2076 tl_cff=tl_cff*pmask(i,j
2077# endif
2078# ifdef WET_DRY_NOT_YET
2079!^ cff=cff*pmask_wet(i,j)
2080!^
2081 tl_cff=tl_cff*pmask_wet(i,j)
2082# endif
2083!^ UFe(i,j)=om_p(i,j)*om_p(i,j)*cff
2084!^
2085 tl_ufe(i,j)=om_p(i,j)*om_p(i,j)*tl_cff
2086!^ VFx(i,j)=on_p(i,j)*on_p(i,j)*cff
2087!^
2088 tl_vfx(i,j)=on_p(i,j)*on_p(i,j)*tl_cff
2089 END DO
2090 END DO
2091!
2092! Add in harmonic viscosity.
2093!
2094 DO j=jstr,jend
2095 DO i=istr,iend
2096 IF (i.ge.istru) THEN
2097!^ cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(UFx(i,j )-UFx(i-1,j))
2098!^
2099 tl_cff1=0.5_r8*(pn(i-1,j)+pn(i,j))* &
2100 & (tl_ufx(i,j )-tl_ufx(i-1,j))
2101!^ cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(UFe(i,j+1)-UFe(i ,j))
2102!^
2103 tl_cff2=0.5_r8*(pm(i-1,j)+pm(i,j))* &
2104 & (tl_ufe(i,j+1)-tl_ufe(i ,j))
2105!^ fac=cff1+cff2
2106!^
2107 tl_fac=tl_cff1+tl_cff2
2108!^ rubar(i,j)=rubar(i,j)+fac
2109!^
2110 tl_rubar(i,j)=tl_rubar(i,j)+tl_fac
2111# if defined DIAGNOSTICS_UV
2112!! DiaU2rhs(i,j,M2hvis)=fac
2113!! DiaU2rhs(i,j,M2xvis)=cff1
2114!! DiaU2rhs(i,j,M2yvis)=cff2
2115# endif
2116 END IF
2117!
2118 IF (j.ge.jstrv) THEN
2119!^ cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(VFx(i+1,j)-VFx(i,j ))
2120!^
2121 tl_cff1=0.5_r8*(pn(i,j-1)+pn(i,j))* &
2122 & (tl_vfx(i+1,j)-tl_vfx(i,j ))
2123!^ cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(VFe(i ,j)-VFe(i,j-1))
2124!^
2125 tl_cff2=0.5_r8*(pm(i,j-1)+pm(i,j))* &
2126 & (tl_vfe(i ,j)-tl_vfe(i,j-1))
2127!^ fac=cff1-cff2
2128!^
2129 tl_fac=tl_cff1-tl_cff2
2130!^ rvbar(i,j)=rvbar(i,j)+fac
2131!^
2132 tl_rvbar(i,j)=tl_rvbar(i,j)+tl_fac
2133# if defined DIAGNOSTICS_UV
2134!! DiaV2rhs(i,j,M2hvis)=fac
2135!! DiaV2rhs(i,j,M2xvis)= cff1
2136!! DiaV2rhs(i,j,M2yvis)=-cff2
2137# endif
2138 END IF
2139 END DO
2140 END DO
2141#endif
2142
2143#ifndef SOLVE3D
2144!
2145!-----------------------------------------------------------------------
2146! Add in bottom stress.
2147!-----------------------------------------------------------------------
2148!
2149 DO j=jstr,jend
2150 DO i=istru,iend
2151!^ fac=bustr(i,j)*om_u(i,j)*on_u(i,j)
2152!^
2153 tl_fac=tl_bustr(i,j)*om_u(i,j)*on_u(i,j)
2154!^ rubar(i,j)=rubar(i,j)-fac
2155!^
2156 tl_rubar(i,j)=tl_rubar(i,j)-tl_fac
2157# ifdef DIAGNOSTICS_UV
2158!! DiaU2rhs(i,j,M2bstr)=-fac
2159# endif
2160 END DO
2161 END DO
2162 DO j=jstrv,jend
2163 DO i=istr,iend
2164!^ fac=bvstr(i,j)*om_v(i,j)*on_v(i,j)
2165!^
2166 tl_fac=tl_bvstr(i,j)*om_v(i,j)*on_v(i,j)
2167!^ rvbar(i,j)=rvbar(i,j)-fac
2168!^
2169 tl_rvbar(i,j)=tl_rvbar(i,j)-tl_fac
2170# ifdef DIAGNOSTICS_UV
2171!! DiaV2rhs(i,j,M2bstr)=-fac
2172# endif
2173 END DO
2174 END DO
2175#else
2176# ifdef DIAGNOSTICS_UV
2177!!
2178!! Initialize the stress term if no bottom friction is defined.
2179!!
2180!! DO j=Jstr,Jend
2181!! DO i=IstrU,Iend
2182!! DiaU2rhs(i,j,M2bstr)=0.0_r8
2183!! END DO
2184!! END DO
2185!! DO j=JstrV,Jend
2186!! DO i=Istr,Iend
2187!! DiaV2rhs(i,j,M2bstr)=0.0_r8
2188!! END DO
2189!! END DO
2190# endif
2191#endif
2192
2193#ifdef SOLVE3D
2194!
2195!-----------------------------------------------------------------------
2196! Coupling between 2D and 3D equations.
2197!-----------------------------------------------------------------------
2198!
2199! Before the predictor step of the first barotropic time step, arrays
2200! "rufrc" and "rvfrc" contain vertical integrals of the 3D RHS terms
2201! for the momentum equations (including surface and bottom stresses,
2202! if so prescribed). During the first barotropic time step, convert
2203! them into forcing terms by subtracting the fast-time "rubar" and
2204! "rvbar" from them.
2205!
2206! These forcing terms are then extrapolated forward in time using
2207! optimized Adams-Bashforth weights, so that the resultant "rufrc"
2208! and "rvfrc" are centered effectively at time n+1/2 in baroclinic
2209! time.
2210!
2211! From now on, these newly computed forcing terms remain unchanged
2212! during the fast time stepping and will be added to "rubar" and
2213! "rvbar" during all subsequent barotropic time steps.
2214!
2215! Thus, the algorithm below is designed for coupling during the 3D
2216! predictor sub-step. The forcing terms "rufrc" and "rvfrc" are
2217! computed as instantaneous values at 3D time index "nstp" first and
2218! then extrapolated half-step forward using AM3-like weights optimized
2219! for maximum stability (with particular care for startup).
2220!
2221 IF (first_2d_step.and.predictor_2d_step) THEN
2222 IF (first_time_step) THEN
2223 cff3=0.0_r8
2224 cff2=0.0_r8
2225 cff1=1.0_r8
2226 ELSE IF (first_time_step+1) THEN
2227 cff3=0.0_r8
2228 cff2=-0.5_r8
2229 cff1=1.5_r8
2230 ELSE
2231 cff3=0.281105_r8
2232 cff2=-0.5_r8-2.0_r8*cff3
2233 cff1=1.5_r8+cff3
2234 END IF
2235!
2236 DO j=jstr,jend
2237 DO i=istru,iend
2238!^ cff=rufrc(i,j)-rubar(i,j)
2239!^
2240 tl_cff=tl_rufrc(i,j)-tl_rubar(i,j)
2241!^ rufrc(i,j)=cff1*cff+ &
2242!^ & cff2*rufrc_bak(i,j,3-nstp)+ &
2243!^ & cff3*rufrc_bak(i,j,nstp )
2244!^
2245 tl_rufrc(i,j)=cff1*tl_cff+ &
2246 & cff2*tl_rufrc_bak(i,j,3-nstp)+ &
2247 & cff3*tl_rufrc_bak(i,j,nstp )
2248!^ rufrc_bak(i,j,nstp)=cff
2249!^
2250 tl_rufrc_bak(i,j,nstp)=tl_cff
2251 END DO
2252 END DO
2253 DO j=jstrv,jend
2254 DO i=istr,iend
2255!^ cff=rvfrc(i,j)-rvbar(i,j)
2256!^
2257 tl_cff=tl_rvfrc(i,j)-tl_rvbar(i,j)
2258!^ rvfrc(i,j)=cff1*cff+ &
2259!^ & cff2*rvfrc_bak(i,j,3-nstp)+ &
2260!^ & cff3*rvfrc_bak(i,j,nstp )
2261!^
2262 tl_rvfrc(i,j)=cff1*tl_cff+ &
2263 & cff2*tl_rvfrc_bak(i,j,3-nstp)+ &
2264 & cff3*tl_rvfrc_bak(i,j,nstp )
2265!^ rvfrc_bak(i,j,nstp)=cff
2266!^
2267 tl_rvfrc_bak(i,j,nstp)=tl_cff
2268 END DO
2269 END DO
2270!
2271! Since coupling requires that the pressure gradient term is computed
2272! using zeta(:,:,kstp) instead of 1/3 toward zeta_new(:,:) as needed
2273! by generalized RK2 scheme, apply compensation to shift pressure
2274! gradient terms from "kstp" to 1/3 toward "knew".
2275!
2276 cff1=0.5_r8*g
2277 cff2=0.333333333333_r8
2278 cff3=1.666666666666_r8
2279
2280 DO j=jstrv-1,jend
2281 DO i=istru-1,iend
2282!^ zwrk(i,j)=cff2*(zeta_new(i,j)-zeta(i,j,kstp))
2283!^
2284 tl_zwrk(i,j)=cff2*(tl_zeta_new(i,j)-tl_zeta(i,j,kstp))
2285# if defined VAR_RHO_2D && defined SOLVE3D
2286!^ rzeta(i,j)=(1.0_r8+rhoS(i,j))*zwrk(i,j)
2287!^
2288 tl_rzeta(i,j)=(1.0_r8+rhos(i,j))*tl_zwrk(i,j)+ &
2289 & tl_rhos(i,j)*zwrk(i,j)
2290!^ rzeta2(i,j)=rzeta(i,j)* &
2291!^ & (cff2*zeta_new(i,j)+ &
2292!^ & cff3*zeta(i,j,kstp))
2293!^
2294 tl_rzeta2(i,j)=tl_rzeta(i,j)* &
2295 & (cff2*zeta_new(i,j)+ &
2296 & cff3*zeta(i,j,kstp))+ &
2297 & rzeta(i,j)* &
2298 & (cff2*tl_zeta_new(i,j)+ &
2299 & cff3*tl_zeta(i,j,kstp))
2300!^ rzetaSA(i,j)=zwrk(i,j)*(rhoS(i,j)-rhoA(i,j))
2301!^
2302 tl_rzetasa(i,j)=tl_zwrk(i,j)* &
2303 & (rhos(i,j)-rhoa(i,j))+ &
2304 & zwrk(i,j)* &
2305 & (tl_rhos(i,j)-tl_rhoa(i,j))
2306# else
2307!^ rzeta(i,j)=zwrk(i,j)
2308!^
2309 tl_rzeta(i,j)=tl_zwrk(i,j)
2310!^ rzeta2(i,j)=zwrk(i,j)* &
2311!^ & (cff2*zeta_new(i,j)+ &
2312!^ & cff3*zeta(i,j,kstp))
2313!^
2314 tl_rzeta2(i,j)=tl_zwrk(i,j)* &
2315 & (cff2*zeta_new(i,j)+ &
2316 & cff3*zeta(i,j,kstp))+ &
2317 & zwrk(i,j)* &
2318 & (cff2*tl_zeta_new(i,j)+ &
2319 & cff3*tl_zeta(i,j,kstp))
2320# endif
2321 END DO
2322 END DO
2323!
2324 DO j=jstr,jend
2325 DO i=istr,iend
2326 IF (i.ge.istru) THEN
2327!^ rubar(i,j)=rubar(i,j)+ &
2328!^ & cff1*on_u(i,j)* &
2329!^ & ((h(i-1,j)+ &
2330!^ & h(i ,j))* &
2331!^ & (rzeta(i-1,j)- &
2332!^ & rzeta(i ,j))+ &
2333# if defined VAR_RHO_2D && defined SOLVE3D
2334!^ & (h(i-1,j)- &
2335!^ & h(i ,j))* &
2336!^ & (rzetaSA(i-1,j)+ &
2337!^ & rzetaSA(i ,j)+ &
2338!^ & cff2*(rhoA(i-1,j)- &
2339!^ & rhoA(i ,j))* &
2340!^ & (zwrk(i-1,j)- &
2341!^ & zwrk(i ,j)))+ &
2342# endif
2343!^ & (rzeta2(i-1,j)- &
2344!^ & rzeta2(i ,j)))
2345!^
2346 tl_rubar(i,j)=tl_rubar(i,j)+ &
2347 & cff1*on_u(i,j)* &
2348 & ((tl_h(i-1,j)+ &
2349 & tl_h(i ,j))* &
2350 & (rzeta(i-1,j)- &
2351 & rzeta(i ,j))+ &
2352 & (h(i-1,j)+ &
2353 & h(i ,j))* &
2354 & (tl_rzeta(i-1,j)- &
2355 & tl_rzeta(i ,j))+ &
2356# if defined VAR_RHO_2D && defined SOLVE3D
2357 & (tl_h(i-1,j)- &
2358 & tl_h(i ,j))* &
2359 & (rzetasa(i-1,j)+ &
2360 & rzetasa(i ,j)+ &
2361 & cff2*(rhoa(i-1,j)- &
2362 & rhoa(i ,j))* &
2363 & (zwrk(i-1,j)- &
2364 & zwrk(i ,j)))+ &
2365 & (h(i-1,j)- &
2366 & h(i ,j))* &
2367 & (tl_rzetasa(i-1,j)+ &
2368 & tl_rzetasa(i ,j)+ &
2369 & cff2*((tl_rhoa(i-1,j)- &
2370 & tl_rhoa(i ,j))* &
2371 & (zwrk(i-1,j)- &
2372 & zwrk(i ,j))+ &
2373 & (rhoa(i-1,j)- &
2374 & rhoa(i ,j))* &
2375 & (tl_zwrk(i-1,j)- &
2376 & tl_zwrk(i ,j))))+ &
2377# endif
2378 & (tl_rzeta2(i-1,j)- &
2379 & tl_rzeta2(i ,j)))
2380# ifdef DIAGNOSTICS_UV
2381!! DiaU2rhs(i,j,M2pgrd)=DiaU2rhs(i,j,M2pgrd)+ &
2382!! & rubar(i,j)
2383# endif
2384 END IF
2385!
2386 IF (j.ge.jstrv) THEN
2387!^ rvbar(i,j)=rvbar(i,j)+ &
2388!^ & cff1*om_v(i,j)* &
2389!^ & ((h(i,j-1)+ &
2390!^ & h(i,j ))* &
2391!^ & (rzeta(i,j-1)- &
2392!^ & rzeta(i,j ))+ &
2393# if defined VAR_RHO_2D && defined SOLVE3D
2394!^ & (h(i,j-1)- &
2395!^ & h(i,j ))* &
2396!^ & (rzetaSA(i,j-1)+ &
2397!^ & rzetaSA(i,j )+ &
2398!^ & cff2*(rhoA(i,j-1)- &
2399!^ & rhoA(i,j ))* &
2400!^ & (zwrk(i,j-1)- &
2401!^ & zwrk(i,j )))+ &
2402# endif
2403!^ & (rzeta2(i,j-1)- &
2404!^ & rzeta2(i,j )))
2405!^
2406 tl_rvbar(i,j)=tl_rvbar(i,j)+ &
2407 & cff1*om_v(i,j)* &
2408 & ((tl_h(i,j-1)+ &
2409 & tl_h(i,j ))* &
2410 & (rzeta(i,j-1)- &
2411 & rzeta(i,j ))+ &
2412 & (h(i,j-1)+ &
2413 & h(i,j ))* &
2414 & (tl_rzeta(i,j-1)- &
2415 & tl_rzeta(i,j ))+ &
2416# if defined VAR_RHO_2D && defined SOLVE3D
2417 & (tl_h(i,j-1)- &
2418 & tl_h(i,j ))* &
2419 & (rzetasa(i,j-1)+ &
2420 & rzetasa(i,j )+ &
2421 & cff2*(rhoa(i,j-1)- &
2422 & rhoa(i,j ))* &
2423 & (zwrk(i,j-1)- &
2424 & zwrk(i,j )))+ &
2425 & (h(i,j-1)- &
2426 & h(i,j ))* &
2427 & (tl_rzetasa(i,j-1)+ &
2428 & tl_rzetasa(i,j )+ &
2429 & cff2*((tl_rhoa(i,j-1)- &
2430 & tl_rhoa(i,j ))* &
2431 & (zwrk(i,j-1)- &
2432 & zwrk(i,j ))+ &
2433 & (rhoa(i,j-1)- &
2434 & rhoa(i,j ))* &
2435 & (tl_zwrk(i,j-1)- &
2436 & tl_zwrk(i,j ))))+ &
2437# endif
2438 & (tl_rzeta2(i,j-1)- &
2439 & tl_rzeta2(i,j )))
2440# ifdef DIAGNOSTICS_UV
2441!! DiaV2rhs(i,j,M2pgrd)=DiaV2rhs(i,j,M2pgrd)+ &
2442!! & rvbar(i,j)
2443# endif
2444 END IF
2445 END DO
2446 END DO
2447 END IF
2448#endif
2449!
2450!=======================================================================
2451! Time step 2D momentum equations.
2452!=======================================================================
2453!
2454! Compute total water column depth.
2455!
2456 IF (first_2d_step.or.(.not.predictor_2d_step)) THEN
2457 DO j=jstrv-1,jend
2458 DO i=istru-1,iend
2459!^ Dstp(i,j)=h(i,j)+zeta(i,j,kstp)
2460!^
2461 tl_dstp(i,j)=tl_h(i,j)+tl_zeta(i,j,kstp)
2462 END DO
2463 END DO
2464 ELSE
2465 DO j=jstrv-1,jend
2466 DO i=istru-1,iend
2467!^ Dstp(i,j)=h(i,j)+zeta(i,j,kbak)
2468!^
2469 tl_dstp(i,j)=tl_h(i,j)+tl_zeta(i,j,kbak)
2470 END DO
2471 END DO
2472 END IF
2473!
2474! During the predictor sub-step, once newly computed "ubar" and "vbar"
2475! become available, interpolate them half-step backward in barotropic
2476! time (i.e., they end up time-centered at n+1/2) in order to use it
2477! during subsequent corrector sub-step.
2478!
2479 IF (predictor_2d_step) THEN
2480 IF (first_2d_step) THEN
2481 cff1=0.5_r8*dtfast(ng)
2482 cff2=0.5_r8
2483 cff3=0.5_r8
2484 cff4=0.0_r8
2485 ELSE
2486 cff1=dtfast(ng)
2487 cff2=0.5_r8-gamma
2488 cff3=0.5_r8+2.0_r8*gamma
2489 cff4=-gamma
2490 ENDIF
2491!
2492 DO j=jstr,jend
2493 DO i=istru,iend
2494 cff=cff1*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2495 fac1=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2496 tl_fac1=-fac1*fac1*(tl_dnew(i,j)+tl_dnew(i-1,j))
2497!^ ubar(i,j,knew)=fac1* &
2498!^ & (ubar(i,j,kbak)* &
2499!^ & (Dstp(i,j)+Dstp(i-1,j))+ &
2500#ifdef SOLVE3D
2501!^ & cff*(rubar(i,j)+rufrc(i,j)))
2502#else
2503!^ & cff*rubar(i,j)+4.0_r8*cff1*sustr(i,j))
2504#endif
2505!^
2506 tl_ubar(i,j,knew)=tl_fac1* &
2507 & (ubar(i,j,kbak)* &
2508 & (dstp(i,j)+dstp(i-1,j))+ &
2509#ifdef SOLVE3D
2510 & cff*(rubar(i,j)+rufrc(i,j)))+ &
2511#else
2512 & cff*rubar(i,j)+4.0_r8*cff1*sustr(i,j))+ &
2513#endif
2514 & fac1* &
2515 & (tl_ubar(i,j,kbak)* &
2516 & (dstp(i,j)+dstp(i-1,j))+ &
2517 & ubar(i,j,kbak)* &
2518 & (tl_dstp(i,j)+tl_dstp(i-1,j))+ &
2519#ifdef SOLVE3D
2520 & cff*(tl_rubar(i,j)+tl_rufrc(i,j)))
2521#else
2522 & cff*tl_rubar(i,j)+ &
2523 & 4.0_r8*cff1*tl_sustr(i,j))
2524#endif
2525#ifdef MASKING
2526!^ ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2527!^
2528 tl_ubar(i,j,knew)=tl_ubar(i,j,knew)*umask(i,j)
2529#endif
2530!^ ubar(i,j,knew)=cff2*ubar(i,j,knew)+ &
2531!^ & cff3*ubar(i,j,kstp)+ &
2532!^ & cff4*ubar(i,j,kbak)
2533!^
2534 tl_ubar(i,j,knew)=cff2*tl_ubar(i,j,knew)+ &
2535 & cff3*tl_ubar(i,j,kstp)+ &
2536 & cff4*tl_ubar(i,j,kbak)
2537#ifdef WET_DRY_NOT_YET
2538!^ cff5=ABS(ABS(umask_wet(i,j))-1.0_r8)
2539!^ cff6=0.5_r8+DSIGN(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2540!^ cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2541!^ ubar(i,j,knew)=ubar(i,j,knew)*cff7
2542!^
2543!^ HGA: TLM code needed here.
2544!^
2545#endif
2546#if defined NESTING && !defined SOLVE3D
2547!^ DU_flux(i,j)=0.5_r8*on_u(i,j)* &
2548!^ & (Dnew(i,j)+Dnew(i-1,j))*ubar(i,j,knew)
2549!^
2550 tl_du_flux(i,j)=0.5_r8*on_u(i,j)* &
2551 & ((dnew(i,j)+dnew(i-1,j))* &
2552 & tl_ubar(i,j,knew)+ &
2553 & (tl_dnew(i,j)+tl_dnew(i-1,j))* &
2554 & ubar(i,j,knew))
2555#endif
2556 END DO
2557 END DO
2558!
2559 DO j=jstrv,jend
2560 DO i=istr,iend
2561 cff=cff1*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2562 fac2=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2563 tl_fac2=-fac2*fac2*(tl_dnew(i,j)+tl_dnew(i,j-1))
2564!^ vbar(i,j,knew)=fac2* &
2565!^ & (vbar(i,j,kbak)* &
2566!^ & (Dstp(i,j)+Dstp(i,j-1))+ &
2567#ifdef SOLVE3D
2568!^ & cff*(rvbar(i,j)+rvfrc(i,j)))
2569#else
2570!^ & cff*rvbar(i,j)+4.0_r8*cff1*svstr(i,j))
2571#endif
2572!^
2573 tl_vbar(i,j,knew)=tl_fac2* &
2574 & (vbar(i,j,kbak)* &
2575 & (dstp(i,j)+dstp(i,j-1))+ &
2576#ifdef SOLVE3D
2577 & cff*(rvbar(i,j)+rvfrc(i,j)))+ &
2578#else
2579 & cff*rvbar(i,j)+4.0_r8*cff1*svstr(i,j))+ &
2580#endif
2581 & fac2* &
2582 & (tl_vbar(i,j,kbak)* &
2583 & (dstp(i,j)+dstp(i,j-1))+ &
2584 & vbar(i,j,kbak)* &
2585 & (tl_dstp(i,j)+tl_dstp(i,j-1))+ &
2586#ifdef SOLVE3D
2587 & cff*(tl_rvbar(i,j)+tl_rvfrc(i,j)))
2588#else
2589 & cff*tl_rvbar(i,j)+ &
2590 & 4.0_r8*cff1*tl_svstr(i,j))
2591#endif
2592#ifdef MASKING
2593!^ vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2594!^
2595 tl_vbar(i,j,knew)=tl_vbar(i,j,knew)*vmask(i,j)
2596#endif
2597!^ vbar(i,j,knew)=cff2*vbar(i,j,knew)+ &
2598!^ & cff3*vbar(i,j,kstp)+ &
2599!^ & cff4*vbar(i,j,kbak)
2600!^
2601 tl_vbar(i,j,knew)=cff2*tl_vbar(i,j,knew)+ &
2602 & cff3*tl_vbar(i,j,kstp)+ &
2603 & cff4*tl_vbar(i,j,kbak)
2604#ifdef WET_DRY_NOT_YET
2605!^ cff5=ABS(ABS(vmask_wet(i,j))-1.0_r8)
2606!^ cff6=0.5_r8+DSIGN(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
2607!^ cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2608!^ vbar(i,j,knew)=vbar(i,j,knew)*cff7
2609!^
2610!^ HGA: TLM code needed here.
2611!^
2612#endif
2613#if defined NESTING && !defined SOLVE3D
2614!^ DV_flux(i,j)=0.5_r8*om_v(i,j)* &
2615!^ & (Dnew(i,j)+Dnew(i,j-1))*vbar(i,j,knew)
2616!^
2617 tl_dv_flux(i,j)=0.5_r8*om_v(i,j)* &
2618 & ((dnew(i,j)+dnew(i,j-1))* &
2619 & tl_vbar(i,j,knew)+ &
2620 & (tl_dnew(i,j)+tl_dnew(i,j-1))* &
2621 & vbar(i,j,knew))
2622#endif
2623 END DO
2624 END DO
2625
2626 ELSE !--> CORRECTOR_2D_STEP
2627
2628 cff1=0.5_r8*dtfast(ng)
2629 DO j=jstr,jend
2630 DO i=istru,iend
2631 cff=cff1*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2632 fac1=1.0_r8/(dnew(i,j)+dnew(i-1,j))
2633 tl_fac1=-fac1*fac1*(dnew(i,j)+dnew(i-1,j))
2634!^ ubar(i,j,knew)=fac1* &
2635!^ & (ubar(i,j,kstp)* &
2636!^ & (Dstp(i,j)+Dstp(i-1,j))+ &
2637#ifdef SOLVE3D
2638!^ & cff*(rubar(i,j)+rufrc(i,j)))
2639#else
2640!^ & cff*rubar(i,j)+4.0_r8*cff1*sustr(i,j))
2641#endif
2642!^
2643 tl_ubar(i,j,knew)=tl_fac1* &
2644 & (ubar(i,j,kstp)* &
2645 & (dstp(i,j)+dstp(i-1,j))+ &
2646#ifdef SOLVE3D
2647 & cff*(rubar(i,j)+rufrc(i,j)))+ &
2648#else
2649 & cff*rubar(i,j)+4.0_r8*cff1*sustr(i,j))+ &
2650#endif
2651 & fac1* &
2652 & (tl_ubar(i,j,kstp)* &
2653 & (dstp(i,j)+dstp(i-1,j))+ &
2654 & ubar(i,j,kstp)* &
2655 & (tl_dstp(i,j)+tl_dstp(i-1,j))+ &
2656#ifdef SOLVE3D
2657 & cff*(tl_rubar(i,j)+tl_rufrc(i,j)))
2658#else
2659 & cff*tl_rubar(i,j)+ &
2660 & 4.0_r8*cff1*tl_sustr(i,j))
2661#endif
2662#ifdef MASKING
2663!^ ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
2664!^
2665 tl_ubar(i,j,knew)=tl_ubar(i,j,knew)*umask(i,j)
2666#endif
2667#ifdef WET_DRY_NOT_YET
2668!^ cff5=ABS(ABS(umask_wet(i,j))-1.0_r8)
2669!^ cff6=0.5_r8+DSIGN(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
2670!^ cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2671!^ ubar(i,j,knew)=ubar(i,j,knew)*cff7
2672!^
2673!^ HGA: TLM code needed here.
2674!^
2675#endif
2676#if defined NESTING && !defined SOLVE3D
2677!^ DU_flux(i,j)=0.5_r8*on_u(i,j)* &
2678!^ & (Dnew(i,j)+Dnew(i-1,j))*ubar(i,j,knew)
2679!^
2680 tl_du_flux(i,j)=0.5_r8*on_u(i,j)* &
2681 & ((dnew(i,j)+dnew(i-1,j))* &
2682 & tl_ubar(i,j,knew)+ &
2683 & (tl_dnew(i,j)+tl_dnew(i-1,j))* &
2684 & ubar(i,j,knew))
2685#endif
2686 END DO
2687 END DO
2688!
2689 DO j=jstrv,jend
2690 DO i=istr,iend
2691 cff=cff1*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2692 fac2=1.0_r8/(dnew(i,j)+dnew(i,j-1))
2693 tl_fac2=-fac2*fac2*(tl_dnew(i,j)+tl_dnew(i,j-1))
2694!^ vbar(i,j,knew)=fac2* &
2695!^ & (vbar(i,j,kstp)* &
2696!^ & (Dstp(i,j)+Dstp(i,j-1))+ &
2697#ifdef SOLVE3D
2698!^ & cff*(rvbar(i,j)+rvfrc(i,j)))
2699#else
2700!^ & cff*rvbar(i,j)+4.0_r8*cff1*svstr(i,j))
2701#endif
2702!^
2703 tl_vbar(i,j,knew)=tl_fac2* &
2704 & (vbar(i,j,kstp)* &
2705 & (dstp(i,j)+dstp(i,j-1))+ &
2706#ifdef SOLVE3D
2707 & cff*(rvbar(i,j)+rvfrc(i,j)))+ &
2708#else
2709 & cff*rvbar(i,j)+4.0_r8*cff1*svstr(i,j))+ &
2710#endif
2711 & fac2* &
2712 & (tl_vbar(i,j,kstp)* &
2713 & (dstp(i,j)+dstp(i,j-1))+ &
2714 & vbar(i,j,kstp)* &
2715 & (tl_dstp(i,j)+tl_dstp(i,j-1))+ &
2716#ifdef SOLVE3D
2717 & cff*(tl_rvbar(i,j)+tl_rvfrc(i,j)))
2718#else
2719 & cff*tl_rvbar(i,j)+ &
2720 & 4.0_r8*cff1*svstr(i,j))
2721#endif
2722#ifdef MASKING
2723!^ vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
2724!^
2725 tl_vbar(i,j,knew)=tl_vbar(i,j,knew)*vmask(i,j)
2726#endif
2727#ifdef WET_DRY_NOT_YET
2728!^ cff5=ABS(ABS(vmask_wet(i,j))-1.0_r8)
2729!^ cff6=0.5_r8+DSIGN(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
2730!^ cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
2731!^ vbar(i,j,knew)=vbar(i,j,knew)*cff7
2732!^
2733!^ HGA: TLM code needed here.
2734!^
2735#endif
2736#if defined NESTING && !defined SOLVE3D
2737!^ DV_flux(i,j)=0.5_r8*om_v(i,j)* &
2738!^ & (Dnew(i,j)+Dnew(i,j-1))*vbar(i,j,knew)
2739!^
2740 tl_dv_flux(i,j)=0.5_r8*om_v(i,j)* &
2741 & ((dnew(i,j)+dnew(i,j-1))* &
2742 & tl_vbar(i,j,knew)+ &
2743 & (tl_dnew(i,j)+tl_dnew(i,j-1))* &
2744 & vbar(i,j,knew))
2745#endif
2746 END DO
2747 END DO
2748 END IF
2749!
2750! Apply lateral boundary conditions.
2751!
2752!^ CALL u2dbc_tile (ng, tile, &
2753!^ & LBi, UBi, LBj, UBj, &
2754!^ & IminS, ImaxS, JminS, JmaxS, &
2755!^ & krhs, kstp, knew, &
2756!^ & ubar, vbar, zeta)
2757!^
2758 CALL tl_u2dbc_tile (ng, tile, &
2759 & lbi, ubi, lbj, ubj, &
2760 & imins, imaxs, jmins, jmaxs, &
2761 & krhs, kstp, knew, &
2762 & ubar, vbar, zeta, &
2763 & tl_ubar, tl_vbar, tl_zeta)
2764!^ CALL v2dbc_tile (ng, tile, &
2765!^ & LBi, UBi, LBj, UBj, &
2766!^ & IminS, ImaxS, JminS, JmaxS, &
2767!^ & krhs, kstp, knew, &
2768!^ & ubar, vbar, zeta)
2769!^
2770 CALL tl_v2dbc_tile (ng, tile, &
2771 & lbi, ubi, lbj, ubj, &
2772 & imins, imaxs, jmins, jmaxs, &
2773 & krhs, kstp, knew, &
2774 & ubar, vbar, zeta, &
2775 & tl_ubar, tl_vbar, tl_zeta)
2776!
2777! Compute integral mass flux across open boundaries and adjust
2778! for volume conservation.
2779!
2780 IF (any(volcons(:,ng))) THEN
2781!^ CALL obc_flux_tile (ng, tile, &
2782!^ & LBi, UBi, LBj, UBj, &
2783!^ & IminS, ImaxS, JminS, JmaxS, &
2784!^ & knew, &
2785#ifdef MASKING
2786!^ & umask, vmask, &
2787#endif
2788!^ & h, om_v, on_u, &
2789!^ & ubar, vbar, zeta)
2790!^
2791 CALL tl_obc_flux_tile (ng, tile, &
2792 & lbi, ubi, lbj, ubj, &
2793 & imins, imaxs, jmins, jmaxs, &
2794 & knew, &
2795#ifdef MASKING
2796 & umask, vmask, &
2797#endif
2798 & h, tl_h, om_v, on_u, &
2799 & ubar, vbar, zeta, &
2800 & tl_ubar, tl_vbar, tl_zeta)
2801 END IF
2802
2803#if defined NESTING && !defined SOLVE3D
2804!
2805! Set barotropic fluxes along physical boundaries.
2806!
2807 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
2808 IF (domain(ng)%Western_Edge(tile)) THEN
2809 DO j=jstr-1,jendr
2810!^ Dnew(Istr-1,j)=h(Istr-1,j)+zeta_new(Istr-1,j)
2811!^
2812 tl_dnew(istr-1,j)=tl_h(istr-1,j)+tl_zeta_new(istr-1,j)
2813 END DO
2814 END IF
2815 END IF
2816 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
2817 IF (domain(ng)%Eastern_Edge(tile)) THEN
2818 DO j=jstr-1,jendr
2819!^ Dnew(Iend+1,j)=h(Iend+1,j)+zeta_new(Iend+1,j)
2820!^
2821 tl_dnew(iend+1,j)=tl_h(iend+1,j)+tl_zeta_new(iend+1,j)
2822 END DO
2823 END IF
2824 END IF
2825 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
2826 IF (domain(ng)%Southern_Edge(tile)) THEN
2827 DO i=istr-1,iendr
2828!^ Dnew(i,Jstr-1)=h(i,Jstr-1)+zeta_new(i,Jstr-1)
2829!^
2830 tl_dnew(i,jstr-1)=tl_h(i,jstr-1)+tl_zeta_new(i,jstr-1)
2831 END DO
2832 END IF
2833 END IF
2834 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
2835 IF (domain(ng)%Northern_Edge(tile)) THEN
2836 DO i=istr-1,iendr
2837!^ Dnew(i,Jend+1)=h(i,Jend+1)+zeta_new(i,Jend+1)
2838!^
2839 tl_dnew(i,jend+1)=tl_h(i,jend+1)+tl_zeta_new(i,jend+1)
2840 END DO
2841 END IF
2842 END IF
2843!
2844 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
2845 IF (domain(ng)%Western_Edge(tile)) THEN
2846 DO j=jstrr,jendr
2847!^ DU_flux(IstrU-1,j)=0.5_r8*on_u(IstrU-1,j)* &
2848!^ & (Dnew(IstrU-1,j)+Dnew(IstrU-2,j))* &
2849!^ & ubar(IstrU-1,j,knew)
2850!^
2851 tl_du_flux(istru-1,j)=0.5_r8*on_u(istru-1,j)* &
2852 & ((dnew(istru-1,j)+ &
2853 & dnew(istru-2,j))* &
2854 & tl_ubar(istru-1,j,knew)+ &
2855 & (tl_dnew(istru-1,j)+ &
2856 & tl_dnew(istru-2,j))* &
2857 & ubar(istru-1,j,knew))
2858 END DO
2859 DO j=jstrv,jend
2860!^ DV_flux(Istr-1,j)=0.5_r8*om_v(Istr-1,j)* &
2861!^ & (Dnew(Istr-1,j)+Dnew(Istr-1,j-1))* &
2862!^ & vbar(Istr-1,j,knew)
2863!^
2864 tl_dv_flux(istr-1,j)=0.5_r8*om_v(istr-1,j)* &
2865 & ((dnew(istr-1,j )+ &
2866 & dnew(istr-1,j-1))* &
2867 & tl_vbar(istr-1,j,knew)+ &
2868 & (tl_dnew(istr-1,j )+ &
2869 & tl_dnew(istr-1,j-1))* &
2870 & vbar(istr-1,j,knew))
2871 END DO
2872 END IF
2873 END IF
2874 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
2875 IF (domain(ng)%Eastern_Edge(tile)) THEN
2876 DO j=jstrr,jendr
2877!^ DU_flux(Iend+1,j)=0.5_r8*on_u(Iend+1,j)* &
2878!^ & (Dnew(Iend+1,j)+Dnew(Iend,j))* &
2879!^ & ubar(Iend+1,j,knew)
2880!^
2881 tl_du_flux(iend+1,j)=0.5_r8*on_u(iend+1,j)* &
2882 & ((dnew(iend+1,j)+ &
2883 & dnew(iend ,j))* &
2884 & tl_ubar(iend+1,j,knew)+ &
2885 & (tl_dnew(iend+1,j)+ &
2886 & tl_dnew(iend ,j))* &
2887 & ubar(iend+1,j,knew))
2888 END DO
2889 DO j=jstrv,jend
2890!^ DV_flux(Iend+1,j)=0.5_r8*om_v(Iend+1,j)* &
2891!^ & (Dnew(Iend+1,j)+Dnew(Iend+1,j-1))* &
2892!^ & vbar(Iend+1,j,knew)
2893!^
2894 tl_dv_flux(iend+1,j)=0.5_r8*om_v(iend+1,j)* &
2895 & ((dnew(iend+1,j )+ &
2896 & dnew(iend+1,j-1))* &
2897 & tl_vbar(iend+1,j,knew)+ &
2898 & (tl_dnew(iend+1,j )+ &
2899 & tl_dnew(iend+1,j-1))* &
2900 & vbar(iend+1,j,knew))
2901 END DO
2902 END IF
2903 END IF
2904 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
2905 IF (domain(ng)%Southern_Edge(tile)) THEN
2906 DO i=istru,iend
2907!^ DU_flux(i,Jstr-1)=0.5_r8*on_u(i,Jstr-1)* &
2908!^ & (Dnew(i,Jstr-1)+Dnew(i-1,Jstr-1))* &
2909!^ & ubar(i,Jstr-1,knew)
2910!^
2911 tl_du_flux(i,jstr-1)=0.5_r8*on_u(i,jstr-1)* &
2912 & ((dnew(i ,jstr-1)+ &
2913 & dnew(i-1,jstr-1))* &
2914 & tl_ubar(i,jstr-1,knew)+ &
2915 & (tl_dnew(i ,jstr-1)+ &
2916 & tl_dnew(i-1,jstr-1))* &
2917 & ubar(i,jstr-1,knew))
2918 END DO
2919 DO i=istrr,iendr
2920!^ DV_flux(i,JstrV-1)=0.5_r8*om_v(i,JstrV-1)* &
2921!^ & (Dnew(i,JstrV-1)+Dnew(i,JstrV-2))* &
2922!^ & vbar(i,JstrV-1,knew)
2923!^
2924 tl_dv_flux(i,jstrv-1)=0.5_r8*om_v(i,jstrv-1)* &
2925 & ((dnew(i,jstrv-1)+ &
2926 & dnew(i,jstrv-2))* &
2927 & tl_vbar(i,jstrv-1,knew)+ &
2928 & (tl_dnew(i,jstrv-1)+ &
2929 & tl_dnew(i,jstrv-2))* &
2930 & vbar(i,jstrv-1,knew))
2931 END DO
2932 END IF
2933 END IF
2934 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
2935 IF (domain(ng)%Northern_Edge(tile)) THEN
2936 DO i=istru,iend
2937!^ DU_flux(i,Jend+1)=0.5_r8*on_u(i,Jend+1)* &
2938!^ & (Dnew(i,Jend+1)+Dnew(i-1,Jend+1))* &
2939!^ & ubar(i,Jend+1,knew)
2940!^
2941 tl_du_flux(i,jend+1)=0.5_r8*on_u(i,jend+1)* &
2942 & ((dnew(i ,jend+1)+ &
2943 & dnew(i-1,jend+1))* &
2944 & tl_ubar(i,jend+1,knew)+ &
2945 & (tl_dnew(i ,jend+1)+ &
2946 & tl_dnew(i-1,jend+1))* &
2947 & ubar(i,jend+1,knew))
2948 END DO
2949 DO i=istrr,iendr
2950!^ DV_flux(i,Jend+1)=0.5_r8*om_v(i,Jend+1)* &
2951!^ & (Dnew(i,Jend+1)+Dnew(i,Jend))* &
2952!^ & vbar(i,Jend+1,knew)
2953!^
2954 tl_dv_flux(i,jend+1)=0.5_r8*om_v(i,jend+1)* &
2955 & ((dnew(i,jend+1)+ &
2956 & dnew(i,jend ))* &
2957 & tl_vbar(i,jend+1,knew)+ &
2958 & (tl_dnew(i,jend+1)+ &
2959 & tl_dnew(i,jend ))* &
2960 & vbar(i,jend+1,knew))
2961 END DO
2962 END IF
2963 END IF
2964#endif
2965!
2966! Apply momentum transport point sources (like river runoff), if any.
2967!
2968! Dsrc(is) = 0, flow across grid cell u-face (positive or negative)
2969! Dsrc(is) = 1, flow across grid cell v-face (positive or negative)
2970!
2971 IF (luvsrc(ng)) THEN
2972 DO is=1,nsrc(ng)
2973 i=sources(ng)%Isrc(is)
2974 j=sources(ng)%Jsrc(is)
2975 IF (((istrr.le.i).and.(i.le.iendr)).and. &
2976 & ((jstrr.le.j).and.(j.le.jendr))) THEN
2977 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
2978 cff=1.0_r8/(on_u(i,j)* &
2979 & 0.5_r8*(dnew(i-1,j)+dnew(i,j)))
2980 tl_cff=-cff*cff*on_u(i,j)* &
2981 & 0.5_r8*(tl_dnew(i-1,j)+tl_dnew(i ,j))
2982!^ ubar(i,j,knew)=SOURCES(ng)%Qbar(is)*cff
2983!^
2984 tl_ubar(i,j,knew)=sources(ng)%tl_Qbar(is)*cff+ &
2985 & sources(ng)%Qbar(is)*tl_cff
2986#ifdef SOLVE3D
2987!^ DU_avg1(i,j)=SOURCES(ng)%Qbar(is)
2988!^
2989 tl_du_avg1(i,j)=sources(ng)%tl_Qbar(is)
2990#endif
2991#if defined NESTING && !defined SOLVE3D
2992!^ DU_flux(i,j)=SOURCES(ng)%Qbar(is)
2993!^
2994 tl_du_flux(i,j)=sources(ng)%tl_Qbar(is)
2995#endif
2996 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
2997 cff=1.0_r8/(om_v(i,j)* &
2998 & 0.5_r8*(dnew(i,j-1)+dnew(i,j)))
2999 tl_cff=-cff*cff*om_v(i,j)* &
3000 & 0.5_r8*(tl_dnew(i,j-1)+tl_dnew(i,j))
3001!^ vbar(i,j,knew)=SOURCES(ng)%Qbar(is)*cff
3002!^
3003 tl_vbar(i,j,knew)=sources(ng)%tl_Qbar(is)*cff+ &
3004 & sources(ng)%Qbar(is)*tl_cff
3005#ifdef SOLVE3D
3006!^ DV_avg1(i,j)=SOURCES(ng)%Qbar(is)
3007!^
3008 tl_dv_avg1(i,j)=sources(ng)%tl_Qbar(is)
3009#endif
3010#if defined NESTING && !defined SOLVE3D
3011!^ DV_flux(i,j)=SOURCES(ng)%Qbar(is)
3012!^
3013 tl_dv_flux(i,j)=sources(ng)%tl_Qbar(is)
3014#endif
3015 END IF
3016 END IF
3017 END DO
3018 END IF
3019
3020#ifdef SOLVE3D
3021!
3022!-----------------------------------------------------------------------
3023! Finalize computation of barotropic mode averages.
3024!-----------------------------------------------------------------------
3025!
3026! This procedure starts with filling in boundary rows of total depths
3027! at the new time step, which is needed to be done only during the
3028! last barotropic time step, Normally, the computation of averages
3029! occurs at the beginning of the next predictor step because "DUon"
3030! and "DVom" are being computed anyway. Strictly speaking, the filling
3031! the boundaries are necessary only in the case of open boundaries,
3032! otherwise, the associated fluxes are all zeros.
3033!
3034 IF ((iif(ng).eq.nfast(ng)).and.(knew.lt.3)) THEN
3035 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
3036 IF (domain(ng)%Western_Edge(tile)) THEN
3037 DO j=jstr-1,jendr
3038!^ Dnew(Istr-1,j)=h(Istr-1,j)+zeta_new(Istr-1,j)
3039!^
3040 tl_dnew(istr-1,j)=tl_h(istr-1,j)+tl_zeta_new(istr-1,j)
3041 END DO
3042 END IF
3043 END IF
3044 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
3045 IF (domain(ng)%Eastern_Edge(tile)) THEN
3046 DO j=jstr-1,jendr
3047!^ Dnew(Iend+1,j)=h(Iend+1,j)+zeta_new(Iend+1,j)
3048!^
3049 tl_dnew(iend+1,j)=tl_h(iend+1,j)+tl_zeta_new(iend+1,j)
3050 END DO
3051 END IF
3052 END IF
3053 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
3054 IF (domain(ng)%Southern_Edge(tile)) THEN
3055 DO i=istr-1,iendr
3056!^ Dnew(i,Jstr-1)=h(i,Jstr-1)+zeta_new(i,Jstr-1)
3057!^
3058 tl_dnew(i,jstr-1)=tl_h(i,jstr-1)+tl_zeta_new(i,jstr-1)
3059 END DO
3060 END IF
3061 END IF
3062 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
3063 IF (domain(ng)%Northern_Edge(tile)) THEN
3064 DO i=istr-1,iendr
3065!^ Dnew(i,Jend+1)=h(i,Jend+1)+zeta_new(i,Jend+1)
3066!^
3067 tl_dnew(i,jend+1)=tl_h(i,jend+1)+tl_zeta_new(i,jend+1)
3068 END DO
3069 END IF
3070 END IF
3071!
3072! At the end of the last 2D time step replace the new free-surface
3073! zeta(:,:,knew) with it fast time-averaged value, Zt_avg1. Recall
3074! this is state variable is the one that communicates with the 3D
3075! kernel. Then, compute time-dependent depths.
3076!
3077 cff=weight(1,iif(ng),ng)
3078 cff1=0.5*cff
3079 DO j=jstrr,jendr
3080 DO i=istrr,iendr
3081!^ Zt_avg1(i,j)=Zt_avg1(i,j)+ &
3082!^ & cff*zeta(i,j,knew)
3083!^
3084 tl_zt_avg1(i,j)=tl_zt_avg1(i,j)+ &
3085 & cff*tl_zeta(i,j,knew)
3086 IF (i.ge.istr) THEN
3087!^ DU_avg1(i,j)=DU_avg1(i,j)+ &
3088!^ & cff1*on_u(i,j)* &
3089!^ & (Dnew(i,j)+Dnew(i-1,j))*ubar(i,j,knew)
3090!^
3091 tl_du_avg1(i,j)=tl_du_avg1(i,j)+ &
3092 & cff1*on_u(i,j)* &
3093 & ((dnew(i,j)+dnew(i-1,j))* &
3094 & tl_ubar(i,j,knew)+ &
3095 & (tl_dnew(i,j)+tl_dnew(i-1,j))* &
3096 & ubar(i,j,knew))
3097 END IF
3098 IF (j.ge.jstr) THEN
3099!^ DV_avg1(i,j)=DV_avg1(i,j)+ &
3100!^ & cff1*om_v(i,j)* &
3101!^ & (Dnew(i,j)+Dnew(i,j-1))*vbar(i,j,knew)
3102!^
3103 tl_dv_avg1(i,j)=tl_dv_avg1(i,j)+ &
3104 & cff1*om_v(i,j)* &
3105 & ((dnew(i,j)+dnew(i,j-1))* &
3106 & tl_vbar(i,j,knew)+ &
3107 & (tl_dnew(i,j)+tl_dnew(i,j-1))* &
3108 & vbar(i,j,knew))
3109 END IF
3110!^ zeta(i,j,knew)=Zt_avg1(i,j)
3111!^
3112 tl_zeta(i,j,knew)=tl_zt_avg1(i,j)
3113 END DO
3114 END DO
3115!^ CALL set_depth (ng, tile, iNLM)
3116!^
3117 CALL tl_set_depth (ng, tile, itlm)
3118
3119# ifdef NESTING
3120!
3121! After all fast time steps are completed, apply boundary conditions
3122! to time averaged fields.
3123!
3124! In nesting applications with refinement grids, we need to exchange
3125! the DU_avg2 and DV_avg2 fluxes boundary information for the case
3126! that a contact point is at a tile partition. Notice that in such
3127! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
3128!
3129 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
3130!^ CALL exchange_r2d_tile (ng, tile, &
3131!^ & LBi, UBi, LBj, UBj, &
3132!^ & Zt_avg1)
3133!^
3134 CALL exchange_r2d_tile (ng, tile, &
3135 & lbi, ubi, lbj, ubj, &
3136 & tl_zt_avg1)
3137!^ CALL exchange_u2d_tile (ng, tile, &
3138!^ & LBi, UBi, LBj, UBj, &
3139!^ & DU_avg1)
3140!^
3141 CALL exchange_u2d_tile (ng, tile, &
3142 & lbi, ubi, lbj, ubj, &
3143 & tl_du_avg1)
3144!^ CALL exchange_v2d_tile (ng, tile, &
3145!^ & LBi, UBi, LBj, UBj, &
3146!^ & DV_avg1)
3147!^
3148 CALL exchange_v2d_tile (ng, tile, &
3149 & lbi, ubi, lbj, ubj, &
3150 & tl_dv_avg1)
3151!^ CALL exchange_u2d_tile (ng, tile, &
3152!^ & LBi, UBi, LBj, UBj, &
3153!^ & DU_avg2)
3154!^
3155 CALL exchange_u2d_tile (ng, tile, &
3156 & lbi, ubi, lbj, ubj, &
3157 & tl_du_avg2)
3158!^ CALL exchange_v2d_tile (ng, tile, &
3159!^ & LBi, UBi, LBj, UBj, &
3160!^ & DV_avg2)
3161!^
3162 CALL exchange_v2d_tile (ng, tile, &
3163 & lbi, ubi, lbj, ubj, &
3164 & tl_dv_avg2)
3165 END IF
3166
3167# ifdef DISTRIBUTE
3168!^ CALL mp_exchange2d (ng, tile, iNLM, 3, &
3169!^ & LBi, UBi, LBj, UBj, &
3170!^ & NghostPoints, &
3171!^ & EWperiodic(ng), NSperiodic(ng), &
3172!^ & Zt_avg1, DU_avg1, DV_avg1)
3173!^
3174 CALL mp_exchange2d (ng, tile, itlm, 3, &
3175 & lbi, ubi, lbj, ubj, &
3176 & nghostpoints, &
3177 & ewperiodic(ng), nsperiodic(ng), &
3178 & tl_zt_avg1, tl_du_avg1, tl_dv_avg1)
3179!^ CALL mp_exchange2d (ng, tile, iNLM, 2, &
3180!^ & LBi, UBi, LBj, UBj, &
3181!^ & NghostPoints, &
3182!^ & EWperiodic(ng), NSperiodic(ng), &
3183!^ & DU_avg2, DV_avg2)
3184!^
3185 CALL mp_exchange2d (ng, tile, itlm, 2, &
3186 & lbi, ubi, lbj, ubj, &
3187 & nghostpoints, &
3188 & ewperiodic(ng), nsperiodic(ng), &
3189 & tl_du_avg2, tl_dv_avg2)
3190# endif
3191# endif
3192 END IF
3193#endif
3194#if defined NESTING && !defined SOLVE3D
3195!
3196! In nesting applications with refinement grids, we need to exchange
3197! the DU_flux and DV_flux fluxes boundary information for the case
3198! that a contact point is at a tile partition. Notice that in such
3199! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
3200!
3201 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
3202!^ CALL exchange_u2d_tile (ng, tile, &
3203!^ & LBi, UBi, LBj, UBj, &
3204!^ & DU_flux)
3205!^
3206 CALL exchange_u2d_tile (ng, tile, &
3207 & lbi, ubi, lbj, ubj, &
3208 & tl_du_flux)
3209!^ CALL exchange_v2d_tile (ng, tile, &
3210!^ & LBi, UBi, LBj, UBj, &
3211!^ & DV_flux)
3212!^
3213 CALL exchange_v2d_tile (ng, tile, &
3214 & lbi, ubi, lbj, ubj, &
3215 & tl_dv_flux)
3216 END IF
3217
3218# ifdef DISTRIBUTE
3219!
3220!^ CALL mp_exchange2d (ng, tile, iNLM, 2, &
3221!^ & LBi, UBi, LBj, UBj, &
3222!^ & NghostPoints, &
3223!^ & EWperiodic(ng), NSperiodic(ng), &
3224!^ & DU_flux, DV_flux)
3225!^
3226 CALL mp_exchange2d (ng, tile, itlm, 2, &
3227 & lbi, ubi, lbj, ubj, &
3228 & nghostpoints, &
3229 & ewperiodic(ng), nsperiodic(ng), &
3230 & tl_du_flux, tl_dv_flux)
3231# endif
3232#endif
3233!
3234! Deallocate local new free-surface.
3235!
3236 deallocate ( tl_zeta_new )
3237
3238#ifdef WET_DRY_NOT_YET
3239!
3240!-----------------------------------------------------------------------
3241! Compute new wet/dry masks.
3242!-----------------------------------------------------------------------
3243!
3244!^ CALL wetdry_tile (ng, tile, &
3245!^ & LBi, UBi, LBj, UBj, &
3246!^ & IminS, ImaxS, JminS, JmaxS, &
3247# ifdef MASKING
3248!^ & pmask, rmask, umask, vmask, &
3249# endif
3250!^ & h, zeta(:,:,knew), &
3251# ifdef SOLVE3D
3252!^ & DU_avg1, DV_avg1, &
3253!^ & rmask_wet_avg, &
3254# endif
3255!^ & pmask_wet, pmask_full, &
3256!^ & rmask_wet, rmask_full, &
3257!^ & umask_wet, umask_full, &
3258!^ & vmask_wet, vmask_full)
3259!^
3260!^ HGA: Need the TLM code here.
3261!^
3262#endif
3263!
3264!-----------------------------------------------------------------------
3265! Exchange boundary information.
3266!-----------------------------------------------------------------------
3267!
3268 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
3269!^ CALL exchange_r2d_tile (ng, tile, &
3270!^ & LBi, UBi, LBj, UBj, &
3271!^ & zeta(:,:,knew))
3272!^
3273 CALL exchange_r2d_tile (ng, tile, &
3274 & lbi, ubi, lbj, ubj, &
3275 & tl_zeta(:,:,knew))
3276!^ CALL exchange_u2d_tile (ng, tile, &
3277!^ & LBi, UBi, LBj, UBj, &
3278!^ & ubar(:,:,knew))
3279!^
3280 CALL exchange_u2d_tile (ng, tile, &
3281 & lbi, ubi, lbj, ubj, &
3282 & tl_ubar(:,:,knew))
3283!^ CALL exchange_v2d_tile (ng, tile, &
3284!^ & LBi, UBi, LBj, UBj, &
3285!^ & vbar(:,:,knew))
3286!^
3287 CALL exchange_v2d_tile (ng, tile, &
3288 & lbi, ubi, lbj, ubj, &
3289 & tl_vbar(:,:,knew))
3290 END IF
3291
3292#ifdef DISTRIBUTE
3293!
3294!^ CALL mp_exchange2d (ng, tile, iNLM, 3, &
3295!^ & LBi, UBi, LBj, UBj, &
3296!^ & NghostPoints, &
3297!^ & EWperiodic(ng), NSperiodic(ng), &
3298!^ & zeta(:,:,knew), &
3299!^ & ubar(:,:,knew), &
3300!^ & vbar(:,:,knew))
3301!^
3302 CALL mp_exchange2d (ng, tile, itlm, 3, &
3303 & lbi, ubi, lbj, ubj, &
3304 & nghostpoints, &
3305 & ewperiodic(ng), nsperiodic(ng), &
3306 & tl_zeta(:,:,knew), &
3307 & tl_ubar(:,:,knew), &
3308 & tl_vbar(:,:,knew))
3309#endif
3310!
3311 RETURN
3312 END SUBROUTINE tl_step2d_tile
3313!
3314 END MODULE tl_step2d_mod
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_r2d_tile
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:253
exchange_2d_mod::exchange_u2d_tile
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:440
exchange_2d_mod::exchange_v2d_tile
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:627
mod_coupling
Definition mod_coupling.F:2
mod_coupling::coupling
type(t_coupling), dimension(:), allocatable coupling
Definition mod_coupling.F:102
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_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::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::itlm
integer, parameter itlm
Definition mod_param.F:663
mod_scalars
Definition mod_scalars.F:2
mod_scalars::luvsrc
logical, dimension(:), allocatable luvsrc
Definition mod_scalars.F:542
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::iwest
integer, parameter iwest
Definition mod_scalars.F:1432
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::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
tl_obc_volcons_mod
Definition tl_obc_volcons.F:2
tl_obc_volcons_mod::tl_obc_flux_tile
subroutine, public tl_obc_flux_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, h, tl_h, om_v, on_u, ubar, vbar, zeta, tl_ubar, tl_vbar, tl_zeta)
Definition tl_obc_volcons.F:76
tl_obc_volcons_mod::tl_set_duv_bc_tile
subroutine, public tl_set_duv_bc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, om_v, on_u, ubar, vbar, tl_ubar, tl_vbar, drhs, duon, dvom, tl_drhs, tl_duon, tl_dvom)
Definition tl_obc_volcons.F:322
tl_set_depth_mod
Definition tl_set_depth.F:2
tl_set_depth_mod::tl_set_depth
subroutine, public tl_set_depth(ng, tile, model)
Definition tl_set_depth.F:41
tl_step2d_mod
Definition tl_step2d_FB.h:3
tl_step2d_mod::tl_step2d
subroutine, public tl_step2d(ng, tile)
Definition tl_step2d_FB.h:80
tl_step2d_mod::tl_step2d_tile
subroutine tl_step2d_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, krhs, kstp, knew, nstp, nnew, pmask, rmask, umask, vmask, pmask_wet, pmask_full, rmask_wet, rmask_full, umask_wet, umask_full, vmask_wet, vmask_full, rmask_wet_avg, fomn, h, tl_h, om_u, om_v, on_u, on_v, pm, pn, dndx, dmde, rdrag, rdrag2, pmon_r, pnom_r, pmon_p, pnom_p, om_r, on_r, om_p, on_p, visc2_p, visc2_r, visc4_p, visc4_r, tl_bed_thick, tl_rustr2d, tl_rvstr2d, tl_rulag2d, tl_rvlag2d, ubar_stokes, tl_ubar_stokes, vbar_stokes, tl_vbar_stokes, eq_tide, tl_eq_tide, sustr, tl_sustr, svstr, tl_svstr, pair, rhoa, tl_rhoa, rhos, tl_rhos, tl_du_avg1, tl_du_avg2, tl_dv_avg1, tl_dv_avg2, tl_zt_avg1, rufrc, tl_rufrc, rvfrc, tl_rvfrc, tl_rufrc_bak, tl_rvfrc_bak, tl_du_flux, tl_dv_flux, ubar, tl_ubar, vbar, tl_vbar, zeta, tl_zeta)
Definition tl_step2d_FB.h:294
tl_u2dbc_mod
Definition tl_u2dbc_im.F:2
tl_u2dbc_mod::tl_u2dbc_tile
subroutine, public tl_u2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta, tl_ubar, tl_vbar, tl_zeta)
Definition tl_u2dbc_im.F:64
tl_v2dbc_mod
Definition tl_v2dbc_im.F:2
tl_v2dbc_mod::tl_v2dbc_tile
subroutine, public tl_v2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta, tl_ubar, tl_vbar, tl_zeta)
Definition tl_v2dbc_im.F:64
tl_zetabc_mod
Definition tl_zetabc.F:2
wetdry_mod
Definition wetdry.F:2
wetdry_mod::wetdry_tile
subroutine wetdry_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, pmask, rmask, umask, vmask, h, zeta, du_avg1, dv_avg1, rmask_wet_avg, pmask_wet, pmask_full, rmask_wet, rmask_full, umask_wet, umask_full, vmask_wet, vmask_full)
Definition wetdry.F:108
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3