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