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