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