ROMS
Loading...
Searching...
No Matches
step2d_FB_LF_AM3.h
Go to the documentation of this file.
1#define PREDICTOR_2D_STEP (knew.eq.3)
2
3 MODULE step2d_mod
4!
5!git $Id$
6!=======================================================================
7! !
8! Solves Nonlinear (NLM) shallow-water primitive equations predictor !
9! (Leap-frog) and corrector (Adams-Moulton) time-stepping engine with !
10! a Forward Backward feedback. !
11! !
12! The kernel formulation is based on Shchepetkin and McWilliams !
13! (2005), equations (2.38)-(2.39) and (2.40)-(2.41). !
14! !
15! Reference: !
16! !
17! Shchepetkin, A.F. and J.C. McWilliams, 2005: The regional oceanic !
18! modeling system (ROMS): a split-explicit, free-surface, !
19! topography-following-coordinate oceanic model, Ocean Modelling, !
20! 9, 347-404, doi:10.1016/j.ocemod.2004.08.002. !
21! !
22!=======================================================================
23!
24 USE mod_param
25 USE mod_parallel
26#ifdef SOLVE3D
27 USE mod_coupling
28#endif
29#ifdef DIAGNOSTICS_UV
30 USE mod_diags
31#endif
32 USE mod_forces
33 USE mod_grid
34#if defined UV_VIS2
35 USE mod_mixing
36#endif
37 USE mod_ncparam
38 USE mod_scalars
39 USE mod_ocean
40#if defined SEDIMENT && defined SED_MORPH && defined SOLVE3D
41 USE mod_sedbed
42#endif
43 USE mod_sources
44 USE mod_stepping
45!
46 USE exchange_2d_mod
47#ifdef DISTRIBUTE
48 USE mp_exchange_mod, ONLY : mp_exchange2d
49#endif
50 USE obc_volcons_mod, ONLY : obc_flux_tile, &
51 & set_duv_bc_tile
52#ifdef SOLVE3D
53 USE set_depth_mod, ONLY : set_depth
54#endif
55 USE u2dbc_mod, ONLY : u2dbc_tile
56 USE v2dbc_mod, ONLY : v2dbc_tile
57#ifdef WET_DRY
58 USE wetdry_mod, ONLY : wetdry_tile
59#endif
60 USE zetabc_mod, ONLY : zetabc_local
61!
62 implicit none
63!
64 PRIVATE
65 PUBLIC :: step2d
66!
67 CONTAINS
68!
69!************************************************************************
70 SUBROUTINE step2d (ng, tile)
71!************************************************************************
72!
73! Imported variable declarations.
74!
75 integer, intent(in) :: ng, tile
76!
77! Local variable declarations.
78!
79 character (len=*), parameter :: MyFile = &
80 & __FILE__
81!
82#include "tile.h"
83!
84#ifdef PROFILE
85 CALL wclock_on (ng, inlm, 9, __line__, myfile)
86#endif
87 CALL step2d_tile (ng, tile, &
88 & lbi, ubi, lbj, ubj, &
89 & imins, imaxs, jmins, jmaxs, &
90 & kstp(ng), knew(ng), &
91#ifdef SOLVE3D
92 & nstp(ng), nnew(ng), &
93#endif
94#ifdef MASKING
95 & grid(ng) % pmask, grid(ng) % rmask, &
96 & grid(ng) % umask, grid(ng) % vmask, &
97#endif
98#ifdef WET_DRY
99 & grid(ng) % pmask_wet, grid(ng) % pmask_full, &
100 & grid(ng) % rmask_wet, grid(ng) % rmask_full, &
101 & grid(ng) % umask_wet, grid(ng) % umask_full, &
102 & grid(ng) % vmask_wet, grid(ng) % vmask_full, &
103# ifdef SOLVE3D
104 & grid(ng) % rmask_wet_avg, &
105# endif
106#endif
107#if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
108 & grid(ng) % fomn, &
109#endif
110 & grid(ng) % h, &
111 & grid(ng) % om_u, grid(ng) % om_v, &
112 & grid(ng) % on_u, grid(ng) % on_v, &
113 & grid(ng) % omn, &
114 & grid(ng) % pm, grid(ng) % pn, &
115#if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
116 & grid(ng) % dndx, grid(ng) % dmde, &
117#endif
118#if defined UV_VIS2 && !defined SOLVE3D
119 & grid(ng) % pmon_r, grid(ng) % pnom_r, &
120 & grid(ng) % pmon_p, grid(ng) % pnom_p, &
121 & grid(ng) % om_r, grid(ng) % on_r, &
122 & grid(ng) % om_p, grid(ng) % on_p, &
123 & mixing(ng) % visc2_p, &
124 & mixing(ng) % visc2_r, &
125#endif
126#if defined SEDIMENT && defined SED_MORPH
127 & sedbed(ng) % bed_thick, &
128#endif
129#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
130 & ocean(ng) % eq_tide, &
131#endif
132#ifndef SOLVE3D
133 & forces(ng) % sustr, forces(ng) % svstr, &
134 & forces(ng) % bustr, forces(ng) % bvstr, &
135# ifdef ATM_PRESS
136 & forces(ng) % Pair, &
137# endif
138#else
139# ifdef VAR_RHO_2D
140 & coupling(ng) % rhoA, coupling(ng) % rhoS, &
141# endif
142 & coupling(ng) % DU_avg1, coupling(ng) % DU_avg2, &
143 & coupling(ng) % DV_avg1, coupling(ng) % DV_avg2, &
144 & coupling(ng) % Zt_avg1, &
145 & coupling(ng) % rufrc, &
146 & coupling(ng) % rvfrc, &
147 & coupling(ng) % rufrc_bak, &
148 & coupling(ng) % rvfrc_bak, &
149#endif
150#ifdef DIAGNOSTICS_UV
151 & diags(ng) % DiaU2wrk, diags(ng) % DiaV2wrk, &
152 & diags(ng) % DiaRUbar, diags(ng) % DiaRVbar, &
153# ifdef SOLVE3D
154 & diags(ng) % DiaU2int, diags(ng) % DiaV2int, &
155 & diags(ng) % DiaRUfrc, diags(ng) % DiaRVfrc, &
156# endif
157#endif
158#if defined NESTING && !defined SOLVE3D
159 & ocean(ng) % DU_flux, ocean(ng) % DV_flux, &
160#endif
161 & ocean(ng) % ubar, ocean(ng) % vbar, &
162 & ocean(ng) % zeta)
163#ifdef PROFILE
164 CALL wclock_off (ng, inlm, 9, __line__, myfile)
165#endif
166!
167 RETURN
168 END SUBROUTINE step2d
169!
170!***********************************************************************
171 SUBROUTINE step2d_tile (ng, tile, &
172 & LBi, UBi, LBj, UBj, &
173 & IminS, ImaxS, JminS, JmaxS, &
174 & kstp, knew, &
175#ifdef SOLVE3D
176 & nstp, nnew, &
177#endif
178#ifdef MASKING
179 & pmask, rmask, umask, vmask, &
180#endif
181#ifdef WET_DRY
182 & pmask_wet, pmask_full, &
183 & rmask_wet, rmask_full, &
184 & umask_wet, umask_full, &
185 & vmask_wet, vmask_full, &
186# ifdef SOLVE3D
187 & rmask_wet_avg, &
188# endif
189#endif
190#if (defined UV_COR && !defined SOLVE3D) || defined step2d_coriolis
191 & fomn, &
192#endif
193 & h, &
194 & om_u, om_v, on_u, on_v, omn, pm, pn, &
195#if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
196 & dndx, dmde, &
197#endif
198#if defined UV_VIS2 && !defined SOLVE3D
199 & pmon_r, pnom_r, pmon_p, pnom_p, &
200 & om_r, on_r, om_p, on_p, &
201 & visc2_p, visc2_r, &
202#endif
203#if defined SEDIMENT && defined SED_MORPH
204 & bed_thick, &
205#endif
206#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
207 & eq_tide, &
208#endif
209#ifndef SOLVE3D
210 & sustr, svstr, bustr, bvstr, &
211# ifdef ATM_PRESS
212 & pair, &
213# endif
214#else
215# ifdef VAR_RHO_2D
216 & rhoa, rhos, &
217# endif
218 & du_avg1, du_avg2, &
219 & dv_avg1, dv_avg2, &
220 & zt_avg1, &
221 & rufrc, rvfrc, &
222 & rufrc_bak, rvfrc_bak, &
223#endif
224#ifdef DIAGNOSTICS_UV
225 & diau2wrk, diav2wrk, &
226 & diarubar, diarvbar, &
227# ifdef SOLVE3D
228 & diau2int, diav2int, &
229 & diarufrc, diarvfrc, &
230# endif
231#endif
232#if defined NESTING && !defined SOLVE3D
233 & du_flux, dv_flux, &
234#endif
235 & ubar, vbar, zeta)
236!***********************************************************************
237!
238! Imported variable declarations.
239!
240 integer, intent(in ) :: ng, tile
241 integer, intent(in ) :: LBi, UBi, LBj, UBj
242 integer, intent(in ) :: IminS, ImaxS, JminS, JmaxS
243 integer, intent(in ) :: kstp, knew
244#ifdef SOLVE3D
245 integer, intent(in ) :: nstp, nnew
246#endif
247!
248#ifdef ASSUMED_SHAPE
249# ifdef MASKING
250 real(r8), intent(in ) :: pmask(LBi:,LBj:)
251 real(r8), intent(in ) :: rmask(LBi:,LBj:)
252 real(r8), intent(in ) :: umask(LBi:,LBj:)
253 real(r8), intent(in ) :: vmask(LBi:,LBj:)
254# endif
255# if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
256 real(r8), intent(in ) :: fomn(LBi:,LBj:)
257# endif
258# if defined SEDIMENT && defined SED_MORPH
259 real(r8), intent(inout) :: h(LBi:,LBj:)
260# else
261 real(r8), intent(in ) :: h(LBi:,LBj:)
262# endif
263 real(r8), intent(in ) :: om_u(LBi:,LBj:)
264 real(r8), intent(in ) :: om_v(LBi:,LBj:)
265 real(r8), intent(in ) :: on_u(LBi:,LBj:)
266 real(r8), intent(in ) :: on_v(LBi:,LBj:)
267 real(r8), intent(in ) :: omn(LBi:,LBj:)
268 real(r8), intent(in ) :: pm(LBi:,LBj:)
269 real(r8), intent(in ) :: pn(LBi:,LBj:)
270# if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
271 real(r8), intent(in ) :: dndx(LBi:,LBj:)
272 real(r8), intent(in ) :: dmde(LBi:,LBj:)
273# endif
274# if defined UV_VIS2 && !defined SOLVE3D
275 real(r8), intent(in ) :: pmon_r(LBi:,LBj:)
276 real(r8), intent(in ) :: pnom_r(LBi:,LBj:)
277 real(r8), intent(in ) :: pmon_p(LBi:,LBj:)
278 real(r8), intent(in ) :: pnom_p(LBi:,LBj:)
279 real(r8), intent(in ) :: om_r(LBi:,LBj:)
280 real(r8), intent(in ) :: on_r(LBi:,LBj:)
281 real(r8), intent(in ) :: om_p(LBi:,LBj:)
282 real(r8), intent(in ) :: on_p(LBi:,LBj:)
283 real(r8), intent(in ) :: visc2_p(LBi:,LBj:)
284 real(r8), intent(in ) :: visc2_r(LBi:,LBj:)
285# endif
286# if defined SEDIMENT && defined SED_MORPH
287 real(r8), intent(in ) :: bed_thick(LBi:,LBj:,:)
288# endif
289# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
290 real(r8), intent(in ) :: eq_tide(LBi:,LBj:)
291# endif
292# ifndef SOLVE3D
293 real(r8), intent(in ) :: sustr(LBi:,LBj:)
294 real(r8), intent(in ) :: svstr(LBi:,LBj:)
295 real(r8), intent(in ) :: bustr(LBi:,LBj:)
296 real(r8), intent(in ) :: bvstr(LBi:,LBj:)
297# ifdef ATM_PRESS
298 real(r8), intent(in ) :: Pair(LBi:,LBj:)
299# endif
300# else
301# ifdef VAR_RHO_2D
302 real(r8), intent(in ) :: rhoA(LBi:,LBj:)
303 real(r8), intent(in ) :: rhoS(LBi:,LBj:)
304# endif
305 real(r8), intent(inout) :: DU_avg1(LBi:,LBj:)
306 real(r8), intent(inout) :: DU_avg2(LBi:,LBj:)
307 real(r8), intent(inout) :: DV_avg1(LBi:,LBj:)
308 real(r8), intent(inout) :: DV_avg2(LBi:,LBj:)
309 real(r8), intent(inout) :: Zt_avg1(LBi:,LBj:)
310 real(r8), intent(inout) :: rufrc(LBi:,LBj:)
311 real(r8), intent(inout) :: rvfrc(LBi:,LBj:)
312 real(r8), intent(inout) :: rufrc_bak(LBi:,LBj:,:)
313 real(r8), intent(inout) :: rvfrc_bak(LBi:,LBj:,:)
314# endif
315# ifdef WET_DRY
316 real(r8), intent(inout) :: pmask_full(LBi:,LBj:)
317 real(r8), intent(inout) :: rmask_full(LBi:,LBj:)
318 real(r8), intent(inout) :: umask_full(LBi:,LBj:)
319 real(r8), intent(inout) :: vmask_full(LBi:,LBj:)
320
321 real(r8), intent(inout) :: pmask_wet(LBi:,LBj:)
322 real(r8), intent(inout) :: rmask_wet(LBi:,LBj:)
323 real(r8), intent(inout) :: umask_wet(LBi:,LBj:)
324 real(r8), intent(inout) :: vmask_wet(LBi:,LBj:)
325# ifdef SOLVE3D
326 real(r8), intent(inout) :: rmask_wet_avg(LBi:,LBj:)
327# endif
328# endif
329# ifdef DIAGNOSTICS_UV
330 real(r8), intent(inout) :: DiaU2wrk(LBi:,LBj:,:)
331 real(r8), intent(inout) :: DiaV2wrk(LBi:,LBj:,:)
332 real(r8), intent(inout) :: DiaRUbar(LBi:,LBj:,:,:)
333 real(r8), intent(inout) :: DiaRVbar(LBi:,LBj:,:,:)
334# ifdef SOLVE3D
335 real(r8), intent(inout) :: DiaU2int(LBi:,LBj:,:)
336 real(r8), intent(inout) :: DiaV2int(LBi:,LBj:,:)
337 real(r8), intent(inout) :: DiaRUfrc(LBi:,LBj:,:,:)
338 real(r8), intent(inout) :: DiaRVfrc(LBi:,LBj:,:,:)
339# endif
340# endif
341 real(r8), intent(inout) :: ubar(LBi:,LBj:,:)
342 real(r8), intent(inout) :: vbar(LBi:,LBj:,:)
343 real(r8), intent(inout) :: zeta(LBi:,LBj:,:)
344# if defined NESTING && !defined SOLVE3D
345 real(r8), intent(out ) :: DU_flux(LBi:,LBj:)
346 real(r8), intent(out ) :: DV_flux(LBi:,LBj:)
347# endif
348
349#else
350
351# ifdef MASKING
352 real(r8), intent(in ) :: pmask(LBi:UBi,LBj:UBj)
353 real(r8), intent(in ) :: rmask(LBi:UBi,LBj:UBj)
354 real(r8), intent(in ) :: umask(LBi:UBi,LBj:UBj)
355 real(r8), intent(in ) :: vmask(LBi:UBi,LBj:UBj)
356# endif
357# if (defined UV_COR && !defined SOLVE3D) || defined STEP2D_CORIOLIS
358 real(r8), intent(in ) :: fomn(LBi:UBi,LBj:UBj)
359# endif
360# if defined SEDIMENT && defined SED_MORPH
361 real(r8), intent(inout) :: h(LBi:UBi,LBj:UBj)
362# else
363 real(r8), intent(in ) :: h(LBi:UBi,LBj:UBj)
364# endif
365 real(r8), intent(in ) :: om_u(LBi:UBi,LBj:UBj)
366 real(r8), intent(in ) :: om_v(LBi:UBi,LBj:UBj)
367 real(r8), intent(in ) :: on_u(LBi:UBi,LBj:UBj)
368 real(r8), intent(in ) :: on_v(LBi:UBi,LBj:UBj)
369 real(r8), intent(in ) :: omn(LBi:UBi,LBj:UBj)
370 real(r8), intent(in ) :: pm(LBi:UBi,LBj:UBj)
371 real(r8), intent(in ) :: pn(LBi:UBi,LBj:UBj)
372# if defined CURVGRID && defined UV_ADV && !defined SOLVE3D
373 real(r8), intent(in ) :: dndx(LBi:UBi,LBj:UBj)
374 real(r8), intent(in ) :: dmde(LBi:UBi,LBj:UBj)
375# endif
376# if defined UV_VIS2 && !defined SOLVE3D
377 real(r8), intent(in ) :: pmon_r(LBi:UBi,LBj:UBj)
378 real(r8), intent(in ) :: pnom_r(LBi:UBi,LBj:UBj)
379 real(r8), intent(in ) :: pmon_p(LBi:UBi,LBj:UBj)
380 real(r8), intent(in ) :: pnom_p(LBi:UBi,LBj:UBj)
381 real(r8), intent(in ) :: om_r(LBi:UBi,LBj:UBj)
382 real(r8), intent(in ) :: on_r(LBi:UBi,LBj:UBj)
383 real(r8), intent(in ) :: om_p(LBi:UBi,LBj:UBj)
384 real(r8), intent(in ) :: on_p(LBi:UBi,LBj:UBj)
385 real(r8), intent(in ) :: visc2_p(LBi:UBi,LBj:UBj)
386 real(r8), intent(in ) :: visc2_r(LBi:UBi,LBj:UBj)
387# endif
388# if defined SEDIMENT && defined SED_MORPH
389 real(r8), intent(in ) :: bed_thick(LBi:UBi,LBj:UBj,1:3)
390# endif
391# if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
392 real(r8), intent(in ) :: eq_tide(LBi:UBi,LBj:UBj)
393# endif
394# ifndef SOLVE3D
395 real(r8), intent(in ) :: sustr(LBi:UBi,LBj:UBj)
396 real(r8), intent(in ) :: svstr(LBi:UBi,LBj:UBj)
397 real(r8), intent(in ) :: bustr(LBi:UBi,LBj:UBj)
398 real(r8), intent(in ) :: bvstr(LBi:UBi,LBj:UBj)
399# ifdef ATM_PRESS
400 real(r8), intent(in ) :: Pair(LBi:UBi,LBj:UBj)
401# endif
402# else
403# ifdef VAR_RHO_2D
404 real(r8), intent(in ) :: rhoA(LBi:UBi,LBj:UBj)
405 real(r8), intent(in ) :: rhoS(LBi:UBi,LBj:UBj)
406# endif
407 real(r8), intent(inout) :: DU_avg1(LBi:UBi,LBj:UBj)
408 real(r8), intent(inout) :: DU_avg2(LBi:UBi,LBj:UBj)
409 real(r8), intent(inout) :: DV_avg1(LBi:UBi,LBj:UBj)
410 real(r8), intent(inout) :: DV_avg2(LBi:UBi,LBj:UBj)
411 real(r8), intent(inout) :: Zt_avg1(LBi:UBi,LBj:UBj)
412 real(r8), intent(inout) :: rufrc(LBi:UBi,LBj:UBj)
413 real(r8), intent(inout) :: rvfrc(LBi:UBi,LBj:UBj)
414 real(r8), intent(inout) :: rufrc_bak(LBi:UBi,LBj:UBj,2)
415 real(r8), intent(inout) :: rvfrc_bak(LBi:UBi,LBj:UBj,2)
416# endif
417# ifdef WET_DRY
418 real(r8), intent(inout) :: pmask_full(LBi:UBi,LBj:UBj)
419 real(r8), intent(inout) :: rmask_full(LBi:UBi,LBj:UBj)
420 real(r8), intent(inout) :: umask_full(LBi:UBi,LBj:UBj)
421 real(r8), intent(inout) :: vmask_full(LBi:UBi,LBj:UBj)
422
423 real(r8), intent(inout) :: pmask_wet(LBi:UBi,LBj:UBj)
424 real(r8), intent(inout) :: rmask_wet(LBi:UBi,LBj:UBj)
425 real(r8), intent(inout) :: umask_wet(LBi:UBi,LBj:UBj)
426 real(r8), intent(inout) :: vmask_wet(LBi:UBi,LBj:UBj)
427# ifdef SOLVE3D
428 real(r8), intent(inout) :: rmask_wet_avg(LBi:UBi,LBj:UBj)
429# endif
430# endif
431# ifdef DIAGNOSTICS_UV
432 real(r8), intent(inout) :: DiaU2wrk(LBi:UBi,LBj:UBj,NDM2d)
433 real(r8), intent(inout) :: DiaV2wrk(LBi:UBi,LBj:UBj,NDM2d)
434 real(r8), intent(inout) :: DiaRUbar(LBi:UBi,LBj:UBj,2,NDM2d-1)
435 real(r8), intent(inout) :: DiaRVbar(LBi:UBi,LBj:UBj,2,NDM2d-1)
436# ifdef SOLVE3D
437 real(r8), intent(inout) :: DiaU2int(LBi:UBi,LBj:UBj,NDM2d)
438 real(r8), intent(inout) :: DiaV2int(LBi:UBi,LBj:UBj,NDM2d)
439 real(r8), intent(inout) :: DiaRUfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
440 real(r8), intent(inout) :: DiaRVfrc(LBi:UBi,LBj:UBj,3,NDM2d-1)
441# endif
442# endif
443 real(r8), intent(inout) :: ubar(LBi:UBi,LBj:UBj,:)
444 real(r8), intent(inout) :: vbar(LBi:UBi,LBj:UBj,:)
445 real(r8), intent(inout) :: zeta(LBi:UBi,LBj:UBj,:)
446# if defined NESTING && !defined SOLVE3D
447 real(r8), intent(out ) :: DU_flux(LBi:UBi,LBj:UBj)
448 real(r8), intent(out ) :: DV_flux(LBi:UBi,LBj:UBj)
449# endif
450#endif
451!
452! Local variable declarations.
453!
454 integer :: i, is, j
455 integer :: krhs, kbak
456#ifdef DIAGNOSTICS_UV
457 integer :: idiag
458#endif
459!
460 real(r8) :: cff, cff1, cff2, cff3, cff4
461#ifdef WET_DRY
462 real(r8) :: cff5, cff6, cff7
463#endif
464 real(r8) :: fac, fac1, fac2
465!
466#if defined UV_C4ADVECTION && !defined SOLVE3D
467 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dgrad
468#endif
469 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dnew
470 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs
471#if defined UV_VIS2 && !defined SOLVE3D
472 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Drhs_p
473#endif
474 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Dstp
475 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DUon
476 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: DVom
477#if defined STEP2D_CORIOLIS || !defined SOLVE3D
478 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFx
479 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFe
480#endif
481#if !defined SOLVE3D
482 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: UFe
483 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: VFx
484#endif
485#if defined UV_C4ADVECTION && !defined SOLVE3D
486 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
487#endif
488 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rubar
489 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rvbar
490 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzeta
491 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzeta2
492#if defined VAR_RHO_2D && defined SOLVE3D
493 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rzetaSA
494#endif
495 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: urhs
496 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: vrhs
497 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zwrk
498#ifdef DIAGNOSTICS_UV
499 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Uwrk
500 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vwrk
501 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaU2rhs
502 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,NDM2d-1) :: DiaV2rhs
503#endif
504!
505 real(r8), allocatable :: zeta_new(:,:)
506!
507! The following stability limits are obtained empirically using 3/4
508! degree Atlantic model configuration. In all these cases barotropic
509! mode time step is about 180...250 seconds, which is much less than
510! the inertial period. The maximum stability coefficients turned out
511! to be slightly different than predicted by linear theory, although
512! all theoretical tendencies agree with the practice. Note the nearly
513! 70% gain in stability compared with LF-TR for appropriate
514! coefficients (linear theory predicts beta=0.166, epsil=0.84).
515!
516!! real(r8), parameter :: gamma=0.0_r8, &
517!! & beta =0.0_r8, epsil=0.0_r8 !--> Cu=0.818
518!! real(r8), parameter :: gamma=1.0_r8/12.0_r8, &
519!! & beta =0.0_r8, epsil=0.0_r8 !--> Cu=0.878
520!! real(r8), parameter :: gamma=1./12., &
521!! beta =0.1_r8, epsil=0.6_r8 !--> Cu=1.050
522 real(r8), parameter :: gamma=0.0_r8, &
523 & beta =0.14_r8, epsil=0.74_r8 !==> Cu=1.341
524
525#include "set_bounds.h"
526!
527!-----------------------------------------------------------------------
528! Timestep vertically integrated (barotropic) equations.
529!-----------------------------------------------------------------------
530!
531! In the code below it is assumed that variables with time index "krhs"
532! are time-centered at step "n" in barotropic time during predictor
533! sub-step and "n+1/2" during corrector.
534!
535 IF (predictor_2d_step) THEN
536 krhs=kstp
537 ELSE
538 krhs=3
539 END IF
540 IF (first_2d_step) THEN
541 kbak=kstp ! "kbak" is used as "from"
542 ELSE ! time index for LF timestep
543 kbak=3-kstp
544 END IF
545!
546!-----------------------------------------------------------------------
547! Preliminary steps.
548!-----------------------------------------------------------------------
549!
550! Compute total depth of the water column and vertically integrated
551! mass fluxes, which are used in computation of free-surface elevation
552! time tendency and advection terms for the barotropic momentum
553! equations.
554!
555#if defined DISTRIBUTE && !defined NESTING
556# define IR_RANGE IstrUm2-1,Iendp2
557# define JR_RANGE JstrVm2-1,Jendp2
558# define IU_RANGE IstrUm1-1,Iendp2
559# define JU_RANGE Jstrm1-1,Jendp2
560# define IV_RANGE Istrm1-1,Iendp2
561# define JV_RANGE JstrVm1-1,Jendp2
562#else
563# define IR_RANGE IstrUm2-1,Iendp2
564# define JR_RANGE JstrVm2-1,Jendp2
565# define IU_RANGE IstrUm2,Iendp2
566# define JU_RANGE JstrVm2-1,Jendp2
567# define IV_RANGE IstrUm2-1,Iendp2
568# define JV_RANGE JstrVm2,Jendp2
569#endif
570
571 DO j=jr_range
572 DO i=ir_range
573 drhs(i,j)=zeta(i,j,krhs)+h(i,j)
574 END DO
575 END DO
576 DO j=ju_range
577 DO i=iu_range
578 cff=0.5_r8*on_u(i,j)
579 cff1=cff*(drhs(i,j)+drhs(i-1,j))
580 duon(i,j)=ubar(i,j,krhs)*cff1
581 END DO
582 END DO
583 DO j=jv_range
584 DO i=iv_range
585 cff=0.5_r8*om_v(i,j)
586 cff1=cff*(drhs(i,j)+drhs(i,j-1))
587 dvom(i,j)=vbar(i,j,krhs)*cff1
588 END DO
589 END DO
590
591#undef IR_RANGE
592#undef IU_RANGE
593#undef IV_RANGE
594#undef JR_RANGE
595#undef JU_RANGE
596#undef JV_RANGE
597
598#if defined DISTRIBUTE && \
599 defined uv_adv && defined uv_c4advection && !defined SOLVE3D
600!
601! In distributed-memory, the I- and J-ranges are different and a
602! special exchange is done here to avoid having three ghost points
603! for high-order numerical stencils. Notice that a private array is
604! passed below to the exchange routine. It also applies periodic
605! boundary conditions, if appropriate and no partitions in I- or
606! J-directions.
607!
608 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
609 CALL exchange_u2d_tile (ng, tile, &
610 & imins, imaxs, jmins, jmaxs, &
611 & duon)
612 CALL exchange_v2d_tile (ng, tile, &
613 & imins, imaxs, jmins, jmaxs, &
614 & dvom)
615 END IF
616 CALL mp_exchange2d (ng, tile, inlm, 2, &
617 & imins, imaxs, jmins, jmaxs, &
618 & nghostpoints, &
619 & ewperiodic(ng), nsperiodic(ng), &
620 & duon, dvom)
621#endif
622!
623! Set vertically integrated mass fluxes DUon and DVom along the open
624! boundaries in such a way that the integral volume is conserved.
625!
626 IF (any(volcons(:,ng))) THEN
627 CALL set_duv_bc_tile (ng, tile, &
628 & lbi, ubi, lbj, ubj, &
629 & imins, imaxs, jmins, jmaxs, &
630 & krhs, &
631#ifdef MASKING
632 & umask, vmask, &
633#endif
634 & om_v, on_u, &
635 & ubar, vbar, &
636 & drhs, duon, dvom)
637 END IF
638
639#ifdef SOLVE3D
640!
641!-----------------------------------------------------------------------
642! Fields averaged over all barotropic time steps.
643!-----------------------------------------------------------------------
644!
645! Notice that the index ranges here are designed to include physical
646! boundaries only. Periodic ghost points and internal mpi computational
647! margins are NOT included.
648!
649! Reset all barotropic mode time-averaged arrays during the first
650! predictor step. At all subsequent time steps, accumulate averages
651! of the first kind using the DELAYED way. For example, "Zt_avg1" is
652! not summed immediately after the corrector step when computed but
653! during the subsequent predictor substep. It allows saving operations
654! because "DUon" and "DVom" are calculated anyway. The last time step
655! has a special code to add all three barotropic variables after the
656! last corrector substep.
657!
658 IF (predictor_2d_step) THEN ! PREDICTOR STEP
659 IF (first_2d_step) THEN
660 DO j=jstrr,jendr
661 DO i=istrr,iendr
662 zt_avg1(i,j)=0.0_r8
663 du_avg1(i,j)=0.0_r8
664 dv_avg1(i,j)=0.0_r8
665 du_avg2(i,j)=0.0_r8
666 dv_avg2(i,j)=0.0_r8
667 END DO
668 END DO
669 ELSE
670 cff=weight(1,iif(ng)-1,ng)
671 DO j=jstrr,jendr
672 DO i=istrr,iendr
673 zt_avg1(i,j)=zt_avg1(i,j)+cff*zeta(i,j,krhs)
674 IF (i.ge.istr) THEN
675 du_avg1(i,j)=du_avg1(i,j)+cff*duon(i,j)
676 END IF
677 IF (j.ge.jstr) THEN
678 dv_avg1(i,j)=dv_avg1(i,j)+cff*dvom(i,j)
679 END IF
680 END DO
681 END DO
682 END IF
683 ELSE ! CORRECTOR STEP
684 cff=weight(2,iif(ng),ng)
685 DO j=jstrr,jendr
686 DO i=istrr,iendr
687 IF (i.ge.istr) THEN
688 du_avg2(i,j)=du_avg2(i,j)+cff*duon(i,j)
689 END IF
690 IF (j.ge.jstr) THEN
691 dv_avg2(i,j)=dv_avg2(i,j)+cff*dvom(i,j)
692 END IF
693 END DO
694 END DO
695 END IF
696#endif
697!
698!-----------------------------------------------------------------------
699! Advance free-surface.
700!-----------------------------------------------------------------------
701!
702! Notice that the new local free-surface is allocated so it can be
703! passed as an argumment to "zetabc" to avoid memory issues.
704!
705 allocate ( zeta_new(imins:imaxs,jmins:jmaxs) )
706 zeta_new = 0.0_r8
707!
708! Compute "zeta_new" at the new time step and interpolate backward for
709! the subsequent computation of barotropic pressure-gradient terms.
710! Notice that during the predictor of the first 2D step in 3D mode,
711! the pressure gradient terms are computed using just zeta(:,:,kstp),
712! i.e., like in the Forward Euler step, rather than the more accurate
713! predictor of generalized RK2. This is to keep it consistent with the
714! computation of pressure gradient in 3D mode, which uses precisely
715! the initial value of "zeta" rather than the value changed by the
716! first barotropic predictor step. Later in this code, just after
717! "rufrc, rvfrc" are finalized, a correction term based on the
718! difference zeta_new(:,:)-zeta(:,:,kstp) to "rubar, rvbar" to make
719! them consistent with generalized RK2 stepping for pressure gradient
720! terms.
721!
722 IF (predictor_2d_step) THEN
723 IF (first_2d_step) THEN ! Modified RK2 time step (with
724 cff=dtfast(ng) ! Forward-Backward feedback with
725#ifdef SOLVE3D
726 cff1=0.0_r8 !==> Forward Euler
727 cff2=1.0_r8
728#else
729 cff1=0.333333333333_r8 ! optimally chosen beta=1/3 and
730 cff2=0.666666666667_r8 ! epsilon=2/3, see below) is used
731#endif
732 cff3=0.0_r8 ! here for the start up.
733 ELSE
734 cff=2.0_r8*dtfast(ng) ! In the code below "zwrk" is
735 cff1=beta ! time-centered at time step "n"
736 cff2=1.0_r8-2.0_r8*beta ! in the case of LF (for all but
737 cff3=beta ! the first time step)
738 END IF
739!
740 DO j=jstrv-1,jend
741 DO i=istru-1,iend
742 fac=cff*pm(i,j)*pn(i,j)
743 zeta_new(i,j)=zeta(i,j,kbak)+ &
744 & fac*(duon(i,j)-duon(i+1,j)+ &
745 & dvom(i,j)-dvom(i,j+1))
746#ifdef MASKING
747 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
748# ifdef WET_DRY
749 zeta_new(i,j)=zeta_new(i,j)+ &
750 & (dcrit(ng)-h(i,j))*(1.0_r8-rmask(i,j))
751# endif
752#endif
753 dnew(i,j)=zeta_new(i,j)+h(i,j)
754
755 zwrk(i,j)=cff1*zeta_new(i,j)+ &
756 & cff2*zeta(i,j,kstp)+ &
757 & cff3*zeta(i,j,kbak)
758#if defined VAR_RHO_2D && defined SOLVE3D
759 rzeta(i,j)=(1.0_r8+rhos(i,j))*zwrk(i,j)
760 rzeta2(i,j)=rzeta(i,j)*zwrk(i,j)
761 rzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
762#else
763 rzeta(i,j)=zwrk(i,j)
764 rzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
765#endif
766 END DO
767 END DO
768 ELSE !--> CORRECTOR STEP
769 IF (first_2d_step) THEN
770 cff =0.333333333333_r8 ! Modified RK2 weighting:
771 cff1=0.333333333333_r8 ! here "zwrk" is time-
772 cff2=0.333333333333_r8 ! centered at "n+1/2".
773 cff3=0.0_r8
774 ELSE
775 cff =1.0_r8-epsil ! zwrk is always time-
776 cff1=(0.5_r8-gamma)*epsil ! centered at n+1/2
777 cff2=(0.5_r8+2.0_r8*gamma)*epsil ! during corrector sub-
778 cff3=-gamma *epsil ! step.
779 END IF
780!
781 DO j=jstrv-1,jend
782 DO i=istru-1,iend
783 fac=dtfast(ng)*pm(i,j)*pn(i,j)
784 zeta_new(i,j)=zeta(i,j,kstp)+ &
785 & fac*(duon(i,j)-duon(i+1,j)+ &
786 & dvom(i,j)-dvom(i,j+1))
787#ifdef MASKING
788 zeta_new(i,j)=zeta_new(i,j)*rmask(i,j)
789#endif
790 dnew(i,j)=zeta_new(i,j)+h(i,j)
791
792 zwrk(i,j)=cff *zeta(i,j,krhs)+ &
793 & cff1*zeta_new(i,j)+ &
794 & cff2*zeta(i,j,kstp)+ &
795 & cff3*zeta(i,j,kbak)
796#if defined VAR_RHO_2D && defined SOLVE3D
797 rzeta(i,j)=(1.0_r8+rhos(i,j))*zwrk(i,j)
798 rzeta2(i,j)=rzeta(i,j)*zwrk(i,j)
799 rzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
800#else
801 rzeta(i,j)=zwrk(i,j)
802 rzeta2(i,j)=zwrk(i,j)*zwrk(i,j)
803#endif
804 END DO
805 END DO
806 END IF
807!
808! Apply mass point sources (volume vertical influx), if any.
809!
810! Dsrc(is) = 2, flow across grid cell w-face (positive or negative)
811!
812 IF (lwsrc(ng)) THEN
813 DO is=1,nsrc(ng)
814 IF (int(sources(ng)%Dsrc(is)).eq.2) THEN
815 i=sources(ng)%Isrc(is)
816 j=sources(ng)%Jsrc(is)
817 IF (((istrr.le.i).and.(i.le.iendr)).and. &
818 & ((jstrr.le.j).and.(j.le.jendr))) THEN
819 zeta_new(i,j)=zeta_new(i,j)+ &
820 & sources(ng)%Qbar(is)* &
821 & pm(i,j)*pn(i,j)*dtfast(ng)
822 END IF
823 END IF
824 END DO
825 END IF
826!
827! Apply boundary conditions to newly computed free-surface "zeta_new"
828! and load into global state array. Notice that "zeta_new" is always
829! centered at time step "m+1", while zeta(:,:,knew) should be centered
830! either at "m+1/2" after predictor step and at "m+1" after corrector.
831! Chosing it to be this way makes it possible avoid storing RHS for
832! zeta, ubar, and vbar between predictor and corrector sub-steps.
833!
834! Here, we use the local "zetabc" since the private array "zeta_new"
835! is passed as an argument to allow computing the lateral boundary
836! conditions on the range IstrU-1:Iend and JstrV-1:Jend, so parallel
837! tile exchanges are avoided.
838!
839 CALL zetabc_local (ng, tile, &
840 & lbi, ubi, lbj, ubj, &
841 & imins, imaxs, jmins, jmaxs, &
842 & kstp, &
843 & zeta, &
844 & zeta_new)
845!
846 IF (predictor_2d_step) THEN
847 IF (first_2d_step) THEN
848 cff1=0.5_r8
849 cff2=0.5_r8
850 cff3=0.0_r8
851 ELSE
852 cff1=0.5_r8-gamma
853 cff2=0.5_r8+2.0_r8*gamma
854 cff3=-gamma
855 END IF
856 DO j=jstrr,jendr
857 DO i=istrr,iendr
858 zeta(i,j,knew)=cff1*zeta_new(i,j)+ &
859 & cff2*zeta(i,j,kstp)+ &
860 & cff3*zeta(i,j,kbak)
861 END DO
862 END DO
863 ELSE
864 DO j=jstrr,jendr
865 DO i=istrr,iendr
866 zeta(i,j,knew)=zeta_new(i,j)
867 END DO
868 END DO
869 END IF
870!
871!=======================================================================
872! Compute right-hand-side for the 2D momentum equations.
873!=======================================================================
874#ifdef SOLVE3D
875!
876! Notice that we are suppressing the computation of momentum advection,
877! Coriolis, and lateral viscosity terms in 3D Applications because
878! these terms are already included in the baroclinic-to-barotropic
879! forcing arrays "rufrc" and "rvfrc". It does not mean we are entirely
880! omitting them, but it is a choice between recomputing them at every
881! barotropic step or keeping them "frozen" during the fast-time
882! stepping.
883# ifdef STEP2D_CORIOLIS
884! However, in some coarse grid applications with larger baroclinic
885! timestep (say, DT around 20 minutes or larger), adding the Coriolis
886! term in the barotropic equations is useful since f*DT is no longer
887! small.
888# endif
889#endif
890!
891!-----------------------------------------------------------------------
892! Compute pressure-gradient terms.
893!-----------------------------------------------------------------------
894!
895! Notice that "rubar" and "rvbar" are computed within the same to allow
896! shared references to array elements (i,j), which increases the
897! computational density by almost a factor of 1.5 resulting in overall
898! more efficient code.
899!
900 cff1=0.5*g
901 cff2=0.333333333333_r8
902#if !defined SOLVE3D && defined ATM_PRESS
903 fac=0.5_r8*100.0_r8/rho0
904#endif
905 DO j=jstr,jend
906 DO i=istr,iend
907 IF (i.ge.istru) THEN
908 rubar(i,j)=cff1*on_u(i,j)* &
909 & ((h(i-1,j)+ &
910 & h(i ,j))* &
911 & (rzeta(i-1,j)- &
912 & rzeta(i ,j))+ &
913#if defined VAR_RHO_2D && defined SOLVE3D
914 & (h(i-1,j)- &
915 & h(i ,j))* &
916 & (rzetasa(i-1,j)+ &
917 & rzetasa(i ,j)+ &
918 & cff2*(rhoa(i-1,j)- &
919 & rhoa(i ,j))* &
920 & (zwrk(i-1,j)- &
921 & zwrk(i ,j)))+ &
922#endif
923 & (rzeta2(i-1,j)- &
924 & rzeta2(i ,j)))
925#if defined ATM_PRESS && !defined SOLVE3D
926 rubar(i,j)=rubar(i,j)- &
927 & fac*on_u(i,j)* &
928 & (h(i-1,j)+h(i,j)+ &
929 & rzeta(i-1,j)+rzeta(i,j))* &
930 & (pair(i,j)-pair(i-1,j))
931#endif
932#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
933 rubar(i,j)=rubar(i,j)- &
934 & cff1*on_u(i,j)* &
935 & (h(i-1,j)+h(i,j)+ &
936 & rzeta(i-1,j)+rzeta(i,j))* &
937 & (eq_tide(i,j)-eq_tide(i-1,j))
938#endif
939#ifdef DIAGNOSTICS_UV
940 diau2rhs(i,j,m2pgrd)=rubar(i,j)
941#endif
942 END IF
943!
944 IF (j.ge.jstrv) THEN
945 rvbar(i,j)=cff1*om_v(i,j)* &
946 & ((h(i,j-1)+ &
947 & h(i,j ))* &
948 & (rzeta(i,j-1)- &
949 & rzeta(i,j ))+ &
950#if defined VAR_RHO_2D && defined SOLVE3D
951 & (h(i,j-1)- &
952 & h(i,j ))* &
953 & (rzetasa(i,j-1)+ &
954 & rzetasa(i,j )+ &
955 & cff2*(rhoa(i,j-1)- &
956 & rhoa(i,j ))* &
957 & (zwrk(i,j-1)- &
958 & zwrk(i,j )))+ &
959#endif
960 & (rzeta2(i,j-1)- &
961 & rzeta2(i,j )))
962#if defined ATM_PRESS && !defined SOLVE3D
963 rvbar(i,j)=rvbar(i,j)- &
964 & fac*om_v(i,j)* &
965 & (h(i,j-1)+h(i,j)+ &
966 & rzeta(i,j-1)+rzeta(i,j))* &
967 & (pair(i,j)-pair(i,j-1))
968#endif
969#if defined TIDE_GENERATING_FORCES && !defined SOLVE3D
970 rvbar(i,j)=rvbar(i,j)- &
971 & cff1*om_v(i,j)* &
972 & (h(i,j-1)+h(i,j)+ &
973 & rzeta(i,j-1)+rzeta(i,j))* &
974 & (eq_tide(i,j)-eq_tide(i,j-1))
975#endif
976#ifdef DIAGNOSTICS_UV
977 diav2rhs(i,j,m2pgrd)=rvbar(i,j)
978#endif
979 END IF
980 END DO
981 END DO
982
983#if defined UV_ADV && !defined SOLVE3D
984!
985!-----------------------------------------------------------------------
986! Add in horizontal advection of momentum.
987!-----------------------------------------------------------------------
988
989# ifdef UV_C2ADVECTION
990!
991! Second-order, centered differences advection fluxes.
992!
993 DO j=jstr,jend
994 DO i=istru-1,iend
995 ufx(i,j)=0.25_r8* &
996 & (duon(i,j)+duon(i+1,j))* &
997 & (ubar(i ,j,krhs)+ &
998 & ubar(i+1,j,krhs))
999 END DO
1000 END DO
1001!
1002 DO j=jstr,jend+1
1003 DO i=istru,iend
1004 ufe(i,j)=0.25_r8* &
1005 & (dvom(i,j)+dvom(i-1,j))* &
1006 & (ubar(i,j ,krhs)+ &
1007 & ubar(i,j-1,krhs))
1008 END DO
1009 END DO
1010!
1011 DO j=jstrv,jend
1012 DO i=istr,iend+1
1013 vfx(i,j)=0.25_r8* &
1014 & (duon(i,j)+duon(i,j-1))* &
1015 & (vbar(i ,j,krhs)+ &
1016 & vbar(i-1,j,krhs))
1017 END DO
1018 END DO
1019!
1020 DO j=jstrv-1,jend
1021 DO i=istr,iend
1022 vfe(i,j)=0.25_r8* &
1023 & (dvom(i,j)+dvom(i,j+1))* &
1024 & (vbar(i,j ,krhs)+ &
1025 & vbar(i,j+1,krhs))
1026 END DO
1027 END DO
1028
1029# elif defined UV_C4ADVECTION
1030!
1031! Fourth-order, centered differences u-momentum advection fluxes.
1032!
1033 DO j=jstr,jend
1034 DO i=istrum1,iendp1
1035 grad(i,j)=ubar(i-1,j,krhs)-2.0_r8*ubar(i,j,krhs)+ &
1036 & ubar(i+1,j,krhs)
1037 dgrad(i,j)=duon(i-1,j)-2.0_r8*duon(i,j)+duon(i+1,j)
1038 END DO
1039 END DO
1040 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1041 IF (domain(ng)%Western_Edge(tile)) THEN
1042 DO j=jstr,jend
1043 grad(istr,j)=grad(istr+1,j)
1044 dgrad(istr,j)=dgrad(istr+1,j)
1045 END DO
1046 END IF
1047 END IF
1048 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1049 IF (domain(ng)%Eastern_Edge(tile)) THEN
1050 DO j=jstr,jend
1051 grad(iend+1,j)=grad(iend,j)
1052 dgrad(iend+1,j)=dgrad(iend,j)
1053 END DO
1054 END IF
1055 END IF
1056! d/dx(Duu/n)
1057 cff=1.0_r8/6.0_r8
1058 DO j=jstr,jend
1059 DO i=istru-1,iend
1060 ufx(i,j)=0.25_r8*(ubar(i ,j,krhs)+ &
1061 & ubar(i+1,j,krhs)- &
1062 & cff*(grad(i,j)+grad(i+1,j)))* &
1063 & (duon(i,j)+duon(i+1,j)- &
1064 & cff*(dgrad(i,j)+dgrad(i+1,j)))
1065 END DO
1066 END DO
1067!
1068 DO j=jstrm1,jendp1
1069 DO i=istru,iend
1070 grad(i,j)=ubar(i,j-1,krhs)-2.0_r8*ubar(i,j,krhs)+ &
1071 & ubar(i,j+1,krhs)
1072 END DO
1073 END DO
1074!
1075 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1076 IF (domain(ng)%Southern_Edge(tile)) THEN
1077 DO i=istru,iend
1078 grad(i,jstr-1)=grad(i,jstr)
1079 END DO
1080 END IF
1081 END IF
1082 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1083 IF (domain(ng)%Northern_Edge(tile)) THEN
1084 DO i=istru,iend
1085 grad(i,jend+1)=grad(i,jend)
1086 END DO
1087 END IF
1088 END IF
1089 DO j=jstr,jend+1
1090 DO i=istru-1,iend
1091 dgrad(i,j)=dvom(i-1,j)-2.0_r8*dvom(i,j)+dvom(i+1,j)
1092 END DO
1093 END DO
1094! d/dy(Duv/m)
1095 cff=1.0_r8/6.0_r8
1096 DO j=jstr,jend+1
1097 DO i=istru,iend
1098 ufe(i,j)=0.25_r8*(ubar(i,j ,krhs)+ &
1099 & ubar(i,j-1,krhs)- &
1100 & cff*(grad(i,j)+grad(i,j-1)))* &
1101 & (dvom(i,j)+dvom(i-1,j)- &
1102 & cff*(dgrad(i,j)+dgrad(i-1,j)))
1103 END DO
1104 END DO
1105!
1106! Fourth-order, centered differences v-momentum advection fluxes.
1107!
1108 DO j=jstrv,jend
1109 DO i=istrm1,iendp1
1110 grad(i,j)=vbar(i-1,j,krhs)-2.0_r8*vbar(i,j,krhs)+ &
1111 & vbar(i+1,j,krhs)
1112 END DO
1113 END DO
1114!
1115 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1116 IF (domain(ng)%Western_Edge(tile)) THEN
1117 DO j=jstrv,jend
1118 grad(istr-1,j)=grad(istr,j)
1119 END DO
1120 END IF
1121 END IF
1122 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1123 IF (domain(ng)%Eastern_Edge(tile)) THEN
1124 DO j=jstrv,jend
1125 grad(iend+1,j)=grad(iend,j)
1126 END DO
1127 END IF
1128 END IF
1129 DO j=jstrv-1,jend
1130 DO i=istr,iend+1
1131 dgrad(i,j)=duon(i,j-1)-2.0_r8*duon(i,j)+duon(i,j+1)
1132 END DO
1133 END DO
1134! d/dx(Duv/n)
1135 cff=1.0_r8/6.0_r8
1136 DO j=jstrv,jend
1137 DO i=istr,iend+1
1138 vfx(i,j)=0.25_r8*(vbar(i ,j,krhs)+ &
1139 & vbar(i-1,j,krhs)- &
1140 & cff*(grad(i,j)+grad(i-1,j)))* &
1141 & (duon(i,j)+duon(i,j-1)- &
1142 & cff*(dgrad(i,j)+dgrad(i,j-1)))
1143 END DO
1144 END DO
1145!
1146 DO j=jstrvm1,jendp1
1147 DO i=istr,iend
1148 grad(i,j)=vbar(i,j-1,krhs)-2.0_r8*vbar(i,j,krhs)+ &
1149 & vbar(i,j+1,krhs)
1150 dgrad(i,j)=dvom(i,j-1)-2.0_r8*dvom(i,j)+dvom(i,j+1)
1151 END DO
1152 END DO
1153 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1154 IF (domain(ng)%Southern_Edge(tile)) THEN
1155 DO i=istr,iend
1156 grad(i,jstr)=grad(i,jstr+1)
1157 dgrad(i,jstr)=dgrad(i,jstr+1)
1158 END DO
1159 END IF
1160 END IF
1161 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1162 IF (domain(ng)%Northern_Edge(tile)) THEN
1163 DO i=istr,iend
1164 grad(i,jend+1)=grad(i,jend)
1165 dgrad(i,jend+1)=dgrad(i,jend)
1166 END DO
1167 END IF
1168 END IF
1169! d/dy(Dvv/m)
1170 cff=1.0_r8/6.0_r8
1171 DO j=jstrv-1,jend
1172 DO i=istr,iend
1173 vfe(i,j)=0.25_r8*(vbar(i,j ,krhs)+ &
1174 & vbar(i,j+1,krhs)- &
1175 & cff*(grad(i,j)+grad(i,j+1)))* &
1176 & (dvom(i,j)+dvom(i,j+1)- &
1177 & cff*(dgrad(i,j)+dgrad(i,j+1)))
1178 END DO
1179 END DO
1180# endif
1181!
1182! Add advection to RHS terms.
1183!
1184 DO j=jstr,jend
1185 DO i=istr,iend
1186 IF (i.ge.istru) THEN
1187 cff1=ufx(i,j)-ufx(i-1,j)
1188 cff2=ufe(i,j+1)-ufe(i,j)
1189 fac=cff1+cff2
1190 rubar(i,j)=rubar(i,j)-fac
1191# if defined DIAGNOSTICS_UV
1192 diau2rhs(i,j,m2xadv)=-cff1
1193 diau2rhs(i,j,m2yadv)=-cff2
1194 diau2rhs(i,j,m2hadv)=-fac
1195# endif
1196 END IF
1197!
1198 IF (j.ge.jstrv) THEN
1199 cff1=vfx(i+1,j)-vfx(i,j)
1200 cff2=vfe(i,j)-vfe(i,j-1)
1201 fac=cff1+cff2
1202 rvbar(i,j)=rvbar(i,j)-fac
1203# if defined DIAGNOSTICS_UV
1204 diav2rhs(i,j,m2xadv)=-cff1
1205 diav2rhs(i,j,m2yadv)=-cff2
1206 diav2rhs(i,j,m2hadv)=-fac
1207# endif
1208 END IF
1209 END DO
1210 END DO
1211#endif
1212
1213#if (defined UV_COR & !defined SOLVE3D) || defined STEP2D_CORIOLIS
1214!
1215!-----------------------------------------------------------------------
1216! Add in Coriolis term.
1217!-----------------------------------------------------------------------
1218!
1219 DO j=jstrv-1,jend
1220 DO i=istru-1,iend
1221 cff=0.5_r8*drhs(i,j)*fomn(i,j)
1222 ufx(i,j)=cff*(vbar(i,j ,krhs)+ &
1223 & vbar(i,j+1,krhs))
1224 vfe(i,j)=cff*(ubar(i ,j,krhs)+ &
1225 & ubar(i+1,j,krhs))
1226 END DO
1227 END DO
1228!
1229 DO j=jstr,jend
1230 DO i=istr,iend
1231 IF (i.ge.istru) THEN
1232 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1233 rubar(i,j)=rubar(i,j)+fac1
1234# if defined DIAGNOSTICS_UV
1235 diau2rhs(i,j,m2fcor)=fac1
1236# endif
1237 END IF
1238!
1239 IF (j.ge.jstrv) THEN
1240 fac2=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1241 rvbar(i,j)=rvbar(i,j)-fac2
1242# if defined DIAGNOSTICS_UV
1243 diav2rhs(i,j,m2fcor)=-fac2
1244# endif
1245 END IF
1246 END DO
1247 END DO
1248#endif
1249
1250#if (defined CURVGRID && defined UV_ADV) && !defined SOLVE3D
1251!
1252!-----------------------------------------------------------------------
1253! Add in curvilinear transformation terms.
1254!-----------------------------------------------------------------------
1255!
1256 DO j=jstrv-1,jend
1257 DO i=istru-1,iend
1258 cff1=0.5_r8*(vbar(i,j ,krhs)+ &
1259 & vbar(i,j+1,krhs))
1260 cff2=0.5_r8*(ubar(i ,j,krhs)+ &
1261 & ubar(i+1,j,krhs))
1262 cff3=cff1*dndx(i,j)
1263 cff4=cff2*dmde(i,j)
1264 cff=drhs(i,j)*(cff3-cff4)
1265 ufx(i,j)=cff*cff1
1266 vfe(i,j)=cff*cff2
1267# if defined DIAGNOSTICS_UV
1268 cff=drhs(i,j)*cff4
1269 uwrk(i,j)=-cff*cff1 ! ubar equation, ETA-term
1270 vwrk(i,j)=-cff*cff2 ! vbar equation, ETA-term
1271# endif
1272 END DO
1273 END DO
1274!
1275 DO j=jstr,jend
1276 DO i=istr,iend
1277 IF (i.ge.istru) THEN
1278 fac1=0.5_r8*(ufx(i,j)+ufx(i-1,j))
1279 rubar(i,j)=rubar(i,j)+fac1
1280# if defined DIAGNOSTICS_UV
1281 fac2=0.5_r8*(uwrk(i,j)+uwrk(i-1,j))
1282 diau2rhs(i,j,m2xadv)=diau2rhs(i,j,m2xadv)+fac1-fac2
1283 diau2rhs(i,j,m2yadv)=diau2rhs(i,j,m2yadv)+fac2
1284 diau2rhs(i,j,m2hadv)=diau2rhs(i,j,m2hadv)+fac1
1285# endif
1286 END IF
1287!
1288 IF (j.ge.jstrv) THEN
1289 fac1=0.5_r8*(vfe(i,j)+vfe(i,j-1))
1290 rvbar(i,j)=rvbar(i,j)-fac1
1291# if defined DIAGNOSTICS_UV
1292 fac2=0.5_r8*(vwrk(i,j)+vwrk(i,j-1))
1293 diav2rhs(i,j,m2xadv)=diav2rhs(i,j,m2xadv)-fac1+fac2
1294 diav2rhs(i,j,m2yadv)=diav2rhs(i,j,m2yadv)-fac2
1295 diav2rhs(i,j,m2hadv)=diav2rhs(i,j,m2hadv)-fac1
1296# endif
1297 END IF
1298 END DO
1299 END DO
1300#endif
1301
1302#if defined UV_VIS2 && !defined SOLVE3D
1303!
1304!-----------------------------------------------------------------------
1305! Add in horizontal harmonic viscosity.
1306!-----------------------------------------------------------------------
1307!
1308! Compute total depth at PSI-points.
1309!
1310 DO j=jstr,jend+1
1311 DO i=istr,iend+1
1312 drhs_p(i,j)=0.25_r8*(drhs(i,j )+drhs(i-1,j )+ &
1313 & drhs(i,j-1)+drhs(i-1,j-1))
1314 END DO
1315 END DO
1316!
1317! Compute flux-components of the horizontal divergence of the stress
1318! tensor (m5/s2) in XI- and ETA-directions.
1319!
1320 DO j=jstrv-1,jend
1321 DO i=istru-1,iend
1322 cff=visc2_r(i,j)*drhs(i,j)*0.5_r8* &
1323 & (pmon_r(i,j)* &
1324 & ((pn(i ,j)+pn(i+1,j))*ubar(i+1,j,krhs)- &
1325 & (pn(i-1,j)+pn(i ,j))*ubar(i ,j,krhs))- &
1326 & pnom_r(i,j)* &
1327 & ((pm(i,j )+pm(i,j+1))*vbar(i,j+1,krhs)- &
1328 & (pm(i,j-1)+pm(i,j ))*vbar(i,j ,krhs)))
1329 ufx(i,j)=on_r(i,j)*on_r(i,j)*cff
1330 vfe(i,j)=om_r(i,j)*om_r(i,j)*cff
1331 END DO
1332 END DO
1333!
1334 DO j=jstr,jend+1
1335 DO i=istr,iend+1
1336 cff=visc2_p(i,j)*drhs_p(i,j)*0.5_r8* &
1337 & (pmon_p(i,j)* &
1338 & ((pn(i ,j-1)+pn(i ,j))*vbar(i ,j,krhs)- &
1339 & (pn(i-1,j-1)+pn(i-1,j))*vbar(i-1,j,krhs))+ &
1340 & pnom_p(i,j)* &
1341 & ((pm(i-1,j )+pm(i,j ))*ubar(i,j ,krhs)- &
1342 & (pm(i-1,j-1)+pm(i,j-1))*ubar(i,j-1,krhs)))
1343# ifdef MASKING
1344 cff=cff*pmask(i,j)
1345# endif
1346# ifdef WET_DRY
1347 cff=cff*pmask_wet(i,j)
1348# endif
1349 ufe(i,j)=om_p(i,j)*om_p(i,j)*cff
1350 vfx(i,j)=on_p(i,j)*on_p(i,j)*cff
1351 END DO
1352 END DO
1353!
1354! Add in harmonic viscosity.
1355!
1356 DO j=jstr,jend
1357 DO i=istr,iend
1358 IF (i.ge.istru) THEN
1359 cff1=0.5_r8*(pn(i-1,j)+pn(i,j))*(ufx(i,j )-ufx(i-1,j))
1360 cff2=0.5_r8*(pm(i-1,j)+pm(i,j))*(ufe(i,j+1)-ufe(i ,j))
1361 fac=cff1+cff2
1362 rubar(i,j)=rubar(i,j)+fac
1363# if defined DIAGNOSTICS_UV
1364 diau2rhs(i,j,m2hvis)=fac
1365 diau2rhs(i,j,m2xvis)=cff1
1366 diau2rhs(i,j,m2yvis)=cff2
1367# endif
1368 END IF
1369!
1370 IF (j.ge.jstrv) THEN
1371 cff1=0.5_r8*(pn(i,j-1)+pn(i,j))*(vfx(i+1,j)-vfx(i,j ))
1372 cff2=0.5_r8*(pm(i,j-1)+pm(i,j))*(vfe(i ,j)-vfe(i,j-1))
1373 fac=cff1-cff2
1374 rvbar(i,j)=rvbar(i,j)+fac
1375# if defined DIAGNOSTICS_UV
1376 diav2rhs(i,j,m2hvis)=fac
1377 diav2rhs(i,j,m2xvis)= cff1
1378 diav2rhs(i,j,m2yvis)=-cff2
1379# endif
1380 END IF
1381 END DO
1382 END DO
1383#endif
1384
1385#ifndef SOLVE3D
1386!
1387!-----------------------------------------------------------------------
1388! Add in bottom stress.
1389!-----------------------------------------------------------------------
1390!
1391 DO j=jstr,jend
1392 DO i=istru,iend
1393 fac=bustr(i,j)*om_u(i,j)*on_u(i,j)
1394 rubar(i,j)=rubar(i,j)-fac
1395# ifdef DIAGNOSTICS_UV
1396 diau2rhs(i,j,m2bstr)=-fac
1397# endif
1398 END DO
1399 END DO
1400 DO j=jstrv,jend
1401 DO i=istr,iend
1402 fac=bvstr(i,j)*om_v(i,j)*on_v(i,j)
1403 rvbar(i,j)=rvbar(i,j)-fac
1404# ifdef DIAGNOSTICS_UV
1405 diav2rhs(i,j,m2bstr)=-fac
1406# endif
1407 END DO
1408 END DO
1409#else
1410# ifdef DIAGNOSTICS_UV
1411!
1412! Initialize the stress term if no bottom friction is defined.
1413!
1414 DO j=jstr,jend
1415 DO i=istru,iend
1416 diau2rhs(i,j,m2bstr)=0.0_r8
1417 END DO
1418 END DO
1419 DO j=jstrv,jend
1420 DO i=istr,iend
1421 diav2rhs(i,j,m2bstr)=0.0_r8
1422 END DO
1423 END DO
1424# endif
1425#endif
1426
1427#ifdef SOLVE3D
1428!
1429!-----------------------------------------------------------------------
1430! Coupling between 2D and 3D equations.
1431!-----------------------------------------------------------------------
1432!
1433! Before the predictor step of the first barotropic time step, arrays
1434! "rufrc" and "rvfrc" contain vertical integrals of the 3D RHS terms
1435! for the momentum equations (including surface and bottom stresses,
1436! if so prescribed). During the first barotropic time step, convert
1437! them into forcing terms by subtracting the fast-time "rubar" and
1438! "rvbar" from them.
1439!
1440! These forcing terms are then extrapolated forward in time using
1441! optimized Adams-Bashforth weights, so that the resultant "rufrc"
1442! and "rvfrc" are centered effectively at time n+1/2 in baroclinic
1443! time.
1444!
1445! From now on, these newly computed forcing terms remain unchanged
1446! during the fast time stepping and will be added to "rubar" and
1447! "rvbar" during all subsequent barotropic time steps.
1448!
1449! Thus, the algorithm below is designed for coupling during the 3D
1450! predictor sub-step. The forcing terms "rufrc" and "rvfrc" are
1451! computed as instantaneous values at 3D time index "nstp" first and
1452! then extrapolated half-step forward using AM3-like weights optimized
1453! for maximum stability (with particular care for startup).
1454!
1455 IF (first_2d_step.and.predictor_2d_step) THEN
1456 IF (first_time_step) THEN
1457 cff3=0.0_r8
1458 cff2=0.0_r8
1459 cff1=1.0_r8
1460 ELSE IF (first_time_step+1) THEN
1461 cff3=0.0_r8
1462 cff2=-0.5_r8
1463 cff1=1.5_r8
1464 ELSE
1465 cff3=0.281105_r8
1466 cff2=-0.5_r8-2.0_r8*cff3
1467 cff1=1.5_r8+cff3
1468 END IF
1469!
1470 DO j=jstr,jend
1471 DO i=istru,iend
1472 cff=rufrc(i,j)-rubar(i,j)
1473 rufrc(i,j)=cff1*cff+ &
1474 & cff2*rufrc_bak(i,j,3-nstp)+ &
1475 & cff3*rufrc_bak(i,j,nstp )
1476 rufrc_bak(i,j,nstp)=cff
1477 END DO
1478 END DO
1479 DO j=jstrv,jend
1480 DO i=istr,iend
1481 cff=rvfrc(i,j)-rvbar(i,j)
1482 rvfrc(i,j)=cff1*cff+ &
1483 & cff2*rvfrc_bak(i,j,3-nstp)+ &
1484 & cff3*rvfrc_bak(i,j,nstp )
1485 rvfrc_bak(i,j,nstp)=cff
1486 END DO
1487 END DO
1488!
1489! Since coupling requires that the pressure gradient term is computed
1490! using zeta(:,:,kstp) instead of 1/3 toward zeta_new(:,:) as needed
1491! by generalized RK2 scheme, apply compensation to shift pressure
1492! gradient terms from "kstp" to 1/3 toward "knew".
1493!
1494 cff1=0.5_r8*g
1495 cff2=0.333333333333_r8
1496 cff3=1.666666666666_r8
1497
1498 DO j=jstrv-1,jend
1499 DO i=istru-1,iend
1500 zwrk(i,j)=cff2*(zeta_new(i,j)-zeta(i,j,kstp))
1501# if defined VAR_RHO_2D && defined SOLVE3D
1502 rzeta(i,j)=(1.0_r8+rhos(i,j))*zwrk(i,j)
1503 rzeta2(i,j)=rzeta(i,j)* &
1504 & (cff2*zeta_new(i,j)+ &
1505 & cff3*zeta(i,j,kstp))
1506 rzetasa(i,j)=zwrk(i,j)*(rhos(i,j)-rhoa(i,j))
1507# else
1508 rzeta(i,j)=zwrk(i,j)
1509 rzeta2(i,j)=zwrk(i,j)* &
1510 & (cff2*zeta_new(i,j)+ &
1511 & cff3*zeta(i,j,kstp))
1512# endif
1513 END DO
1514 END DO
1515!
1516 DO j=jstr,jend
1517 DO i=istr,iend
1518 IF (i.ge.istru) THEN
1519 rubar(i,j)=rubar(i,j)+ &
1520 & cff1*on_u(i,j)* &
1521 & ((h(i-1,j)+ &
1522 & h(i ,j))* &
1523 & (rzeta(i-1,j)- &
1524 & rzeta(i ,j))+ &
1525# if defined VAR_RHO_2D && defined SOLVE3D
1526 & (h(i-1,j)- &
1527 & h(i ,j))* &
1528 & (rzetasa(i-1,j)+ &
1529 & rzetasa(i ,j)+ &
1530 & cff2*(rhoa(i-1,j)- &
1531 & rhoa(i ,j))* &
1532 & (zwrk(i-1,j)- &
1533 & zwrk(i ,j)))+ &
1534# endif
1535 & (rzeta2(i-1,j)- &
1536 & rzeta2(i ,j)))
1537# ifdef DIAGNOSTICS_UV
1538 diau2rhs(i,j,m2pgrd)=diau2rhs(i,j,m2pgrd)+ &
1539 & rubar(i,j)
1540# endif
1541 END IF
1542!
1543 IF (j.ge.jstrv) THEN
1544 rvbar(i,j)=rvbar(i,j)+ &
1545 & cff1*om_v(i,j)* &
1546 & ((h(i,j-1)+ &
1547 & h(i,j ))* &
1548 & (rzeta(i,j-1)- &
1549 & rzeta(i,j ))+ &
1550# if defined VAR_RHO_2D && defined SOLVE3D
1551 & (h(i,j-1)- &
1552 & h(i,j ))* &
1553 & (rzetasa(i,j-1)+ &
1554 & rzetasa(i,j )+ &
1555 & cff2*(rhoa(i,j-1)- &
1556 & rhoa(i,j ))* &
1557 & (zwrk(i,j-1)- &
1558 & zwrk(i,j )))+ &
1559# endif
1560 & (rzeta2(i,j-1)- &
1561 & rzeta2(i,j )))
1562# ifdef DIAGNOSTICS_UV
1563 diav2rhs(i,j,m2pgrd)=diav2rhs(i,j,m2pgrd)+ &
1564 & rvbar(i,j)
1565# endif
1566 END IF
1567 END DO
1568 END DO
1569 END IF
1570#endif
1571!
1572!=======================================================================
1573! Time step 2D momentum equations.
1574!=======================================================================
1575!
1576! Compute total water column depth.
1577!
1578 IF (first_2d_step.or.(.not.predictor_2d_step)) THEN
1579 DO j=jstrv-1,jend
1580 DO i=istru-1,iend
1581 dstp(i,j)=h(i,j)+zeta(i,j,kstp)
1582 END DO
1583 END DO
1584 ELSE
1585 DO j=jstrv-1,jend
1586 DO i=istru-1,iend
1587 dstp(i,j)=h(i,j)+zeta(i,j,kbak)
1588 END DO
1589 END DO
1590 END IF
1591!
1592! During the predictor sub-step, once newly computed "ubar" and "vbar"
1593! become available, interpolate them half-step backward in barotropic
1594! time (i.e., they end up time-centered at n+1/2) in order to use it
1595! during subsequent corrector sub-step.
1596!
1597 IF (predictor_2d_step) THEN
1598 IF (first_2d_step) THEN
1599 cff1=0.5_r8*dtfast(ng)
1600 cff2=0.5_r8
1601 cff3=0.5_r8
1602 cff4=0.0_r8
1603 ELSE
1604 cff1=dtfast(ng)
1605 cff2=0.5_r8-gamma
1606 cff3=0.5_r8+2.0_r8*gamma
1607 cff4=-gamma
1608 ENDIF
1609!
1610 DO j=jstr,jend
1611 DO i=istru,iend
1612 cff=cff1*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
1613 fac1=1.0_r8/(dnew(i,j)+dnew(i-1,j))
1614 ubar(i,j,knew)=fac1* &
1615 & (ubar(i,j,kbak)* &
1616 & (dstp(i,j)+dstp(i-1,j))+ &
1617#ifdef SOLVE3D
1618 & cff*(rubar(i,j)+rufrc(i,j)))
1619#else
1620 & cff*rubar(i,j)+4.0_r8*cff1*sustr(i,j))
1621#endif
1622#ifdef MASKING
1623 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
1624#endif
1625 ubar(i,j,knew)=cff2*ubar(i,j,knew)+ &
1626 & cff3*ubar(i,j,kstp)+ &
1627 & cff4*ubar(i,j,kbak)
1628#ifdef WET_DRY
1629 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
1630 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
1631 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
1632 ubar(i,j,knew)=ubar(i,j,knew)*cff7
1633#endif
1634#if defined NESTING && !defined SOLVE3D
1635 du_flux(i,j)=0.5_r8*on_u(i,j)* &
1636 & (dnew(i,j)+dnew(i-1,j))*ubar(i,j,knew)
1637#endif
1638 END DO
1639 END DO
1640!
1641 DO j=jstrv,jend
1642 DO i=istr,iend
1643 cff=cff1*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
1644 fac2=1.0_r8/(dnew(i,j)+dnew(i,j-1))
1645 vbar(i,j,knew)=fac2* &
1646 & (vbar(i,j,kbak)* &
1647 & (dstp(i,j)+dstp(i,j-1))+ &
1648#ifdef SOLVE3D
1649 & cff*(rvbar(i,j)+rvfrc(i,j)))
1650#else
1651 & cff*rvbar(i,j)+4.0_r8*cff1*svstr(i,j))
1652#endif
1653#ifdef MASKING
1654 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
1655#endif
1656 vbar(i,j,knew)=cff2*vbar(i,j,knew)+ &
1657 & cff3*vbar(i,j,kstp)+ &
1658 & cff4*vbar(i,j,kbak)
1659#ifdef WET_DRY
1660 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
1661 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
1662 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
1663 vbar(i,j,knew)=vbar(i,j,knew)*cff7
1664#endif
1665#if defined NESTING && !defined SOLVE3D
1666 dv_flux(i,j)=0.5_r8*om_v(i,j)* &
1667 & (dnew(i,j)+dnew(i,j-1))*vbar(i,j,knew)
1668#endif
1669 END DO
1670 END DO
1671
1672 ELSE !--> CORRECTOR_2D_STEP
1673
1674 cff1=0.5_r8*dtfast(ng)
1675 DO j=jstr,jend
1676 DO i=istru,iend
1677 cff=cff1*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
1678 fac1=1.0_r8/(dnew(i,j)+dnew(i-1,j))
1679 ubar(i,j,knew)=fac1* &
1680 & (ubar(i,j,kstp)* &
1681 & (dstp(i,j)+dstp(i-1,j))+ &
1682#ifdef SOLVE3D
1683 & cff*(rubar(i,j)+rufrc(i,j)))
1684#else
1685 & cff*rubar(i,j)+4.0_r8*cff1*sustr(i,j))
1686#endif
1687#ifdef MASKING
1688 ubar(i,j,knew)=ubar(i,j,knew)*umask(i,j)
1689#endif
1690#ifdef WET_DRY
1691 cff5=abs(abs(umask_wet(i,j))-1.0_r8)
1692 cff6=0.5_r8+dsign(0.5_r8,ubar(i,j,knew))*umask_wet(i,j)
1693 cff7=0.5_r8*umask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
1694 ubar(i,j,knew)=ubar(i,j,knew)*cff7
1695#endif
1696#if defined NESTING && !defined SOLVE3D
1697 du_flux(i,j)=0.5_r8*on_u(i,j)* &
1698 & (dnew(i,j)+dnew(i-1,j))*ubar(i,j,knew)
1699#endif
1700 END DO
1701 END DO
1702!
1703 DO j=jstrv,jend
1704 DO i=istr,iend
1705 cff=cff1*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
1706 fac2=1.0_r8/(dnew(i,j)+dnew(i,j-1))
1707 vbar(i,j,knew)=fac2* &
1708 & (vbar(i,j,kstp)* &
1709 & (dstp(i,j)+dstp(i,j-1))+ &
1710#ifdef SOLVE3D
1711 & cff*(rvbar(i,j)+rvfrc(i,j)))
1712#else
1713 & cff*rvbar(i,j)+4.0_r8*cff1*svstr(i,j))
1714#endif
1715#ifdef MASKING
1716 vbar(i,j,knew)=vbar(i,j,knew)*vmask(i,j)
1717#endif
1718#ifdef WET_DRY
1719 cff5=abs(abs(vmask_wet(i,j))-1.0_r8)
1720 cff6=0.5_r8+dsign(0.5_r8,vbar(i,j,knew))*vmask_wet(i,j)
1721 cff7=0.5_r8*vmask_wet(i,j)*cff5+cff6*(1.0_r8-cff5)
1722 vbar(i,j,knew)=vbar(i,j,knew)*cff7
1723#endif
1724#if defined NESTING && !defined SOLVE3D
1725 dv_flux(i,j)=0.5_r8*om_v(i,j)* &
1726 & (dnew(i,j)+dnew(i,j-1))*vbar(i,j,knew)
1727#endif
1728 END DO
1729 END DO
1730 END IF
1731!
1732! Apply lateral boundary conditions.
1733!
1734 CALL u2dbc_tile (ng, tile, &
1735 & lbi, ubi, lbj, ubj, &
1736 & imins, imaxs, jmins, jmaxs, &
1737 & krhs, kstp, knew, &
1738 & ubar, vbar, zeta)
1739 CALL v2dbc_tile (ng, tile, &
1740 & lbi, ubi, lbj, ubj, &
1741 & imins, imaxs, jmins, jmaxs, &
1742 & krhs, kstp, knew, &
1743 & ubar, vbar, zeta)
1744!
1745! Compute integral mass flux across open boundaries and adjust
1746! for volume conservation.
1747!
1748 IF (any(volcons(:,ng))) THEN
1749 CALL obc_flux_tile (ng, tile, &
1750 & lbi, ubi, lbj, ubj, &
1751 & imins, imaxs, jmins, jmaxs, &
1752 & knew, &
1753#ifdef MASKING
1754 & umask, vmask, &
1755#endif
1756 & h, om_v, on_u, &
1757 & ubar, vbar, zeta)
1758 END IF
1759
1760#if defined NESTING && !defined SOLVE3D
1761!
1762! Set barotropic fluxes along physical boundaries.
1763!
1764 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1765 IF (domain(ng)%Western_Edge(tile)) THEN
1766 DO j=jstr-1,jendr
1767 dnew(istr-1,j)=h(istr-1,j)+zeta_new(istr-1,j)
1768 END DO
1769 END IF
1770 END IF
1771 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1772 IF (domain(ng)%Eastern_Edge(tile)) THEN
1773 DO j=jstr-1,jendr
1774 dnew(iend+1,j)=h(iend+1,j)+zeta_new(iend+1,j)
1775 END DO
1776 END IF
1777 END IF
1778 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1779 IF (domain(ng)%Southern_Edge(tile)) THEN
1780 DO i=istr-1,iendr
1781 dnew(i,jstr-1)=h(i,jstr-1)+zeta_new(i,jstr-1)
1782 END DO
1783 END IF
1784 END IF
1785 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1786 IF (domain(ng)%Northern_Edge(tile)) THEN
1787 DO i=istr-1,iendr
1788 dnew(i,jend+1)=h(i,jend+1)+zeta_new(i,jend+1)
1789 END DO
1790 END IF
1791 END IF
1792!
1793 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1794 IF (domain(ng)%Western_Edge(tile)) THEN
1795 DO j=jstrr,jendr
1796 du_flux(istru-1,j)=0.5_r8*on_u(istru-1,j)* &
1797 & (dnew(istru-1,j)+dnew(istru-2,j))* &
1798 & ubar(istru-1,j,knew)
1799 END DO
1800 DO j=jstrv,jend
1801 dv_flux(istr-1,j)=0.5_r8*om_v(istr-1,j)* &
1802 & (dnew(istr-1,j)+dnew(istr-1,j-1))* &
1803 & vbar(istr-1,j,knew)
1804 END DO
1805 END IF
1806 END IF
1807 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1808 IF (domain(ng)%Eastern_Edge(tile)) THEN
1809 DO j=jstrr,jendr
1810 du_flux(iend+1,j)=0.5_r8*on_u(iend+1,j)* &
1811 & (dnew(iend+1,j)+dnew(iend,j))* &
1812 & ubar(iend+1,j,knew)
1813 END DO
1814 DO j=jstrv,jend
1815 dv_flux(iend+1,j)=0.5_r8*om_v(iend+1,j)* &
1816 & (dnew(iend+1,j)+dnew(iend+1,j-1))* &
1817 & vbar(iend+1,j,knew)
1818 END DO
1819 END IF
1820 END IF
1821 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1822 IF (domain(ng)%Southern_Edge(tile)) THEN
1823 DO i=istru,iend
1824 du_flux(i,jstr-1)=0.5_r8*on_u(i,jstr-1)* &
1825 & (dnew(i,jstr-1)+dnew(i-1,jstr-1))* &
1826 & ubar(i,jstr-1,knew)
1827 END DO
1828 DO i=istrr,iendr
1829 dv_flux(i,jstrv-1)=0.5_r8*om_v(i,jstrv-1)* &
1830 & (dnew(i,jstrv-1)+dnew(i,jstrv-2))* &
1831 & vbar(i,jstrv-1,knew)
1832 END DO
1833 END IF
1834 END IF
1835 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1836 IF (domain(ng)%Northern_Edge(tile)) THEN
1837 DO i=istru,iend
1838 du_flux(i,jend+1)=0.5_r8*on_u(i,jend+1)* &
1839 & (dnew(i,jend+1)+dnew(i-1,jend+1))* &
1840 & ubar(i,jend+1,knew)
1841 END DO
1842 DO i=istrr,iendr
1843 dv_flux(i,jend+1)=0.5_r8*om_v(i,jend+1)* &
1844 & (dnew(i,jend+1)+dnew(i,jend))* &
1845 & vbar(i,jend+1,knew)
1846 END DO
1847 END IF
1848 END IF
1849#endif
1850!
1851! Apply momentum transport point sources (like river runoff), if any.
1852!
1853! Dsrc(is) = 0, flow across grid cell u-face (positive or negative)
1854! Dsrc(is) = 1, flow across grid cell v-face (positive or negative)
1855!
1856 IF (luvsrc(ng)) THEN
1857 DO is=1,nsrc(ng)
1858 i=sources(ng)%Isrc(is)
1859 j=sources(ng)%Jsrc(is)
1860 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1861 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1862 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
1863 cff=1.0_r8/(on_u(i,j)* &
1864 & 0.5_r8*(dnew(i-1,j)+dnew(i,j)))
1865 ubar(i,j,knew)=sources(ng)%Qbar(is)*cff
1866#ifdef SOLVE3D
1867 du_avg1(i,j)=sources(ng)%Qbar(is)
1868#endif
1869#if defined NESTING && !defined SOLVE3D
1870 du_flux(i,j)=sources(ng)%Qbar(is)
1871#endif
1872 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
1873 cff=1.0_r8/(om_v(i,j)* &
1874 & 0.5_r8*(dnew(i,j-1)+dnew(i,j)))
1875 vbar(i,j,knew)=sources(ng)%Qbar(is)*cff
1876#ifdef SOLVE3D
1877 dv_avg1(i,j)=sources(ng)%Qbar(is)
1878#endif
1879#if defined NESTING && !defined SOLVE3D
1880 dv_flux(i,j)=sources(ng)%Qbar(is)
1881#endif
1882 END IF
1883 END IF
1884 END DO
1885 END IF
1886
1887#ifdef SOLVE3D
1888!
1889!-----------------------------------------------------------------------
1890! Finalize computation of barotropic mode averages.
1891!-----------------------------------------------------------------------
1892!
1893! This procedure starts with filling in boundary rows of total depths
1894! at the new time step, which is needed to be done only during the
1895! last barotropic time step, Normally, the computation of averages
1896! occurs at the beginning of the next predictor step because "DUon"
1897! and "DVom" are being computed anyway. Strictly speaking, the filling
1898! the boundaries are necessary only in the case of open boundaries,
1899! otherwise, the associated fluxes are all zeros.
1900!
1901 IF ((iif(ng).eq.nfast(ng)).and.(knew.lt.3)) THEN
1902 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1903 IF (domain(ng)%Western_Edge(tile)) THEN
1904 DO j=jstr-1,jendr
1905 dnew(istr-1,j)=h(istr-1,j)+zeta_new(istr-1,j)
1906 END DO
1907 END IF
1908 END IF
1909 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1910 IF (domain(ng)%Eastern_Edge(tile)) THEN
1911 DO j=jstr-1,jendr
1912 dnew(iend+1,j)=h(iend+1,j)+zeta_new(iend+1,j)
1913 END DO
1914 END IF
1915 END IF
1916 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1917 IF (domain(ng)%Southern_Edge(tile)) THEN
1918 DO i=istr-1,iendr
1919 dnew(i,jstr-1)=h(i,jstr-1)+zeta_new(i,jstr-1)
1920 END DO
1921 END IF
1922 END IF
1923 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1924 IF (domain(ng)%Northern_Edge(tile)) THEN
1925 DO i=istr-1,iendr
1926 dnew(i,jend+1)=h(i,jend+1)+zeta_new(i,jend+1)
1927 END DO
1928 END IF
1929 END IF
1930!
1931! At the end of the last 2D time step replace the new free-surface
1932! zeta(:,:,knew) with it fast time-averaged value, Zt_avg1. Recall
1933! this is state variable is the one that communicates with the 3D
1934! kernel. Then, compute time-dependent depths.
1935!
1936 cff=weight(1,iif(ng),ng)
1937 cff1=0.5*cff
1938 DO j=jstrr,jendr
1939 DO i=istrr,iendr
1940 zt_avg1(i,j)=zt_avg1(i,j)+ &
1941 & cff*zeta(i,j,knew)
1942 IF (i.ge.istr) THEN
1943 du_avg1(i,j)=du_avg1(i,j)+ &
1944 & cff1*on_u(i,j)* &
1945 & (dnew(i,j)+dnew(i-1,j))*ubar(i,j,knew)
1946 END IF
1947 IF (j.ge.jstr) THEN
1948 dv_avg1(i,j)=dv_avg1(i,j)+ &
1949 & cff1*om_v(i,j)* &
1950 & (dnew(i,j)+dnew(i,j-1))*vbar(i,j,knew)
1951 END IF
1952 zeta(i,j,knew)=zt_avg1(i,j)
1953 END DO
1954 END DO
1955 CALL set_depth (ng, tile, inlm)
1956
1957# ifdef NESTING
1958!
1959! After all fast time steps are completed, apply boundary conditions
1960! to time averaged fields.
1961!
1962! In nesting applications with refinement grids, we need to exchange
1963! the DU_avg2 and DV_avg2 fluxes boundary information for the case
1964! that a contact point is at a tile partition. Notice that in such
1965! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
1966!
1967 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1968 CALL exchange_r2d_tile (ng, tile, &
1969 & lbi, ubi, lbj, ubj, &
1970 & zt_avg1)
1971 CALL exchange_u2d_tile (ng, tile, &
1972 & lbi, ubi, lbj, ubj, &
1973 & du_avg1)
1974 CALL exchange_v2d_tile (ng, tile, &
1975 & lbi, ubi, lbj, ubj, &
1976 & dv_avg1)
1977 CALL exchange_u2d_tile (ng, tile, &
1978 & lbi, ubi, lbj, ubj, &
1979 & du_avg2)
1980 CALL exchange_v2d_tile (ng, tile, &
1981 & lbi, ubi, lbj, ubj, &
1982 & dv_avg2)
1983 END IF
1984
1985# ifdef DISTRIBUTE
1986 CALL mp_exchange2d (ng, tile, inlm, 3, &
1987 & lbi, ubi, lbj, ubj, &
1988 & nghostpoints, &
1989 & ewperiodic(ng), nsperiodic(ng), &
1990 & zt_avg1, du_avg1, dv_avg1)
1991 CALL mp_exchange2d (ng, tile, inlm, 2, &
1992 & lbi, ubi, lbj, ubj, &
1993 & nghostpoints, &
1994 & ewperiodic(ng), nsperiodic(ng), &
1995 & du_avg2, dv_avg2)
1996# endif
1997# endif
1998 END IF
1999#endif
2000#if defined NESTING && !defined SOLVE3D
2001!
2002! In nesting applications with refinement grids, we need to exchange
2003! the DU_flux and DV_flux fluxes boundary information for the case
2004! that a contact point is at a tile partition. Notice that in such
2005! cases, we need i+1 and j+1 values for spatial/temporal interpolation.
2006!
2007 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
2008 CALL exchange_u2d_tile (ng, tile, &
2009 & lbi, ubi, lbj, ubj, &
2010 & du_flux)
2011 CALL exchange_v2d_tile (ng, tile, &
2012 & lbi, ubi, lbj, ubj, &
2013 & dv_flux)
2014 END IF
2015
2016# ifdef DISTRIBUTE
2017!
2018 CALL mp_exchange2d (ng, tile, inlm, 2, &
2019 & lbi, ubi, lbj, ubj, &
2020 & nghostpoints, &
2021 & ewperiodic(ng), nsperiodic(ng), &
2022 & du_flux, dv_flux)
2023# endif
2024#endif
2025!
2026! Deallocate local new free-surface.
2027!
2028 deallocate ( zeta_new )
2029
2030#ifdef WET_DRY
2031!
2032!-----------------------------------------------------------------------
2033! Compute new wet/dry masks.
2034!-----------------------------------------------------------------------
2035!
2036 CALL wetdry_tile (ng, tile, &
2037 & lbi, ubi, lbj, ubj, &
2038 & imins, imaxs, jmins, jmaxs, &
2039# ifdef MASKING
2040 & pmask, rmask, umask, vmask, &
2041# endif
2042 & h, zeta(:,:,knew), &
2043# ifdef SOLVE3D
2044 & du_avg1, dv_avg1, &
2045 & rmask_wet_avg, &
2046# endif
2047 & pmask_wet, pmask_full, &
2048 & rmask_wet, rmask_full, &
2049 & umask_wet, umask_full, &
2050 & vmask_wet, vmask_full)
2051#endif
2052!
2053!-----------------------------------------------------------------------
2054! Exchange boundary information.
2055!-----------------------------------------------------------------------
2056!
2057 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
2058 CALL exchange_r2d_tile (ng, tile, &
2059 & lbi, ubi, lbj, ubj, &
2060 & zeta(:,:,knew))
2061 CALL exchange_u2d_tile (ng, tile, &
2062 & lbi, ubi, lbj, ubj, &
2063 & ubar(:,:,knew))
2064 CALL exchange_v2d_tile (ng, tile, &
2065 & lbi, ubi, lbj, ubj, &
2066 & vbar(:,:,knew))
2067 END IF
2068
2069#ifdef DISTRIBUTE
2070!
2071 CALL mp_exchange2d (ng, tile, inlm, 3, &
2072 & lbi, ubi, lbj, ubj, &
2073 & nghostpoints, &
2074 & ewperiodic(ng), nsperiodic(ng), &
2075 & zeta(:,:,knew), &
2076 & ubar(:,:,knew), &
2077 & vbar(:,:,knew))
2078#endif
2079!
2080 RETURN
2081 END SUBROUTINE step2d_tile
2082!
2083 END MODULE step2d_mod
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_r2d_tile
subroutine exchange_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:253
exchange_2d_mod::exchange_u2d_tile
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:440
exchange_2d_mod::exchange_v2d_tile
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:627
mod_coupling
Definition mod_coupling.F:2
mod_coupling::coupling
type(t_coupling), dimension(:), allocatable coupling
Definition mod_coupling.F:102
mod_diags
Definition mod_diags.F:2
mod_diags::diags
type(t_diags), dimension(:), allocatable diags
Definition mod_diags.F:100
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ncparam
Definition mod_ncparam.F:2
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_parallel
Definition mod_parallel.F:2
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
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::m2fcor
integer m2fcor
Definition mod_scalars.F:150
mod_scalars::dcrit
real(r8), dimension(:), allocatable dcrit
Definition mod_scalars.F:782
mod_scalars::m2pgrd
integer m2pgrd
Definition mod_scalars.F:167
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::m2xadv
integer m2xadv
Definition mod_scalars.F:152
mod_scalars::m2yadv
integer m2yadv
Definition mod_scalars.F:153
mod_scalars::predictor_2d_step
logical, dimension(:), allocatable predictor_2d_step
Definition mod_scalars.F:245
mod_scalars::weight
real(dp), dimension(:,:,:), allocatable weight
Definition mod_scalars.F:809
mod_scalars::volcons
logical, dimension(:,:), allocatable volcons
Definition mod_scalars.F:515
mod_scalars::nfast
integer, dimension(:), allocatable nfast
Definition mod_scalars.F:259
mod_scalars::lwsrc
logical, dimension(:), allocatable lwsrc
Definition mod_scalars.F:543
mod_scalars::m2hvis
integer m2hvis
Definition mod_scalars.F:168
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::m2yvis
integer m2yvis
Definition mod_scalars.F:170
mod_scalars::dtfast
real(dp), dimension(:), allocatable dtfast
Definition mod_scalars.F:270
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::m2hadv
integer m2hadv
Definition mod_scalars.F:151
mod_scalars::rho0
real(dp) rho0
Definition mod_scalars.F:787
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::m2xvis
integer m2xvis
Definition mod_scalars.F:169
mod_scalars::m2bstr
integer m2bstr
Definition mod_scalars.F:172
mod_scalars::iif
integer, dimension(:), allocatable iif
Definition mod_scalars.F:250
mod_sedbed
Definition sedbed_mod.h:1
mod_sedbed::sedbed
type(t_sedbed), dimension(:), allocatable sedbed
Definition sedbed_mod.h:157
mod_sources
Definition mod_sources.F:2
mod_sources::sources
type(t_sources), dimension(:), allocatable sources
Definition mod_sources.F:90
mod_sources::nsrc
integer, dimension(:), allocatable nsrc
Definition mod_sources.F:97
mod_stepping
Definition mod_stepping.F:2
mod_stepping::kstp
integer, dimension(:), allocatable kstp
Definition mod_stepping.F:66
mod_stepping::knew
integer, dimension(:), allocatable knew
Definition mod_stepping.F:64
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
obc_volcons_mod
Definition obc_volcons.F:3
obc_volcons_mod::set_duv_bc_tile
subroutine, public set_duv_bc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, om_v, on_u, ubar, vbar, drhs, duon, dvom)
Definition obc_volcons.F:245
obc_volcons_mod::obc_flux_tile
subroutine, public obc_flux_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, kinp, umask, vmask, h, om_v, on_u, ubar, vbar, zeta)
Definition obc_volcons.F:69
set_depth_mod
Definition set_depth.F:2
set_depth_mod::set_depth
subroutine, public set_depth(ng, tile, model)
Definition set_depth.F:34
step2d_mod
Definition step2d_FB.h:3
step2d_mod::step2d_tile
subroutine 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, 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, bed_thick, eq_tide, sustr, svstr,
Definition step2d_FB.h:219
step2d_mod::step2d
subroutine, public step2d(ng, tile)
Definition step2d_FB.h:75
u2dbc_mod
Definition u2dbc_im.F:2
u2dbc_mod::u2dbc_tile
subroutine, public u2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta)
Definition u2dbc_im.F:56
v2dbc_mod
Definition v2dbc_im.F:3
v2dbc_mod::v2dbc_tile
subroutine, public v2dbc_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, krhs, kstp, kout, ubar, vbar, zeta)
Definition v2dbc_im.F:57
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
zetabc_mod
Definition zetabc.F:2
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