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