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