ROMS
Loading...
Searching...
No Matches
ad_step3d_uv.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE ad_step3d_uv_mod
3
4#if defined ADJOINT && defined SOLVE3D
5!
6!git $Id$
7!================================================== Hernan G. Arango ===
8! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
9! Licensed under a MIT/X style license !
10! See License_ROMS.md !
11!=======================================================================
12! !
13! This subroutine time-steps the adjoint horizontal momentum !
14! equations. The vertical viscosity terms are time-stepped using !
15! an implicit algorithm. !
16! !
17! BASIC STATE variables needed: u, v, ru, rv, Akv, !
18! Hz, Huon, Hvom, z_r, z_w, !
19! DU_avg1, DU_avg2, DV_avg1, DV_avg2, !
20! Qsrc !
21! !
22!=======================================================================
23!
24 USE mod_param
25 USE mod_coupling
26# ifdef DIAGNOSTICS_UV
27!! USE mod_diags
28# endif
29 USE mod_forces
30 USE mod_grid
31 USE mod_mixing
32 USE mod_ocean
33 USE mod_scalars
34 USE mod_sources
35!
36 USE ad_exchange_2d_mod
37 USE ad_exchange_3d_mod
38 USE exchange_3d_mod
39# ifdef DISTRIBUTE
40 USE mp_exchange_mod, ONLY : ad_mp_exchange2d, ad_mp_exchange3d
41# endif
42 USE ad_u3dbc_mod, ONLY : ad_u3dbc_tile
43 USE ad_v3dbc_mod, ONLY : ad_v3dbc_tile
44# ifdef UV_DESTAGGERED
45 USE uv_var_change_mod, ONLY : ad_uv_c2a_grid
46# endif
47!
48 implicit none
49!
50 PRIVATE
51 PUBLIC :: ad_step3d_uv
52!
53 CONTAINS
54!
55!***********************************************************************
56 SUBROUTINE ad_step3d_uv (ng, tile)
57!***********************************************************************
58!
59 USE mod_stepping
60!
61! Imported variable declarations.
62!
63 integer, intent(in) :: ng, tile
64!
65! Local variable declarations.
66!
67 character (len=*), parameter :: myfile = &
68 & __FILE__
69!
70# include "tile.h"
71!
72# ifdef PROFILE
73 CALL wclock_on (ng, iadm, 34, __line__, myfile)
74# endif
75 CALL ad_step3d_uv_tile (ng, tile, &
76 & lbi, ubi, lbj, ubj, &
77 & imins, imaxs, jmins, jmaxs, &
78 & knew(ng), nrhs(ng), nstp(ng), nnew(ng), &
79# ifdef MASKING
80 & grid(ng) % umask, &
81 & grid(ng) % vmask, &
82# endif
83# ifdef WET_DRY_NOT_YET
84 & grid(ng) % umask_wet, &
85 & grid(ng) % vmask_wet, &
86# endif
87 & grid(ng) % om_v, &
88 & grid(ng) % on_u, &
89 & grid(ng) % pm, &
90 & grid(ng) % pn, &
91 & grid(ng) % Hz, &
92 & grid(ng) % ad_Hz, &
93 & grid(ng) % z_r, &
94 & grid(ng) % ad_z_r, &
95 & grid(ng) % z_w, &
96 & grid(ng) % ad_z_w, &
97 & mixing(ng) % Akv, &
98 & mixing(ng) % ad_Akv, &
99 & coupling(ng) % DU_avg1, &
100 & coupling(ng) % ad_DU_avg1, &
101 & coupling(ng) % DV_avg1, &
102 & coupling(ng) % ad_DV_avg1, &
103 & coupling(ng) % DU_avg2, &
104 & coupling(ng) % ad_DU_avg2, &
105 & coupling(ng) % DV_avg2, &
106 & coupling(ng) % ad_DV_avg2, &
107 & ocean(ng) % ad_ru, &
108 & ocean(ng) % ad_rv, &
109# ifdef DIAGNOSTICS_UV
110!! & DIAGS(ng) % DiaU2wrk, &
111!! & DIAGS(ng) % DiaV2wrk, &
112!! & DIAGS(ng) % DiaU2int, &
113!! & DIAGS(ng) % DiaV2int, &
114!! & DIAGS(ng) % DiaU3wrk, &
115!! & DIAGS(ng) % DiaV3wrk, &
116!! & DIAGS(ng) % DiaRU, &
117!! & DIAGS(ng) % DiaRV, &
118# endif
119 & ocean(ng) % u, &
120 & ocean(ng) % ad_u, &
121 & ocean(ng) % v, &
122 & ocean(ng) % ad_v, &
123# ifdef AD_OUTPUT_STATE
124 & ocean(ng) % ad_u_sol, &
125 & ocean(ng) % ad_v_sol, &
126# endif
127 & ocean(ng) % ad_ubar, &
128 & ocean(ng) % ad_vbar, &
129 & ocean(ng) % ad_ubar_sol, &
130 & ocean(ng) % ad_vbar_sol, &
131# ifdef NEARSHORE_MELLOR
132 & ocean(ng) % ubar_stokes, &
133 & ocean(ng) % tl_ubar_stokes, &
134 & ocean(ng) % vbar_stokes, &
135 & ocean(ng) % tl_vbar_stokes, &
136 & ocean(ng) % u_stokes, &
137 & ocean(ng) % tl_u_stokes, &
138 & ocean(ng) % v_stokes, &
139 & ocean(ng) % tl_v_stokes, &
140# endif
141 & grid(ng) % Huon, &
142 & grid(ng) % ad_Huon, &
143 & grid(ng) % Hvom, &
144 & grid(ng) % ad_Hvom)
145# ifdef PROFILE
146 CALL wclock_off (ng, iadm, 34, __line__, myfile)
147# endif
148!
149 RETURN
150 END SUBROUTINE ad_step3d_uv
151!
152!***********************************************************************
153 SUBROUTINE ad_step3d_uv_tile (ng, tile, &
154 & LBi, UBi, LBj, UBj, &
155 & IminS, ImaxS, JminS, JmaxS, &
156 & knew, nrhs, nstp, nnew, &
157# ifdef MASKING
158 & umask, vmask, &
159# endif
160# ifdef WET_DRY_NOT_YET
161 & umask_wet, vmask_wet, &
162# endif
163 & om_v, on_u, pm, pn, &
164 & Hz, ad_Hz, &
165 & z_r, ad_z_r, &
166 & z_w, ad_z_w, &
167 & Akv, ad_Akv, &
168 & DU_avg1, ad_DU_avg1, &
169 & DV_avg1, ad_DV_avg1, &
170 & DU_avg2, ad_DU_avg2, &
171 & DV_avg2, ad_DV_avg2, &
172 & ad_ru, ad_rv, &
173# ifdef DIAGNOSTICS_UV
174!! & DiaU2wrk, DiaV2wrk, &
175!! & DiaU2int, DiaV2int, &
176!! & DiaU3wrk, DiaV3wrk, &
177!! & DiaRU, DiaRV, &
178# endif
179 & u, ad_u, &
180 & v, ad_v, &
181# ifdef AD_OUTPUT_STATE
182 & ad_u_sol, ad_v_sol, &
183# endif
184 & ad_ubar, ad_vbar, &
185 & ad_ubar_sol, ad_vbar_sol, &
186# ifdef NEARSHORE_MELLOR
187 & ubar_stokes, ad_ubar_stokes, &
188 & vbar_stokes, ad_vbar_stokes, &
189 & u_stokes, ad_u_stokes, &
190 & v_stokes, ad_v_stokes, &
191# endif
192 & Huon, ad_Huon, &
193 & Hvom, ad_Hvom)
194!***********************************************************************
195!
196! Imported variable declarations.
197!
198 integer, intent(in) :: ng, tile
199 integer, intent(in) :: LBi, UBi, LBj, UBj
200 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
201 integer, intent(in) :: knew, nrhs, nstp, nnew
202!
203# ifdef ASSUMED_SHAPE
204# ifdef MASKING
205 real(r8), intent(in) :: umask(LBi:,LBj:)
206 real(r8), intent(in) :: vmask(LBi:,LBj:)
207# endif
208# ifdef WET_DRY_NOT_YET
209 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
210 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
211# endif
212 real(r8), intent(in) :: om_v(LBi:,LBj:)
213 real(r8), intent(in) :: on_u(LBi:,LBj:)
214 real(r8), intent(in) :: pm(LBi:,LBj:)
215 real(r8), intent(in) :: pn(LBi:,LBj:)
216 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
217 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
218 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
219 real(r8), intent(in) :: Akv(LBi:,LBj:,0:)
220 real(r8), intent(in) :: DU_avg1(LBi:,LBj:)
221 real(r8), intent(in) :: DV_avg1(LBi:,LBj:)
222 real(r8), intent(in) :: DU_avg2(LBi:,LBj:)
223 real(r8), intent(in) :: DV_avg2(LBi:,LBj:)
224 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
225 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
226# ifdef NEARSHORE_MELLOR
227 real(r8), intent(in) :: ubar_stokes(LBi:,LBj:)
228 real(r8), intent(in) :: vbar_stokes(LBi:,LBj:)
229 real(r8), intent(in) :: u_stokes(LBi:,LBj:,:)
230 real(r8), intent(in) :: v_stokes(LBi:,LBj:,:)
231# endif
232# ifdef DIAGNOSTICS_UV
233!! real(r8), intent(inout) :: DiaU2wrk(LBi:,LBj:,:)
234!! real(r8), intent(inout) :: DiaV2wrk(LBi:,LBj:,:)
235!! real(r8), intent(inout) :: DiaU2int(LBi:,LBj:,:)
236!! real(r8), intent(inout) :: DiaV2int(LBi:,LBj:,:)
237!! real(r8), intent(inout) :: DiaU3wrk(LBi:,LBj:,:,:)
238!! real(r8), intent(inout) :: DiaV3wrk(LBi:,LBj:,:,:)
239!! real(r8), intent(inout) :: DiaRU(LBi:,LBj:,:,:,:)
240!! real(r8), intent(inout) :: DiaRV(LBi:,LBj:,:,:,:)
241# endif
242 real(r8), intent(inout) :: Huon(LBi:,LBj:,:)
243 real(r8), intent(inout) :: Hvom(LBi:,LBj:,:)
244 real(r8), intent(inout) :: ad_Hz(LBi:,LBj:,:)
245 real(r8), intent(inout) :: ad_z_r(LBi:,LBj:,:)
246 real(r8), intent(inout) :: ad_z_w(LBi:,LBj:,0:)
247 real(r8), intent(inout) :: ad_Akv(LBi:,LBj:,0:)
248 real(r8), intent(inout) :: ad_DU_avg1(LBi:,LBj:)
249 real(r8), intent(inout) :: ad_DV_avg1(LBi:,LBj:)
250 real(r8), intent(inout) :: ad_DU_avg2(LBi:,LBj:)
251 real(r8), intent(inout) :: ad_DV_avg2(LBi:,LBj:)
252 real(r8), intent(inout) :: ad_ru(LBi:,LBj:,0:,:)
253 real(r8), intent(inout) :: ad_rv(LBi:,LBj:,0:,:)
254 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
255 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
256# ifdef AD_OUTPUT_STATE
257 real(r8), intent(inout) :: ad_u_sol(LBi:,LBj:,:)
258 real(r8), intent(inout) :: ad_v_sol(LBi:,LBj:,:)
259# endif
260# ifdef NEARSHORE_MELLOR
261 real(r8), intent(inout) :: ad_ubar_stokes(LBi:,LBj:)
262 real(r8), intent(inout) :: ad_vbar_stokes(LBi:,LBj:)
263 real(r8), intent(inout) :: ad_u_stokes(LBi:,LBj:,:)
264 real(r8), intent(inout) :: ad_v_stokes(LBi:,LBj:,:)
265# endif
266 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
267 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
268 real(r8), intent(inout) :: ad_Huon(LBi:,LBj:,:)
269 real(r8), intent(inout) :: ad_Hvom(LBi:,LBj:,:)
270
271 real(r8), intent(out) :: ad_ubar_sol(LBi:,LBj:)
272 real(r8), intent(out) :: ad_vbar_sol(LBi:,LBj:)
273
274# else
275
276# ifdef MASKING
277 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
278 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
279# endif
280# ifdef WET_DRY_NOT_YET
281 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
282 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
283# endif
284 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
285 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
286 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
287 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
288 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
289 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
290 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
291 real(r8), intent(in) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
292 real(r8), intent(in) :: DU_avg1(LBi:UBi,LBj:UBj)
293 real(r8), intent(in) :: DV_avg1(LBi:UBi,LBj:UBj)
294 real(r8), intent(in) :: DU_avg2(LBi:UBi,LBj:UBj)
295 real(r8), intent(in) :: DV_avg2(LBi:UBi,LBj:UBj)
296 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
297 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
298# ifdef NEARSHORE_MELLOR
299 real(r8), intent(in) :: ubar_stokes(LBi:UBi,LBj:UBj)
300 real(r8), intent(in) :: vbar_stokes(LBi:UBi,LBj:UBj)
301 real(r8), intent(in) :: u_stokes(LBi:UBi,LBj:UBj,N(ng))
302 real(r8), intent(in) :: v_stokes(LBi:UBi,LBj:UBj,N(ng))
303# endif
304# ifdef DIAGNOSTICS_UV
305!! real(r8), intent(inout) :: DiaU2wrk(LBi:UBi,LBj:UBj,NDM2d)
306!! real(r8), intent(inout) :: DiaV2wrk(LBi:UBi,LBj:UBj,NDM2d)
307!! real(r8), intent(inout) :: DiaU2int(LBi:UBi,LBj:UBj,NDM2d)
308!! real(r8), intent(inout) :: DiaV2int(LBi:UBi,LBj:UBj,NDM2d)
309!! real(r8), intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
310!! real(r8), intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
311!! real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
312!! real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
313# endif
314 real(r8), intent(inout) :: Huon(LBi:UBi,LBj:UBj,N(ng))
315 real(r8), intent(inout) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
316 real(r8), intent(inout) :: ad_Hz(LBi:UBi,LBj:UBj,N(ng))
317 real(r8), intent(inout) :: ad_z_r(LBi:UBi,LBj:UBj,N(ng))
318 real(r8), intent(inout) :: ad_z_w(LBi:UBi,LBj:UBj,0:N(ng))
319 real(r8), intent(inout) :: ad_Akv(LBi:UBi,LBj:UBj,0:N(ng))
320 real(r8), intent(inout) :: ad_DU_avg1(LBi:UBi,LBj:UBj)
321 real(r8), intent(inout) :: ad_DV_avg1(LBi:UBi,LBj:UBj)
322 real(r8), intent(inout) :: ad_DU_avg2(LBi:UBi,LBj:UBj)
323 real(r8), intent(inout) :: ad_DV_avg2(LBi:UBi,LBj:UBj)
324 real(r8), intent(inout) :: ad_ru(LBi:UBi,LBj:UBj,0:N(ng),2)
325 real(r8), intent(inout) :: ad_rv(LBi:UBi,LBj:UBj,0:N(ng),2)
326 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
327 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
328# ifdef AD_OUTPUT_STATE
329 real(r8), intent(inout) :: ad_u_sol(LBi:UBi,LBj:UBj,N(ng))
330 real(r8), intent(inout) :: ad_v_sol(LBi:UBi,LBj:UBj,N(ng))
331# endif
332# ifdef NEARSHORE_MELLOR
333 real(r8), intent(in) :: ad_ubar_stokes(LBi:UBi,LBj:UBj)
334 real(r8), intent(in) :: ad_vbar_stokes(LBi:UBi,LBj:UBj)
335 real(r8), intent(inout) :: ad_u_stokes(LBi:UBi,LBj:UBj,N(ng))
336 real(r8), intent(inout) :: ad_v_stokes(LBi:UBi,LBj:UBj,N(ng))
337# endif
338 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,3)
339 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,3)
340 real(r8), intent(inout) :: ad_Huon(LBi:UBi,LBj:UBj,N(ng))
341 real(r8), intent(inout) :: ad_Hvom(LBi:UBi,LBj:UBj,N(ng))
342
343 real(r8), intent(out) :: ad_ubar_sol(LBi:UBi,LBj:UBj)
344 real(r8), intent(out) :: ad_vbar_sol(LBi:UBi,LBj:UBj)
345# endif
346!
347! Local variable declarations.
348!
349 integer :: i, idiag, is, j, k
350!
351 real(r8) :: cff, cff1, cff2
352 real(r8) :: ad_cff, ad_cff1, ad_cff2
353 real(r8) :: adfac, adfac1, adfac2
354!
355 real(r8), dimension(IminS:ImaxS) :: CF1
356 real(r8), dimension(IminS:ImaxS) :: FC1
357# ifdef NEARSHORE_MELLOR
358 real(r8), dimension(IminS:ImaxS) :: CFs1
359 real(r8), dimension(IminS:ImaxS) :: DCs1
360# endif
361 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: AK
362 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
363 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
364 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
365 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC1
366 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
367# ifdef NEARSHORE_MELLOR
368 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CFs
369 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DCs
370# endif
371 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hzk
372 real(r8), dimension(IminS:ImaxS,N(ng)) :: oHz
373# ifdef DIAGNOSTICS_UV
374!! real(r8), dimension(IminS:ImaxS,0:N(ng)) :: wrk
375!! real(r8), dimension(IminS:ImaxS,1:NDM2d-1) :: Dwrk
376# endif
377 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_AK
378 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_BC
379 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_CF
380 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_DC
381 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_FC
382# ifdef NEARSHORE_MELLOR
383 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_CFs
384 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: ad_DCs
385# endif
386 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_Hzk
387 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_oHz
388!
389# include "set_bounds.h"
390!
391!-----------------------------------------------------------------------
392! Initialize adjoint private variables.
393!-----------------------------------------------------------------------
394!
395 ad_cff=0.0_r8
396 ad_cff1=0.0_r8
397 DO k=0,n(ng)
398 DO i=imins,imaxs
399 ad_ak(i,k)=0.0_r8
400 ad_bc(i,k)=0.0_r8
401 ad_cf(i,k)=0.0_r8
402 ad_dc(i,k)=0.0_r8
403# ifdef NEARSHORE_MELLOR
404 ad_cfs(i,k)=0.0_r8
405 ad_dcs(i,k)=0.0_r8
406# endif
407 ad_fc(i,k)=0.0_r8
408 END DO
409 END DO
410 DO k=1,n(ng)
411 DO i=imins,imaxs
412 ad_hzk(i,k)=0.0_r8
413 ad_ohz(i,k)=0.0_r8
414 END DO
415 END DO
416
417# ifdef UV_DESTAGGERED
418!
419!-----------------------------------------------------------------------
420! Compute adjoint 3D momentum (ad_ua, ad_va) at RHO-points (A-Grid)
421! for output purposes and data assimilation where the observations
422! and state vector is located at the cell-center.
423!-----------------------------------------------------------------------
424!
425!> CALL tl_uv_C2A_grid (ng, tile, iTLM, nnew)
426!>
427 CALL ad_uv_c2a_grid (ng, tile, itlm, nnew)
428# endif
429!
430!-----------------------------------------------------------------------
431! Exchange boundary data.
432!-----------------------------------------------------------------------
433!
434# ifdef DISTRIBUTE
435!^ CALL mp_exchange2d (ng, tile, iTLM, 4, &
436!^ & LBi, UBi, LBj, UBj, &
437!^ & NghostPoints, &
438!^ & EWperiodic(ng), NSperiodic(ng), &
439!^ & tl_ubar(:,:,1), tl_vbar(:,:,1), &
440!^ & tl_ubar(:,:,2), tl_vbar(:,:,2))
441!^
442 CALL ad_mp_exchange2d (ng, tile, iadm, 4, &
443 & lbi, ubi, lbj, ubj, &
444 & nghostpoints, &
445 & ewperiodic(ng), nsperiodic(ng), &
446 & ad_ubar(:,:,1), ad_vbar(:,:,1), &
447 & ad_ubar(:,:,2), ad_vbar(:,:,2))
448
449!^ CALL mp_exchange3d (ng, tile, iTLM, 4, &
450!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
451!^ & NghostPoints, &
452!^ & EWperiodic(ng), NSperiodic(ng), &
453!^ & tl_u(:,:,:,nnew), tl_v(:,:,:,nnew), &
454!^ & tl_Huon, tl_Hvom)
455!^
456 CALL ad_mp_exchange3d (ng, tile, iadm, 4, &
457 & lbi, ubi, lbj, ubj, 1, n(ng), &
458 & nghostpoints, &
459 & ewperiodic(ng), nsperiodic(ng), &
460 & ad_u(:,:,:,nnew), ad_v(:,:,:,nnew), &
461 & ad_huon, ad_hvom)
462!
463# endif
464
465 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
466 DO k=1,2
467!^ CALL exchange_v2d_tile (ng, tile, &
468!^ & LBi, UBi, LBj, UBj, &
469!^ & tl_vbar(:,:,k))
470!^
471 CALL ad_exchange_v2d_tile (ng, tile, &
472 & lbi, ubi, lbj, ubj, &
473 & ad_vbar(:,:,k))
474!^ CALL exchange_u2d_tile (ng, tile, &
475!^ & LBi, UBi, LBj, UBj, &
476!^ & tl_ubar(:,:,k))
477!^
478 CALL ad_exchange_u2d_tile (ng, tile, &
479 & lbi, ubi, lbj, ubj, &
480 & ad_ubar(:,:,k))
481 END DO
482!
483 CALL exchange_v3d_tile (ng, tile, &
484 & lbi, ubi, lbj, ubj, 1, n(ng), &
485 & hvom)
486!^ CALL exchange_v3d_tile (ng, tile, &
487!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
488!^ & tl_Hvom)
489!^
490 CALL ad_exchange_v3d_tile (ng, tile, &
491 & lbi, ubi, lbj, ubj, 1, n(ng), &
492 & ad_hvom)
493
494 CALL exchange_u3d_tile (ng, tile, &
495 & lbi, ubi, lbj, ubj, 1, n(ng), &
496 & huon)
497!^ CALL exchange_u3d_tile (ng, tile, &
498!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
499!^ & tl_Huon)
500!^
501 CALL ad_exchange_u3d_tile (ng, tile, &
502 & lbi, ubi, lbj, ubj, 1, n(ng), &
503 & ad_huon)
504
505!^ CALL exchange_v3d_tile (ng, tile, &
506!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
507!^ & tl_v(:,:,:,nnew))
508!^
509 CALL ad_exchange_v3d_tile (ng, tile, &
510 & lbi, ubi, lbj, ubj, 1, n(ng), &
511 & ad_v(:,:,:,nnew))
512!^ CALL exchange_u3d_tile (ng, tile, &
513!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
514!^ & tl_u(:,:,:,nnew))
515!^
516 CALL ad_exchange_u3d_tile (ng, tile, &
517 & lbi, ubi, lbj, ubj, 1, n(ng), &
518 & ad_u(:,:,:,nnew))
519 END IF
520!
521!-----------------------------------------------------------------------
522! Couple 2D and 3D momentum equations.
523!-----------------------------------------------------------------------
524!
525! Couple adjoint velocity component in the ETA-direction.
526!
527 j_loop1 : DO j=jstrt,jendt
528 IF (j.ge.jstr) THEN
529!
530! First, compute BASIC STATE quantities. This section was computed
531! three times in the original code. This can avoided by saving the
532! intermediate values in scratch arrays DC1, CF1, and FC1 (HGA).
533!
534 DO i=istrt,iendt
535 dc(i,0)=0.0_r8
536 cf(i,0)=0.0_r8
537# ifdef NEARSHORE_MELLOR
538 cfs(i,0)=0.0_r8
539# endif
540 fc(i,0)=0.0_r8
541 END DO
542 DO k=1,n(ng)
543 DO i=istrt,iendt
544 cff=0.5_r8*om_v(i,j)
545 dc(i,k)=cff*(hz(i,j,k)+hz(i,j-1,k))
546 dc(i,0)=dc(i,0)+dc(i,k)
547 cf(i,0)=cf(i,0)+ &
548 & dc(i,k)*v(i,j,k,nnew)
549# ifdef NEARSHORE_MELLOR
550 cfs(i,0)=cfs(i,0)+ &
551 & dc(i,k)*v_stokes(i,j,k)
552# endif
553 END DO
554 END DO
555 DO i=istrt,iendt
556 dc1(i,0)=dc(i,0) ! intermediate
557# ifdef NEARSHORE_MELLOR
558 cff2=dc(i,0)*vbar_stokes(i,j)
559# endif
560 dc(i,0)=1.0_r8/dc(i,0) ! recursive
561 cf1(i)=cf(i,0) ! intermediate
562 cf(i,0)=dc(i,0)*(cf(i,0)-dv_avg1(i,j)) ! recursive
563# ifdef NEARSHORE_MELLOR
564 cfs1(i)=cfs(i,0)
565 cfs(i,0)=dc(i,0)*(cfs(i,0)-cff2) ! recursive
566# endif
567 END DO
568 DO k=n(ng),1,-1
569 DO i=istrt,iendt
570!^ Hvom(i,j,k)=0.5_r8*(Hvom(i,j,k)+v(i,j,k,nnew)*DC(i,k))
571# ifdef NEARSHORE_MELLOR
572!^ Hvom(i,j,k)=Hvom(i,j,k)+0.5_r8*v_stokes(i,j,k)*DC(i,k)
573# endif
574!^ FC(i,0)=FC(i,0)+Hvom(i,j,k)
575!^
576 fc(i,0)=fc(i,0)+ &
577 & 0.5_r8*(hvom(i,j,k)+v(i,j,k,nnew)*dc(i,k))
578# ifdef NEARSHORE_MELLOR
579 fc(i,0)=fc(i,0)+ &
580 & 0.5_r8*v_stokes(i,j,k)*dc(i,k)
581# endif
582 END DO
583 END DO
584 DO i=istrt,iendt
585 fc1(i)=fc(i,0) ! intermediate
586 fc(i,0)=dc(i,0)*(fc(i,0)-dv_avg2(i,j)) ! recursive
587 END DO
588!
589! Compute correct mass flux, Hz*v/m.
590!
591 DO k=1,n(ng)
592 DO i=istrt,iendt
593!^ tl_Hvom(i,j,k)=tl_Hvom(i,j,k)- &
594!^ & tl_DC(i,k)*FC(i,0)- &
595!^ & DC(i,k)*tl_FC(i,0)
596!^
597 ad_dc(i,k)=ad_dc(i,k)-ad_hvom(i,j,k)*fc(i,0)
598 ad_fc(i,0)=ad_fc(i,0)-ad_hvom(i,j,k)*dc(i,k)
599 END DO
600 END DO
601 DO i=istrt,iendt
602!^ tl_FC(i,0)=tl_DC(i,0)*(FC1(i)-DV_avg2(i,j))+ &
603!^ & DC(i,0)*(tl_FC(i,0)-tl_DV_avg2(i,j))
604!^
605 adfac=dc(i,0)*ad_fc(i,0)
606 ad_dc(i,0)=ad_dc(i,0)+ad_fc(i,0)*(fc1(i)-dv_avg2(i,j))
607 ad_dv_avg2(i,j)=ad_dv_avg2(i,j)-adfac
608 ad_fc(i,0)=adfac
609 END DO
610!^ DO k=N(ng),1,-1
611!^
612 DO k=1,n(ng)
613 DO i=istrt,iendt
614# ifdef DIAGNOSTICS_UV
615!! DiaV3wrk(i,j,k,M3rate)=v(i,j,k,nnew)- &
616!! & DiaV3wrk(i,j,k,M3rate)
617# endif
618!^ tl_FC(i,0)=tl_FC(i,0)+tl_Hvom(i,j,k)
619!^
620 ad_hvom(i,j,k)=ad_hvom(i,j,k)+ad_fc(i,0)
621# ifdef NEARSHORE_MELLOR
622!^ tl_Hvom(i,j,k)=tl_Hvom(i,j,k)+ &
623!^ & 0.5_r8*(tl_v_stokes(i,j,k)*DC(i,k)+ &
624!^ & v_stokes(i,j,k)*tl_DC(i,k))
625!^
626 adfac=0.5_r8*ad_hvom(i,j,k)
627 ad_dc(i,k)=ad_dc(i,k)+v_stokes(i,j,k)*adfac
628 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)+dc(i,k)*adfac
629# endif
630!^ tl_Hvom(i,j,k)=0.5_r8*(tl_Hvom(i,j,k)+ &
631!^ & tl_v(i,j,k,nnew)*DC(i,k)+ &
632!^ & v(i,j,k,nnew)*tl_DC(i,k))
633!^
634 adfac=0.5_r8*ad_hvom(i,j,k)
635 ad_dc(i,k)=ad_dc(i,k)+v(i,j,k,nnew)*adfac
636 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)+dc(i,k)*adfac
637 ad_hvom(i,j,k)=adfac
638 END DO
639 END DO
640!
641! Replace only BOUNDARY POINTS incorrect vertical mean with more
642! accurate barotropic component, vbar=DV_avg1/(D*om_v). Recall that,
643! D=CF(:,0). The replacement is avoided if the boundary edge is
644! periodic. The J-loop is pipelined, so we need to do a special
645! test for the southern and northern domain edges.
646!
647 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
648 IF (j.eq.mm(ng)+1) THEN ! northern boundary
649 DO k=1,n(ng) ! J-loop pipelined
650 DO i=istr,iend
651# ifdef NEARSHORE_MELLOR
652# ifdef WET_DRY_NOT_YET
653!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)* &
654!^ & vmask_wet(i,j)
655!^
656 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)* &
657 & vmask_wet(i,j)
658# endif
659# ifdef MASKING
660!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)* &
661!^ & vmask(i,j)
662!^
663 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)* &
664 & vmask(i,j)
665# endif
666!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)- &
667!^ & tl_CFs(i,0)
668!^
669 ad_cfs(i,0)=ad_cfs(i,0)-ad_v_stokes(i,j,k)
670# endif
671# ifdef WET_DRY_NOT_YET
672!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)* &
673!^ & vmask_wet(i,j)
674!^
675 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)* &
676 & vmask_wet(i,j)
677# endif
678# ifdef MASKING
679!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)* &
680!^ & vmask(i,j)
681!^
682 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)* &
683 & vmask(i,j)
684# endif
685!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)- &
686!^ & tl_CF(i,0)
687!^
688 ad_cf(i,0)=ad_cf(i,0)-ad_v(i,j,k,nnew)
689 END DO
690 END DO
691 END IF
692 END IF
693!
694 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
695 IF (j.eq.1) THEN ! southern boundary
696 DO k=1,n(ng) ! J-loop pipelined
697 DO i=istr,iend
698# ifdef NEARSHORE_MELLOR
699# ifdef WET_DRY_NOT_YET
700!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)* &
701!^ & vmask_wet(i,j)
702!^
703 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)* &
704 & vmask_wet(i,j)
705# endif
706# ifdef MASKING
707!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)* &
708!^ & vmask(i,j)
709!^
710 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)* &
711 & vmask(i,j)
712# endif
713!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)- &
714!^ & tl_CFs(i,0)
715!^
716 ad_cfs(i,0)=ad_cfs(i,0)-ad_v_stokes(i,j,k)
717# endif
718# ifdef WET_DRY_NOT_YET
719!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)* &
720!^ & vmask_wet(i,j)
721!^
722 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)* &
723 & vmask_wet(i,j)
724# endif
725# ifdef MASKING
726!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)* &
727!^ & vmask(i,j)
728!^
729 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)* &
730 & vmask(i,j)
731# endif
732!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)- &
733!^ & tl_CF(i,0)
734!^
735 ad_cf(i,0)=ad_cf(i,0)-ad_v(i,j,k,nnew)
736 END DO
737 END DO
738 END IF
739 END IF
740!
741 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
742 IF (domain(ng)%Eastern_Edge(tile)) THEN
743 DO k=1,n(ng)
744# ifdef NEARSHORE_MELLOR
745# ifdef WET_DRY_NOT_YET
746!^ tl_v_stokes(Iend+1,j,k)=tl_v_stokes(Iend+1,j,k)* &
747!^ & vmask_wet(Iend+1,j)
748!^
749 ad_v_stokes(iend+1,j,k)=ad_v_stokes(iend+1,j,k)* &
750 & vmask_wet(iend+1,j)
751# endif
752# ifdef MASKING
753!^ tl_v_stokes(Iend+1,j,k)=tl_v_stokes(Iend+1,j,k)* &
754!^ & vmask(Iend+1,j)
755!^
756 ad_v_stokes(iend+1,j,k)=ad_v_stokes(iend+1,j,k)* &
757 & vmask(iend+1,j)
758# endif
759!^ tl_v_stokes(Iend+1,j,k)=tl_v_stokes(Iend+1,j,k)- &
760!^ & tl_CFs(Iend+1,0)
761!^
762 ad_cfs(iend+1,0)=ad_cfs(iend+1,0)- &
763 & ad_v_stokes(iend+1,j,k)
764# endif
765# ifdef WET_DRY_NOT_YET
766!^ tl_v(Iend+1,j,k,nnew)=tl_v(Iend+1,j,k,nnew)* &
767!^ & vmask_wet(Iend+1,j)
768!^
769 ad_v(iend+1,j,k,nnew)=ad_v(iend+1,j,k,nnew)* &
770 & vmask_wet(iend+1,j)
771# endif
772# ifdef MASKING
773!^ tl_v(Iend+1,j,k,nnew)=tl_v(Iend+1,j,k,nnew)* &
774!^ & vmask(Iend+1,j)
775!^
776 ad_v(iend+1,j,k,nnew)=ad_v(iend+1,j,k,nnew)* &
777 & vmask(iend+1,j)
778# endif
779!^ tl_v(Iend+1,j,k,nnew)=tl_v(Iend+1,j,k,nnew)- &
780!^ & tl_CF(Iend+1,0)
781!^
782 ad_cf(iend+1,0)=ad_cf(iend+1,0)- &
783 & ad_v(iend+1,j,k,nnew)
784 END DO
785 END IF
786 END IF
787!
788 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
789 IF (domain(ng)%Western_Edge(tile)) THEN
790 DO k=1,n(ng)
791# ifdef NEARSHORE_MELLOR
792# ifdef WET_DRY_NOT_YET
793!^ tl_v_stokes(Istr-1,j,k)=tl_v_stokes(Istr-1,j,k)* &
794!^ & vmask_wet(Istr-1,j)
795!^
796 ad_v_stokes(istr-1,j,k)=ad_v_stokes(istr-1,j,k)* &
797 & vmask_wet(istr-1,j)
798# endif
799# ifdef MASKING
800!^ tl_v_stokes(Istr-1,j,k)=tl_v_stokes(Istr-1,j,k)* &
801!^ & vmask(Istr-1,j)
802!^
803 ad_v_stokes(istr-1,j,k)=ad_v_stokes(istr-1,j,k)* &
804 & vmask(istr-1,j)
805# endif
806!^ tl_v_stokes(Istr-1,j,k)=tl_v_stokes(Istr-1,j,k)- &
807!^ & tl_CFs(Istr-1,0)
808!^
809 ad_cfs(istr-1,0)=ad_cfs(istr-1,0)- &
810 & ad_v_stokes(istr-1,j,k)
811# endif
812# ifdef WET_DRY_NOT_YET
813!^ tl_v(Istr-1,j,k,nnew)=tl_v(Istr-1,j,k,nnew)* &
814!^ & vmask_wet(Istr-1,j)
815!^
816 ad_v(istr-1,j,k,nnew)=ad_v(istr-1,j,k,nnew)* &
817 & vmask_wet(istr-1,j)
818# endif
819# ifdef MASKING
820!^ tl_v(Istr-1,j,k,nnew)=tl_v(Istr-1,j,k,nnew)* &
821!^ & vmask(Istr-1,j)
822!^
823 ad_v(istr-1,j,k,nnew)=ad_v(istr-1,j,k,nnew)* &
824 & vmask(istr-1,j)
825# endif
826!^ tl_v(Istr-1,j,k,nnew)=tl_v(Istr-1,j,k,nnew)- &
827!^ & tl_CF(Istr-1,0)
828!^
829 ad_cf(istr-1,0)=ad_cf(istr-1,0)- &
830 & ad_v(istr-1,j,k,nnew)
831 END DO
832 END IF
833 END IF
834
835# ifdef DIAGNOSTICS_UV
836!!
837!! Convert the units of the fast-time integrated change in mass flux
838!! diagnostic values to change in velocity (m s-1).
839!!
840!! DO idiag=1,NDM2d-1
841!! DO i=IstrT,IendT
842# ifdef MASKING
843!! DiaV2wrk(i,j,idiag)=DiaV2wrk(i,j,idiag)*vmask(i,j)
844# endif
845!! DiaV2wrk(i,j,idiag)=DC(i,0)*DiaV2wrk(i,j,idiag)
846!! END DO
847!! END DO
848# endif
849!
850! Save adjoint solution for time-step iic(ng)-1 as the sum of time
851! records 1 and 2.
852!
853 DO i=istrt,iendt
854 ad_vbar_sol(i,j)=ad_vbar(i,j,1)+ad_vbar(i,j,2)
855 END DO
856!
857! Compute adjoint thicknesses of V-boxes DC(i,1:N), total depth of the
858! water column DC(i,0), and incorrect vertical mean CF(i,0). Notice
859! that barotropic component is replaced with its fast-time averaged
860! values.
861!
862 DO i=istrt,iendt
863# ifdef DIAGNOSTICS_UV
864!! DiaV2int(i,j,M2rate)=vbar(i,j,1) ! HGA check
865!! DiaV2wrk(i,j,M2rate)=vbar(i,j,1)-DiaV2int(i,j,M2rate)
866!! DiaV2int(i,j,M2rate)=vbar(i,j,1)*DC1(i,0)
867!! DiaV2wrk(i,j,M2rate)=vbar(i,j,1)- &
868!! & DiaV2int(i,j,M2rate)*DC(i,0)
869# endif
870!^ tl_vbar(i,j,2)=tl_vbar(i,j,1)
871!^
872 ad_vbar(i,j,1)=ad_vbar(i,j,1)+ad_vbar(i,j,2)
873 ad_vbar(i,j,2)=0.0_r8
874!^ tl_vbar(i,j,1)=tl_DC(i,0)*DV_avg1(i,j)+ &
875!^ & DC(i,0)*tl_DV_avg1(i,j)
876!^
877 ad_dc(i,0)=ad_dc(i,0)+ad_vbar(i,j,1)*dv_avg1(i,j)
878 ad_dv_avg1(i,j)=ad_dv_avg1(i,j)+ad_vbar(i,j,1)*dc(i,0)
879 ad_vbar(i,j,1)=0.0_r8
880# ifdef NEARSHORE_MELLOR
881!^ tl_CFs(i,0)=tl_DC(i,0)*(CFs1(i)-cff2)+ &
882!^ & DC(i,0)*(tl_CFs(i,0)-tl_cff2)
883!^
884 adfac=dc(i,0)*ad_cfs(i,0)
885 ad_dc(i,0)=ad_dc(i,0)+(cfs1(i)-cff2)*ad_cfs(i,0)
886 ad_cff2=ad_cff2-adfac
887 ad_cfs(i,0)=ad_cfs(i,0)+adfac
888 ad_cfs(i,0)=0.0_r8
889# endif
890!^ tl_CF(i,0)=tl_DC(i,0)*(CF1(i)-DV_avg1(i,j))+ &
891!^ & DC(i,0)*(tl_CF(i,0)-tl_DV_avg1(i,j))
892!^
893 adfac=dc(i,0)*ad_cf(i,0)
894 ad_dc(i,0)=ad_dc(i,0)+ad_cf(i,0)*(cf1(i)-dv_avg1(i,j))
895 ad_dv_avg1(i,j)=ad_dv_avg1(i,j)-adfac
896 ad_cf(i,0)=adfac
897!^ tl_DC(i,0)=-tl_DC(i,0/(DC1(i,0)*DC1(i,0))
898!^
899 ad_dc(i,0)=-ad_dc(i,0)/(dc1(i,0)*dc1(i,0))
900# ifdef NEARSHORE_MELLOR
901!^ tl_cff2=tl_DC(i,0)*vbar_stokes(i,j)+ &
902!^ & DC(i,0)*tl_vbar_stokes(i,j))
903!^
904 ad_vbar_stokes(i,j)=ad_vbar_stokes(i,j)+dc(i,0)*ad_cff2
905 ad_dc(i,0)=ad_dc(i,0)+vbar_stokes(i,j)*ad_cff2
906 ad_cff2=0.0_r8
907# endif
908 END DO
909 DO i=istrt,iendt
910 dc(i,0)=0.0_r8
911 cf(i,0)=0.0_r8
912# ifdef NEARSHORE_MELLOR
913 cfs(i,0)=0.0_r8
914# endif
915 fc(i,0)=0.0_r8
916 END DO
917 DO k=1,n(ng)
918 DO i=istrt,iendt
919 cff=0.5_r8*om_v(i,j)
920 dc(i,k)=cff*(hz(i,j,k)+hz(i,j-1,k))
921 dc(i,0)=dc(i,0)+dc(i,k)
922 cf(i,0)=cf(i,0)+ &
923 & dc(i,k)*v(i,j,k,nnew)
924# ifdef NEARSHORE_MELLOR
925 cfs(i,0)=cfs(i,0)+dc(i,k)*v_stokes(i,j,k)
926!^ tl_CFs(i,0)=tl_CFs(i,0)+ &
927!^ & tl_DC(i,k)*v_stokes(i,j,k)+ &
928!^ & DC(i,k)*tl_v_stokes(i,j,k)
929!^
930 ad_dc(i,k)=ad_dc(i,k)+v_stokes(i,j,k)*ad_cfs(i,0)
931 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)+dc(i,k)*ad_cfs(i,0)
932# endif
933!^ tl_CF(i,0)=tl_CF(i,0)+ &
934!^ & tl_DC(i,k)*v(i,j,k,nnew)+ &
935!^ & DC(i,k)*tl_v(i,j,k,nnew)
936!^
937 ad_dc(i,k)=ad_dc(i,k)+ad_cf(i,0)*v(i,j,k,nnew)
938 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)+ad_cf(i,0)*dc(i,k)
939!^ tl_DC(i,0)=tl_DC(i,0)+tl_DC(i,k)
940!^
941 ad_dc(i,k)=ad_dc(i,k)+ad_dc(i,0)
942!^ tl_DC(i,k)=cff*(tl_Hz(i,j,k)+tl_Hz(i,j-1,k))
943!^
944 adfac=cff*ad_dc(i,k)
945 ad_hz(i,j-1,k)=ad_hz(i,j-1,k)+adfac
946 ad_hz(i,j ,k)=ad_hz(i,j ,k)+adfac
947 ad_dc(i,k)=0.0_r8
948 END DO
949 END DO
950 DO i=istrt,iendt
951!^ tl_FC(i,0)=0.0_r8
952!^
953 ad_fc(i,0)=0.0_r8
954# ifdef NEARSHORE_MELLOR
955!^ tl_CFs(i,0)=0.0_r8
956!^
957 ad_cfs(i,0)=0.0_r8
958# endif
959!^ tl_CF(i,0)=0.0_r8
960!^
961 ad_cf(i,0)=0.0_r8
962!^ tl_DC(i,0)=0.0_r8
963!^
964 ad_dc(i,0)=0.0_r8
965 END DO
966 END IF
967!
968! Couple adjoint velocity component in the XI-direction. First,
969! compute BASIC STATE quantities. This section was computed three
970! times in the original code. This can avoided by saving the
971! intermediate values in scratch arrays DC1, CF1, and FC1 (HGA).
972!
973 DO i=istrp,iendt
974 dc(i,0)=0.0_r8
975 cf(i,0)=0.0_r8
976# ifdef NEARSHORE_MELLOR
977 cfs(i,0)=0.0_r8
978# endif
979 fc(i,0)=0.0_r8
980 END DO
981 DO k=1,n(ng)
982 DO i=istrp,iendt
983 cff=0.5_r8*on_u(i,j)
984 dc(i,k)=cff*(hz(i,j,k)+hz(i-1,j,k))
985 dc(i,0)=dc(i,0)+dc(i,k)
986 cf(i,0)=cf(i,0)+ &
987 & dc(i,k)*u(i,j,k,nnew)
988# ifdef NEARSHORE_MELLOR
989 cfs(i,0)=cfs(i,0)+ &
990 & dc(i,k)*u_stokes(i,j,k)
991# endif
992 END DO
993 END DO
994 DO i=istrp,iendt
995 dc1(i,0)=dc(i,0) ! intermediate
996# ifdef NEARSHORE_MELLOR
997 cff2=dc(i,0)*ubar_stokes(i,j)
998# endif
999 dc(i,0)=1.0_r8/dc(i,0) ! recursive
1000 cf1(i)=cf(i,0) ! intermediate
1001 cf(i,0)=dc(i,0)*(cf(i,0)-du_avg1(i,j)) ! recursive
1002# ifdef NEARSHORE_MELLOR
1003 cfs1(i)=cfs(i,0) ! intermediate
1004 cfs(i,0)=dc(i,0)*(cfs(i,0)-cff2) ! recursive
1005# endif
1006 END DO
1007 DO k=n(ng),1,-1
1008 DO i=istrp,iendt
1009!^ Huon(i,j,k)=0.5_r8*(Huon(i,j,k)+u(i,j,k,nnew)*DC(i,k))
1010# ifdef NEARSHORE_MELLOR
1011!^ Huon(i,j,k)=Huon(i,j,k)+0.5_r8*u_stokes(i,j,k)*DC(i,k)
1012# endif
1013!^ FC(i,0)=FC(i,0)+Huon(i,j,k)
1014!^
1015 fc(i,0)=fc(i,0)+ &
1016 & 0.5_r8*(huon(i,j,k)+u(i,j,k,nnew)*dc(i,k))
1017# ifdef NEARSHORE_MELLOR
1018 fc(i,0)=fc(i,0)+ &
1019 & 0.5_r8*u_stokes(i,j,k)*dc(i,k)
1020# endif
1021 END DO
1022 END DO
1023 DO i=istrp,iendt
1024 fc1(i)=fc(i,0) ! intermediate
1025 fc(i,0)=dc(i,0)*(fc(i,0)-du_avg2(i,j)) ! recursive
1026 END DO
1027!
1028! Compute correct adjoint mass flux, Hz*u/n.
1029!
1030 DO k=1,n(ng)
1031 DO i=istrp,iendt
1032!^ tl_Huon(i,j,k)=tl_Huon(i,j,k)- &
1033!^ & tl_DC(i,k)*FC(i,0)- &
1034!^ & DC(i,k)*tl_FC(i,0)
1035!^
1036 ad_dc(i,k)=ad_dc(i,k)-ad_huon(i,j,k)*fc(i,0)
1037 ad_fc(i,0)=ad_fc(i,0)-ad_huon(i,j,k)*dc(i,k)
1038 END DO
1039 END DO
1040 DO i=istrp,iendt
1041!^ tl_FC(i,0)=tl_DC(i,0)*(FC1(i)-DU_avg2(i,j))+ &
1042!^ & DC(i,0)*(tl_FC(i,0)-tl_DU_avg2(i,j))
1043!^
1044 adfac=dc(i,0)*ad_fc(i,0)
1045 ad_dc(i,0)=ad_dc(i,0)+ad_fc(i,0)*(fc1(i)-du_avg2(i,j))
1046 ad_du_avg2(i,j)=ad_du_avg2(i,j)-adfac
1047 ad_fc(i,0)=adfac
1048 END DO
1049!^ DO k=N(ng),1,-1
1050!^
1051 DO k=1,n(ng)
1052 DO i=istrp,iendt
1053# ifdef DIAGNOSTICS_UV
1054!! DiaU3wrk(i,j,k,M3rate)=u(i,j,k,nnew)-DiaU3wrk(i,j,k,M3rate)
1055# endif
1056!^ tl_FC(i,0)=tl_FC(i,0)+tl_Huon(i,j,k)
1057!^
1058 ad_huon(i,j,k)=ad_huon(i,j,k)+ad_fc(i,0)
1059# ifdef NEARSHORE_MELLOR
1060!^ tl_Huon(i,j,k)=tl_Huon(i,j,k)+ &
1061!^ & 0.5_r8*(tl_u_stokes(i,j,k)*DC(i,k)+ &
1062!^ & u_stokes(i,j,k)*tl_DC(i,k))
1063!^
1064 adfac=0.5_r8*ad_huon(i,j,k)
1065 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)+dc(i,k)*adfac
1066 ad_dc(i,k)=ad_dc(i,k)+u_stokes(i,j,k)*adfac
1067# endif
1068!^ tl_Huon(i,j,k)=0.5_r8*(tl_Huon(i,j,k)+ &
1069!^ & tl_u(i,j,k,nnew)*DC(i,k)+ &
1070!^ & u(i,j,k,nnew)*tl_DC(i,k))
1071!^
1072 adfac=0.5_r8*ad_huon(i,j,k)
1073 ad_dc(i,k)=ad_dc(i,k)+u(i,j,k,nnew)*adfac
1074 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)+dc(i,k)*adfac
1075 ad_huon(i,j,k)=adfac
1076 END DO
1077 END DO
1078!
1079! Replace only BOUNDARY POINTS incorrect vertical mean with more
1080! accurate barotropic component, ubar=DU_avg1/(D*on_u). Recall that,
1081! D=CF(:,0). The replacement is avoided if the boundary edge is
1082! periodic. The J-loop is pipelined, so we need to do a special
1083! test for the southern and northern domain edges.
1084!
1085 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1086 IF (j.eq.mm(ng)+1) THEN ! northern boundary
1087 DO k=1,n(ng) ! J-loop piplined
1088 DO i=istru,iend
1089# ifdef NEARSHORE_MELLOR
1090# ifdef WET_DRY_NOT_YET
1091!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)* &
1092!^ & umask_wet(i,j)
1093!^
1094 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)* &
1095 & umask_wet(i,j)
1096# endif
1097# ifdef MASKING
1098!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)* &
1099!^ & umask(i,j)
1100!^
1101 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)* &
1102 & umask(i,j)
1103# endif
1104!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)- &
1105!^ & tl_CFs(i,0)
1106!^
1107 ad_cfs(i,0)=ad_cfs(i,0)-ad_u_stokes(i,j,k)
1108# endif
1109# ifdef WET_DRY_NOT_YET
1110!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)* &
1111!^ & umask_wet(i,j)
1112!^
1113 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)* &
1114 & umask_wet(i,j)
1115# endif
1116# ifdef MASKING
1117!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)* &
1118!^ & umask(i,j)
1119!^
1120 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)* &
1121 & umask(i,j)
1122# endif
1123!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)-tl_CF(i,0)
1124!^
1125 ad_cf(i,0)=ad_cf(i,0)-ad_u(i,j,k,nnew)
1126 END DO
1127 END DO
1128 END IF
1129 END IF
1130!
1131 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1132 IF (j.eq.0) THEN ! southern boundary
1133 DO k=1,n(ng) ! J-loop pipelined
1134 DO i=istru,iend
1135# ifdef NEARSHORE_MELLOR
1136# ifdef WET_DRY_NOT_YET
1137!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)* &
1138!^ & umask_wet(i,j)
1139!^
1140 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)* &
1141 & umask_wet(i,j)
1142# endif
1143# ifdef MASKING
1144!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)* &
1145!^ & umask(i,j)
1146!^
1147 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)* &
1148 & umask(i,j)
1149# endif
1150!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)- &
1151!^ & tl_CFs(i,0)
1152!^
1153 ad_cfs(i,0)=ad_cfs(i,0)-ad_u_stokes(i,j,k)
1154# endif
1155# ifdef WET_DRY_NOT_YET
1156!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)* &
1157!^ & umask_wet(i,j)
1158!^
1159 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)* &
1160 & umask_wet(i,j)
1161# endif
1162# ifdef MASKING
1163!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)* &
1164!^ & umask(i,j)
1165!^
1166 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)* &
1167 & umask(i,j)
1168# endif
1169!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)-tl_CF(i,0)
1170!^
1171 ad_cf(i,0)=ad_cf(i,0)-ad_u(i,j,k,nnew)
1172 END DO
1173 END DO
1174 END IF
1175 END IF
1176!
1177 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1178 IF (domain(ng)%Eastern_Edge(tile)) THEN
1179 DO k=1,n(ng)
1180# ifdef NEARSHORE_MELLOR
1181# ifdef WET_DRY_NOT_YET
1182!^ tl_u_stokes(Iend+1,j,k)=tl_u_stokes(Iend+1,j,k)* &
1183!^ & umask_wet(Iend+1,j)
1184!^
1185 ad_u_stokes(iend+1,j,k)=ad_u_stokes(iend+1,j,k)* &
1186 & umask_wet(iend+1,j)
1187# endif
1188# ifdef MASKING
1189!^ tl_u_stokes(Iend+1,j,k)=tl_u_stokes(Iend+1,j,k)* &
1190!^ & umask(Iend+1,j)
1191!^
1192 ad_u_stokes(iend+1,j,k)=ad_u_stokes(iend+1,j,k)* &
1193 & umask(iend+1,j)
1194# endif
1195!^ tl_u_stokes(Iend+1,j,k)=tl_u_stokes(Iend+1,j,k)- &
1196!^ & tl_CFs(Iend+1,0)
1197!^
1198 ad_cfs(iend+1,0)=ad_cfs(iend+1,0)- &
1199 & ad_u_stokes(iend+1,j,k)
1200# endif
1201# ifdef WET_DRY_NOT_YET
1202!^ tl_u(Iend+1,j,k,nnew)=tl_u(Iend+1,j,k,nnew)* &
1203!^ & umask_wet(Iend+1,j)
1204!^
1205 ad_u(iend+1,j,k,nnew)=ad_u(iend+1,j,k,nnew)* &
1206 & umask_wet(iend+1,j)
1207# endif
1208# ifdef MASKING
1209!^ tl_u(Iend+1,j,k,nnew)=tl_u(Iend+1,j,k,nnew)* &
1210!^ & umask(Iend+1,j)
1211!^
1212 ad_u(iend+1,j,k,nnew)=ad_u(iend+1,j,k,nnew)* &
1213 & umask(iend+1,j)
1214# endif
1215!^ tl_u(Iend+1,j,k,nnew)=tl_u(Iend+1,j,k,nnew)- &
1216!^ & tl_CF(Iend+1,0)
1217!^
1218 ad_cf(iend+1,0)=ad_cf(iend+1,0)- &
1219 & ad_u(iend+1,j,k,nnew)
1220 END DO
1221 END IF
1222 END IF
1223!
1224 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1225 IF (domain(ng)%Western_Edge(tile)) THEN
1226 DO k=1,n(ng)
1227# ifdef NEARSHORE_MELLOR
1228# ifdef WET_DRY_NOT_YET
1229!^ tl_u_stokes(Istr,j,k)=tl_u_stokes(Istr,j,k)* &
1230!^ & umask_wet(Istr,j)
1231!^
1232 ad_u_stokes(istr,j,k)=ad_u_stokes(istr,j,k)* &
1233 & umask_wet(istr,j)
1234# endif
1235# ifdef MASKING
1236!^ tl_u_stokes(Istr,j,k)=tl_u_stokes(Istr,j,k)* &
1237!^ & umask(Istr,j)
1238!^
1239 ad_u_stokes(istr,j,k)=ad_u_stokes(istr,j,k)* &
1240 & umask(istr,j)
1241# endif
1242!^ tl_u_stokes(Istr,j,k)=tl_u_stokes(Istr,j,k)- &
1243!^ & tl_CFs(Istr,0)
1244!^
1245 ad_cfs(istr,0)=ad_cfs(istr,0)- &
1246 & ad_u_stokes(istr,j,k)
1247# endif
1248# ifdef WET_DRY_NOT_YET
1249!^ tl_u(Istr,j,k,nnew)=tl_u(Istr,j,k,nnew)* &
1250!^ & umask_wet(Istr,j)
1251!^
1252 ad_u(istr,j,k,nnew)=ad_u(istr,j,k,nnew)* &
1253 & umask_wet(istr,j)
1254# endif
1255# ifdef MASKING
1256!^ tl_u(Istr,j,k,nnew)=tl_u(Istr,j,k,nnew)* &
1257!^ & umask(Istr,j)
1258!^
1259 ad_u(istr,j,k,nnew)=ad_u(istr,j,k,nnew)* &
1260 & umask(istr,j)
1261# endif
1262!^ tl_u(Istr,j,k,nnew)=tl_u(Istr,j,k,nnew)- &
1263!^ & tl_CF(Istr,0)
1264!^
1265 ad_cf(istr,0)=ad_cf(istr,0)- &
1266 & ad_u(istr,j,k,nnew)
1267 END DO
1268 END IF
1269 END IF
1270
1271# ifdef DIAGNOSTICS_UV
1272!!
1273!! Convert the units of the fast-time integrated change in mass flux
1274!! diagnostic values to change in velocity (m s-1).
1275!!
1276!! DO idiag=1,NDM2d-1
1277!! DO i=IstrP,IendT
1278# ifdef MASKING
1279!! DiaU2wrk(i,j,idiag)=DiaU2wrk(i,j,idiag)*umask(i,j)
1280# endif
1281!! DiaU2wrk(i,j,idiag)=DC(i,0)*DiaU2wrk(i,j,idiag)
1282!! END DO
1283!! END DO
1284# endif
1285!
1286! Save adjoint solution for time-step iic(ng)-1 as the sum of time
1287! records 1 and 2.
1288!
1289 DO i=istrp,iendt
1290 ad_ubar_sol(i,j)=ad_ubar(i,j,1)+ad_ubar(i,j,2)
1291 END DO
1292!
1293! Compute thicknesses of U-boxes DC(i,1:N), total depth of the water
1294! column DC(i,0), and incorrect vertical mean CF(i,0). Notice that
1295! barotropic component is replaced with its fast-time averaged
1296! values.
1297!
1298 DO i=istrp,iendt
1299# ifdef DIAGNOSTICS_UV
1300!! DiaU2int(i,j,M2rate)=ubar(i,j,1)*DC1(i,0)
1301!! DiaU2wrk(i,j,M2rate)=ubar(i,j,1)-DiaU2int(i,j,M2rate)*DC(i,0)
1302# endif
1303!^ tl_ubar(i,j,2)=tl_ubar(i,j,1)
1304!^
1305 ad_ubar(i,j,1)=ad_ubar(i,j,1)+ad_ubar(i,j,2)
1306 ad_ubar(i,j,2)=0.0_r8
1307!^ tl_ubar(i,j,1)=tl_DC(i,0)*DU_avg1(i,j)+ &
1308!^ & DC(i,0)*tl_DU_avg1(i,j)
1309!^
1310 ad_dc(i,0)=ad_dc(i,0)+ad_ubar(i,j,1)*du_avg1(i,j)
1311 ad_du_avg1(i,j)=ad_du_avg1(i,j)+ad_ubar(i,j,1)*dc(i,0)
1312 ad_ubar(i,j,1)=0.0_r8
1313# ifdef NEARSHORE_MELLOR
1314!^ tl_CFs(i,0)=tl_DC(i,0)*(CFs1(i)-cff2)+ &
1315!^ & DC(i,0)*(tl_CFs(i,0)-tl_cff2)
1316!^
1317 adfac=dc(i,0)*ad_cfs(i,0)
1318 ad_dc(i,0)=ad_dc(i,0)+(cfs1(i)-cff2)*ad_cfs(i,0)
1319 ad_cfs(i,0)=ad_cfs(i,0)+adfac
1320 ad_cff2=tl_cff2-adfac
1321 ad_cfs(i,0)=0.0_r8
1322# endif
1323!^ tl_CF(i,0)=tl_DC(i,0)*(CF1(i)-DU_avg1(i,j))+ &
1324!^ & DC(i,0)*(tl_CF(i,0)-tl_DU_avg1(i,j))
1325!^
1326 adfac=dc(i,0)*ad_cf(i,0)
1327 ad_dc(i,0)=ad_dc(i,0)+ad_cf(i,0)*(cf1(i)-du_avg1(i,j))
1328 ad_du_avg1(i,j)=ad_du_avg1(i,j)-adfac
1329 ad_cf(i,0)=adfac
1330!^ tl_DC(i,0)=-tl_DC(i,0)/(DC1(i,0)*DC1(i,0))
1331!^
1332 ad_dc(i,0)=-ad_dc(i,0)/(dc1(i,0)*dc1(i,0))
1333# ifdef NEARSHORE_MELLOR
1334!^ tl_cff2=tl_DC(i,0)*ubar_stokes(i,j)+ &
1335!^ & DC(i,0)*tl_ubar_stokes(i,j)
1336!^
1337 ad_dc(i,0)=ad_dc(i,0)+ubar_stokes(i,j)*ad_cff2
1338 ad_ubar_stokes(i,j)=ad_ubar_stokes(i,j)+dc(i,0)*ad_cff2
1339 ad_cff2=0.0_r8
1340# endif
1341 END DO
1342 DO i=istrp,iendt
1343 dc(i,0)=0.0_r8
1344 cf(i,0)=0.0_r8
1345# ifdef NEARSHORE_MELLOR
1346 cfs(i,0)=0.0_r8
1347# endif
1348 fc(i,0)=0.0_r8
1349 END DO
1350 DO k=1,n(ng)
1351 DO i=istrp,iendt
1352 cff=0.5_r8*on_u(i,j)
1353 dc(i,k)=0.5_r8*(hz(i,j,k)+hz(i-1,j,k))*on_u(i,j)
1354 dc(i,0)=dc(i,0)+dc(i,k)
1355 cf(i,0)=cf(i,0)+ &
1356 & dc(i,k)*u(i,j,k,nnew)
1357# ifdef NEARSHORE_MELLOR
1358 cfs(i,0)=cfs(i,0)+ &
1359 & dc(i,k)*u_stokes(i,j,k)
1360!^ tl_CFs(i,0)=tl_CFs(i,0)+ &
1361!^ & tl_DC(i,k)*u_stokes(i,j,k)+ &
1362!^ & DC(i,k)*tl_u_stokes(i,j,k)
1363!^
1364 ad_dc(i,k)=ad_dc(i,k)+u_stokes(i,j,k)*ad_cfs(i,0)
1365 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)+dc(i,k)*ad_cfs(i,0)
1366# endif
1367!^ tl_CF(i,0)=tl_CF(i,0)+ &
1368!^ & tl_DC(i,k)*u(i,j,k,nnew)+ &
1369!^ & DC(i,k)*tl_u(i,j,k,nnew)
1370!^
1371 ad_dc(i,k)=ad_dc(i,k)+ad_cf(i,0)*u(i,j,k,nnew)
1372 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)+ad_cf(i,0)*dc(i,k)
1373!^ tl_DC(i,0)=tl_DC(i,0)+tl_DC(i,k)
1374!^
1375 ad_dc(i,k)=ad_dc(i,k)+ad_dc(i,0)
1376!^ tl_DC(i,k)=cff*(tl_Hz(i,j,k)+tl_Hz(i-1,j,k))
1377!^
1378 adfac=cff*ad_dc(i,k)
1379 ad_hz(i-1,j,k)=ad_hz(i-1,j,k)+adfac
1380 ad_hz(i ,j,k)=ad_hz(i ,j,k)+adfac
1381 ad_dc(i,k)=0.0_r8
1382 END DO
1383 END DO
1384 DO i=istrp,iendt
1385!^ tl_FC(i,0)=0.0_r8
1386!^
1387 ad_fc(i,0)=0.0_r8
1388# ifdef NEARSHORE_MELLOR
1389!^ tl_CFs(i,0)=0.0_r8
1390!^
1391 ad_cfs(i,0)=0.0_r8
1392# endif
1393!^ tl_CF(i,0)=0.0_r8
1394!^
1395 ad_cf(i,0)=0.0_r8
1396!^ tl_DC(i,0)=0.0_r8
1397!^
1398 ad_dc(i,0)=0.0_r8
1399 END DO
1400 END DO j_loop1
1401!
1402!-----------------------------------------------------------------------
1403! Apply adjoint momentum transport point sources (like river runoff),
1404! if any.
1405!-----------------------------------------------------------------------
1406!
1407 IF (luvsrc(ng)) THEN
1408 DO is=1,nsrc(ng)
1409 i=sources(ng)%Isrc(is)
1410 j=sources(ng)%Jsrc(is)
1411 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1412 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1413 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
1414 DO k=1,n(ng)
1415 cff1=1.0_r8/(on_u(i,j)* &
1416 & 0.5_r8*(z_w(i-1,j,k)-z_w(i-1,j,k-1)+ &
1417 & z_w(i ,j,k)-z_w(i ,j,k-1)))
1418!^ tl_u(i,j,k,nnew)=SOURCES(ng)%tl_Qsrc(is,k)*cff1+ &
1419!^ & SOURCES(ng)%Qsrc(is,k)*tl_cff1
1420!^
1421 sources(ng)%ad_Qsrc(is,k)=sources(ng)%ad_Qsrc(is,k)+ &
1422 & cff1*ad_u(i,j,k,nnew)
1423 ad_cff1=ad_cff1+ &
1424 & sources(ng)%Qsrc(is,k)*ad_u(i,j,k,nnew)
1425 ad_u(i,j,k,nnew)=0.0_r8
1426!^ tl_cff1=-cff1*cff1*on_u(i,j)* &
1427!^ & 0.5_r8*(tl_z_w(i-1,j,k)-tl_z_w(i-1,j,k-1)+ &
1428!^ & tl_z_w(i ,j,k)-tl_z_w(i ,j,k-1))
1429!^
1430 adfac=-0.5_r8*cff1*cff1*on_u(i,j)*ad_cff1
1431 ad_z_w(i-1,j,k-1)=ad_z_w(i-1,j,k-1)-adfac
1432 ad_z_w(i ,j,k-1)=ad_z_w(i ,j,k-1)-adfac
1433 ad_z_w(i-1,j,k )=ad_z_w(i-1,j,k )+adfac
1434 ad_z_w(i ,j,k )=ad_z_w(i ,j,k )+adfac
1435 ad_cff1=0.0_r8
1436 END DO
1437 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
1438 DO k=1,n(ng)
1439 cff1=1.0_r8/(om_v(i,j)* &
1440 & 0.5_r8*(z_w(i,j-1,k)-z_w(i,j-1,k-1)+ &
1441 & z_w(i,j ,k)-z_w(i,j ,k-1)))
1442!^ tl_v(i,j,k,nnew)=SOURCES(ng)%tl_Qsrc(is,k)*cff1+ &
1443!^ & SOURCES(ng)%Qsrc(is,k)*tl_cff1
1444!^
1445 sources(ng)%ad_Qsrc(is,k)=sources(ng)%ad_Qsrc(is,k)+ &
1446 & cff1*ad_v(i,j,k,nnew)
1447 ad_cff1=ad_cff1+ &
1448 & sources(ng)%Qsrc(is,k)*ad_v(i,j,k,nnew)
1449 ad_v(i,j,k,nnew)=0.0_r8
1450!^ tl_cff1=-cff1*cff1*om_v(i,j)* &
1451!^ & 0.5_r8*(tl_z_w(i,j-1,k)-tl_z_w(i,j-1,k-1)+ &
1452!^ & tl_z_w(i,j ,k)-tl_z_w(i,j ,k-1))
1453!^
1454 adfac=-0.5_r8*cff1*cff1*om_v(i,j)*ad_cff1
1455 ad_z_w(i,j-1,k-1)=ad_z_w(i,j-1,k-1)-adfac
1456 ad_z_w(i,j ,k-1)=ad_z_w(i,j ,k-1)-adfac
1457 ad_z_w(i,j-1,k )=ad_z_w(i,j-1,k )+adfac
1458 ad_z_w(i,j ,k )=ad_z_w(i,j ,k )+adfac
1459 ad_cff1=0.0_r8
1460 END DO
1461 END IF
1462 END IF
1463 END DO
1464 END IF
1465!
1466!-----------------------------------------------------------------------
1467! Set lateral adjoint boundary conditions.
1468!-----------------------------------------------------------------------
1469!
1470!^ CALL tl_v3dbc_tile (ng, tile, &
1471!^ & LBi, UBi, LBj, UBj, N(ng), &
1472!^ & IminS, ImaxS, JminS, JmaxS, &
1473!^ & nstp, nnew, &
1474!^ & v)
1475!^
1476 CALL ad_v3dbc_tile (ng, tile, &
1477 & lbi, ubi, lbj, ubj, n(ng), &
1478 & imins, imaxs, jmins, jmaxs, &
1479 & nstp, nnew, &
1480 & ad_v)
1481!^ CALL tl_u3dbc_tile (ng, tile, &
1482!^ & LBi, UBi, LBj, UBj, N(ng), &
1483!^ & IminS, ImaxS, JminS, JmaxS, &
1484!^ & nstp, nnew, &
1485!^ & tl_u)
1486!^
1487 CALL ad_u3dbc_tile (ng, tile, &
1488 & lbi, ubi, lbj, ubj, n(ng), &
1489 & imins, imaxs, jmins, jmaxs, &
1490 & nstp, nnew, &
1491 & ad_u)
1492!
1493!-----------------------------------------------------------------------
1494! Time step momentum equation in the ETA-direction.
1495!-----------------------------------------------------------------------
1496!
1497 j_loop2 : DO j=jstr,jend
1498 IF (j.ge.jstrv) THEN
1499!
1500! Compute BASIC STATE quantities.
1501!
1502 DO i=istr,iend
1503 ak(i,0)=0.5_r8*(akv(i,j-1,0)+ &
1504 & akv(i,j ,0))
1505 DO k=1,n(ng)
1506 ak(i,k)=0.5_r8*(akv(i,j-1,k)+ &
1507 & akv(i,j ,k))
1508 hzk(i,k)=0.5_r8*(hz(i,j-1,k)+ &
1509 & hz(i,j ,k))
1510# if defined SPLINES_VVISC || defined DIAGNOSTICS_UV
1511 ohz(i,k)=1.0_r8/hzk(i,k)
1512# endif
1513 END DO
1514 END DO
1515!
1516! Couple and update new adjoint solution.
1517!
1518# if defined DIAGNOSTICS_UV && defined MASKING
1519!! DO k=1,N(ng)
1520!! DO i=Istr,Iend
1521!! DO idiag=1,NDM3d
1522!! DiaV3wrk(i,j,k,idiag)=DiaV3wrk(i,j,k,idiag)*vmask(i,j)
1523!! END DO
1524!! END DO
1525!! END DO
1526# endif
1527
1528 DO k=1,n(ng)
1529 DO i=istr,iend
1530# ifdef DIAGNOSTICS_UV
1531# ifdef NEARSHORE_MELLOR
1532!! DiaV3wrk(i,j,k,M3hrad)=DiaV3wrk(i,j,k,M3hrad)- &
1533!! & Dwrk(i,M2hrad)
1534# endif
1535# ifdef UV_ADV
1536!! DiaV3wrk(i,j,k,M3hadv)=DiaV3wrk(i,j,k,M3hadv)- &
1537!! & Dwrk(i,M2hadv)
1538!! DiaV3wrk(i,j,k,M3yadv)=DiaV3wrk(i,j,k,M3yadv)- &
1539!! & Dwrk(i,M2yadv)
1540!! DiaV3wrk(i,j,k,M3xadv)=DiaV3wrk(i,j,k,M3xadv)- &
1541!! & Dwrk(i,M2xadv)
1542# endif
1543# if defined UV_VIS2 || defined UV_VIS4
1544!! DiaV3wrk(i,j,k,M3hvis)=DiaV3wrk(i,j,k,M3hvis)- &
1545!! & Dwrk(i,M2hvis)
1546!! DiaV3wrk(i,j,k,M3yvis)=DiaV3wrk(i,j,k,M3yvis)- &
1547!! & Dwrk(i,M2yvis)
1548!! DiaV3wrk(i,j,k,M3xvis)=DiaV3wrk(i,j,k,M3xvis)- &
1549!! & Dwrk(i,M2xvis)
1550# endif
1551# ifdef UV_COR
1552!! DiaV3wrk(i,j,k,M3fcor)=DiaV3wrk(i,j,k,M3fcor)- &
1553!! & Dwrk(i,M2fcor)
1554# endif
1555!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)- &
1556!! & Dwrk(i,M2bstr)
1557!! DiaV3wrk(i,j,k,M3pgrd)=DiaV3wrk(i,j,k,M3pgrd)- &
1558!! & Dwrk(i,M2pgrd)
1559# endif
1560# ifdef NEARSHORE_MELLOR
1561# ifdef WET_DRY_NOT_YET
1562!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*vmask_wet(i,j)
1563!^
1564 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)*vmask_wet(i,j)
1565# endif
1566# ifdef MASKING
1567!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*vmask(i,j)
1568!^
1569 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)*vmask(i,j)
1570# endif
1571!^ tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)-tl_DCs(i,0)
1572!^
1573 ad_dcs(i,0)=ad_dcs(i,0)-ad_v_stokes(i,j,k)
1574# endif
1575# ifdef WET_DRY_NOT_YET
1576!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*vmask_wet(i,j)
1577!^
1578 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)*vmask_wet(i,j)
1579# endif
1580# ifdef MASKING
1581!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*vmask(i,j)
1582!^
1583 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)*vmask(i,j)
1584# endif
1585!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)-tl_DC(i,0)
1586!^
1587 ad_dc(i,0)=ad_dc(i,0)-ad_v(i,j,k,nnew)
1588 END DO
1589 END DO
1590!
1591! Replace INTERIOR POINTS incorrect vertical mean with more accurate
1592! barotropic component, vbar=DV_avg1/(D*om_v). Recall that, D=CF(:,0).
1593!
1594 DO i=istr,iend
1595 cf(i,0)=hzk(i,1)
1596 dc(i,0)=v(i,j,1,nnew)*hzk(i,1)
1597 END DO
1598 DO k=2,n(ng)
1599 DO i=istr,iend
1600 cf(i,0)=cf(i,0)+hzk(i,k)
1601 dc(i,0)=dc(i,0)+v(i,j,k,nnew)*hzk(i,k)
1602 END DO
1603 END DO
1604 DO i=istr,iend
1605 dc1(i,0)=dc(i,0) ! intermediate
1606 cff1=1.0_r8/(cf(i,0)*om_v(i,j))
1607 dc(i,0)=(dc(i,0)*om_v(i,j)-dv_avg1(i,j))*cff1 ! recursive
1608# ifdef NEARSHORE_MELLOR
1609 dcs1(i)=dcs(i,0) ! intermediate
1610 cff2=1.0_r8/cf(i,0)
1611 dcs(i,0)=dcs(i,0)*cff2-vbar_stokes(i,j) ! recursive
1612# endif
1613# ifdef DIAGNOSTICS_UV
1614!! Dwrk(i,M2bstr)=(Dwrk(i,M2bstr)*om_v(i,j)- &
1615!! & DiaV2wrk(i,j,M2bstr)-DiaV2wrk(i,j,M2sstr))* &
1616!! & cff1
1617!! DO idiag=1,M2pgrd
1618!! Dwrk(i,idiag)=(Dwrk(i,idiag)*om_v(i,j)- &
1619!! & DiaV2wrk(i,j,idiag))*cff1
1620!! END DO
1621# endif
1622# ifdef NEARSHORE_MELLOR
1623!^ tl_DCs(i,0)=tl_DCs(i,0)*cff2+ &
1624!^ & DCs1(i,0)*tl_cff2-tl_vbar_stokes(i,j)
1625!^
1626 ad_vbar_stokes(i,j)=ad_vbar_stokes(i,j)-ad_dcs(i,0)
1627 ad_cff2=ad_cff2+dcs1(i,0)*ad_dcs(i,0)
1628 ad_dcs(i,0)=cff2*ad_dcs(i,0)
1629!^ tl_cff2=-cff2*cff2*tl_CF(i,0)
1630!^
1631 ad_cf(i,0)=ad_cf(i,0)-cff2*cff2*tl_cff2
1632 ad_cff2=0.0_r8
1633# endif
1634!^ tl_DC(i,0)=(tl_DC(i,0)*om_v(i,j)-tl_DV_avg1(i,j))*cff1+ &
1635!^ & (DC1(i,0)*om_v(i,j)-DV_avg1(i,j))*tl_cff1
1636!^
1637 adfac=cff1*ad_dc(i,0)
1638 ad_dv_avg1(i,j)=ad_dv_avg1(i,j)-adfac
1639 ad_cff1=ad_cff1+ &
1640 & (dc1(i,0)*om_v(i,j)-dv_avg1(i,j))*ad_dc(i,0)
1641 ad_dc(i,0)=om_v(i,j)*adfac
1642!^ tl_cff1=-cff1*cff1*tl_CF(i,0)*om_v(i,j)
1643!^
1644 ad_cf(i,0)=ad_cf(i,0)-cff1*cff1*om_v(i,j)*ad_cff1
1645 ad_cff1=0.0_r8
1646 END DO
1647 DO i=istr,iend
1648 cf(i,0)=hzk(i,1)
1649 dc(i,0)=v(i,j,1,nnew)*hzk(i,1)
1650 END DO
1651 DO k=2,n(ng)
1652 DO i=istr,iend
1653 cf(i,0)=cf(i,0)+hzk(i,k)
1654 dc(i,0)=dc(i,0)+v(i,j,k,nnew)*hzk(i,k)
1655# ifdef DIAGNOSTICS_UV
1656# ifdef NEARSHORE_MELLOR
1657!! Dwrk(i,M2hrad)=Dwrk(i,M2hrad)+ &
1658!! & DiaV3wrk(i,j,k,M3hrad)*Hzk(i,k)
1659# endif
1660# ifdef UV_ADV
1661!! Dwrk(i,M2hadv)=Dwrk(i,M2hadv)+ &
1662!! & DiaV3wrk(i,j,k,M3hadv)*Hzk(i,k)
1663!! Dwrk(i,M2yadv)=Dwrk(i,M2yadv)+ &
1664!! & DiaV3wrk(i,j,k,M3yadv)*Hzk(i,k)
1665!! Dwrk(i,M2xadv)=Dwrk(i,M2xadv)+ &
1666!! & DiaV3wrk(i,j,k,M3xadv)*Hzk(i,k)
1667# endif
1668# if defined UV_VIS2 || defined UV_VIS4
1669!! Dwrk(i,M2hvis)=Dwrk(i,M2hvis)+ &
1670!! & DiaV3wrk(i,j,k,M3hvis)*Hzk(i,k)
1671!! Dwrk(i,M2yvis)=Dwrk(i,M2yvis)+ &
1672!! & DiaV3wrk(i,j,k,M3yvis)*Hzk(i,k)
1673!! Dwrk(i,M2xvis)=Dwrk(i,M2xvis)+ &
1674!! & DiaV3wrk(i,j,k,M3xvis)*Hzk(i,k)
1675# endif
1676# ifdef UV_COR
1677!! Dwrk(i,M2fcor)=Dwrk(i,M2fcor)+ &
1678!! & DiaV3wrk(i,j,k,M3fcor)*Hzk(i,k)
1679# endif
1680!! Dwrk(i,M2bstr)=Dwrk(i,M2bstr)+ &
1681!! & DiaV3wrk(i,j,k,M3vvis)*Hzk(i,k)
1682!! Dwrk(i,M2pgrd)=Dwrk(i,M2pgrd)+ &
1683!! & DiaV3wrk(i,j,k,M3pgrd)*Hzk(i,k)
1684# endif
1685# ifdef NEARSHORE_MELLOR
1686!^ tl_DCs(i,0)=tl_DCs(i,0)+ &
1687!^ & tl_v_stokes(i,j,k)*Hzk(i,k)+ &
1688!^ & v_stokes(i,j,k)*tl_Hzk(i,k)
1689!^
1690 ad_v_stokes(i,j,k)=ad_v_stokes(i,j,k)+hzk(i,k)*ad_dcs(i,0)
1691 ad_hzk(i,k)=ad_hzk(i,k)+v_stokes(i,j,k)*ad_dcs(i,0)
1692# endif
1693!^ tl_DC(i,0)=tl_DC(i,0)+ &
1694!^ & tl_v(i,j,k,nnew)*Hzk(i,k)+ &
1695!^ & v(i,j,k,nnew)*tl_Hzk(i,k)
1696!^
1697 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)+ad_dc(i,0)*hzk(i,k)
1698 ad_hzk(i,k)=ad_hzk(i,k)+v(i,j,k,nnew)*ad_dc(i,0)
1699!^ tl_CF(i,0)=tl_CF(i,0)+tl_Hzk(i,k)
1700!^
1701 ad_hzk(i,k)=ad_hzk(i,k)+ad_cf(i,0)
1702 END DO
1703 END DO
1704 DO i=istr,iend
1705 cf(i,0)=hzk(i,1)
1706 dc(i,0)=v(i,j,1,nnew)*hzk(i,1)
1707# ifdef DIAGNOSTICS_UV
1708# ifdef NEARSHORE_MELLOR
1709!! Dwrk(i,M2hrad)=DiaV3wrk(i,j,1,M3hrad)*Hzk(i,1)
1710# endif
1711# ifdef UV_ADV
1712!! Dwrk(i,M2hadv)=DiaV3wrk(i,j,1,M3hadv)*Hzk(i,1)
1713!! Dwrk(i,M2yadv)=DiaV3wrk(i,j,1,M3yadv)*Hzk(i,1)
1714!! Dwrk(i,M2xadv)=DiaV3wrk(i,j,1,M3xadv)*Hzk(i,1)
1715# endif
1716# if defined UV_VIS2 || defined UV_VIS4
1717!! Dwrk(i,M2hvis)=DiaV3wrk(i,j,1,M3hvis)*Hzk(i,1)
1718!! Dwrk(i,M2yvis)=DiaV3wrk(i,j,1,M3yvis)*Hzk(i,1)
1719!! Dwrk(i,M2xvis)=DiaV3wrk(i,j,1,M3xvis)*Hzk(i,1)
1720# endif
1721# ifdef UV_COR
1722!! Dwrk(i,M2fcor)=DiaV3wrk(i,j,1,M3fcor)*Hzk(i,1)
1723# endif
1724!! Dwrk(i,M2bstr)=DiaV3wrk(i,j,1,M3vvis)*Hzk(i,1)
1725!! Dwrk(i,M2pgrd)=DiaV3wrk(i,j,1,M3pgrd)*Hzk(i,1)
1726# endif
1727# ifdef NEARSHORE_MELLOR
1728!^ tl_DCs(i,0)=tl_v_stokes(i,j,1)*Hzk(i,1)+ &
1729!^ & v_stokes(i,j,1)*tl_Hzk(i,1)
1730!^
1731 ad_v_stokes(i,j,1)=ad_v_stokes(i,j,1)+hzk(i,1)*ad_dcs(i,0)
1732 ad_hzk(i,1)=ad_hzk(i,1)+v_stokes(i,j,1)*ad_dcs(i,0)
1733 ad_dcs(i,0)=0.0_r8
1734# endif
1735!^ tl_DC(i,0)=tl_v(i,j,1,nnew)*Hzk(i,1)+ &
1736!^ & v(i,j,1,nnew)*tl_Hzk(i,1)
1737!^
1738 ad_v(i,j,1,nnew)=ad_v(i,j,1,nnew)+ad_dc(i,0)*hzk(i,1)
1739 ad_hzk(i,1)=ad_hzk(i,1)+v(i,j,1,nnew)*ad_dc(i,0)
1740 ad_dc(i,0)=0.0_r8
1741!^ tl_CF(i,0)=tl_Hzk(i,1)
1742!^
1743 ad_hzk(i,1)=ad_hzk(i,1)+ad_cf(i,0)
1744 ad_cf(i,0)=0.0_r8
1745 END DO
1746
1747# if defined SPLINES_VVISC
1748!
1749! Apply spline algorithm to BASIC STATE "v" which should be
1750! in units of velocity (m/s). DC will be used in the tangent
1751! linear spline algorithm below. Solve BASIC STATE tridiagonal
1752! system.
1753!
1754 cff1=1.0_r8/6.0_r8
1755 DO k=1,n(ng)-1
1756 DO i=istr,iend
1757 fc(i,k)=cff1*hzk(i,k )-dt(ng)*ak(i,k-1)*ohz(i,k )
1758 cf(i,k)=cff1*hzk(i,k+1)-dt(ng)*ak(i,k+1)*ohz(i,k+1)
1759 END DO
1760 END DO
1761 DO i=istr,iend
1762 cf(i,0)=0.0_r8
1763 dc(i,0)=0.0_r8
1764 END DO
1765!
1766! LU decomposition and forward substitution.
1767!
1768 cff1=1.0_r8/3.0_r8
1769 DO k=1,n(ng)-1
1770 DO i=istr,iend
1771 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
1772 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
1773 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
1774 cf(i,k)=cff*cf(i,k)
1775 dc(i,k)=cff*(v(i,j,k+1,nnew)-v(i,j,k,nnew)- &
1776 & fc(i,k)*dc(i,k-1))
1777 END DO
1778 END DO
1779!
1780! Backward substitution. Save DC for the adjoint below.
1781! DC is scaled later by AK.
1782!
1783 DO i=istr,iend
1784 dc(i,n(ng))=0.0_r8
1785 END DO
1786 DO k=n(ng)-1,1,-1
1787 DO i=istr,iend
1788 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1789 END DO
1790 END DO
1791!
1792! Multiply DC (BASIC STATE spline gradients) by AK and save it in DC1.
1793!
1794 DO k=0,n(ng)
1795 DO i=istr,iend
1796 dc1(i,k)=dc(i,k)*ak(i,k)
1797 END DO
1798 END DO
1799!^ DO k=1,N(ng)
1800!^
1801 DO k=n(ng),1,-1
1802 DO i=istr,iend
1803# ifdef DIAGNOSTICS_UV
1804!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+cff
1805# endif
1806!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)+tl_cff
1807!^
1808 ad_cff=ad_cff+ad_v(i,j,k,nnew)
1809!^ tl_cff=dt(ng)*(tl_oHz(i,k)*(DC(i,k)-DC(i,k-1))+ &
1810!^ & oHz(i,k)*(tl_DC(i,k)-tl_DC(i,k-1)))
1811!^ use DC1
1812 adfac=dt(ng)*ad_cff
1813 adfac1=adfac*ohz(i,k)
1814 ad_dc(i,k-1)=ad_dc(i,k-1)-adfac1
1815 ad_dc(i,k )=ad_dc(i,k )+adfac1
1816 ad_ohz(i,k)=ad_ohz(i,k)+(dc1(i,k)-dc1(i,k-1))*adfac
1817 ad_cff=0.0_r8
1818!^ tl_DC(i,k)=tl_DC(i,k)*AK(i,k)+DC(i,k)*tl_AK(i,k)
1819!^ use DC
1820 ad_ak(i,k)=ad_ak(i,k)+dc(i,k)*ad_dc(i,k)
1821 ad_dc(i,k)=ak(i,k)*ad_dc(i,k)
1822 END DO
1823 END DO
1824!
1825! Adjoint back substitution.
1826!
1827 DO k=1,n(ng)-1
1828 DO i=istr,iend
1829!^ tl_DC(i,k)=tl_DC(i,k)-CF(i,k)*tl_DC(i,k+1)
1830!^
1831 ad_dc(i,k+1)=ad_dc(i,k+1)-cf(i,k)*ad_dc(i,k)
1832 END DO
1833 END DO
1834 DO i=istr,iend
1835!^ tl_DC(i,N(ng))=0.0_r8
1836!^
1837 ad_dc(i,n(ng))=0.0_r8
1838 END DO
1839!
1840! Adjoint LU decomposition and forward substitution.
1841!
1842 cff1=1.0_r8/3.0_r8
1843 DO k=n(ng)-1,1,-1
1844 DO i=istr,iend
1845 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
1846 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
1847 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
1848!^ tl_DC(i,k)=cff*(tl_v(i,j,k+1,nnew)- &
1849!^ & tl_v(i,j,k ,nnew)- &
1850!^ & (tl_FC(i,k)*DC(i,k-1)+ &
1851!^ & tl_BC(i,k)*DC(i,k )+ &
1852!^ & tl_CF(i,k)*DC(i,k+1))- &
1853!^ & FC(i,k)*tl_DC(i,k-1))
1854!^
1855 adfac=cff*ad_dc(i,k)
1856 ad_dc(i,k-1)=ad_dc(i,k-1)-fc(i,k)*adfac
1857 ad_cf(i,k)=ad_cf(i,k)-dc(i,k+1)*adfac
1858 ad_bc(i,k)=ad_bc(i,k)-dc(i,k )*adfac
1859 ad_fc(i,k)=ad_fc(i,k)-dc(i,k-1)*adfac
1860 ad_v(i,j,k, nnew)=ad_v(i,j,k ,nnew)-adfac
1861 ad_v(i,j,k+1,nnew)=ad_v(i,j,k+1,nnew)+adfac
1862 ad_dc(i,k)=0.0_r8
1863!^ tl_BC(i,k)=cff1*(tl_Hzk(i,k)+tl_Hzk(i,k+1))+ &
1864!^ & dt(ng)*(tl_AK(i,k)*(oHz(i,k)+oHz(i,k+1))+ &
1865!^ & AK(i,k)*(tl_oHz(i,k)+tl_oHz(i,k+1)))
1866!^
1867 adfac=dt(ng)*ad_bc(i,k)
1868 adfac1=adfac*ak(i,k)
1869 adfac2=cff1*ad_bc(i,k)
1870 ad_ohz(i,k )=ad_ohz(i,k )+adfac1
1871 ad_ohz(i,k+1)=ad_ohz(i,k+1)+adfac1
1872 ad_ak(i,k)=ad_ak(i,k)+(ohz(i,k)+ohz(i,k+1))*adfac
1873 ad_hzk(i,k )=ad_hzk(i,k )+adfac2
1874 ad_hzk(i,k+1)=ad_hzk(i,k+1)+adfac2
1875 ad_bc(i,k)=0.0_r8
1876 END DO
1877 END DO
1878!
1879! Use conservative, parabolic spline reconstruction of adjoint
1880! vertical viscosity derivatives. Then, time step vertical
1881! viscosity term implicitly by solving a tridiagonal system.
1882!
1883 DO i=istr,iend
1884!^ tl_DC(i,0)=0.0_r8
1885!^
1886 ad_dc(i,0)=0.0_r8
1887!^ tl_CF(i,0)=0.0_r8
1888!^
1889 ad_cf(i,0)=0.0_r8
1890 END DO
1891 cff1=1.0_r8/6.0_r8
1892 DO k=1,n(ng)-1
1893 DO i=istr,iend
1894!^ tl_CF(i,k)=cff1*tl_Hzk(i,k+1)- &
1895!^ & dt(ng)*(tl_AK(i,k+1)*oHz(i,k+1)+ &
1896!^ & AK(i,k+1)*tl_oHz(i,k+1))
1897!^
1898 adfac=dt(ng)*ad_cf(i,k)
1899 ad_ohz(i,k+1)=ad_ohz(i,k+1)-ak(i,k+1)*adfac
1900 ad_ak(i,k+1)=ad_ak(i,k+1)-ohz(i,k+1)*adfac
1901 ad_hzk(i,k+1)=ad_hzk(i,k+1)+cff1*ad_cf(i,k)
1902 ad_cf(i,k)=0.0_r8
1903!^ tl_FC(i,k)=cff1*tl_Hzk(i,k )-
1904!^ & dt(ng)*(tl_AK(i,k-1)*oHz(i,k )+ &
1905!^ & AK(i,k-1)*tl_oHz(i,k ))
1906!^
1907 adfac=dt(ng)*ad_fc(i,k)
1908 ad_ohz(i,k)=ad_ohz(i,k)-ak(i,k-1)*adfac
1909 ad_ak(i,k-1)=ad_ak(i,k-1)-ohz(i,k)*adfac
1910 ad_hzk(i,k)=ad_hzk(i,k)+cff1*ad_fc(i,k)
1911 ad_fc(i,k)=0.0_r8
1912 END DO
1913 END DO
1914# else
1915!
1916! Compute BASIC STATE off-diagonal coefficients [lambda*dt*Akv/Hz]
1917! for the implicit vertical viscosity term at horizontal V-points
1918! and vertical W-points.
1919!
1920! NOTE: The original code solves the tridiagonal system A*v=r where
1921! A is a matrix and v and r are vectors. We need to solve the
1922! tangent linear of this system which is A*tl_v+tl_A*v=tl_r.
1923! Here, tl_A*v and tl_r are known, so we must solve for tl_v
1924! by inverting A*tl_v=tl_r-tl_A*v.
1925!
1926 cff=-lambda*dt(ng)/0.5_r8
1927 DO k=1,n(ng)-1
1928 DO i=istr,iend
1929 cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
1930 & z_r(i,j,k )-z_r(i,j-1,k ))
1931 fc(i,k)=cff*cff1*ak(i,k)
1932 END DO
1933 END DO
1934 DO i=istr,iend
1935 fc(i,0)=0.0_r8
1936 fc(i,n(ng))=0.0_r8
1937 END DO
1938 DO k=1,n(ng)
1939 DO i=istr,iend
1940 bc(i,k)=hzk(i,k)-fc(i,k)-fc(i,k-1)
1941 END DO
1942 END DO
1943 DO i=istr,iend
1944 cff=1.0_r8/bc(i,1)
1945 cf(i,1)=cff*fc(i,1)
1946 END DO
1947 DO k=2,n(ng)-1
1948 DO i=istr,iend
1949 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
1950 cf(i,k)=cff*fc(i,k)
1951 END DO
1952 END DO
1953!
1954! Compute new adjoint solution by back substitution.
1955! (DC is a tangent linear variable here).
1956!
1957!^ DO k=N(ng)-1,1,-1
1958!^
1959 DO k=1,n(ng)-1
1960 DO i=istr,iend
1961# ifdef DIAGNOSTICS_UV
1962!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+ &
1963!! & v(i,j,k,nnew)-wrk(i,k)
1964# endif
1965!^ tl_v(i,j,k,nnew)=DC(i,k)
1966!^
1967 ad_dc(i,k)=ad_dc(i,k)+ad_v(i,j,k,nnew)
1968 ad_v(i,j,k,nnew)=0.0_r8
1969!^ DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
1970!^
1971 ad_dc(i,k+1)=-cf(i,k)*ad_dc(i,k)
1972# ifdef DIAGNOSTICS_UV
1973!! wrk(i,k)=v(i,j,k,nnew)*oHz(i,k)
1974# endif
1975 END DO
1976 END DO
1977 DO i=istr,iend
1978# ifdef DIAGNOSTICS_UV
1979!! DiaV3wrk(i,j,N(ng),M3vvis)=DiaV3wrk(i,j,N(ng),M3vvis)+ &
1980!! & v(i,j,N(ng),nnew)-wrk(i,N(ng))
1981# endif
1982!^ tl_v(i,j,N(ng),nnew)=DC(i,N(ng))
1983!^
1984 ad_dc(i,n(ng))=ad_dc(i,n(ng))+ad_v(i,j,n(ng),nnew)
1985 ad_v(i,j,n(ng),nnew)=0.0_r8
1986!^ DC(i,N(ng))=(DC(i,N(ng))-FC(i,N(ng)-1)*DC(i,N(ng)-1))/ &
1987!^ & (BC(i,N(ng))-FC(i,N(ng)-1)*CF(i,N(ng)-1))
1988!^
1989 adfac=ad_dc(i,n(ng))/ &
1990 & (bc(i,n(ng))-fc(i,n(ng)-1)*cf(i,n(ng)-1))
1991 ad_dc(i,n(ng)-1)=ad_dc(i,n(ng)-1)-fc(i,n(ng)-1)*adfac
1992 ad_dc(i,n(ng) )=adfac
1993 END DO
1994!
1995! Solve the adjoint tridiagonal system.
1996! (DC is a tangent linear variable here).
1997!
1998!^ DO k=2,N(ng)-1
1999!^
2000 DO k=n(ng)-1,2,-1
2001 DO i=istr,iend
2002 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
2003!^ DC(i,k)=cff*(DC(i,k)-FC(i,k-1)*DC(i,k-1))
2004!^
2005 adfac=cff*ad_dc(i,k)
2006 ad_dc(i,k-1)=ad_dc(i,k-1)-fc(i,k-1)*adfac
2007 ad_dc(i,k )=adfac
2008 END DO
2009 END DO
2010 DO i=istr,iend
2011 cff=1.0_r8/bc(i,1)
2012!^ DC(i,1)=cff*DC(i,1)
2013!^
2014 ad_dc(i,1)=cff*ad_dc(i,1)
2015 END DO
2016 DO i=istr,iend
2017!^ DC(i,N(ng))=tl_v(i,j,N(ng),nnew)- &
2018!^ & (tl_FC(i,N(ng)-1)*v(i,j,N(ng)-1,nnew)+ &
2019!^ & tl_BC(i,N(ng) )*v(i,j,N(ng) ,nnew))
2020!^
2021 ad_bc(i,n(ng) )=-v(i,j,n(ng) ,nnew)*ad_dc(i,n(ng))
2022 ad_fc(i,n(ng)-1)=-v(i,j,n(ng)-1,nnew)*ad_dc(i,n(ng))
2023 ad_v(i,j,n(ng),nnew)=ad_dc(i,n(ng))
2024 ad_dc(i,n(ng))=0.0_r8
2025!^ DC(i,1)=tl_v(i,j,1,nnew)- &
2026!^ & (tl_BC(i,1)*v(i,j,1,nnew)+ &
2027!^ & tl_FC(i,1)*v(i,j,2,nnew))
2028!^
2029 ad_fc(i,1)=-v(i,j,2,nnew)*ad_dc(i,1)
2030 ad_bc(i,1)=-v(i,j,1,nnew)*ad_dc(i,1)
2031 ad_v(i,j,1,nnew)=ad_dc(i,1)
2032 ad_dc(i,1)=0.0_r8
2033 END DO
2034 DO k=2,n(ng)-1
2035 DO i=istr,iend
2036!^ DC(i,k)=tl_v(i,j,k,nnew)- &
2037!^ & (tl_FC(i,k-1)*v(i,j,k-1,nnew)+ &
2038!^ & tl_BC(i,k )*v(i,j,k ,nnew)+ &
2039!^ & tl_FC(i,k )*v(i,j,k+1,nnew))
2040!^
2041 ad_fc(i,k-1)=ad_fc(i,k-1)-v(i,j,k-1,nnew)*ad_dc(i,k)
2042 ad_fc(i,k )=ad_fc(i,k )-v(i,j,k+1,nnew)*ad_dc(i,k)
2043 ad_bc(i,k )=ad_bc(i,k )-v(i,j,k ,nnew)*ad_dc(i,k)
2044 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)+ad_dc(i,k)
2045 ad_dc(i,k)=0.0_r8
2046 END DO
2047 END DO
2048 DO k=1,n(ng)
2049 DO i=istr,iend
2050!^ tl_BC(i,k)=tl_Hzk(i,k)-tl_FC(i,k)-tl_FC(i,k-1)
2051!^
2052 ad_hzk(i,k)=ad_hzk(i,k)+ad_bc(i,k)
2053 ad_fc(i,k-1)=ad_fc(i,k-1)-ad_bc(i,k)
2054 ad_fc(i,k )=ad_fc(i,k )-ad_bc(i,k)
2055 ad_bc(i,k)=0.0_r8
2056 END DO
2057 END DO
2058!
2059! Compute adjoint off-diagonal coefficients [lambda*dt*Akv/Hz] for
2060! the implicit vertical viscosity term at horizontal V-points and
2061! vertical W-points.
2062!
2063 DO i=istr,iend
2064!^ tl_FC(i,N(ng))=0.0_r8
2065!^
2066 ad_fc(i,n(ng))=0.0_r8
2067!^ tl_FC(i,0)=0.0_r8
2068!^
2069 ad_fc(i,0)=0.0_r8
2070 END DO
2071 cff=-lambda*dt(ng)/0.5_r8
2072 DO k=1,n(ng)-1
2073 DO i=istr,iend
2074 cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
2075 & z_r(i,j,k )-z_r(i,j-1,k ))
2076!^ tl_FC(i,k)=cff*(tl_cff1*AK(i,k)+cff1*tl_AK(i,k))
2077!^
2078 adfac=cff*ad_fc(i,k)
2079 ad_ak(i,k)=ad_ak(i,k)+cff1*adfac
2080 ad_cff1=ad_cff1+ak(i,k)*adfac
2081 ad_fc(i,k)=0.0_r8
2082!^ tl_cff1=-cff1*cff1*(tl_z_r(i,j,k+1)+tl_z_r(i,j-1,k+1)- &
2083!^ & tl_z_r(i,j,k )-tl_z_r(i,j-1,k ))
2084!^
2085 adfac=-cff1*cff1*ad_cff1
2086 ad_z_r(i,j-1,k )=ad_z_r(i,j-1,k )-adfac
2087 ad_z_r(i,j ,k )=ad_z_r(i,j ,k )-adfac
2088 ad_z_r(i,j-1,k+1)=ad_z_r(i,j-1,k+1)+adfac
2089 ad_z_r(i,j ,k+1)=ad_z_r(i,j ,k+1)+adfac
2090 ad_cff1=0.0_r8
2091 END DO
2092 END DO
2093# endif
2094!
2095! Time step adjoint right-hand-side terms.
2096!
2097 IF (iic(ng).eq.ntfirst(ng)) THEN
2098 cff=0.25_r8*dt(ng)
2099 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
2100 cff=0.25_r8*dt(ng)*3.0_r8/2.0_r8
2101 ELSE
2102 cff=0.25_r8*dt(ng)*23.0_r8/12.0_r8
2103 END IF
2104 DO i=istr,iend
2105 dc(i,0)=cff*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
2106 END DO
2107!
2108! The BASIC STATE "v" used below must be in transport units, but "v"
2109! retrived is in m/s so we multiply by "Hzk".
2110!
2111 DO k=1,n(ng)
2112 DO i=istr,iend
2113# ifdef DIAGNOSTICS_UV
2114!! DiaV3wrk(i,j,k,M3rate)=DiaV3wrk(i,j,k,M3rate)*oHz(i,k)
2115# ifdef BODYFORCE
2116!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+ &
2117!! & DC(i,0)*DiaRV(i,j,k,nrhs,M3vvis)* &
2118!! & oHz(i,k)
2119# endif
2120!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)*oHz(i,k)
2121# if defined UV_VIS2 || defined UV_VIS4
2122!! DiaV3wrk(i,j,k,M3hvis)=DiaV3wrk(i,j,k,M3hvis)*oHz(i,k)
2123!! DiaV3wrk(i,j,k,M3yvis)=DiaV3wrk(i,j,k,M3yvis)*oHz(i,k)
2124!! DiaV3wrk(i,j,k,M3xvis)=DiaV3wrk(i,j,k,M3xvis)*oHz(i,k)
2125# endif
2126!! DO idiag=1,M3pgrd
2127!! DiaV3wrk(i,j,k,idiag)=(DiaV3wrk(i,j,k,idiag)+ &
2128!! & DC(i,0)* &
2129!! & DiaRV(i,j,k,nrhs,idiag))* &
2130!! & oHz(i,k)
2131!! END DO
2132# endif
2133# ifdef SPLINES_VVISC
2134!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*oHz(i,k)+ &
2135!^ & (v(i,j,k,nnew)*Hzk(i,k))*tl_oHz(i,k)
2136!^
2137 ad_ohz(i,k)=ad_ohz(i,k)+ &
2138 & (v(i,j,k,nnew)*hzk(i,k))*ad_v(i,j,k,nnew)
2139 ad_v(i,j,k,nnew)=ad_v(i,j,k,nnew)*ohz(i,k)
2140# endif
2141!^ tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)+ &
2142!^ & DC(i,0)*tl_rv(i,j,k,nrhs)
2143!^
2144 ad_rv(i,j,k,nrhs)=ad_rv(i,j,k,nrhs)+ &
2145 & dc(i,0)*ad_v(i,j,k,nnew)
2146 END DO
2147 END DO
2148 DO i=istr,iend
2149 DO k=1,n(ng)
2150# ifdef SPLINES_VVISC
2151 ohz(i,k)=1.0_r8/hzk(i,k)
2152!^ tl_oHz(i,k)=-oHz(i,k)*oHz(i,k)*tl_Hzk(i,k)
2153!^
2154 ad_hzk(i,k)=ad_hzk(i,k)-ohz(i,k)*ohz(i,k)*ad_ohz(i,k)
2155 ad_ohz(i,k)=0.0_r8
2156# endif
2157!^ tl_Hzk(i,k)=0.5_r8*(tl_Hz(i,j-1,k)+ &
2158!^ & tl_Hz(i,j ,k))
2159!^
2160 adfac=0.5_r8*ad_hzk(i,k)
2161 ad_hz(i,j-1,k)=ad_hz(i,j-1,k)+adfac
2162 ad_hz(i,j ,k)=ad_hz(i,j ,k)+adfac
2163 ad_hzk(i,k)=0.0_r8
2164!^ tl_AK(i,k)=0.5_r8*(tl_Akv(i,j-1,k)+ &
2165!^ & tl_Akv(i,j ,k))
2166!^
2167 adfac=0.5_r8*ad_ak(i,k)
2168 ad_akv(i,j-1,k)=ad_akv(i,j-1,k)+adfac
2169 ad_akv(i,j ,k)=ad_akv(i,j ,k)+adfac
2170 ad_ak(i,k)=0.0_r8
2171 END DO
2172!^ tl_AK(i,0)=0.5_r8*(tl_Akv(i,j-1,0)+ &
2173!^ & tl_Akv(i,j ,0))
2174!^
2175 adfac=0.5_r8*ad_ak(i,0)
2176 ad_akv(i,j-1,0)=ad_akv(i,j-1,0)+adfac
2177 ad_akv(i,j ,0)=ad_akv(i,j ,0)+adfac
2178 ad_ak(i,0)=0.0_r8
2179 END DO
2180 END IF
2181!
2182!-----------------------------------------------------------------------
2183! Time step momentum equation in the XI-direction.
2184!-----------------------------------------------------------------------
2185!
2186! Compute BASIC STATE quatities.
2187!
2188 DO i=istru,iend
2189 ak(i,0)=0.5_r8*(akv(i-1,j,0)+ &
2190 & akv(i ,j,0))
2191 DO k=1,n(ng)
2192 ak(i,k)=0.5_r8*(akv(i-1,j,k)+ &
2193 & akv(i ,j,k))
2194 hzk(i,k)=0.5_r8*(hz(i-1,j,k)+ &
2195 & hz(i ,j,k))
2196# ifdef SPLINES_VVISC
2197 ohz(i,k)=1.0_r8/hzk(i,k)
2198# endif
2199 END DO
2200 END DO
2201!
2202! Couple and update new adjoint solution.
2203!
2204# if defined DIAGNOSTICS_UV && defined MASKING
2205!! DO k=1,N(ng)
2206!! DO i=IstrU,Iend
2207!! DO idiag=1,NDM3d
2208!! DiaU3wrk(i,j,k,idiag)=DiaU3wrk(i,j,k,idiag)*umask(i,j)
2209!! END DO
2210!! END DO
2211!! END DO
2212# endif
2213
2214 DO k=1,n(ng)
2215 DO i=istru,iend
2216# ifdef DIAGNOSTICS_UV
2217# ifdef NEARSHORE_MELLOR
2218!! DiaU3wrk(i,j,k,M3hrad)=DiaU3wrk(i,j,k,M3hrad)- &
2219!! & Dwrk(i,M2hrad)
2220# endif
2221# ifdef UV_ADV
2222!! DiaU3wrk(i,j,k,M3hadv)=DiaU3wrk(i,j,k,M3hadv)- &
2223!! & Dwrk(i,M2hadv)
2224!! DiaU3wrk(i,j,k,M3yadv)=DiaU3wrk(i,j,k,M3yadv)- &
2225!! & Dwrk(i,M2yadv)
2226!! DiaU3wrk(i,j,k,M3xadv)=DiaU3wrk(i,j,k,M3xadv)- &
2227!! & Dwrk(i,M2xadv)
2228# endif
2229# if defined UV_VIS2 || defined UV_VIS4
2230!! DiaU3wrk(i,j,k,M3hvis)=DiaU3wrk(i,j,k,M3hvis)- &
2231!! & Dwrk(i,M2hvis)
2232!! DiaU3wrk(i,j,k,M3yvis)=DiaU3wrk(i,j,k,M3yvis)- &
2233!! & Dwrk(i,M2yvis)
2234!! DiaU3wrk(i,j,k,M3xvis)=DiaU3wrk(i,j,k,M3xvis)- &
2235!! & Dwrk(i,M2xvis)
2236# endif
2237# ifdef UV_COR
2238!! DiaU3wrk(i,j,k,M3fcor)=DiaU3wrk(i,j,k,M3fcor)- &
2239!! & Dwrk(i,M2fcor)
2240# endif
2241!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)- &
2242!! & Dwrk(i,M2bstr)
2243!! DiaU3wrk(i,j,k,M3pgrd)=DiaU3wrk(i,j,k,M3pgrd)- &
2244!! & Dwrk(i,M2pgrd)
2245# endif
2246# ifdef NEARSHORE_MELLOR
2247# ifdef WET_DRY_NOT_YET
2248!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*umask_wet(i,j)
2249!^
2250 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)*umask_wet(i,j)
2251# endif
2252# ifdef MASKING
2253!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*umask(i,j)
2254!^
2255 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)*umask(i,j)
2256# endif
2257!^ tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)-tl_DCs(i,0)
2258!^
2259 ad_dcs(i,0)=ad_dcs(i,0)-ad_u_stokes(i,j,k)
2260# endif
2261# ifdef WET_DRY_NOT_YET
2262!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*umask_wet(i,j)
2263!^
2264 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)*umask_wet(i,j)
2265# endif
2266# ifdef MASKING
2267!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*umask(i,j)
2268!^
2269 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)*umask(i,j)
2270# endif
2271!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)-tl_DC(i,0)
2272!^
2273 ad_dc(i,0)=ad_dc(i,0)-ad_u(i,j,k,nnew)
2274 END DO
2275 END DO
2276!
2277! Replace INTERIOR POINTS incorrect vertical mean with more accurate
2278! barotropic component, ubar=DU_avg1/(D*on_u). Recall that, D=CF(:,0).
2279!
2280 DO i=istru,iend
2281 cf(i,0)=hzk(i,1)
2282 dc(i,0)=u(i,j,1,nnew)*hzk(i,1)
2283 END DO
2284 DO k=2,n(ng)
2285 DO i=istru,iend
2286 cf(i,0)=cf(i,0)+hzk(i,k)
2287 dc(i,0)=dc(i,0)+u(i,j,k,nnew)*hzk(i,k)
2288 END DO
2289 END DO
2290 DO i=istru,iend
2291 dc1(i,0)=dc(i,0) ! intermediate
2292 cff1=1.0/(cf(i,0)*on_u(i,j))
2293 dc(i,0)=(dc(i,0)*on_u(i,j)-du_avg1(i,j))*cff1 ! recursive
2294# ifdef NEARSHORE_MELLOR
2295 dcs1(i)=dcs(i,0) ! intermediate
2296 cff2=1.0_r8/cf(i,0)
2297 dcs(i,0)=dcs(i,0)*cff2-ubar_stokes(i,j) ! recursive
2298# endif
2299# ifdef DIAGNOSTICS_UV
2300!! Dwrk(i,M2bstr)=(Dwrk(i,M2bstr)*on_u(i,j)- &
2301!! & DiaU2wrk(i,j,M2bstr)-DiaU2wrk(i,j,M2sstr))* &
2302!! & cff1
2303!! DO idiag=1,M2pgrd
2304!! Dwrk(i,idiag)=(Dwrk(i,idiag)*on_u(i,j)- &
2305!! & DiaU2wrk(i,j,idiag))*cff1
2306!! END DO
2307# endif
2308# ifdef NEARSHORE_MELLOR
2309!^ tl_DCs(i,0)=tl_DCs(i,0)*cff2+ &
2310!^ & DCs1(i)*tl_cff2-tl_ubar_stokes(i,j)
2311!^
2312 ad_cff2=ad_cff2+dcs1(i)*ad_dcs(i,0)
2313 ad_ubar_stokes(i,j)=ad_ubar_stokes(i,j)-ad_dcs(i,0)
2314 ad_dcs(i,0)=cff2*ad_dcs(i,0)
2315!^ tl_cff2=-cff2*cff2*tl_CF(i,0)
2316!^
2317 ad_cf(i,0)=ad_cf(i,0)-cff2*cff2*ad_cff2
2318 ad_cff2=0.0_r8
2319# endif
2320!^ tl_DC(i,0)=(tl_DC(i,0)*on_u(i,j)-tl_DU_avg1(i,j))*cff1+ &
2321!^ & (DC1(i,0)*on_u(i,j)-DU_avg1(i,j))*tl_cff1
2322!^
2323 adfac=cff1*ad_dc(i,0)
2324 ad_du_avg1(i,j)=ad_du_avg1(i,j)-adfac
2325 ad_cff1=ad_cff1+ &
2326 & (dc1(i,0)*on_u(i,j)-du_avg1(i,j))*ad_dc(i,0)
2327 ad_dc(i,0)=on_u(i,j)*adfac
2328!^ tl_cff1=-cff1*cff1*tl_CF(i,0)*on_u(i,j)
2329!^
2330 ad_cf(i,0)=ad_cf(i,0)-cff1*cff1*on_u(i,j)*ad_cff1
2331 ad_cff1=0.0_r8
2332 END DO
2333 DO i=istru,iend
2334 cf(i,0)=hzk(i,1)
2335 dc(i,0)=u(i,j,1,nnew)*hzk(i,1)
2336 END DO
2337 DO k=2,n(ng)
2338 DO i=istru,iend
2339 cf(i,0)=cf(i,0)+hzk(i,k)
2340 dc(i,0)=dc(i,0)+u(i,j,k,nnew)*hzk(i,k)
2341# ifdef NEARSHORE_MELLOR
2342 dcs(i,0)=dcs(i,0)+u_stokes(i,j,k)*hzk(i,k)
2343# endif
2344# ifdef DIAGNOSTICS_UV
2345# ifdef NEARSHORE_MELLOR
2346!! Dwrk(i,M2hrad)=Dwrk(i,M2hrad)+ &
2347!! & DiaU3wrk(i,j,k,M3hrad)*Hzk(i,k)
2348# endif
2349# ifdef UV_ADV
2350!! Dwrk(i,M2hadv)=Dwrk(i,M2hadv)+ &
2351!! & DiaU3wrk(i,j,k,M3hadv)*Hzk(i,k)
2352!! Dwrk(i,M2yadv)=Dwrk(i,M2yadv)+ &
2353!! & DiaU3wrk(i,j,k,M3yadv)*Hzk(i,k)
2354!! Dwrk(i,M2xadv)=Dwrk(i,M2xadv)+ &
2355!! & DiaU3wrk(i,j,k,M3xadv)*Hzk(i,k)
2356# endif
2357# if defined UV_VIS2 || defined UV_VIS4
2358!! Dwrk(i,M2hvis)=Dwrk(i,M2hvis)+ &
2359!! & DiaU3wrk(i,j,k,M3hvis)*Hzk(i,k)
2360!! Dwrk(i,M2yvis)=Dwrk(i,M2yvis)+ &
2361!! & DiaU3wrk(i,j,k,M3yvis)*Hzk(i,k)
2362!! Dwrk(i,M2xvis)=Dwrk(i,M2xvis)+ &
2363!! & DiaU3wrk(i,j,k,M3xvis)*Hzk(i,k)
2364# endif
2365# ifdef UV_COR
2366!! Dwrk(i,M2fcor)=Dwrk(i,M2fcor)+ &
2367!! & DiaU3wrk(i,j,k,M3fcor)*Hzk(i,k)
2368# endif
2369!! Dwrk(i,M2bstr)=Dwrk(i,M2bstr)+ &
2370!! & DiaU3wrk(i,j,k,M3vvis)*Hzk(i,k)
2371!! Dwrk(i,M2pgrd)=Dwrk(i,M2pgrd)+ &
2372!! & DiaU3wrk(i,j,k,M3pgrd)*Hzk(i,k)
2373# endif
2374# ifdef NEARSHORE_MELLOR
2375!^ tl_DCs(i,0)=tl_DCs(i,0)+ &
2376!^ & tl_u_stokes(i,j,k)*Hzk(i,k)+ &
2377!^ & u_stokes(i,j,k)*tl_Hzk(i,k)
2378!^
2379 ad_u_stokes(i,j,k)=ad_u_stokes(i,j,k)+hzk(i,k)*ad_dcs(i,0)
2380 ad_hzk(i,k)=ad_hzk(i,k)+u_stokes(i,j,k)*ad_dcs(i,0)
2381# endif
2382!^ tl_DC(i,0)=tl_DC(i,0)+ &
2383!^ & tl_u(i,j,k,nnew)*Hzk(i,k)+ &
2384!^ & u(i,j,k,nnew)*tl_Hzk(i,k)
2385!^
2386 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)+ad_dc(i,0)*hzk(i,k)
2387 ad_hzk(i,k)=ad_hzk(i,k)+u(i,j,k,nnew)*ad_dc(i,0)
2388!^ tl_CF(i,0)=tl_CF(i,0)+tl_Hzk(i,k)
2389!^
2390 ad_hzk(i,k)=ad_hzk(i,k)+ad_cf(i,0)
2391 END DO
2392 END DO
2393 DO i=istru,iend
2394 cf(i,0)=hzk(i,1)
2395 dc(i,0)=u(i,j,1,nnew)*hzk(i,1)
2396# ifdef NEARSHORE_MELLOR
2397 dcs(i,0)=u_stokes(i,j,1)*hzk(i,1)
2398# endif
2399# ifdef DIAGNOSTICS_UV
2400# ifdef NEARSHORE_MELLOR
2401!! Dwrk(i,M2hrad)=DiaU3wrk(i,j,1,M3hrad)*Hzk(i,1)
2402# endif
2403# ifdef UV_ADV
2404!! Dwrk(i,M2hadv)=DiaU3wrk(i,j,1,M3hadv)*Hzk(i,1)
2405!! Dwrk(i,M2yadv)=DiaU3wrk(i,j,1,M3yadv)*Hzk(i,1)
2406!! Dwrk(i,M2xadv)=DiaU3wrk(i,j,1,M3xadv)*Hzk(i,1)
2407# endif
2408# if defined UV_VIS2 || defined UV_VIS4
2409!! Dwrk(i,M2hvis)=DiaU3wrk(i,j,1,M3hvis)*Hzk(i,1)
2410!! Dwrk(i,M2yvis)=DiaU3wrk(i,j,1,M3yvis)*Hzk(i,1)
2411!! Dwrk(i,M2xvis)=DiaU3wrk(i,j,1,M3xvis)*Hzk(i,1)
2412# endif
2413# ifdef UV_COR
2414!! Dwrk(i,M2fcor)=DiaU3wrk(i,j,1,M3fcor)*Hzk(i,1)
2415# endif
2416!! Dwrk(i,M2bstr)=DiaU3wrk(i,j,1,M3vvis)*Hzk(i,1)
2417!! Dwrk(i,M2pgrd)=DiaU3wrk(i,j,1,M3pgrd)*Hzk(i,1)
2418# endif
2419# ifdef NEARSHORE_MELLOR
2420!^ tl_DCs(i,0)=tl_u_stokes(i,j,1)*Hzk(i,1)+ &
2421!^ & u_stokes(i,j,1)*tl_Hzk(i,1)
2422!^
2423 ad_u_stokes(i,j,1)=ad_u_stokes(i,j,1)+hzk(i,1)*ad_dcs(i,0)
2424 ad_hzk(i,1)=ad_hzk(i,1)+u_stokes(i,j,1)*ad_dcs(i,0)
2425 ad_dcs(i,0)=0.0_r8
2426# endif
2427!^ tl_DC(i,0)=tl_u(i,j,1,nnew)*Hzk(i,1)+ &
2428!^ & u(i,j,1,nnew)*tl_Hzk(i,1)
2429!^
2430 ad_u(i,j,1,nnew)=ad_u(i,j,1,nnew)+ad_dc(i,0)*hzk(i,1)
2431 ad_hzk(i,1)=ad_hzk(i,1)+u(i,j,1,nnew)*ad_dc(i,0)
2432 ad_dc(i,0)=0.0_r8
2433!^ tl_CF(i,0)=tl_Hzk(i,1)
2434!^
2435 ad_hzk(i,1)=ad_hzk(i,1)+ad_cf(i,0)
2436 ad_cf(i,0)=0.0_r8
2437 END DO
2438
2439# if defined SPLINES_VVISC
2440!
2441! Apply spline algorithm to BASIC STATE "u" which should be
2442! in units of velocity (m/s). DC will be used in the tangent
2443! linear spline algorithm below. Solve BASIC STATE tridiagonal
2444! system.
2445!
2446 cff1=1.0_r8/6.0_r8
2447 DO k=1,n(ng)-1
2448 DO i=istru,iend
2449 fc(i,k)=cff1*hzk(i,k )-dt(ng)*ak(i,k-1)*ohz(i,k )
2450 cf(i,k)=cff1*hzk(i,k+1)-dt(ng)*ak(i,k+1)*ohz(i,k+1)
2451 END DO
2452 END DO
2453 DO i=istru,iend
2454 cf(i,0)=0.0_r8
2455 dc(i,0)=0.0_r8
2456 END DO
2457!
2458! LU decomposition and forward substitution.
2459!
2460 cff1=1.0_r8/3.0_r8
2461 DO k=1,n(ng)-1
2462 DO i=istru,iend
2463 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
2464 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
2465 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
2466 cf(i,k)=cff*cf(i,k)
2467 dc(i,k)=cff*(u(i,j,k+1,nnew)-u(i,j,k,nnew)- &
2468 & fc(i,k)*dc(i,k-1))
2469 END DO
2470 END DO
2471!
2472! Backward substitution. Save DC for the adjoint below.
2473! DC is scaled later by AK.
2474!
2475 DO i=istru,iend
2476 dc(i,n(ng))=0.0_r8
2477 END DO
2478 DO k=n(ng)-1,1,-1
2479 DO i=istru,iend
2480 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
2481 END DO
2482 END DO
2483!
2484! Multiply DC (BASIC STATE spline gradients) by AK and save it in DC1.
2485!
2486 DO k=0,n(ng)
2487 DO i=istru,iend
2488 dc1(i,k)=dc(i,k)*ak(i,k)
2489 END DO
2490 END DO
2491!^ DO k=1,N(ng)
2492!^
2493 DO k=n(ng),1,-1
2494 DO i=istru,iend
2495# ifdef DIAGNOSTICS_UV
2496!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+cff
2497# endif
2498!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)+tl_cff
2499!^
2500 ad_cff=ad_cff+ad_u(i,j,k,nnew)
2501!^ tl_cff=dt(ng)*(tl_oHz(i,k)*(DC(i,k)-DC(i,k-1))+ &
2502!^ & oHz(i,k)*(tl_DC(i,k)-tl_DC(i,k-1)))
2503!^ use DC1
2504 adfac=dt(ng)*ad_cff
2505 adfac1=adfac*ohz(i,k)
2506 ad_dc(i,k-1)=ad_dc(i,k-1)-adfac1
2507 ad_dc(i,k )=ad_dc(i,k )+adfac1
2508 ad_ohz(i,k)=ad_ohz(i,k)+(dc1(i,k)-dc1(i,k-1))*adfac
2509 ad_cff=0.0_r8
2510!^ tl_DC(i,k)=tl_DC(i,k)*AK(i,k)+DC(i,k)*tl_AK(i,k)
2511!^ use DC
2512 ad_ak(i,k)=ad_ak(i,k)+dc(i,k)*ad_dc(i,k)
2513 ad_dc(i,k)=ak(i,k)*ad_dc(i,k)
2514 END DO
2515 END DO
2516!
2517! Adjoint back substitution.
2518!
2519 DO k=1,n(ng)-1
2520 DO i=istru,iend
2521!^ tl_DC(i,k)=tl_DC(i,k)-CF(i,k)*tl_DC(i,k+1)
2522!^
2523 ad_dc(i,k+1)=ad_dc(i,k+1)-cf(i,k)*ad_dc(i,k)
2524 END DO
2525 END DO
2526 DO i=istru,iend
2527!^ tl_DC(i,N(ng))=0.0_r8
2528!^
2529 ad_dc(i,n(ng))=0.0_r8
2530 END DO
2531!
2532! Adjoint LU decomposition and forward substitution.
2533!
2534 cff1=1.0_r8/3.0_r8
2535 DO k=n(ng)-1,1,-1
2536 DO i=istru,iend
2537 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
2538 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
2539 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
2540!^ tl_DC(i,k)=cff*(tl_u(i,j,k+1,nnew)- &
2541!^ & tl_u(i,j,k ,nnew)- &
2542!^ & (tl_FC(i,k)*DC(i,k-1)+ &
2543!^ & tl_BC(i,k)*DC(i,k )+ &
2544!^ & tl_CF(i,k)*DC(i,k+1))- &
2545!^ & FC(i,k)*tl_DC(i,k-1))
2546!^
2547 adfac=cff*ad_dc(i,k)
2548 ad_dc(i,k-1)=ad_dc(i,k-1)-fc(i,k)*adfac
2549 ad_cf(i,k)=ad_cf(i,k)-dc(i,k+1)*adfac
2550 ad_bc(i,k)=ad_bc(i,k)-dc(i,k )*adfac
2551 ad_fc(i,k)=ad_fc(i,k)-dc(i,k-1)*adfac
2552 ad_u(i,j,k ,nnew)=ad_u(i,j,k ,nnew)-adfac
2553 ad_u(i,j,k+1,nnew)=ad_u(i,j,k+1,nnew)+adfac
2554 ad_dc(i,k)=0.0_r8
2555!^ tl_BC(i,k)=cff1*(tl_Hzk(i,k)+tl_Hzk(i,k+1))+ &
2556!^ & dt(ng)*(tl_AK(i,k)*(oHz(i,k)+oHz(i,k+1))+ &
2557!^ & AK(i,k)*(tl_oHz(i,k)+tl_oHz(i,k+1)))
2558!^
2559 adfac=dt(ng)*ad_bc(i,k)
2560 adfac1=adfac*ak(i,k)
2561 adfac2=cff1*ad_bc(i,k)
2562 ad_ohz(i,k )=ad_ohz(i,k )+adfac1
2563 ad_ohz(i,k+1)=ad_ohz(i,k+1)+adfac1
2564 ad_ak(i,k)=ad_ak(i,k)+(ohz(i,k)+ohz(i,k+1))*adfac
2565 ad_hzk(i,k )=ad_hzk(i,k )+adfac2
2566 ad_hzk(i,k+1)=ad_hzk(i,k+1)+adfac2
2567 ad_bc(i,k)=0.0_r8
2568 END DO
2569 END DO
2570!
2571! Use conservative, parabolic spline reconstruction of adjoint
2572! vertical viscosity derivatives. Then, time step vertical
2573! viscosity term implicitly by solving a tridiagonal system.
2574!
2575 DO i=istru,iend
2576!^ tl_DC(i,0)=0.0_r8
2577!^
2578 ad_dc(i,0)=0.0_r8
2579!^ tl_CF(i,0)=0.0_r8
2580!^
2581 ad_cf(i,0)=0.0_r8
2582 END DO
2583 cff1=1.0_r8/6.0_r8
2584 DO k=1,n(ng)-1
2585 DO i=istru,iend
2586!^ tl_CF(i,k)=cff1*tl_Hzk(i,k+1)- &
2587!^ & dt(ng)*(tl_AK(i,k+1)*oHz(i,k+1)+ &
2588!^ & AK(i,k+1)*tl_oHz(i,k+1))
2589!^
2590 adfac=dt(ng)*ad_cf(i,k)
2591 ad_ohz(i,k+1)=ad_ohz(i,k+1)-ak(i,k+1)*adfac
2592 ad_ak(i,k+1)=ad_ak(i,k+1)-ohz(i,k+1)*adfac
2593 ad_hzk(i,k+1)=ad_hzk(i,k+1)+cff1*ad_cf(i,k)
2594 ad_cf(i,k)=0.0_r8
2595!^ tl_FC(i,k)=cff1*tl_Hzk(i,k )- &
2596!^ & dt(ng)*(tl_AK(i,k-1)*oHz(i,k )+ &
2597!^ & AK(i,k-1)*tl_oHz(i,k ))
2598!^
2599 adfac=dt(ng)*ad_fc(i,k)
2600 ad_ohz(i,k)=ad_ohz(i,k)-ak(i,k-1)*adfac
2601 ad_ak(i,k-1)=ad_ak(i,k-1)-ohz(i,k)*adfac
2602 ad_hzk(i,k)=ad_hzk(i,k)+cff1*ad_fc(i,k)
2603 ad_fc(i,k)=0.0_r8
2604 END DO
2605 END DO
2606# else
2607!
2608! Compute BASIC STATE off-diagonal coefficients [lambda*dt*Akv/Hz]
2609! for the implicit vertical viscosity term at horizontal U-points
2610! and vertical W-points.
2611!
2612! NOTE: The original code solves the tridiagonal system A*u=r where
2613! A is a matrix and u and r are vectors. We need to solve the
2614! tangent linear of this system which is A*tl_u+tl_A*u=tl_r.
2615! Here, tl_A*u and tl_r are known, so we must solve for tl_u
2616! by inverting A*tl_u=tl_r-tl_A*u.
2617!
2618 cff=-lambda*dt(ng)/0.5_r8
2619 DO k=1,n(ng)-1
2620 DO i=istru,iend
2621 cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
2622 & z_r(i,j,k )-z_r(i-1,j,k ))
2623 fc(i,k)=cff*cff1*ak(i,k)
2624 END DO
2625 END DO
2626 DO i=istru,iend
2627 fc(i,0)=0.0_r8
2628 fc(i,n(ng))=0.0_r8
2629 END DO
2630 DO k=1,n(ng)
2631 DO i=istru,iend
2632 bc(i,k)=hzk(i,k)-fc(i,k)-fc(i,k-1)
2633 END DO
2634 END DO
2635 DO i=istru,iend
2636 cff=1.0_r8/bc(i,1)
2637 cf(i,1)=cff*fc(i,1)
2638 END DO
2639 DO k=2,n(ng)-1
2640 DO i=istru,iend
2641 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
2642 cf(i,k)=cff*fc(i,k)
2643 END DO
2644 END DO
2645!
2646! Compute new adjoint solution by back substitution.
2647! (DC is a tangent linear variable here).
2648!
2649!^ DO k=N(ng)-1,1,-1
2650!^
2651 DO k=1,n(ng)-1
2652 DO i=istru,iend
2653# ifdef DIAGNOSTICS_UV
2654!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+ &
2655!! & u(i,j,k,nnew)-wrk(i,k)
2656# endif
2657!^ tl_u(i,j,k,nnew)=DC(i,k)
2658!^
2659 ad_dc(i,k)=ad_dc(i,k)+ad_u(i,j,k,nnew)
2660 ad_u(i,j,k,nnew)=0.0_r8
2661!^ DC(i,k)=DC(i,k)-CF(i,k)*DC(i,k+1)
2662!^
2663 ad_dc(i,k+1)=-cf(i,k)*ad_dc(i,k)
2664# ifdef DIAGNOSTICS_UV
2665!! wrk(i,k)=u(i,j,k,nnew)*oHz(i,k)
2666# endif
2667 END DO
2668 END DO
2669 DO i=istru,iend
2670# ifdef DIAGNOSTICS_UV
2671!! DiaU3wrk(i,j,N(ng),M3vvis)=DiaU3wrk(i,j,N(ng),M3vvis)+ &
2672!! & u(i,j,N(ng),nnew)-wrk(i,N(ng))
2673# endif
2674!^ tl_u(i,j,N(ng),nnew)=DC(i,N(ng))
2675!^
2676 ad_dc(i,n(ng))=ad_dc(i,n(ng))+ad_u(i,j,n(ng),nnew)
2677 ad_u(i,j,n(ng),nnew)=0.0_r8
2678!^ DC(i,N(ng))=(DC(i,N(ng))-FC(i,N(ng)-1)*DC(i,N(ng)-1))/ &
2679!^ & (BC(i,N(ng))-FC(i,N(ng)-1)*CF(i,N(ng)-1))
2680!^
2681 adfac=ad_dc(i,n(ng))/ &
2682 & (bc(i,n(ng))-fc(i,n(ng)-1)*cf(i,n(ng)-1))
2683 ad_dc(i,n(ng)-1)=ad_dc(i,n(ng)-1)-fc(i,n(ng)-1)*adfac
2684 ad_dc(i,n(ng) )=adfac
2685 END DO
2686!
2687! Solve the adjoint tridiagonal system.
2688! (DC is a tangent linear variable here).
2689!
2690!^ DO k=2,N(ng)-1
2691!^
2692 DO k=n(ng)-1,2,-1
2693 DO i=istru,iend
2694 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
2695!^ DC(i,k)=cff*(DC(i,k)-FC(i,k-1)*DC(i,k-1))
2696!^
2697 adfac=cff*ad_dc(i,k)
2698 ad_dc(i,k-1)=ad_dc(i,k-1)-fc(i,k-1)*adfac
2699 ad_dc(i,k )=adfac
2700 END DO
2701 END DO
2702 DO i=istru,iend
2703 cff=1.0_r8/bc(i,1)
2704!^ DC(i,1)=cff*DC(i,1)
2705!^
2706 ad_dc(i,1)=cff*ad_dc(i,1)
2707 END DO
2708 DO i=istru,iend
2709!^ DC(i,N(ng))=tl_u(i,j,N(ng),nnew)- &
2710!^ & (tl_FC(i,N(ng)-1)*u(i,j,N(ng)-1,nnew)+ &
2711!^ & tl_BC(i,N(ng) )*u(i,j,N(ng) ,nnew))
2712!^
2713 ad_bc(i,n(ng) )=-u(i,j,n(ng) ,nnew)*ad_dc(i,n(ng))
2714 ad_fc(i,n(ng)-1)=-u(i,j,n(ng)-1,nnew)*ad_dc(i,n(ng))
2715 ad_u(i,j,n(ng),nnew)=ad_dc(i,n(ng))
2716 ad_dc(i,n(ng))=0.0_r8
2717!^ DC(i,1)=tl_u(i,j,1,nnew)- &
2718!^ & (tl_BC(i,1)*u(i,j,1,nnew)+ &
2719!^ & tl_FC(i,1)*u(i,j,2,nnew))
2720!^
2721 ad_fc(i,1)=-u(i,j,2,nnew)*ad_dc(i,1)
2722 ad_bc(i,1)=-u(i,j,1,nnew)*ad_dc(i,1)
2723 ad_u(i,j,1,nnew)=ad_dc(i,1)
2724 ad_dc(i,1)=0.0_r8
2725 END DO
2726 DO k=2,n(ng)-1
2727 DO i=istru,iend
2728!^ DC(i,k)=tl_u(i,j,k,nnew)- &
2729!^ & (tl_FC(i,k-1)*u(i,j,k-1,nnew)+ &
2730!^ & tl_BC(i,k )*u(i,j,k ,nnew)+ &
2731!^ & tl_FC(i,k )*u(i,j,k+1,nnew))
2732!^
2733 ad_fc(i,k-1)=ad_fc(i,k-1)-u(i,j,k-1,nnew)*ad_dc(i,k)
2734 ad_fc(i,k )=ad_fc(i,k )-u(i,j,k+1,nnew)*ad_dc(i,k)
2735 ad_bc(i,k )=ad_bc(i,k )-u(i,j,k ,nnew)*ad_dc(i,k)
2736 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)+ad_dc(i,k)
2737 ad_dc(i,k)=0.0_r8
2738 END DO
2739 END DO
2740 DO k=1,n(ng)
2741 DO i=istru,iend
2742!^ tl_BC(i,k)=tl_Hzk(i,k)-tl_FC(i,k)-tl_FC(i,k-1)
2743!^
2744 ad_hzk(i,k)=ad_hzk(i,k)+ad_bc(i,k)
2745 ad_fc(i,k-1)=ad_fc(i,k-1)-ad_bc(i,k)
2746 ad_fc(i,k )=ad_fc(i,k )-ad_bc(i,k)
2747 ad_bc(i,k)=0.0_r8
2748 END DO
2749 END DO
2750!
2751! Compute adjoint off-diagonal coefficients [lambda*dt*Akv/Hz] for
2752! the implicit vertical viscosity term at horizontal U-points and
2753! vertical W-points.
2754!
2755 DO i=istru,iend
2756!^ tl_FC(i,N(ng))=0.0_r8
2757!^
2758 ad_fc(i,n)=0.0_r8
2759!^ tl_FC(i,0)=0.0_r8
2760!^
2761 ad_fc(i,0)=0.0_r8
2762 END DO
2763 cff=-lambda*dt(ng)/0.5_r8
2764 DO k=1,n(ng)-1
2765 DO i=istru,iend
2766 cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
2767 & z_r(i,j,k )-z_r(i-1,j,k ))
2768!^ tl_FC(i,k)=cff*(tl_cff1*AK(i,k)+cff1*tl_AK(i,k))
2769!^
2770 adfac=cff*ad_fc(i,k)
2771 ad_ak(i,k)=ad_ak(i,k)+cff1*adfac
2772 ad_cff1=ad_cff1+ak(i,k)*adfac
2773 ad_fc(i,k)=0.0_r8
2774!^ tl_cff1=-cff1*cff1*(tl_z_r(i,j,k+1)+tl_z_r(i-1,j,k+1)- &
2775!^ & tl_z_r(i,j,k )-tl_z_r(i-1,j,k ))
2776!^
2777 adfac=-cff1*cff1*ad_cff1
2778 ad_z_r(i-1,j,k )=ad_z_r(i-1,j,k )-adfac
2779 ad_z_r(i ,j,k )=ad_z_r(i ,j,k )-adfac
2780 ad_z_r(i-1,j,k+1)=ad_z_r(i-1,j,k+1)+adfac
2781 ad_z_r(i ,j,k+1)=ad_z_r(i ,j,k+1)+adfac
2782 ad_cff1=0.0_r8
2783 END DO
2784 END DO
2785# endif
2786!
2787! Time step adjoint right-hand-side terms.
2788!
2789 IF (iic(ng).eq.ntfirst(ng)) THEN
2790 cff=0.25_r8*dt(ng)
2791 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
2792 cff=0.25_r8*dt(ng)*3.0_r8/2.0_r8
2793 ELSE
2794 cff=0.25_r8*dt(ng)*23.0_r8/12.0_r8
2795 END IF
2796 DO i=istru,iend
2797 dc(i,0)=cff*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
2798 END DO
2799!
2800! The BASIC STATE "u" used below must be in transport units, but "u"
2801! retrived is in m/s so we multiply by "Hzk".
2802!
2803 DO k=1,n(ng)
2804 DO i=istru,iend
2805# ifdef DIAGNOSTICS_UV
2806!! DiaU3wrk(i,j,k,M3rate)=DiaU3wrk(i,j,k,M3rate)*oHz(i,k)
2807# ifdef BODYFORCE
2808!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+ &
2809!! & DC(i,0)*DiaRU(i,j,k,nrhs,M3vvis)* &
2810!! & oHz(i,k)
2811# endif
2812!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)*oHz(i,k)
2813# if defined UV_VIS2 || defined UV_VIS4
2814!! DiaU3wrk(i,j,k,M3hvis)=DiaU3wrk(i,j,k,M3hvis)*oHz(i,k)
2815!! DiaU3wrk(i,j,k,M3yvis)=DiaU3wrk(i,j,k,M3yvis)*oHz(i,k)
2816!! DiaU3wrk(i,j,k,M3xvis)=DiaU3wrk(i,j,k,M3xvis)*oHz(i,k)
2817# endif
2818!! DO idiag=1,M3pgrd
2819!! DiaU3wrk(i,j,k,idiag)=(DiaU3wrk(i,j,k,idiag)+ &
2820!! & DC(i,0)*DiaRU(i,j,k,nrhs,idiag))* &
2821!! & oHz(i,k)
2822!! END DO
2823# endif
2824# ifdef SPLINES_VVISC
2825!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*oHz(i,k)+ &
2826!^ & (u(i,j,k,nnew)*Hzk(i,k))*tl_oHz(i,k)
2827!^
2828 ad_ohz(i,k)=ad_ohz(i,k)+ &
2829 & (u(i,j,k,nnew)*hzk(i,k))*ad_u(i,j,k,nnew)
2830 ad_u(i,j,k,nnew)=ad_u(i,j,k,nnew)*ohz(i,k)
2831# endif
2832!^ tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)+ &
2833!^ & DC(i,0)*tl_ru(i,j,k,nrhs)
2834!^
2835 ad_ru(i,j,k,nrhs)=ad_ru(i,j,k,nrhs)+ &
2836 & dc(i,0)*ad_u(i,j,k,nnew)
2837 END DO
2838 END DO
2839 DO i=istru,iend
2840 DO k=1,n(ng)
2841# if defined SPLINES_VVISC || defined DIAGNOSTICS_UV
2842 ohz(i,k)=1.0_r8/hzk(i,k)
2843!^ tl_oHz(i,k)=-oHz(i,k)*oHz(i,k)*tl_Hzk(i,k)
2844!^
2845 ad_hzk(i,k)=ad_hzk(i,k)-ohz(i,k)*ohz(i,k)*ad_ohz(i,k)
2846 ad_ohz(i,k)=0.0_r8
2847# endif
2848!^ tl_Hzk(i,k)=0.5_r8*(tl_Hz(i-1,j,k)+ &
2849!^ & tl_Hz(i ,j,k))
2850!^
2851 adfac=0.5_r8*ad_hzk(i,k)
2852 ad_hz(i-1,j,k)=ad_hz(i-1,j,k)+adfac
2853 ad_hz(i ,j,k)=ad_hz(i ,j,k)+adfac
2854 ad_hzk(i,k)=0.0_r8
2855!^ tl_AK(i,k)=0.5_r8*(tl_Akv(i-1,j,k)+ &
2856!^ & tl_Akv(i ,j,k))
2857!^
2858 adfac=0.5_r8*ad_ak(i,k)
2859 ad_akv(i-1,j,k)=ad_akv(i-1,j,k)+adfac
2860 ad_akv(i ,j,k)=ad_akv(i ,j,k)+adfac
2861 ad_ak(i,k)=0.0_r8
2862 END DO
2863!^ tl_AK(i,0)=0.5_r8*(tl_Akv(i-1,j,0)+ &
2864!^ & tl_Akv(i ,j,0))
2865!^
2866 adfac=0.5_r8*ad_ak(i,0)
2867 ad_akv(i-1,j,0)=ad_akv(i-1,j,0)+adfac
2868 ad_akv(i ,j,0)=ad_akv(i ,j,0)+adfac
2869 ad_ak(i,0)=0.0_r8
2870 END DO
2871 END DO j_loop2
2872
2873# ifdef AD_OUTPUT_STATE
2874!
2875!-----------------------------------------------------------------------
2876! Initialize full 3D momentum adjoint output arrays. Due to the exact
2877! discrete adjoint and the predictor/corrector time stepping scheme
2878! with multiple time level schemes, this routine contains the first
2879! piece of the adjoint solution at time level "nnew". The other
2880! piece is computed in "ad_pre_ste3d".
2881!-----------------------------------------------------------------------
2882!
2883 DO j=jstrb,jendb
2884 IF (j.ge.jstrm) THEN
2885 DO k=1,n(ng)
2886 DO i=istrb,iendb
2887 ad_v_sol(i,j,k)=ad_v(i,j,k,nnew)
2888 END DO
2889 END DO
2890 END IF
2891 DO k=1,n(ng)
2892 DO i=istrm,iendb
2893 ad_u_sol(i,j,k)=ad_u(i,j,k,nnew)
2894 END DO
2895 END DO
2896 END DO
2897# endif
2898!
2899 RETURN
2900 END SUBROUTINE ad_step3d_uv_tile
2901#endif
2902 END MODULE ad_step3d_uv_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_3d_mod
Definition ad_exchange_3d.F:2
ad_exchange_3d_mod::ad_exchange_v3d_tile
subroutine ad_exchange_v3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, ad_a)
Definition ad_exchange_3d.F:1283
ad_exchange_3d_mod::ad_exchange_u3d_tile
subroutine ad_exchange_u3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, ad_a)
Definition ad_exchange_3d.F:872
ad_step3d_uv_mod
Definition ad_step3d_uv.F:2
ad_step3d_uv_mod::ad_step3d_uv_tile
subroutine ad_step3d_uv_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, knew, nrhs, nstp, nnew, umask, vmask, umask_wet, vmask_wet, om_v, on_u, pm, pn, hz, ad_hz, z_r, ad_z_r, z_w, ad_z_w, akv, ad_akv, du_avg1, ad_du_avg1, dv_avg1, ad_dv_avg1, du_avg2, ad_du_avg2, dv_avg2, ad_dv_avg2, ad_ru, ad_rv, u, ad_u, v, ad_v, ad_ubar, ad_vbar, ad_ubar_sol, ad_vbar_sol, ubar_stokes, ad_ubar_stokes, vbar_stokes, ad_vbar_stokes, u_stokes, ad_u_stokes, v_stokes, ad_v_stokes, huon, ad_huon, hvom, ad_hvom)
Definition ad_step3d_uv.F:194
ad_step3d_uv_mod::ad_step3d_uv
subroutine, public ad_step3d_uv(ng, tile)
Definition ad_step3d_uv.F:57
ad_u3dbc_mod
Definition ad_u3dbc_im.F:2
ad_u3dbc_mod::ad_u3dbc_tile
subroutine, public ad_u3dbc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nstp, nout, ad_u)
Definition ad_u3dbc_im.F:57
ad_v3dbc_mod
Definition ad_v3dbc_im.F:2
ad_v3dbc_mod::ad_v3dbc_tile
subroutine, public ad_v3dbc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nstp, nout, ad_v)
Definition ad_v3dbc_im.F:57
exchange_3d_mod
Definition exchange_3d.F:2
exchange_3d_mod::exchange_u3d_tile
subroutine exchange_u3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:495
exchange_3d_mod::exchange_v3d_tile
subroutine exchange_v3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:707
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_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
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_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
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::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::lambda
real(r8) lambda
Definition mod_scalars.F:750
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::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::ntfirst
integer, dimension(:), allocatable ntfirst
Definition mod_scalars.F:377
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
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::knew
integer, dimension(:), allocatable knew
Definition mod_stepping.F:64
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
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_exchange3d
subroutine ad_mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, ad_a, ad_b, ad_c, ad_d)
Definition mp_exchange.F:5451
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
uv_var_change_mod
Definition uv_var_change.F:2
uv_var_change_mod::ad_uv_c2a_grid
subroutine, public ad_uv_c2a_grid(ng, tile, model, ninp)
Definition uv_var_change.F:715
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