ROMS
Loading...
Searching...
No Matches
tl_step3d_uv.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE tl_step3d_uv_mod
3#if defined TANGENT && defined SOLVE3D
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine time-steps the tangent linear horizontal momentum !
13! equations. The vertical viscosity terms are time-stepped using !
14! an implicit algorithm. !
15! !
16! BASIC STATE variables needed: u, v, ru, rv, Akv, !
17! Hz, Huon, Hvom, z_r, z_w, !
18! DU_avg1, DU_avg2, DV_avg1, DV_avg2, !
19! Qsrc !
20! !
21!=======================================================================
22!
23 USE mod_param
24 USE mod_coupling
25# ifdef DIAGNOSTICS_UV
26!! USE mod_diags
27# endif
28 USE mod_forces
29 USE mod_grid
30 USE mod_mixing
31 USE mod_ncparam
32 USE mod_ocean
33 USE mod_scalars
34 USE mod_sources
35!
36 USE exchange_2d_mod
37 USE exchange_3d_mod
38# ifdef DISTRIBUTE
39 USE mp_exchange_mod, ONLY : mp_exchange2d, mp_exchange3d
40# endif
41 USE tl_u3dbc_mod, ONLY : tl_u3dbc_tile
42 USE tl_v3dbc_mod, ONLY : tl_v3dbc_tile
43# ifdef UV_DESTAGGERED
44 USE uv_var_change_mod, ONLY : tl_uv_c2a_grid
45# endif
46!
47 implicit none
48!
49 PRIVATE
50 PUBLIC :: tl_step3d_uv
51!
52 CONTAINS
53!
54!***********************************************************************
55 SUBROUTINE tl_step3d_uv (ng, tile)
56!***********************************************************************
57!
58 USE mod_stepping
59!
60! Imported variable declarations.
61!
62 integer, intent(in) :: ng, tile
63!
64! Local variable declarations.
65!
66 character (len=*), parameter :: myfile = &
67 & __FILE__
68!
69# include "tile.h"
70!
71# ifdef PROFILE
72 CALL wclock_on (ng, itlm, 34, __line__, myfile)
73# endif
74 CALL tl_step3d_uv_tile (ng, tile, &
75 & lbi, ubi, lbj, ubj, &
76 & imins, imaxs, jmins, jmaxs, &
77 & nrhs(ng), nstp(ng), nnew(ng), &
78# ifdef MASKING
79 & grid(ng) % umask, &
80 & grid(ng) % vmask, &
81# endif
82# ifdef WET_DRY_NOT_YET
83 & grid(ng) % umask_wet, &
84 & grid(ng) % vmask_wet, &
85# endif
86 & grid(ng) % om_v, &
87 & grid(ng) % on_u, &
88 & grid(ng) % pm, &
89 & grid(ng) % pn, &
90 & grid(ng) % Hz, &
91 & grid(ng) % tl_Hz, &
92 & grid(ng) % z_r, &
93 & grid(ng) % tl_z_r, &
94 & grid(ng) % z_w, &
95 & grid(ng) % tl_z_w, &
96 & mixing(ng) % Akv, &
97 & mixing(ng) % tl_Akv, &
98 & coupling(ng) % DU_avg1, &
99 & coupling(ng) % tl_DU_avg1, &
100 & coupling(ng) % DV_avg1, &
101 & coupling(ng) % tl_DV_avg1, &
102 & coupling(ng) % DU_avg2, &
103 & coupling(ng) % tl_DU_avg2, &
104 & coupling(ng) % DV_avg2, &
105 & coupling(ng) % tl_DV_avg2, &
106 & ocean(ng) % tl_ru, &
107 & ocean(ng) % tl_rv, &
108# ifdef DIAGNOSTICS_UV
109!! & DIAGS(ng) % DiaU2wrk, &
110!! & DIAGS(ng) % DiaV2wrk, &
111!! & DIAGS(ng) % DiaU2int, &
112!! & DIAGS(ng) % DiaV2int, &
113!! & DIAGS(ng) % DiaU3wrk, &
114!! & DIAGS(ng) % DiaV3wrk, &
115!! & DIAGS(ng) % DiaRU, &
116!! & DIAGS(ng) % DiaRV, &
117# endif
118 & ocean(ng) % u, &
119 & ocean(ng) % tl_u, &
120 & ocean(ng) % v, &
121 & ocean(ng) % tl_v, &
122 & ocean(ng) % tl_ubar, &
123 & ocean(ng) % tl_vbar, &
124# ifdef NEARSHORE_MELLOR
125 & ocean(ng) % ubar_stokes, &
126 & ocean(ng) % tl_ubar_stokes, &
127 & ocean(ng) % vbar_stokes, &
128 & ocean(ng) % tl_vbar_stokes, &
129 & ocean(ng) % u_stokes, &
130 & ocean(ng) % tl_u_stokes, &
131 & ocean(ng) % v_stokes, &
132 & ocean(ng) % tl_v_stokes, &
133# endif
134 & grid(ng) % Huon, &
135 & grid(ng) % tl_Huon, &
136 & grid(ng) % Hvom, &
137 & grid(ng) % tl_Hvom)
138# ifdef PROFILE
139 CALL wclock_off (ng, itlm, 34, __line__, myfile)
140# endif
141!
142 RETURN
143 END SUBROUTINE tl_step3d_uv
144!
145!***********************************************************************
146 SUBROUTINE tl_step3d_uv_tile (ng, tile, &
147 & LBi, UBi, LBj, UBj, &
148 & IminS, ImaxS, JminS, JmaxS, &
149 & nrhs, nstp, nnew, &
150# ifdef MASKING
151 & umask, vmask, &
152# endif
153# ifdef WET_DRY_NOT_YET
154 & umask_wet, vmask_wet, &
155# endif
156 & om_v, on_u, pm, pn, &
157 & Hz, tl_Hz, &
158 & z_r, tl_z_r, &
159 & z_w, tl_z_w, &
160 & Akv, tl_Akv, &
161 & DU_avg1, tl_DU_avg1, &
162 & DV_avg1, tl_DV_avg1, &
163 & DU_avg2, tl_DU_avg2, &
164 & DV_avg2, tl_DV_avg2, &
165 & tl_ru, tl_rv, &
166# ifdef DIAGNOSTICS_UV
167!! & DiaU2wrk, DiaV2wrk, &
168!! & DiaU2int, DiaV2int, &
169!! & DiaU3wrk, DiaV3wrk, &
170!! & DiaRU, DiaRV, &
171# endif
172 & u, tl_u, &
173 & v, tl_v, &
174 & tl_ubar, tl_vbar, &
175# ifdef NEARSHORE_MELLOR
176 & ubar_stokes, tl_ubar_stokes, &
177 & vbar_stokes, tl_vbar_stokes, &
178 & u_stokes, tl_u_stokes, &
179 & v_stokes, tl_v_stokes, &
180# endif
181 & Huon, tl_Huon, &
182 & Hvom, tl_Hvom)
183!***********************************************************************
184!
185! Imported variable declarations.
186!
187 integer, intent(in) :: ng, tile
188 integer, intent(in) :: LBi, UBi, LBj, UBj
189 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
190 integer, intent(in) :: nrhs, nstp, nnew
191!
192# ifdef ASSUMED_SHAPE
193# ifdef MASKING
194 real(r8), intent(in) :: umask(LBi:,LBj:)
195 real(r8), intent(in) :: vmask(LBi:,LBj:)
196# endif
197# ifdef WET_DRY_NOT_YET
198 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
199 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
200# endif
201 real(r8), intent(in) :: om_v(LBi:,LBj:)
202 real(r8), intent(in) :: on_u(LBi:,LBj:)
203 real(r8), intent(in) :: pm(LBi:,LBj:)
204 real(r8), intent(in) :: pn(LBi:,LBj:)
205 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
206 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
207 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
208 real(r8), intent(in) :: Akv(LBi:,LBj:,0:)
209 real(r8), intent(in) :: DU_avg1(LBi:,LBj:)
210 real(r8), intent(in) :: DV_avg1(LBi:,LBj:)
211 real(r8), intent(in) :: DU_avg2(LBi:,LBj:)
212 real(r8), intent(in) :: DV_avg2(LBi:,LBj:)
213 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
214 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
215# ifdef NEARSHORE_MELLOR
216 real(r8), intent(in) :: ubar_stokes(LBi:,LBj:)
217 real(r8), intent(in) :: vbar_stokes(LBi:,LBj:)
218 real(r8), intent(in) :: tl_ubar_stokes(LBi:,LBj:)
219 real(r8), intent(in) :: tl_vbar_stokes(LBi:,LBj:)
220# endif
221 real(r8), intent(in) :: tl_Hz(LBi:,LBj:,:)
222 real(r8), intent(in) :: tl_z_r(LBi:,LBj:,:)
223 real(r8), intent(in) :: tl_z_w(LBi:,LBj:,0:)
224 real(r8), intent(in) :: tl_Akv(LBi:,LBj:,0:)
225 real(r8), intent(in) :: tl_DU_avg1(LBi:,LBj:)
226 real(r8), intent(in) :: tl_DV_avg1(LBi:,LBj:)
227 real(r8), intent(in) :: tl_DU_avg2(LBi:,LBj:)
228 real(r8), intent(in) :: tl_DV_avg2(LBi:,LBj:)
229 real(r8), intent(in) :: tl_ru(LBi:,LBj:,0:,:)
230 real(r8), intent(in) :: tl_rv(LBi:,LBj:,0:,:)
231
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) :: tl_u(LBi:,LBj:,:,:)
245 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
246# ifdef NEARSHORE_MELLOR
247 real(r8), intent(inout) :: u_stokes(LBi:,LBj:,:)
248 real(r8), intent(inout) :: v_stokes(LBi:,LBj:,:)
249 real(r8), intent(inout) :: tl_u_stokes(LBi:,LBj:,:)
250 real(r8), intent(inout) :: tl_v_stokes(LBi:,LBj:,:)
251# endif
252 real(r8), intent(out) :: tl_ubar(LBi:,LBj:,:)
253 real(r8), intent(out) :: tl_vbar(LBi:,LBj:,:)
254 real(r8), intent(out) :: tl_Huon(LBi:,LBj:,:)
255 real(r8), intent(out) :: tl_Hvom(LBi:,LBj:,:)
256
257# else
258
259# ifdef MASKING
260 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
261 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
262# endif
263# ifdef WET_DRY_NOT_YET
264 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
265 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
266# endif
267 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
268 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
269 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
270 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
271 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
272 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
273 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
274 real(r8), intent(in) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
275 real(r8), intent(in) :: DU_avg1(LBi:UBi,LBj:UBj)
276 real(r8), intent(in) :: DV_avg1(LBi:UBi,LBj:UBj)
277 real(r8), intent(in) :: DU_avg2(LBi:UBi,LBj:UBj)
278 real(r8), intent(in) :: DV_avg2(LBi:UBi,LBj:UBj)
279 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
280 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
281# ifdef NEARSHORE_MELLOR
282 real(r8), intent(in) :: ubar_stokes(LBi:UBi,LBj:UBj)
283 real(r8), intent(in) :: vbar_stokes(LBi:UBi,LBj:UBj)
284 real(r8), intent(in) :: tl_ubar_stokes(LBi:UBi,LBj:UBj)
285 real(r8), intent(in) :: tl_vbar_stokes(LBi:UBi,LBj:UBj)
286# endif
287
288 real(r8), intent(in) :: tl_Hz(LBi:UBi,LBj:UBj,N(ng))
289 real(r8), intent(in) :: tl_z_r(LBi:UBi,LBj:UBj,N(ng))
290 real(r8), intent(in) :: tl_z_w(LBi:UBi,LBj:UBj,0:N(ng))
291 real(r8), intent(in) :: tl_Akv(LBi:UBi,LBj:UBj,0:N(ng))
292 real(r8), intent(in) :: tl_DU_avg1(LBi:UBi,LBj:UBj)
293 real(r8), intent(in) :: tl_DV_avg1(LBi:UBi,LBj:UBj)
294 real(r8), intent(in) :: tl_DU_avg2(LBi:UBi,LBj:UBj)
295 real(r8), intent(in) :: tl_DV_avg2(LBi:UBi,LBj:UBj)
296 real(r8), intent(in) :: tl_ru(LBi:UBi,LBj:UBj,0:N(ng),2)
297 real(r8), intent(in) :: tl_rv(LBi:UBi,LBj:UBj,0:N(ng),2)
298# ifdef DIAGNOSTICS_UV
299!! real(r8), intent(inout) :: DiaU2wrk(LBi:UBi,LBj:UBj,NDM2d)
300!! real(r8), intent(inout) :: DiaV2wrk(LBi:UBi,LBj:UBj,NDM2d)
301!! real(r8), intent(inout) :: DiaU2int(LBi:UBi,LBj:UBj,NDM2d)
302!! real(r8), intent(inout) :: DiaV2int(LBi:UBi,LBj:UBj,NDM2d)
303!! real(r8), intent(inout) :: DiaU3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
304!! real(r8), intent(inout) :: DiaV3wrk(LBi:UBi,LBj:UBj,N(ng),NDM3d)
305!! real(r8), intent(inout) :: DiaRU(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
306!! real(r8), intent(inout) :: DiaRV(LBi:UBi,LBj:UBj,N(ng),2,NDrhs)
307# endif
308 real(r8), intent(inout) :: Huon(LBi:UBi,LBj:UBj,N(ng))
309 real(r8), intent(inout) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
310 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
311 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
312# ifdef NEARSHORE_MELLOR
313 real(r8), intent(inout) :: u_stokes(LBi:UBi,LBj:UBj,N(ng))
314 real(r8), intent(inout) :: v_stokes(LBi:UBi,LBj:UBj,N(ng))
315 real(r8), intent(inout) :: tl_u_stokes(LBi:UBi,LBj:UBj,N(ng))
316 real(r8), intent(inout) :: tl_v_stokes(LBi:UBi,LBj:UBj,N(ng))
317# endif
318 real(r8), intent(out) :: tl_ubar(LBi:UBi,LBj:UBj,3)
319 real(r8), intent(out) :: tl_vbar(LBi:UBi,LBj:UBj,3)
320 real(r8), intent(out) :: tl_Huon(LBi:UBi,LBj:UBj,N(ng))
321 real(r8), intent(out) :: tl_Hvom(LBi:UBi,LBj:UBj,N(ng))
322# endif
323!
324! Local variable declarations.
325!
326 integer :: i, idiag, is, j, k
327!
328 real(r8) :: cff, cff1, cff2
329 real(r8) :: tl_cff, tl_cff1, tl_cff2
330!
331 real(r8), dimension(IminS:ImaxS) :: CF1
332 real(r8), dimension(IminS:ImaxS) :: FC1
333# ifdef NEARSHORE_MELLOR
334 real(r8), dimension(IminS:ImaxS) :: CFs1
335 real(r8), dimension(IminS:ImaxS) :: DCs1
336# endif
337 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: AK
338 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BC
339 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
340 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC
341 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC1
342 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
343# ifdef NEARSHORE_MELLOR
344 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CFs
345 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DCs
346# endif
347 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hzk
348 real(r8), dimension(IminS:ImaxS,N(ng)) :: oHz
349# ifdef DIAGNOSTICS_UV
350!! real(r8), dimension(IminS:ImaxS,0:N(ng)) :: wrk
351!! real(r8), dimension(IminS:ImaxS,1:NDM2d-1) :: Dwrk
352# endif
353 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_AK
354 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_BC
355 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_CF
356 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_DC
357 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_FC
358# ifdef NEARSHORE_MELLOR
359 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_CFs
360 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: tl_DCs
361# endif
362 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_Hzk
363 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_oHz
364!
365# include "set_bounds.h"
366!
367!-----------------------------------------------------------------------
368! Time step momentum equation in the XI-direction.
369!-----------------------------------------------------------------------
370!
371 DO j=jstr,jend
372 DO i=istru,iend
373 ak(i,0)=0.5_r8*(akv(i-1,j,0)+ &
374 & akv(i ,j,0))
375 tl_ak(i,0)=0.5_r8*(tl_akv(i-1,j,0)+ &
376 & tl_akv(i ,j,0))
377 DO k=1,n(ng)
378 ak(i,k)=0.5_r8*(akv(i-1,j,k)+ &
379 & akv(i ,j,k))
380 tl_ak(i,k)=0.5_r8*(tl_akv(i-1,j,k)+ &
381 & tl_akv(i ,j,k))
382 hzk(i,k)=0.5_r8*(hz(i-1,j,k)+ &
383 & hz(i ,j,k))
384 tl_hzk(i,k)=0.5_r8*(tl_hz(i-1,j,k)+ &
385 & tl_hz(i ,j,k))
386# if defined SPLINES_VVISC || defined DIAGNOSTICS_UV
387 ohz(i,k)=1.0_r8/hzk(i,k)
388 tl_ohz(i,k)=-ohz(i,k)*ohz(i,k)*tl_hzk(i,k)
389# endif
390 END DO
391 END DO
392!
393! Time step right-hand-side terms.
394!
395# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
396 IF (iic(ng).eq.ntfirst(ng).and.soinitial(ng)) THEN
397# else
398 IF (iic(ng).eq.ntfirst(ng)) THEN
399# endif
400 cff=0.25_r8*dt(ng)
401# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
402 ELSE IF (iic(ng).eq.(ntfirst(ng)+1).and.soinitial(ng)) THEN
403# else
404 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
405# endif
406 cff=0.25_r8*dt(ng)*3.0_r8/2.0_r8
407 ELSE
408 cff=0.25_r8*dt(ng)*23.0_r8/12.0_r8
409 END IF
410 DO i=istru,iend
411 dc(i,0)=cff*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
412 END DO
413!
414! The BASIC STATE "u" used below must be in transport units, but "u"
415! retrived is in m/s so we multiply by "Hzk".
416!
417 DO k=1,n(ng)
418 DO i=istru,iend
419!^ u(i,j,k,nnew)=u(i,j,k,nnew)+ &
420!^ & DC(i,0)*ru(i,j,k,nrhs)
421!^
422 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)+ &
423 & dc(i,0)*tl_ru(i,j,k,nrhs)
424# ifdef SPLINES_VVISC
425!^ u(i,j,k,nnew)=u(i,j,k,nnew)*oHz(i,k)
426!^
427 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*ohz(i,k)+ &
428 & (u(i,j,k,nnew)*hzk(i,k))*tl_ohz(i,k)
429# endif
430# ifdef DIAGNOSTICS_UV
431!! DO idiag=1,M3pgrd
432!! DiaU3wrk(i,j,k,idiag)=(DiaU3wrk(i,j,k,idiag)+ &
433!! & DC(i,0)*DiaRU(i,j,k,nrhs,idiag))* &
434!! & oHz(i,k)
435!! END DO
436# if defined UV_VIS2 || defined UV_VIS4
437!! DiaU3wrk(i,j,k,M3xvis)=DiaU3wrk(i,j,k,M3xvis)*oHz(i,k)
438!! DiaU3wrk(i,j,k,M3yvis)=DiaU3wrk(i,j,k,M3yvis)*oHz(i,k)
439!! DiaU3wrk(i,j,k,M3hvis)=DiaU3wrk(i,j,k,M3hvis)*oHz(i,k)
440# endif
441!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)*oHz(i,k)
442# ifdef BODYFORCE
443!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+ &
444!! & DC(i,0)*DiaRU(i,j,k,nrhs,M3vvis)* &
445!! & oHz(i,k)
446# endif
447!! DiaU3wrk(i,j,k,M3rate)=DiaU3wrk(i,j,k,M3rate)*oHz(i,k)
448# endif
449 END DO
450 END DO
451
452# ifdef SPLINES_VVISC
453!
454! Apply spline algorithim to BASIC STATE "u" which should be
455! in units of velocity (m/s). DC will be used in the tangent
456! linear spline algorithm below. Solve BASIC STATE tridiagonal
457! system.
458!
459 cff1=1.0_r8/6.0_r8
460 DO k=1,n(ng)-1
461 DO i=istru,iend
462 fc(i,k)=cff1*hzk(i,k )-dt(ng)*ak(i,k-1)*ohz(i,k )
463 cf(i,k)=cff1*hzk(i,k+1)-dt(ng)*ak(i,k+1)*ohz(i,k+1)
464 END DO
465 END DO
466 DO i=istru,iend
467 cf(i,0)=0.0_r8
468 dc(i,0)=0.0_r8
469 END DO
470!
471! LU decomposition and forward substitution.
472!
473 cff1=1.0_r8/3.0_r8
474 DO k=1,n(ng)-1
475 DO i=istru,iend
476 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
477 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
478 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
479 cf(i,k)=cff*cf(i,k)
480 dc(i,k)=cff*(u(i,j,k+1,nnew)-u(i,j,k,nnew)- &
481 & fc(i,k)*dc(i,k-1))
482 END DO
483 END DO
484!
485! Backward substitution. Save DC for the tangent linearization.
486! DC is scaled later by AK.
487!
488 DO i=istru,iend
489 dc(i,n(ng))=0.0_r8
490 END DO
491 DO k=n(ng)-1,1,-1
492 DO i=istru,iend
493 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
494 END DO
495 END DO
496!
497! Use conservative, parabolic spline reconstruction of tangent
498! linear vertical viscosity derivatives. Then, time step vertical
499! viscosity term implicitly by solving a tridiagonal system.
500!
501 cff1=1.0_r8/6.0_r8
502 DO k=1,n(ng)-1
503 DO i=istru,iend
504 fc(i,k)=cff1*hzk(i,k )- &
505 & dt(ng)*ak(i,k-1)*ohz(i,k )
506 tl_fc(i,k)=cff1*tl_hzk(i,k )- &
507 & dt(ng)*(tl_ak(i,k-1)*ohz(i,k )+ &
508 & ak(i,k-1)*tl_ohz(i,k ))
509 cf(i,k)=cff1*hzk(i,k+1)- &
510 & dt(ng)*ak(i,k+1)*ohz(i,k+1)
511 tl_cf(i,k)=cff1*tl_hzk(i,k+1)- &
512 & dt(ng)*(tl_ak(i,k+1)*ohz(i,k+1)+ &
513 & ak(i,k+1)*tl_ohz(i,k+1))
514 END DO
515 END DO
516 DO i=istru,iend
517 cf(i,0)=0.0_r8
518 tl_cf(i,0)=0.0_r8
519 tl_dc(i,0)=0.0_r8
520 END DO
521!
522! Tangent linear LU decomposition and forward substitution.
523!
524 cff1=1.0_r8/3.0_r8
525 DO k=1,n(ng)-1
526 DO i=istru,iend
527 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
528 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
529 tl_bc(i,k)=cff1*(tl_hzk(i,k)+tl_hzk(i,k+1))+ &
530 & dt(ng)*(tl_ak(i,k)*(ohz(i,k)+ohz(i,k+1))+ &
531 & ak(i,k)*(tl_ohz(i,k)+tl_ohz(i,k+1)))
532 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
533 cf(i,k)=cff*cf(i,k)
534 tl_dc(i,k)=cff*(tl_u(i,j,k+1,nnew)- &
535 & tl_u(i,j,k ,nnew)- &
536 & (tl_fc(i,k)*dc(i,k-1)+ &
537 & tl_bc(i,k)*dc(i,k )+ &
538 & tl_cf(i,k)*dc(i,k+1))- &
539 & fc(i,k)*tl_dc(i,k-1))
540 END DO
541 END DO
542!
543! Tangent linear backward substitution.
544!
545 DO i=istru,iend
546 tl_dc(i,n(ng))=0.0_r8
547 END DO
548 DO k=n(ng)-1,1,-1
549 DO i=istru,iend
550 tl_dc(i,k)=tl_dc(i,k)-cf(i,k)*tl_dc(i,k+1)
551 END DO
552 END DO
553!
554! Compute tl_DC before multiplying BASIC STATE spline gradients
555! DC by AK.
556!
557 DO k=1,n(ng)
558 DO i=istru,iend
559 tl_dc(i,k)=tl_dc(i,k)*ak(i,k)+dc(i,k)*tl_ak(i,k)
560 dc(i,k)=dc(i,k)*ak(i,k)
561!^ cff=dt(ng)*oHz(i,k)*(DC(i,k)-DC(i,k-1))
562!^
563 tl_cff=dt(ng)*(tl_ohz(i,k)*(dc(i,k)-dc(i,k-1))+ &
564 & ohz(i,k)*(tl_dc(i,k)-tl_dc(i,k-1)))
565!^ u(i,j,k,nnew)=u(i,j,k,nnew)+cff
566!^
567 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)+tl_cff
568# ifdef DIAGNOSTICS_UV
569!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+cff
570# endif
571 END DO
572 END DO
573# else
574!
575! Compute off-diagonal coefficients [lambda*dt*Akv/Hz] for the
576! implicit vertical viscosity term at horizontal U-points and
577! vertical W-points.
578!
579! NOTE: The original code solves the tridiagonal system A*u=r where
580! A is a matrix and u and r are vectors. We need to solve the
581! tangent linear of this system which is A*tl_u+tl_A*u=tl_r.
582! Here, tl_A*u and tl_r are known, so we must solve for tl_u
583! by inverting A*tl_u=tl_r-tl_A*u.
584!
585 cff=-lambda*dt(ng)/0.5_r8
586 DO k=1,n(ng)-1
587 DO i=istru,iend
588 cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
589 & z_r(i,j,k )-z_r(i-1,j,k ))
590 tl_cff1=-cff1*cff1*(tl_z_r(i,j,k+1)+tl_z_r(i-1,j,k+1)- &
591 & tl_z_r(i,j,k )-tl_z_r(i-1,j,k ))
592 fc(i,k)=cff*cff1*ak(i,k)
593 tl_fc(i,k)=cff*(tl_cff1*ak(i,k)+cff1*tl_ak(i,k))
594 END DO
595 END DO
596 DO i=istru,iend
597 fc(i,0)=0.0_r8
598 tl_fc(i,0)=0.0_r8
599 fc(i,n(ng))=0.0_r8
600 tl_fc(i,n(ng))=0.0_r8
601 END DO
602!
603! Solve the tangent linear tridiagonal system.
604! (DC is a tangent linear variable here).
605!
606 DO k=1,n(ng)
607 DO i=istru,iend
608 bc(i,k)=hzk(i,k)-fc(i,k)-fc(i,k-1)
609 tl_bc(i,k)=tl_hzk(i,k)-tl_fc(i,k)-tl_fc(i,k-1)
610 END DO
611 END DO
612 DO k=2,n(ng)-1
613 DO i=istru,iend
614 dc(i,k)=tl_u(i,j,k,nnew)- &
615 & (tl_fc(i,k-1)*u(i,j,k-1,nnew)+ &
616 & tl_bc(i,k )*u(i,j,k ,nnew)+ &
617 & tl_fc(i,k )*u(i,j,k+1,nnew))
618 END DO
619 END DO
620 DO i=istru,iend
621 dc(i,1)=tl_u(i,j,1,nnew)- &
622 & (tl_bc(i,1)*u(i,j,1,nnew)+ &
623 & tl_fc(i,1)*u(i,j,2,nnew))
624 dc(i,n(ng))=tl_u(i,j,n(ng),nnew)- &
625 & (tl_fc(i,n(ng)-1)*u(i,j,n(ng)-1,nnew)+ &
626 & tl_bc(i,n(ng) )*u(i,j,n(ng) ,nnew))
627 END DO
628 DO i=istru,iend
629 cff=1.0_r8/bc(i,1)
630 cf(i,1)=cff*fc(i,1)
631 dc(i,1)=cff*dc(i,1)
632 END DO
633 DO k=2,n(ng)-1
634 DO i=istru,iend
635 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
636 cf(i,k)=cff*fc(i,k)
637 dc(i,k)=cff*(dc(i,k)-fc(i,k-1)*dc(i,k-1))
638 END DO
639 END DO
640!
641! Compute new solution by back substitution.
642! (DC is a tangent linear variable here).
643!
644 DO i=istru,iend
645# ifdef DIAGNOSTICS_UV
646!! wrk(i,N(ng))=u(i,j,N(ng),nnew)*oHz(i,N(ng))
647# endif
648 dc(i,n(ng))=(dc(i,n(ng))-fc(i,n(ng)-1)*dc(i,n(ng)-1))/ &
649 & (bc(i,n(ng))-fc(i,n(ng)-1)*cf(i,n(ng)-1))
650!^ u(i,j,N(ng),nnew)=DC(i,N(ng))
651!^
652 tl_u(i,j,n(ng),nnew)=dc(i,n(ng))
653# ifdef DIAGNOSTICS_UV
654!! DiaU3wrk(i,j,N(ng),M3vvis)=DiaU3wrk(i,j,N(ng),M3vvis)+ &
655!! & u(i,j,N(ng),nnew)-wrk(i,N(ng))
656# endif
657 END DO
658 DO k=n(ng)-1,1,-1
659 DO i=istru,iend
660# ifdef DIAGNOSTICS_UV
661!! wrk(i,k)=u(i,j,k,nnew)*oHz(i,k)
662# endif
663 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
664!^ u(i,j,k,nnew)=DC(i,k)
665!^
666 tl_u(i,j,k,nnew)=dc(i,k)
667# ifdef DIAGNOSTICS_UV
668!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)+ &
669!! & u(i,j,k,nnew)-wrk(i,k)
670# endif
671 END DO
672 END DO
673# endif
674!
675! Replace INTERIOR POINTS incorrect vertical mean with more accurate
676! barotropic component, ubar=DU_avg1/(D*on_u). Recall that, D=CF(:,0).
677!
678 DO i=istru,iend
679 cf(i,0)=hzk(i,1)
680 tl_cf(i,0)=tl_hzk(i,1)
681 dc(i,0)=u(i,j,1,nnew)*hzk(i,1)
682 tl_dc(i,0)=tl_u(i,j,1,nnew)*hzk(i,1)+ &
683 & u(i,j,1,nnew)*tl_hzk(i,1)
684# ifdef NEARSHORE_MELLOR
685 dcs(i,0)=u_stokes(i,j,1)*hzk(i,1)
686 tl_dcs(i,0)=tl_u_stokes(i,j,1)*hzk(i,1)+ &
687 & u_stokes(i,j,1)*tl_hzk(i,1)
688# endif
689# ifdef DIAGNOSTICS_UV
690!! Dwrk(i,M2pgrd)=DiaU3wrk(i,j,1,M3pgrd)*Hzk(i,1)
691!! Dwrk(i,M2bstr)=DiaU3wrk(i,j,1,M3vvis)*Hzk(i,1)
692# ifdef UV_COR
693!! Dwrk(i,M2fcor)=DiaU3wrk(i,j,1,M3fcor)*Hzk(i,1)
694# endif
695# if defined UV_VIS2 || defined UV_VIS4
696!! Dwrk(i,M2xvis)=DiaU3wrk(i,j,1,M3xvis)*Hzk(i,1)
697!! Dwrk(i,M2yvis)=DiaU3wrk(i,j,1,M3yvis)*Hzk(i,1)
698!! Dwrk(i,M2hvis)=DiaU3wrk(i,j,1,M3hvis)*Hzk(i,1)
699# endif
700# ifdef UV_ADV
701!! Dwrk(i,M2xadv)=DiaU3wrk(i,j,1,M3xadv)*Hzk(i,1)
702!! Dwrk(i,M2yadv)=DiaU3wrk(i,j,1,M3yadv)*Hzk(i,1)
703!! Dwrk(i,M2hadv)=DiaU3wrk(i,j,1,M3hadv)*Hzk(i,1)
704# endif
705# ifdef NEARSHORE_MELLOR
706!! Dwrk(i,M2hrad)=DiaU3wrk(i,j,1,M3hrad)*Hzk(i,1)
707# endif
708# endif
709 END DO
710 DO k=2,n(ng)
711 DO i=istru,iend
712 cf(i,0)=cf(i,0)+hzk(i,k)
713 tl_cf(i,0)=tl_cf(i,0)+tl_hzk(i,k)
714 dc(i,0)=dc(i,0)+u(i,j,k,nnew)*hzk(i,k)
715 tl_dc(i,0)=tl_dc(i,0)+ &
716 & tl_u(i,j,k,nnew)*hzk(i,k)+ &
717 & u(i,j,k,nnew)*tl_hzk(i,k)
718# ifdef NEARSHORE_MELLOR
719 dcs(i,0)=dcs(i,0)+u_stokes(i,j,k)*hzk(i,k)
720 tl_dcs(i,0)=tl_dcs(i,0)+ &
721 & tl_u_stokes(i,j,k)*hzk(i,k)+ &
722 & u_stokes(i,j,k)*tl_hzk(i,k)
723# endif
724# ifdef DIAGNOSTICS_UV
725!! Dwrk(i,M2pgrd)=Dwrk(i,M2pgrd)+ &
726!! & DiaU3wrk(i,j,k,M3pgrd)*Hzk(i,k)
727!! Dwrk(i,M2bstr)=Dwrk(i,M2bstr)+ &
728!! & DiaU3wrk(i,j,k,M3vvis)*Hzk(i,k)
729# ifdef UV_COR
730!! Dwrk(i,M2fcor)=Dwrk(i,M2fcor)+ &
731!! & DiaU3wrk(i,j,k,M3fcor)*Hzk(i,k)
732# endif
733# if defined UV_VIS2 || defined UV_VIS4
734!! Dwrk(i,M2xvis)=Dwrk(i,M2xvis)+ &
735!! & DiaU3wrk(i,j,k,M3xvis)*Hzk(i,k)
736!! Dwrk(i,M2yvis)=Dwrk(i,M2yvis)+ &
737!! & DiaU3wrk(i,j,k,M3yvis)*Hzk(i,k)
738!! Dwrk(i,M2hvis)=Dwrk(i,M2hvis)+ &
739!! & DiaU3wrk(i,j,k,M3hvis)*Hzk(i,k)
740# endif
741# ifdef UV_ADV
742!! Dwrk(i,M2xadv)=Dwrk(i,M2xadv)+ &
743!! & DiaU3wrk(i,j,k,M3xadv)*Hzk(i,k)
744!! Dwrk(i,M2yadv)=Dwrk(i,M2yadv)+ &
745!! & DiaU3wrk(i,j,k,M3yadv)*Hzk(i,k)
746!! Dwrk(i,M2hadv)=Dwrk(i,M2hadv)+ &
747!! & DiaU3wrk(i,j,k,M3hadv)*Hzk(i,k)
748# endif
749# ifdef NEARSHORE_MELLOR
750!! Dwrk(i,M2hrad)=Dwrk(i,M2hrad)+ &
751!! & DiaU3wrk(i,j,k,M3hrad)*Hzk(i,k)
752# endif
753# endif
754 END DO
755 END DO
756 DO i=istru,iend
757 dc1(i,0)=dc(i,0) ! intermediate
758 cff1=1.0_r8/(cf(i,0)*on_u(i,j))
759 tl_cff1=-cff1*cff1*tl_cf(i,0)*on_u(i,j)
760 dc(i,0)=(dc(i,0)*on_u(i,j)-du_avg1(i,j))*cff1 ! recursive
761 tl_dc(i,0)=(tl_dc(i,0)*on_u(i,j)-tl_du_avg1(i,j))*cff1+ &
762 & (dc1(i,0)*on_u(i,j)-du_avg1(i,j))*tl_cff1
763# ifdef NEARSHORE_MELLOR
764 dcs1(i)=dcs(i,0) ! intermediate
765 cff2=1.0_r8/cf(i,0)
766 tl_cff2=-cff2*cff2*tl_cf(i,0)
767 dcs(i,0)=dcs(i,0)*cff2-ubar_stokes(i,j) ! recursive
768 tl_dcs(i,0)=tl_dcs(i,0)*cff2+ &
769 & dcs1(i)*tl_cff2-tl_ubar_stokes(i,j)
770# endif
771# ifdef DIAGNOSTICS_UV
772!! DO idiag=1,M2pgrd
773!! Dwrk(i,idiag)=(Dwrk(i,idiag)*on_u(i,j)- &
774!! & DiaU2wrk(i,j,idiag))*cff1
775!! END DO
776!! Dwrk(i,M2bstr)=(Dwrk(i,M2bstr)*on_u(i,j)- &
777!! & DiaU2wrk(i,j,M2bstr)-DiaU2wrk(i,j,M2sstr))* &
778!! & cff1
779# endif
780 END DO
781!
782! Couple and update new solution.
783!
784 DO k=1,n(ng)
785 DO i=istru,iend
786!^ u(i,j,k,nnew)=u(i,j,k,nnew)-DC(i,0)
787!^
788 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)-tl_dc(i,0)
789# ifdef MASKING
790!^ u(i,j,k,nnew)=u(i,j,k,nnew)*umask(i,j)
791!^
792 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*umask(i,j)
793# endif
794# ifdef WET_DRY_NOT_YET
795!^ u(i,j,k,nnew)=u(i,j,k,nnew)*umask_wet(i,j)
796!^
797 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)*umask_wet(i,j)
798# endif
799# ifdef NEARSHORE_MELLOR
800!^ u_stokes(i,j,k)=u_stokes(i,j,k)-DCs(i,0)
801!^
802 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)-tl_dcs(i,0)
803# ifdef MASKING
804!^ u_stokes(i,j,k)=u_stokes(i,j,k)*umask(i,j)
805!^
806 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*umask(i,j)
807# endif
808# ifdef WET_DRY_NOT_YET
809!^ u_stokes(i,j,k)=u_stokes(i,j,k)*umask_wet(i,j)
810!^
811 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)*umask_wet(i,j)
812# endif
813# endif
814# ifdef DIAGNOSTICS_UV
815!! DiaU3wrk(i,j,k,M3pgrd)=DiaU3wrk(i,j,k,M3pgrd)- &
816!! & Dwrk(i,M2pgrd)
817!! DiaU3wrk(i,j,k,M3vvis)=DiaU3wrk(i,j,k,M3vvis)- &
818!! & Dwrk(i,M2bstr)
819# ifdef UV_COR
820!! DiaU3wrk(i,j,k,M3fcor)=DiaU3wrk(i,j,k,M3fcor)- &
821!! & Dwrk(i,M2fcor)
822# endif
823# if defined UV_VIS2 || defined UV_VIS4
824!! DiaU3wrk(i,j,k,M3xvis)=DiaU3wrk(i,j,k,M3xvis)- &
825!! & Dwrk(i,M2xvis)
826!! DiaU3wrk(i,j,k,M3yvis)=DiaU3wrk(i,j,k,M3yvis)- &
827!! & Dwrk(i,M2yvis)
828!! DiaU3wrk(i,j,k,M3hvis)=DiaU3wrk(i,j,k,M3hvis)- &
829!! & Dwrk(i,M2hvis)
830# endif
831# ifdef UV_ADV
832!! DiaU3wrk(i,j,k,M3xadv)=DiaU3wrk(i,j,k,M3xadv)- &
833!! & Dwrk(i,M2xadv)
834!! DiaU3wrk(i,j,k,M3yadv)=DiaU3wrk(i,j,k,M3yadv)- &
835!! & Dwrk(i,M2yadv)
836!! DiaU3wrk(i,j,k,M3hadv)=DiaU3wrk(i,j,k,M3hadv)- &
837!! & Dwrk(i,M2hadv)
838# endif
839# ifdef NEARSHORE_MELLOR
840!! DiaU3wrk(i,j,k,M3hrad)=DiaU3wrk(i,j,k,M3hrad)- &
841!! & Dwrk(i,M2hrad)
842# endif
843# endif
844 END DO
845 END DO
846
847# if defined DIAGNOSTICS_UV && defined MASKING
848!! DO k=1,N(ng)
849!! DO i=IstrU,Iend
850!! DO idiag=1,NDM3d
851!! DiaU3wrk(i,j,k,idiag)=DiaU3wrk(i,j,k,idiag)*umask(i,j)
852!! END DO
853!! END DO
854!! END DO
855# endif
856!
857!-----------------------------------------------------------------------
858! Time step momentum equation in the ETA-direction.
859!-----------------------------------------------------------------------
860!
861 IF (j.ge.jstrv) THEN
862 DO i=istr,iend
863 ak(i,0)=0.5_r8*(akv(i,j-1,0)+ &
864 & akv(i,j ,0))
865 tl_ak(i,0)=0.5_r8*(tl_akv(i,j-1,0)+ &
866 & tl_akv(i,j ,0))
867 DO k=1,n(ng)
868 ak(i,k)=0.5_r8*(akv(i,j-1,k)+ &
869 & akv(i,j ,k))
870 tl_ak(i,k)=0.5_r8*(tl_akv(i,j-1,k)+ &
871 & tl_akv(i,j ,k))
872 hzk(i,k)=0.5_r8*(hz(i,j-1,k)+ &
873 & hz(i,j ,k))
874 tl_hzk(i,k)=0.5_r8*(tl_hz(i,j-1,k)+ &
875 & tl_hz(i,j ,k))
876# if defined SPLINES_VVISC || defined DIAGNOSTICS_UV
877 ohz(i,k)=1.0_r8/hzk(i,k)
878 tl_ohz(i,k)=-ohz(i,k)*ohz(i,k)*tl_hzk(i,k)
879# endif
880 END DO
881 END DO
882!
883! Time step right-hand-side terms.
884!
885# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
886 IF (iic(ng).eq.ntfirst(ng).and.soinitial(ng)) THEN
887# else
888 IF (iic(ng).eq.ntfirst(ng)) THEN
889# endif
890 cff=0.25_r8*dt(ng)
891# if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
892 ELSE IF (iic(ng).eq.(ntfirst(ng)+1).and.soinitial(ng)) THEN
893# else
894 ELSE IF (iic(ng).eq.(ntfirst(ng)+1)) THEN
895# endif
896 cff=0.25_r8*dt(ng)*3.0_r8/2.0_r8
897 ELSE
898 cff=0.25_r8*dt(ng)*23.0_r8/12.0_r8
899 END IF
900 DO i=istr,iend
901 dc(i,0)=cff*(pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))
902 END DO
903!
904! The BASIC STATE "v" used below must be in transport units, but "v"
905! retrived is in m/s so we multiply by "Hzk".
906!
907 DO k=1,n(ng)
908 DO i=istr,iend
909!^ v(i,j,k,nnew)=v(i,j,k,nnew)+DC(i,0)*rv(i,j,k,nrhs)
910!^
911 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)+ &
912 & dc(i,0)*tl_rv(i,j,k,nrhs)
913# ifdef SPLINES_VVISC
914!^ v(i,j,k,nnew)=v(i,j,k,nnew)*oHz(i,k)
915!^
916 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*ohz(i,k)+ &
917 & (v(i,j,k,nnew)*hzk(i,k))*tl_ohz(i,k)
918# endif
919# ifdef DIAGNOSTICS_UV
920!! DO idiag=1,M3pgrd
921!! DiaV3wrk(i,j,k,idiag)=(DiaV3wrk(i,j,k,idiag)+ &
922!! & DC(i,0)* &
923!! & DiaRV(i,j,k,nrhs,idiag))* &
924!! & oHz(i,k)
925!! END DO
926# if defined UV_VIS2 || defined UV_VIS4
927!! DiaV3wrk(i,j,k,M3xvis)=DiaV3wrk(i,j,k,M3xvis)*oHz(i,k)
928!! DiaV3wrk(i,j,k,M3yvis)=DiaV3wrk(i,j,k,M3yvis)*oHz(i,k)
929!! DiaV3wrk(i,j,k,M3hvis)=DiaV3wrk(i,j,k,M3hvis)*oHz(i,k)
930# endif
931!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)*oHz(i,k)
932# ifdef BODYFORCE
933!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+ &
934!! & DC(i,0)*DiaRV(i,j,k,nrhs,M3vvis)* &
935!! & oHz(i,k)
936# endif
937!! DiaV3wrk(i,j,k,M3rate)=DiaV3wrk(i,j,k,M3rate)*oHz(i,k)
938# endif
939 END DO
940 END DO
941
942# ifdef SPLINES_VVISC
943!
944! Apply spline algorithim to BASIC STATE "v" which should be
945! in units of velocity (m/s). DC will be used in the tangent
946! linear spline algorithm below. Solve BASIC STATE tridiagonal
947! system.
948!
949 cff1=1.0_r8/6.0_r8
950 DO k=1,n(ng)-1
951 DO i=istr,iend
952 fc(i,k)=cff1*hzk(i,k )-dt(ng)*ak(i,k-1)*ohz(i,k )
953 cf(i,k)=cff1*hzk(i,k+1)-dt(ng)*ak(i,k+1)*ohz(i,k+1)
954 END DO
955 END DO
956 DO i=istr,iend
957 cf(i,0)=0.0_r8
958 dc(i,0)=0.0_r8
959 END DO
960!
961! LU decomposition and forward substitution.
962!
963 cff1=1.0_r8/3.0_r8
964 DO k=1,n(ng)-1
965 DO i=istr,iend
966 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
967 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
968 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
969 cf(i,k)=cff*cf(i,k)
970 dc(i,k)=cff*(v(i,j,k+1,nnew)-v(i,j,k,nnew)- &
971 & fc(i,k)*dc(i,k-1))
972 END DO
973 END DO
974!
975! Backward substitution. Save DC for the tangent linearization.
976! DC is scaled later by AK.
977!
978 DO i=istr,iend
979 dc(i,n(ng))=0.0_r8
980 END DO
981 DO k=n(ng)-1,1,-1
982 DO i=istr,iend
983 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
984 END DO
985 END DO
986!
987! Use conservative, parabolic spline reconstruction of tangent
988! linear vertical viscosity derivatives. Then, time step vertical
989! viscosity term implicitly by solving a tridiagonal system.
990!
991 cff1=1.0_r8/6.0_r8
992 DO k=1,n(ng)-1
993 DO i=istr,iend
994 fc(i,k)=cff1*hzk(i,k )- &
995 & dt(ng)*ak(i,k-1)*ohz(i,k )
996 tl_fc(i,k)=cff1*tl_hzk(i,k )- &
997 & dt(ng)*(tl_ak(i,k-1)*ohz(i,k )+ &
998 & ak(i,k-1)*tl_ohz(i,k ))
999 cf(i,k)=cff1*hzk(i,k+1)- &
1000 & dt(ng)*ak(i,k+1)*ohz(i,k+1)
1001 tl_cf(i,k)=cff1*tl_hzk(i,k+1)- &
1002 & dt(ng)*(tl_ak(i,k+1)*ohz(i,k+1)+ &
1003 & ak(i,k+1)*tl_ohz(i,k+1))
1004 END DO
1005 END DO
1006 DO i=istr,iend
1007 cf(i,0)=0.0_r8
1008 tl_cf(i,0)=0.0_r8
1009 tl_dc(i,0)=0.0_r8
1010 END DO
1011!
1012! Tangent linear LU decomposition and forward substitution.
1013!
1014 cff1=1.0_r8/3.0_r8
1015 DO k=1,n(ng)-1
1016 DO i=istr,iend
1017 bc(i,k)=cff1*(hzk(i,k)+hzk(i,k+1))+ &
1018 & dt(ng)*ak(i,k)*(ohz(i,k)+ohz(i,k+1))
1019 tl_bc(i,k)=cff1*(tl_hzk(i,k)+tl_hzk(i,k+1))+ &
1020 & dt(ng)*(tl_ak(i,k)*(ohz(i,k)+ohz(i,k+1))+ &
1021 & ak(i,k)*(tl_ohz(i,k)+tl_ohz(i,k+1)))
1022 cff=1.0_r8/(bc(i,k)-fc(i,k)*cf(i,k-1))
1023 cf(i,k)=cff*cf(i,k)
1024 tl_dc(i,k)=cff*(tl_v(i,j,k+1,nnew)- &
1025 & tl_v(i,j,k ,nnew)- &
1026 & (tl_fc(i,k)*dc(i,k-1)+ &
1027 & tl_bc(i,k)*dc(i,k )+ &
1028 & tl_cf(i,k)*dc(i,k+1))- &
1029 & fc(i,k)*tl_dc(i,k-1))
1030 END DO
1031 END DO
1032!
1033! Tangent linear backward substitution.
1034!
1035 DO i=istr,iend
1036 tl_dc(i,n(ng))=0.0_r8
1037 END DO
1038 DO k=n(ng)-1,1,-1
1039 DO i=istr,iend
1040 tl_dc(i,k)=tl_dc(i,k)-cf(i,k)*tl_dc(i,k+1)
1041 END DO
1042 END DO
1043!
1044! Compute tl_DC before multiplying BASIC STATE spline gradients
1045! DC by AK.
1046!
1047 DO k=1,n(ng)
1048 DO i=istr,iend
1049 tl_dc(i,k)=tl_dc(i,k)*ak(i,k)+dc(i,k)*tl_ak(i,k)
1050 dc(i,k)=dc(i,k)*ak(i,k)
1051!^ cff=dt(ng)*oHz(i,k)*(DC(i,k)-DC(i,k-1))
1052!^
1053 tl_cff=dt(ng)*(tl_ohz(i,k)*(dc(i,k)-dc(i,k-1))+ &
1054 & ohz(i,k)*(tl_dc(i,k)-tl_dc(i,k-1)))
1055!^ v(i,j,k,nnew)=v(i,j,k,nnew)+cff
1056!^
1057 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)+tl_cff
1058# ifdef DIAGNOSTICS_UV
1059!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+cff
1060# endif
1061 END DO
1062 END DO
1063# else
1064!
1065! Compute off-diagonal coefficients [lambda*dt*Akv/Hz] for the
1066! implicit vertical viscosity term at horizontal V-points and
1067! vertical W-points.
1068!
1069! NOTE: The original code solves the tridiagonal system A*v=r where
1070! A is a matrix and v and r are vectors. We need to solve the
1071! tangent linear of this system which is A*tl_v+tl_A*v=tl_r.
1072! Here, tl_A*v and tl_r are known, so we must solve for tl_v
1073! by inverting A*tl_v=tl_r-tl_A*v.
1074!
1075 cff=-lambda*dt(ng)/0.5_r8
1076 DO k=1,n(ng)-1
1077 DO i=istr,iend
1078 cff1=1.0_r8/(z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
1079 & z_r(i,j,k )-z_r(i,j-1,k ))
1080 tl_cff1=-cff1*cff1*(tl_z_r(i,j,k+1)+tl_z_r(i,j-1,k+1)- &
1081 & tl_z_r(i,j,k )-tl_z_r(i,j-1,k ))
1082 fc(i,k)=cff*cff1*ak(i,k)
1083 tl_fc(i,k)=cff*(tl_cff1*ak(i,k)+cff1*tl_ak(i,k))
1084 END DO
1085 END DO
1086 DO i=istr,iend
1087 fc(i,0)=0.0_r8
1088 tl_fc(i,0)=0.0_r8
1089 fc(i,n(ng))=0.0_r8
1090 tl_fc(i,n(ng))=0.0_r8
1091 END DO
1092!
1093! Solve the tangent linear tridiagonal system.
1094!
1095 DO k=1,n(ng)
1096 DO i=istr,iend
1097 bc(i,k)=hzk(i,k)-fc(i,k)-fc(i,k-1)
1098 tl_bc(i,k)=tl_hzk(i,k)-tl_fc(i,k)-tl_fc(i,k-1)
1099 END DO
1100 END DO
1101 DO k=2,n(ng)-1
1102 DO i=istr,iend
1103 dc(i,k)=tl_v(i,j,k,nnew)- &
1104 & (tl_fc(i,k-1)*v(i,j,k-1,nnew)+ &
1105 & tl_bc(i,k )*v(i,j,k ,nnew)+ &
1106 & tl_fc(i,k )*v(i,j,k+1,nnew))
1107 END DO
1108 END DO
1109 DO i=istr,iend
1110 dc(i,1)=tl_v(i,j,1,nnew)- &
1111 & (tl_bc(i,1)*v(i,j,1,nnew)+ &
1112 & tl_fc(i,1)*v(i,j,2,nnew))
1113 dc(i,n(ng))=tl_v(i,j,n(ng),nnew)- &
1114 & (tl_fc(i,n(ng)-1)*v(i,j,n(ng)-1,nnew)+ &
1115 & tl_bc(i,n(ng) )*v(i,j,n(ng) ,nnew))
1116 END DO
1117 DO i=istr,iend
1118 cff=1.0_r8/bc(i,1)
1119 cf(i,1)=cff*fc(i,1)
1120 dc(i,1)=cff*dc(i,1)
1121 END DO
1122 DO k=2,n(ng)-1
1123 DO i=istr,iend
1124 cff=1.0_r8/(bc(i,k)-fc(i,k-1)*cf(i,k-1))
1125 cf(i,k)=cff*fc(i,k)
1126 dc(i,k)=cff*(dc(i,k)-fc(i,k-1)*dc(i,k-1))
1127 END DO
1128 END DO
1129!
1130! Compute new solution by back substitution.
1131! (DC is a tangent linear variable here).
1132!
1133 DO i=istr,iend
1134# ifdef DIAGNOSTICS_UV
1135!! wrk(i,N(ng))=v(i,j,N(ng),nnew)*oHz(i,N(ng))
1136# endif
1137 dc(i,n(ng))=(dc(i,n(ng))-fc(i,n(ng)-1)*dc(i,n(ng)-1))/ &
1138 & (bc(i,n(ng))-fc(i,n(ng)-1)*cf(i,n(ng)-1))
1139!^ v(i,j,N(ng),nnew)=DC(i,N(ng))
1140!^
1141 tl_v(i,j,n(ng),nnew)=dc(i,n(ng))
1142# ifdef DIAGNOSTICS_UV
1143!! DiaV3wrk(i,j,N(ng),M3vvis)=DiaV3wrk(i,j,N(ng),M3vvis)+ &
1144!! & v(i,j,N(ng),nnew)-wrk(i,N(ng))
1145# endif
1146 END DO
1147 DO k=n(ng)-1,1,-1
1148 DO i=istr,iend
1149# ifdef DIAGNOSTICS_UV
1150!! wrk(i,k)=v(i,j,k,nnew)*oHz(i,k)
1151# endif
1152 dc(i,k)=dc(i,k)-cf(i,k)*dc(i,k+1)
1153!^ v(i,j,k,nnew)=DC(i,k)
1154!^
1155 tl_v(i,j,k,nnew)=dc(i,k)
1156# ifdef DIAGNOSTICS_UV
1157!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)+ &
1158!! & v(i,j,k,nnew)-wrk(i,k)
1159# endif
1160 END DO
1161 END DO
1162# endif
1163!
1164! Replace INTERIOR POINTS incorrect vertical mean with more accurate
1165! barotropic component, vbar=DV_avg1/(D*om_v). Recall that, D=CF(:,0).
1166!
1167 DO i=istr,iend
1168 cf(i,0)=hzk(i,1)
1169 tl_cf(i,0)=tl_hzk(i,1)
1170 dc(i,0)=v(i,j,1,nnew)*hzk(i,1)
1171 tl_dc(i,0)=tl_v(i,j,1,nnew)*hzk(i,1)+ &
1172 & v(i,j,1,nnew)*tl_hzk(i,1)
1173# ifdef NEARSHORE_MELLOR
1174 dcs(i,0)=v_stokes(i,j,1)*hzk(i,1)
1175 tl_dcs(i,0)=tl_v_stokes(i,j,1)*hzk(i,1)+ &
1176 & v_stokes(i,j,1)*tl_hzk(i,1)
1177# endif
1178# ifdef DIAGNOSTICS_UV
1179!! Dwrk(i,M2pgrd)=DiaV3wrk(i,j,1,M3pgrd)*Hzk(i,1)
1180!! Dwrk(i,M2bstr)=DiaV3wrk(i,j,1,M3vvis)*Hzk(i,1)
1181# ifdef UV_COR
1182!! Dwrk(i,M2fcor)=DiaV3wrk(i,j,1,M3fcor)*Hzk(i,1)
1183# endif
1184# if defined UV_VIS2 || defined UV_VIS4
1185!! Dwrk(i,M2xvis)=DiaV3wrk(i,j,1,M3xvis)*Hzk(i,1)
1186!! Dwrk(i,M2yvis)=DiaV3wrk(i,j,1,M3yvis)*Hzk(i,1)
1187!! Dwrk(i,M2hvis)=DiaV3wrk(i,j,1,M3hvis)*Hzk(i,1)
1188# endif
1189# ifdef UV_ADV
1190!! Dwrk(i,M2xadv)=DiaV3wrk(i,j,1,M3xadv)*Hzk(i,1)
1191!! Dwrk(i,M2yadv)=DiaV3wrk(i,j,1,M3yadv)*Hzk(i,1)
1192!! Dwrk(i,M2hadv)=DiaV3wrk(i,j,1,M3hadv)*Hzk(i,1)
1193# endif
1194# ifdef NEARSHORE_MELLOR
1195!! Dwrk(i,M2hrad)=DiaV3wrk(i,j,1,M3hrad)*Hzk(i,1)
1196# endif
1197# endif
1198 END DO
1199 DO k=2,n(ng)
1200 DO i=istr,iend
1201 cf(i,0)=cf(i,0)+hzk(i,k)
1202 tl_cf(i,0)=tl_cf(i,0)+tl_hzk(i,k)
1203 dc(i,0)=dc(i,0)+v(i,j,k,nnew)*hzk(i,k)
1204 tl_dc(i,0)=tl_dc(i,0)+ &
1205 & tl_v(i,j,k,nnew)*hzk(i,k)+ &
1206 & v(i,j,k,nnew)*tl_hzk(i,k)
1207# ifdef NEARSHORE_MELLOR
1208 dcs(i,0)=dcs(i,0)+v_stokes(i,j,k)*hzk(i,k)
1209 tl_dcs(i,0)=tl_dcs(i,0)+ &
1210 & tl_v_stokes(i,j,k)*hzk(i,k)+ &
1211 & v_stokes(i,j,k)*tl_hzk(i,k)
1212# endif
1213# ifdef DIAGNOSTICS_UV
1214!! Dwrk(i,M2pgrd)=Dwrk(i,M2pgrd)+ &
1215!! & DiaV3wrk(i,j,k,M3pgrd)*Hzk(i,k)
1216!! Dwrk(i,M2bstr)=Dwrk(i,M2bstr)+ &
1217!! & DiaV3wrk(i,j,k,M3vvis)*Hzk(i,k)
1218# ifdef UV_COR
1219!! Dwrk(i,M2fcor)=Dwrk(i,M2fcor)+ &
1220!! & DiaV3wrk(i,j,k,M3fcor)*Hzk(i,k)
1221# endif
1222# if defined UV_VIS2 || defined UV_VIS4
1223!! Dwrk(i,M2xvis)=Dwrk(i,M2xvis)+ &
1224!! & DiaV3wrk(i,j,k,M3xvis)*Hzk(i,k)
1225!! Dwrk(i,M2yvis)=Dwrk(i,M2yvis)+ &
1226!! & DiaV3wrk(i,j,k,M3yvis)*Hzk(i,k)
1227!! Dwrk(i,M2hvis)=Dwrk(i,M2hvis)+ &
1228!! & DiaV3wrk(i,j,k,M3hvis)*Hzk(i,k)
1229# endif
1230# ifdef UV_ADV
1231!! Dwrk(i,M2xadv)=Dwrk(i,M2xadv)+ &
1232!! & DiaV3wrk(i,j,k,M3xadv)*Hzk(i,k)
1233!! Dwrk(i,M2yadv)=Dwrk(i,M2yadv)+ &
1234!! & DiaV3wrk(i,j,k,M3yadv)*Hzk(i,k)
1235!! Dwrk(i,M2hadv)=Dwrk(i,M2hadv)+ &
1236!! & DiaV3wrk(i,j,k,M3hadv)*Hzk(i,k)
1237# endif
1238# ifdef NEARSHORE_MELLOR
1239!! Dwrk(i,M2hrad)=Dwrk(i,M2hrad)+ &
1240!! & DiaV3wrk(i,j,k,M3hrad)*Hzk(i,k)
1241# endif
1242# endif
1243 END DO
1244 END DO
1245 DO i=istr,iend
1246 dc1(i,0)=dc(i,0) ! intermediate
1247 cff1=1.0_r8/(cf(i,0)*om_v(i,j))
1248 tl_cff1=-cff1*cff1*tl_cf(i,0)*om_v(i,j)
1249 dc(i,0)=(dc(i,0)*om_v(i,j)-dv_avg1(i,j))*cff1 ! recursive
1250 tl_dc(i,0)=(tl_dc(i,0)*om_v(i,j)-tl_dv_avg1(i,j))*cff1+ &
1251 & (dc1(i,0)*om_v(i,j)-dv_avg1(i,j))*tl_cff1
1252# ifdef NEARSHORE_MELLOR
1253 dcs1(i)=dcs(i,0) ! intermediate
1254 cff2=1.0_r8/cf(i,0)
1255 tl_cff2=-cff2*cff2*tl_cf(i,0)
1256 dcs(i,0)=dcs(i,0)*cff2-vbar_stokes(i,j) ! recursive
1257 tl_dcs(i,0)=tl_dcs(i,0)*cff2+ &
1258 & dcs1(i,0)*tl_cff2-tl_vbar_stokes(i,j)
1259# endif
1260# ifdef DIAGNOSTICS_UV
1261!! DO idiag=1,M2pgrd
1262!! Dwrk(i,idiag)=(Dwrk(i,idiag)*om_v(i,j)- &
1263!! & DiaV2wrk(i,j,idiag))*cff1
1264!! END DO
1265!! Dwrk(i,M2bstr)=(Dwrk(i,M2bstr)*om_v(i,j)- &
1266!! & DiaV2wrk(i,j,M2bstr)-DiaV2wrk(i,j,M2sstr))* &
1267!! & cff1
1268# endif
1269 END DO
1270!
1271! Couple and update new solution.
1272!
1273 DO k=1,n(ng)
1274 DO i=istr,iend
1275!^ v(i,j,k,nnew)=v(i,j,k,nnew)-DC(i,0)
1276!^
1277 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)-tl_dc(i,0)
1278# ifdef MASKING
1279!^ v(i,j,k,nnew)=v(i,j,k,nnew)*vmask(i,j)
1280!^
1281 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*vmask(i,j)
1282# endif
1283# ifdef WET_DRY_NOT_YET
1284!^ v(i,j,k,nnew)=v(i,j,k,nnew)*vmask_wet(i,j)
1285!^
1286 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)*vmask_wet(i,j)
1287# endif
1288# ifdef NEARSHORE_MELLOR
1289!^ v_stokes(i,j,k)=v_stokes(i,j,k)-DCs(i,0)
1290!^
1291 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)-tl_dcs(i,0)
1292# ifdef MASKING
1293!^ v_stokes(i,j,k)=v_stokes(i,j,k)*vmask(i,j)
1294!^
1295 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*vmask(i,j)
1296# endif
1297# ifdef WET_DRY_NOT_YET
1298!^ v_stokes(i,j,k)=v_stokes(i,j,k)*vmask_wet(i,j)
1299!^
1300 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)*vmask_wet(i,j)
1301# endif
1302# endif
1303# ifdef DIAGNOSTICS_UV
1304!! DiaV3wrk(i,j,k,M3pgrd)=DiaV3wrk(i,j,k,M3pgrd)- &
1305!! & Dwrk(i,M2pgrd)
1306!! DiaV3wrk(i,j,k,M3vvis)=DiaV3wrk(i,j,k,M3vvis)- &
1307!! & Dwrk(i,M2bstr)
1308# ifdef UV_COR
1309!! DiaV3wrk(i,j,k,M3fcor)=DiaV3wrk(i,j,k,M3fcor)- &
1310!! & Dwrk(i,M2fcor)
1311# endif
1312# if defined UV_VIS2 || defined UV_VIS4
1313!! DiaV3wrk(i,j,k,M3xvis)=DiaV3wrk(i,j,k,M3xvis)- &
1314!! & Dwrk(i,M2xvis)
1315!! DiaV3wrk(i,j,k,M3yvis)=DiaV3wrk(i,j,k,M3yvis)- &
1316!! & Dwrk(i,M2yvis)
1317!! DiaV3wrk(i,j,k,M3hvis)=DiaV3wrk(i,j,k,M3hvis)- &
1318!! & Dwrk(i,M2hvis)
1319# endif
1320# ifdef UV_ADV
1321!! DiaV3wrk(i,j,k,M3xadv)=DiaV3wrk(i,j,k,M3xadv)- &
1322!! & Dwrk(i,M2xadv)
1323!! DiaV3wrk(i,j,k,M3yadv)=DiaV3wrk(i,j,k,M3yadv)- &
1324!! & Dwrk(i,M2yadv)
1325!! DiaV3wrk(i,j,k,M3hadv)=DiaV3wrk(i,j,k,M3hadv)- &
1326!! & Dwrk(i,M2hadv)
1327# endif
1328# ifdef NEARSHORE_MELLOR
1329!! DiaV3wrk(i,j,k,M3hrad)=DiaV3wrk(i,j,k,M3hrad)- &
1330!! & Dwrk(i,M2hrad)
1331# endif
1332# endif
1333 END DO
1334 END DO
1335
1336# if defined DIAGNOSTICS_UV && defined MASKING
1337!! DO k=1,N(ng)
1338!! DO i=Istr,Iend
1339!! DO idiag=1,NDM3d
1340!! DiaV3wrk(i,j,k,idiag)=DiaV3wrk(i,j,k,idiag)*vmask(i,j)
1341!! END DO
1342!! END DO
1343!! END DO
1344# endif
1345 END IF
1346 END DO
1347!
1348!-----------------------------------------------------------------------
1349! Set lateral boundary conditions.
1350!-----------------------------------------------------------------------
1351!
1352!^ CALL u3dbc_tile (ng, tile, &
1353!^ & LBi, UBi, LBj, UBj, N(ng), &
1354!^ & IminS, ImaxS, JminS, JmaxS, &
1355!^ & nstp, nnew, &
1356!^ & u)
1357!^
1358 CALL tl_u3dbc_tile (ng, tile, &
1359 & lbi, ubi, lbj, ubj, n(ng), &
1360 & imins, imaxs, jmins, jmaxs, &
1361 & nstp, nnew, &
1362 & tl_u)
1363!^ CALL v3dbc_tile (ng, tile, &
1364!^ & LBi, UBi, LBj, UBj, N(ng), &
1365!^ & IminS, ImaxS, JminS, JmaxS, &
1366!^ & nstp, nnew, &
1367!^ & v)
1368!^
1369 CALL tl_v3dbc_tile (ng, tile, &
1370 & lbi, ubi, lbj, ubj, n(ng), &
1371 & imins, imaxs, jmins, jmaxs, &
1372 & nstp, nnew, &
1373 & tl_v)
1374!
1375!-----------------------------------------------------------------------
1376! Apply momentum transport point sources (like river runoff), if any.
1377!-----------------------------------------------------------------------
1378!
1379 IF (luvsrc(ng)) THEN
1380 DO is=1,nsrc(ng)
1381 i=sources(ng)%Isrc(is)
1382 j=sources(ng)%Jsrc(is)
1383 IF (((istrr.le.i).and.(i.le.iendr)).and. &
1384 & ((jstrr.le.j).and.(j.le.jendr))) THEN
1385 IF (int(sources(ng)%Dsrc(is)).eq.0) THEN
1386 DO k=1,n(ng)
1387 cff1=1.0_r8/(on_u(i,j)* &
1388 & 0.5_r8*(z_w(i-1,j,k)-z_w(i-1,j,k-1)+ &
1389 & z_w(i ,j,k)-z_w(i ,j,k-1)))
1390 tl_cff1=-cff1*cff1*on_u(i,j)* &
1391 & 0.5_r8*(tl_z_w(i-1,j,k)-tl_z_w(i-1,j,k-1)+ &
1392 & tl_z_w(i ,j,k)-tl_z_w(i ,j,k-1))
1393!^ u(i,j,k,nnew)=SOURCES(ng)%Qsrc(is,k)*cff1
1394!^
1395 tl_u(i,j,k,nnew)=sources(ng)%tl_Qsrc(is,k)*cff1+ &
1396 & sources(ng)%Qsrc(is,k)*tl_cff1
1397 END DO
1398 ELSE IF (int(sources(ng)%Dsrc(is)).eq.1) THEN
1399 DO k=1,n(ng)
1400 cff1=1.0_r8/(om_v(i,j)* &
1401 & 0.5_r8*(z_w(i,j-1,k)-z_w(i,j-1,k-1)+ &
1402 & z_w(i,j ,k)-z_w(i,j ,k-1)))
1403 tl_cff1=-cff1*cff1*om_v(i,j)* &
1404 & 0.5_r8*(tl_z_w(i,j-1,k)-tl_z_w(i,j-1,k-1)+ &
1405 & tl_z_w(i,j ,k)-tl_z_w(i,j ,k-1))
1406!^ v(i,j,k,nnew)=SOURCES(ng)%Qsrc(is,k)*cff1
1407!^
1408 tl_v(i,j,k,nnew)=sources(ng)%tl_Qsrc(is,k)*cff1+ &
1409 & sources(ng)%Qsrc(is,k)*tl_cff1
1410 END DO
1411 END IF
1412 END IF
1413 END DO
1414 END IF
1415!
1416!-----------------------------------------------------------------------
1417! Couple 2D and 3D momentum equations.
1418!-----------------------------------------------------------------------
1419!
1420! Couple velocity component in the XI-direction.
1421!
1422 DO j=jstrt,jendt
1423 DO i=istrp,iendt
1424 dc(i,0)=0.0_r8
1425 tl_dc(i,0)=0.0_r8
1426 cf(i,0)=0.0_r8
1427 tl_cf(i,0)=0.0_r8
1428# ifdef NEARSHORE_MELLOR
1429 cfs(i,0)=0.0_r8
1430 tl_cfs(i,0)=0.0_r8
1431# endif
1432 fc(i,0)=0.0_r8
1433 tl_fc(i,0)=0.0_r8
1434 END DO
1435!
1436! Compute thicknesses of U-boxes DC(i,1:N), total depth of the water
1437! column DC(i,0), and incorrect vertical mean CF(i,0). Notice that
1438! barotropic component is replaced with its fast-time averaged
1439! values.
1440!
1441 DO k=1,n(ng)
1442 DO i=istrp,iendt
1443 cff=0.5_r8*on_u(i,j)
1444 dc(i,k)=cff*(hz(i,j,k)+hz(i-1,j,k))
1445 tl_dc(i,k)=cff*(tl_hz(i,j,k)+tl_hz(i-1,j,k))
1446 dc(i,0)=dc(i,0)+dc(i,k)
1447 tl_dc(i,0)=tl_dc(i,0)+tl_dc(i,k)
1448 cf(i,0)=cf(i,0)+ &
1449 & dc(i,k)*u(i,j,k,nnew)
1450 tl_cf(i,0)=tl_cf(i,0)+ &
1451 & tl_dc(i,k)*u(i,j,k,nnew)+ &
1452 & dc(i,k)*tl_u(i,j,k,nnew)
1453# ifdef NEARSHORE_MELLOR
1454 cfs(i,0)=cfs(i,0)+ &
1455 & dc(i,k)*u_stokes(i,j,k)
1456 tl_cfs(i,0)=tl_cfs(i,0)+ &
1457 & tl_dc(i,k)*u_stokes(i,j,k)+ &
1458 & dc(i,k)*tl_u_stokes(i,j,k)
1459# endif
1460 END DO
1461 END DO
1462 DO i=istrp,iendt
1463 dc1(i,0)=dc(i,0) ! intermediate
1464# ifdef NEARSHORE_MELLOR
1465 cff2=dc(i,0)*ubar_stokes(i,j)
1466 tl_cff2=tl_dc(i,0)*ubar_stokes(i,j)+ &
1467 & dc(i,0)*tl_ubar_stokes(i,j)
1468# endif
1469 dc(i,0)=1.0_r8/dc(i,0) ! recursive
1470 tl_dc(i,0)=-tl_dc(i,0)/(dc1(i,0)*dc1(i,0))
1471 cf1(i)=cf(i,0) ! intermediate
1472 cf(i,0)=dc(i,0)*(cf(i,0)-du_avg1(i,j)) ! recursive
1473 tl_cf(i,0)=tl_dc(i,0)*(cf1(i)-du_avg1(i,j))+ &
1474 & dc(i,0)*(tl_cf(i,0)-tl_du_avg1(i,j))
1475# ifdef NEARSHORE_MELLOR
1476 cfs1(i)=cfs(i,0) ! intermediate
1477 cfs(i,0)=dc(i,0)*(cfs(i,0)-cff2) ! recursive
1478 tl_cfs(i,0)=tl_dc(i,0)*(cfs1(i)-cff2)+ &
1479 & dc(i,0)*(tl_cfs(i,0)-tl_cff2)
1480# endif
1481!^ ubar(i,j,1)=DC(i,0)*DU_avg1(i,j)
1482!^
1483 tl_ubar(i,j,1)=tl_dc(i,0)*du_avg1(i,j)+ &
1484 & dc(i,0)*tl_du_avg1(i,j)
1485!^ ubar(i,j,2)=ubar(i,j,1)
1486!^
1487 tl_ubar(i,j,2)=tl_ubar(i,j,1)
1488# ifdef DIAGNOSTICS_UV
1489!! DiaU2wrk(i,j,M2rate)=ubar(i,j,1)-DiaU2int(i,j,M2rate)*DC(i,0)
1490!! DiaU2int(i,j,M2rate)=ubar(i,j,1)*DC1(i,0)
1491# endif
1492 END DO
1493# ifdef DIAGNOSTICS_UV
1494!!
1495!! Convert the units of the fast-time integrated change in mass flux
1496!! diagnostic values to change in velocity (m s-1).
1497!!
1498!! DO idiag=1,NDM2d-1
1499!! DO i=IstrP,IendT
1500!! DiaU2wrk(i,j,idiag)=DC(i,0)*DiaU2wrk(i,j,idiag)
1501# ifdef MASKING
1502!! DiaU2wrk(i,j,idiag)=DiaU2wrk(i,j,idiag)*umask(i,j)
1503# endif
1504!! END DO
1505!! END DO
1506# endif
1507!
1508! Replace only BOUNDARY POINTS incorrect vertical mean with more
1509! accurate barotropic component, ubar=DU_avg1/(D*on_u). Recall that,
1510! D=CF(:,0).
1511!
1512! NOTE: Only the BOUNDARY POINTS need to be replaced. Avoid redundant
1513! update in the interior again for computational purposes which
1514! will not affect the nonlinear code. However, the adjoint
1515! code is wrong because the interior solution is corrected
1516! twice. The replacement is avoided if the boundary edge is
1517! periodic. The J-loop is pipelined, so we need to do a special
1518! test for the southern and northern domain edges.
1519!
1520 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1521 IF (domain(ng)%Western_Edge(tile)) THEN
1522 DO k=1,n(ng)
1523!^ u(Istr,j,k,nnew)=u(Istr,j,k,nnew)-CF(Istr,0)
1524!^
1525 tl_u(istr,j,k,nnew)=tl_u(istr,j,k,nnew)- &
1526 & tl_cf(istr,0)
1527# ifdef MASKING
1528!^ u(Istr,j,k,nnew)=u(Istr,j,k,nnew)* &
1529!^ & umask(Istr,j)
1530!^
1531 tl_u(istr,j,k,nnew)=tl_u(istr,j,k,nnew)* &
1532 & umask(istr,j)
1533# endif
1534# ifdef WET_DRY_NOT_YET
1535!^ u(Istr,j,k,nnew)=u(Istr,j,k,nnew)* &
1536!^ & umask_wet(Istr,j)
1537!^
1538 tl_u(istr,j,k,nnew)=tl_u(istr,j,k,nnew)* &
1539 & umask_wet(istr,j)
1540# endif
1541# ifdef NEARSHORE_MELLOR
1542!^ u_stokes(Istr,j,k)=u_stokes(Istr,j,k)-CFs(Istr,0)
1543!^
1544 tl_u_stokes(istr,j,k)=tl_u_stokes(istr,j,k)- &
1545 & tl_cfs(istr,0)
1546# ifdef MASKING
1547!^ u_stokes(Istr,j,k)=u_stokes(Istr,j,k)* &
1548!^ & umask(Istr,j)
1549!^
1550 tl_u_stokes(istr,j,k)=tl_u_stokes(istr,j,k)* &
1551 & umask(istr,j)
1552# endif
1553# ifdef WET_DRY_NOT_YET
1554!^ u_stokes(Istr,j,k)=u_stokes(Istr,j,k)* &
1555!^ & umask_wet(Istr,j)
1556!^
1557 tl_u_stokes(istr,j,k)=tl_u_stokes(istr,j,k)* &
1558 & umask_wet(istr,j)
1559# endif
1560# endif
1561 END DO
1562 END IF
1563 END IF
1564!
1565 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1566 IF (domain(ng)%Eastern_Edge(tile)) THEN
1567 DO k=1,n(ng)
1568!^ u(Iend+1,j,k,nnew)=u(Iend+1,j,k,nnew)-CF(Iend+1,0)
1569!^
1570 tl_u(iend+1,j,k,nnew)=tl_u(iend+1,j,k,nnew)- &
1571 & tl_cf(iend+1,0)
1572# ifdef MASKING
1573!^ u(Iend+1,j,k,nnew)=u(Iend+1,j,k,nnew)* &
1574!^ & umask(Iend+1,j)
1575!^
1576 tl_u(iend+1,j,k,nnew)=tl_u(iend+1,j,k,nnew)* &
1577 & umask(iend+1,j)
1578# endif
1579# ifdef WET_DRY_NOT_YET
1580!^ u(Iend+1,j,k,nnew)=u(Iend+1,j,k,nnew)* &
1581!^ & umask_wet(Iend+1,j)
1582!^
1583 tl_u(iend+1,j,k,nnew)=tl_u(iend+1,j,k,nnew)* &
1584 & umask_wet(iend+1,j)
1585# endif
1586# ifdef NEARSHORE_MELLOR
1587!^ u_stokes(Iend+1,j,k)=u_stokes(Iend+1,j,k)-CFs(Iend+1,0)
1588!^
1589 tl_u_stokes(iend+1,j,k)=tl_u_stokes(iend+1,j,k)- &
1590 & tl_cfs(iend+1,0)
1591# ifdef MASKING
1592!^ u_stokes(Iend+1,j,k)=u_stokes(Iend+1,j,k)* &
1593!^ & umask(Iend+1,j)
1594!^
1595 tl_u_stokes(iend+1,j,k)=tl_u_stokes(iend+1,j,k)* &
1596 & umask(iend+1,j)
1597# endif
1598# ifdef WET_DRY_NOT_YET
1599!^ u_stokes(Iend+1,j,k)=u_stokes(Iend+1,j,k)* &
1600!^ & umask_wet(Iend+1,j)
1601!^
1602 tl_u_stokes(iend+1,j,k)=tl_u_stokes(iend+1,j,k)* &
1603 & umask_wet(iend+1,j)
1604# endif
1605# endif
1606 END DO
1607 END IF
1608 END IF
1609!
1610 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1611 IF (j.eq.0) THEN ! southern boundary
1612 DO k=1,n(ng) ! J-loop pipelined
1613 DO i=istru,iend
1614!^ u(i,j,k,nnew)=u(i,j,k,nnew)-CF(i,0)
1615!^
1616 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)- &
1617 & tl_cf(i,0)
1618# ifdef MASKING
1619!^ u(i,j,k,nnew)=u(i,j,k,nnew)*umask(i,j)
1620!^
1621 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)* &
1622 & umask(i,j)
1623# endif
1624# ifdef WET_DRY_NOT_YET
1625!^ u(i,j,k,nnew)=u(i,j,k,nnew)* &
1626!^ & umask_wet(i,j)
1627!^
1628 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)* &
1629 & umask_wet(i,j)
1630# endif
1631# ifdef NEARSHORE_MELLOR
1632!^ u_stokes(i,j,k)=u_stokes(i,j,k)-CFs(i,0)
1633!^
1634 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)- &
1635 & tl_cfs(i,0)
1636# ifdef MASKING
1637!^ u_stokes(i,j,k)=u_stokes(i,j,k)* &
1638!^ & umask(i,j)
1639!^
1640 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)* &
1641 & umask(i,j)
1642# endif
1643# ifdef WET_DRY_NOT_YET
1644!^ u_stokes(i,j,k)=u_stokes(i,j,k)* &
1645!^ & umask_wet(i,j)
1646!^
1647 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)* &
1648 & umask_wet(i,j)
1649# endif
1650# endif
1651 END DO
1652 END DO
1653 END IF
1654 END IF
1655!
1656 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1657 IF (j.eq.mm(ng)+1) THEN ! northern boundary
1658 DO k=1,n(ng) ! J-loop pipelined
1659 DO i=istru,iend
1660!^ u(i,j,k,nnew)=u(i,j,k,nnew)-CF(i,0)
1661!^
1662 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)- &
1663 & tl_cf(i,0)
1664# ifdef MASKING
1665!^ u(i,j,k,nnew)=u(i,j,k,nnew)*umask(i,j)
1666!^
1667 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)* &
1668 & umask(i,j)
1669# endif
1670# ifdef WET_DRY_NOT_YET
1671!^ u(i,j,k,nnew)=u(i,j,k,nnew)*umask_wet(i,j)
1672!^
1673 tl_u(i,j,k,nnew)=tl_u(i,j,k,nnew)* &
1674 & umask_wet(i,j)
1675# endif
1676# ifdef NEARSHORE_MELLOR
1677!^ u_stokes(i,j,k)=u_stokes(i,j,k)-CFs(i,0)
1678!^
1679 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)- &
1680 & tl_cfs(i,0)
1681# ifdef MASKING
1682!^ u_stokes(i,j,k)=u_stokes(i,j,k)* &
1683!^ & umask(i,j)
1684!^
1685 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)* &
1686 & umask(i,j)
1687# endif
1688# ifdef WET_DRY_NOT_YET
1689!^ u_stokes(i,j,k)=u_stokes(i,j,k)* &
1690!^ & umask_wet(i,j)
1691!^
1692 tl_u_stokes(i,j,k)=tl_u_stokes(i,j,k)* &
1693 & umask_wet(i,j)
1694# endif
1695# endif
1696 END DO
1697 END DO
1698 END IF
1699 END IF
1700!
1701! Compute correct mass flux, Hz*u/n.
1702!
1703 DO k=n(ng),1,-1
1704 DO i=istrp,iendt
1705 huon(i,j,k)=0.5_r8*(huon(i,j,k)+u(i,j,k,nnew)*dc(i,k))
1706 tl_huon(i,j,k)=0.5_r8*(tl_huon(i,j,k)+ &
1707 & tl_u(i,j,k,nnew)*dc(i,k)+ &
1708 & u(i,j,k,nnew)*tl_dc(i,k))
1709# ifdef NEARSHORE_MELLOR
1710 huon(i,j,k)=huon(i,j,k)+0.5_r8*u_stokes(i,j,k)*dc(i,k)
1711 tl_huon(i,j,k)=tl_huon(i,j,k)+ &
1712 & 0.5_r8*(tl_u_stokes(i,j,k)*dc(i,k)+ &
1713 & u_stokes(i,j,k)*tl_dc(i,k))
1714# endif
1715 fc(i,0)=fc(i,0)+huon(i,j,k)
1716 tl_fc(i,0)=tl_fc(i,0)+tl_huon(i,j,k)
1717# ifdef DIAGNOSTICS_UV
1718!! DiaU3wrk(i,j,k,M3rate)=u(i,j,k,nnew)-DiaU3wrk(i,j,k,M3rate)
1719# endif
1720 END DO
1721 END DO
1722 DO i=istrp,iendt
1723 fc1(i)=fc(i,0) ! intermediate
1724 fc(i,0)=dc(i,0)*(fc(i,0)-du_avg2(i,j)) ! recursive
1725 tl_fc(i,0)=tl_dc(i,0)*(fc1(i)-du_avg2(i,j))+ &
1726 & dc(i,0)*(tl_fc(i,0)-tl_du_avg2(i,j))
1727 END DO
1728 DO k=1,n(ng)
1729 DO i=istrp,iendt
1730 huon(i,j,k)=huon(i,j,k)-dc(i,k)*fc(i,0)
1731 tl_huon(i,j,k)=tl_huon(i,j,k)- &
1732 & tl_dc(i,k)*fc(i,0)- &
1733 & dc(i,k)*tl_fc(i,0)
1734 END DO
1735 END DO
1736!
1737! Couple velocity component in the ETA-direction.
1738!
1739 IF (j.ge.jstr) THEN
1740 DO i=istrt,iendt
1741 dc(i,0)=0.0_r8
1742 tl_dc(i,0)=0.0_r8
1743 cf(i,0)=0.0_r8
1744 tl_cf(i,0)=0.0_r8
1745# ifdef NEARSHORE_MELLOR
1746 cfs(i,0)=0.0_r8
1747 tl_cfs(i,0)=0.0_r8
1748# endif
1749 fc(i,0)=0.0_r8
1750 tl_fc(i,0)=0.0_r8
1751 END DO
1752!
1753! Compute thicknesses of V-boxes DC(i,1:N), total depth of the water
1754! column DC(i,0), and incorrect vertical mean CF(i,0). Notice that
1755! barotropic component is replaced with its fast-time averaged
1756! values.
1757!
1758 DO k=1,n(ng)
1759 DO i=istrt,iendt
1760 cff=0.5_r8*om_v(i,j)
1761 dc(i,k)=cff*(hz(i,j,k)+hz(i,j-1,k))
1762 tl_dc(i,k)=cff*(tl_hz(i,j,k)+tl_hz(i,j-1,k))
1763 dc(i,0)=dc(i,0)+dc(i,k)
1764 tl_dc(i,0)=tl_dc(i,0)+tl_dc(i,k)
1765 cf(i,0)=cf(i,0)+ &
1766 & dc(i,k)*v(i,j,k,nnew)
1767 tl_cf(i,0)=tl_cf(i,0)+ &
1768 & tl_dc(i,k)*v(i,j,k,nnew)+ &
1769 & dc(i,k)*tl_v(i,j,k,nnew)
1770# ifdef NEARSHORE_MELLOR
1771 cfs(i,0)=cfs(i,0)+ &
1772 & dc(i,k)*v_stokes(i,j,k)
1773 tl_cfs(i,0)=tl_cfs(i,0)+ &
1774 & tl_dc(i,k)*v_stokes(i,j,k)+ &
1775 & dc(i,k)*tl_v_stokes(i,j,k)
1776# endif
1777 END DO
1778 END DO
1779 DO i=istrt,iendt
1780 dc1(i,0)=dc(i,0) ! intermediate
1781# ifdef NEARSHORE_MELLOR
1782 cff2=dc(i,0)*vbar_stokes(i,j)
1783 tl_cff2=tl_dc(i,0)*vbar_stokes(i,j)+ &
1784 & dc(i,0)*tl_vbar_stokes(i,j)
1785# endif
1786 dc(i,0)=1.0_r8/dc(i,0) ! recursive
1787 tl_dc(i,0)=-tl_dc(i,0)/(dc1(i,0)*dc1(i,0))
1788 cf1(i)=cf(i,0) ! intermediate
1789 cf(i,0)=dc(i,0)*(cf(i,0)-dv_avg1(i,j)) ! recursive
1790 tl_cf(i,0)=tl_dc(i,0)*(cf1(i)-dv_avg1(i,j))+ &
1791 & dc(i,0)*(tl_cf(i,0)-tl_dv_avg1(i,j))
1792# ifdef NEARSHORE_MELLOR
1793 cfs1(i)=cfs(i,0)
1794 cfs(i,0)=dc(i,0)*(cfs(i,0)-cff2) ! recursive
1795 tl_cfs(i,0)=tl_dc(i,0)*(cfs1(i)-cff2)+ &
1796 & dc(i,0)*(tl_cfs(i,0)-tl_cff2)
1797# endif
1798!^ vbar(i,j,1)=DC(i,0)*DV_avg1(i,j)
1799!^
1800 tl_vbar(i,j,1)=tl_dc(i,0)*dv_avg1(i,j)+ &
1801 & dc(i,0)*tl_dv_avg1(i,j)
1802!^ vbar(i,j,2)=vbar(i,j,1)
1803!^
1804 tl_vbar(i,j,2)=tl_vbar(i,j,1)
1805# ifdef DIAGNOSTICS_UV
1806!! DiaV2wrk(i,j,M2rate)=vbar(i,j,1)- &
1807!! & DiaV2int(i,j,M2rate)*DC(i,0)
1808!! DiaV2int(i,j,M2rate)=vbar(i,j,1)*DC1(i,0)
1809!! DiaV2wrk(i,j,M2rate)=vbar(i,j,1)-DiaV2int(i,j,M2rate)
1810!! DiaV2int(i,j,M2rate)=vbar(i,j,1)
1811# endif
1812 END DO
1813# ifdef DIAGNOSTICS_UV
1814!!
1815!! Convert the units of the fast-time integrated change in mass flux
1816!! diagnostic values to change in velocity (m s-1).
1817!!
1818!! DO idiag=1,NDM2d-1
1819!! DO i=IstrT,IendT
1820!! DiaV2wrk(i,j,idiag)=DC(i,0)*DiaV2wrk(i,j,idiag)
1821# ifdef MASKING
1822!! DiaV2wrk(i,j,idiag)=DiaV2wrk(i,j,idiag)*vmask(i,j)
1823# endif
1824!! END DO
1825!! END DO
1826# endif
1827!
1828! Replace only BOUNDARY POINTS incorrect vertical mean with more
1829! accurate barotropic component, vbar=DV_avg1/(D*om_v). Recall that,
1830! D=CF(:,0).
1831!
1832! NOTE: Only the BOUNDARY POINTS need to be replaced. Avoid redundant
1833! update in the interior again for computational purposes which
1834! will not affect the nonlinear code. However, the adjoint
1835! code is wrong because the interior solution is corrected
1836! twice. The replacement is avoided if the boundary edge is
1837! periodic. The J-loop is pipelined, so we need to do a special
1838! test for the southern and northern domain edges.
1839!
1840 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1841 IF (domain(ng)%Western_Edge(tile)) THEN
1842 DO k=1,n(ng)
1843!^ v(Istr-1,j,k,nnew)=v(Istr-1,j,k,nnew)-CF(Istr-1,0)
1844!^
1845 tl_v(istr-1,j,k,nnew)=tl_v(istr-1,j,k,nnew)- &
1846 & tl_cf(istr-1,0)
1847# ifdef MASKING
1848!^ v(Istr-1,j,k,nnew)=v(Istr-1,j,k,nnew)* &
1849!^ & vmask(Istr-1,j)
1850!^
1851 tl_v(istr-1,j,k,nnew)=tl_v(istr-1,j,k,nnew)* &
1852 & vmask(istr-1,j)
1853# endif
1854# ifdef WET_DRY_NOT_YET
1855!^ v(Istr-1,j,k,nnew)=v(Istr-1,j,k,nnew)* &
1856!^ & vmask_wet(Istr-1,j)
1857!^
1858 tl_v(istr-1,j,k,nnew)=tl_v(istr-1,j,k,nnew)* &
1859 & vmask_wet(istr-1,j)
1860# endif
1861# ifdef NEARSHORE_MELLOR
1862!^ v_stokes(Istr-1,j,k)=v_stokes(Istr-1,j,k)- &
1863!^ & CFs(Istr-1,0)
1864!^
1865 tl_v_stokes(istr-1,j,k)=tl_v_stokes(istr-1,j,k)- &
1866 & tl_cfs(istr-1,0)
1867# ifdef MASKING
1868!^ v_stokes(Istr-1,j,k)=v_stokes(Istr-1,j,k)* &
1869!^ & vmask(Istr-1,j)
1870!^
1871 tl_v_stokes(istr-1,j,k)=tl_v_stokes(istr-1,j,k)* &
1872 & vmask(istr-1,j)
1873# endif
1874# ifdef WET_DRY_NOT_YET
1875!^ v_stokes(Istr-1,j,k)=v_stokes(Istr-1,j,k)* &
1876!^ & vmask_wet(Istr-1,j)
1877!^
1878 tl_v_stokes(istr-1,j,k)=tl_v_stokes(istr-1,j,k)* &
1879 & vmask_wet(istr-1,j)
1880# endif
1881# endif
1882 END DO
1883 END IF
1884 END IF
1885!
1886 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1887 IF (domain(ng)%Eastern_Edge(tile)) THEN
1888 DO k=1,n(ng)
1889!^ v(Iend+1,j,k,nnew)=v(Iend+1,j,k,nnew)-CF(Iend+1,0)
1890!^
1891 tl_v(iend+1,j,k,nnew)=tl_v(iend+1,j,k,nnew)- &
1892 & tl_cf(iend+1,0)
1893# ifdef MASKING
1894!^ v(Iend+1,j,k,nnew)=v(Iend+1,j,k,nnew)* &
1895!^ & vmask(Iend+1,j)
1896!^
1897 tl_v(iend+1,j,k,nnew)=tl_v(iend+1,j,k,nnew)* &
1898 & vmask(iend+1,j)
1899# endif
1900# ifdef WET_DRY_NOT_YET
1901!^ v(Iend+1,j,k,nnew)=v(Iend+1,j,k,nnew)* &
1902!^ & vmask_wet(Iend+1,j)
1903!^
1904 tl_v(iend+1,j,k,nnew)=tl_v(iend+1,j,k,nnew)* &
1905 & vmask_wet(iend+1,j)
1906# endif
1907# ifdef NEARSHORE_MELLOR
1908!^ v_stokes(Iend+1,j,k)=v_stokes(Iend+1,j,k)- &
1909!^ & CFs(Iend+1,0)
1910!^
1911 tl_v_stokes(iend+1,j,k)=tl_v_stokes(iend+1,j,k)- &
1912 & tl_cfs(iend+1,0)
1913# ifdef MASKING
1914!^ v_stokes(Iend+1,j,k)=v_stokes(Iend+1,j,k)* &
1915!^ & vmask(Iend+1,j)
1916!^
1917 tl_v_stokes(iend+1,j,k)=tl_v_stokes(iend+1,j,k)* &
1918 & vmask(iend+1,j)
1919# endif
1920# ifdef WET_DRY_NOT_YET
1921!^ v_stokes(Iend+1,j,k)=v_stokes(Iend+1,j,k)* &
1922!^ & vmask_wet(Iend+1,j)
1923!^
1924 tl_v_stokes(iend+1,j,k)=tl_v_stokes(iend+1,j,k)* &
1925 & vmask_wet(iend+1,j)
1926# endif
1927# endif
1928 END DO
1929 END IF
1930 END IF
1931!
1932 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1933 IF (j.eq.1) THEN ! southern boundary
1934 DO k=1,n(ng) ! J-loop pipelined
1935 DO i=istr,iend
1936!^ v(i,j,k,nnew)=v(i,j,k,nnew)-CF(i,0)
1937!^
1938 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)- &
1939 & tl_cf(i,0)
1940# ifdef MASKING
1941!^ v(i,j,k,nnew)=v(i,j,k,nnew)*
1942!^ & vmask(i,j)
1943!^
1944 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)* &
1945 & vmask(i,j)
1946# endif
1947# ifdef WET_DRY_NOT_YET
1948!^ v(i,j,k,nnew)=v(i,j,k,nnew)* &
1949!^ & vmask_wet(i,j)
1950!^
1951 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)* &
1952 & vmask_wet(i,j)
1953# endif
1954# ifdef NEARSHORE_MELLOR
1955!^ v_stokes(i,j,k)=v_stokes(i,j,k)-CFs(i,0)
1956!^
1957 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)- &
1958 & tl_cfs(i,0)
1959# ifdef MASKING
1960!^ v_stokes(i,j,k)=v_stokes(i,j,k)*vmask(i,j)
1961!^
1962 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)* &
1963 & vmask(i,j)
1964# endif
1965# ifdef WET_DRY_NOT_YET
1966!^ v_stokes(i,j,k)=v_stokes(i,j,k)* &
1967!^ & vmask_wet(i,j)
1968!^
1969 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)* &
1970 & vmask_wet(i,j)
1971# endif
1972# endif
1973 END DO
1974 END DO
1975 END IF
1976 END IF
1977!
1978 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1979 IF (j.eq.mm(ng)+1) THEN ! northern boundary
1980 DO k=1,n(ng) ! J-loop pipelined
1981 DO i=istr,iend
1982!^ v(i,j,k,nnew)=v(i,j,k,nnew)-CF(i,0)
1983!^
1984 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)- &
1985 & tl_cf(i,0)
1986# ifdef MASKING
1987!^ v(i,j,k,nnew)=v(i,j,k,nnew)* &
1988!^ & vmask(i,j)
1989!^
1990 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)* &
1991 & vmask(i,j)
1992# endif
1993# ifdef WET_DRY_NOT_YET
1994!^ v(i,j,k,nnew)=v(i,j,k,nnew)* &
1995!^ & vmask_wet(i,j)
1996!^
1997 tl_v(i,j,k,nnew)=tl_v(i,j,k,nnew)* &
1998 & vmask_wet(i,j)
1999# endif
2000# ifdef NEARSHORE_MELLOR
2001!^ v_stokes(i,j,k)=v_stokes(i,j,k)-CFs(i,0)
2002!^
2003 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)- &
2004 & tl_cfs(i,0)
2005# ifdef MASKING
2006!^ v_stokes(i,j,k)=v_stokes(i,j,k)* &
2007!^ & vmask(i,j)
2008!^
2009 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)* &
2010 & vmask(i,j)
2011# endif
2012# ifdef WET_DRY_NOT_YET
2013!^ v_stokes(i,j,k)=v_stokes(i,j,k)* &
2014!^ & vmask_wet(i,j)
2015!^
2016 tl_v_stokes(i,j,k)=tl_v_stokes(i,j,k)* &
2017 & vmask_wet(i,j)
2018# endif
2019# endif
2020 END DO
2021 END DO
2022 END IF
2023 END IF
2024!
2025! Compute correct mass flux, Hz*v/m.
2026!
2027 DO k=n(ng),1,-1
2028 DO i=istrt,iendt
2029 hvom(i,j,k)=0.5_r8*(hvom(i,j,k)+v(i,j,k,nnew)*dc(i,k))
2030 tl_hvom(i,j,k)=0.5_r8*(tl_hvom(i,j,k)+ &
2031 & tl_v(i,j,k,nnew)*dc(i,k)+ &
2032 & v(i,j,k,nnew)*tl_dc(i,k))
2033# ifdef NEARSHORE_MELLOR
2034 hvom(i,j,k)=hvom(i,j,k)+0.5_r8*v_stokes(i,j,k)*dc(i,k)
2035 tl_hvom(i,j,k)=tl_hvom(i,j,k)+ &
2036 & 0.5_r8*(tl_v_stokes(i,j,k)*dc(i,k)+ &
2037 & v_stokes(i,j,k)*tl_dc(i,k))
2038# endif
2039 fc(i,0)=fc(i,0)+hvom(i,j,k)
2040 tl_fc(i,0)=tl_fc(i,0)+tl_hvom(i,j,k)
2041# ifdef DIAGNOSTICS_UV
2042!! DiaV3wrk(i,j,k,M3rate)=v(i,j,k,nnew)- &
2043!! & DiaV3wrk(i,j,k,M3rate)
2044# endif
2045 END DO
2046 END DO
2047 DO i=istrt,iendt
2048 fc1(i)=fc(i,0) ! intermediate
2049 fc(i,0)=dc(i,0)*(fc(i,0)-dv_avg2(i,j)) ! recursive
2050 tl_fc(i,0)=tl_dc(i,0)*(fc1(i)-dv_avg2(i,j))+ &
2051 & dc(i,0)*(tl_fc(i,0)-tl_dv_avg2(i,j))
2052 END DO
2053 DO k=1,n(ng)
2054 DO i=istrt,iendt
2055 hvom(i,j,k)=hvom(i,j,k)-dc(i,k)*fc(i,0)
2056 tl_hvom(i,j,k)=tl_hvom(i,j,k)- &
2057 & tl_dc(i,k)*fc(i,0)- &
2058 & dc(i,k)*tl_fc(i,0)
2059 END DO
2060 END DO
2061 END IF
2062 END DO
2063!
2064!-----------------------------------------------------------------------
2065! Exchange boundary data.
2066!-----------------------------------------------------------------------
2067!
2068 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
2069!^ CALL exchange_u3d_tile (ng, tile, &
2070!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
2071!^ & u(:,:,:,nnew))
2072!^
2073 CALL exchange_u3d_tile (ng, tile, &
2074 & lbi, ubi, lbj, ubj, 1, n(ng), &
2075 & tl_u(:,:,:,nnew))
2076!^ CALL exchange_v3d_tile (ng, tile, &
2077!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
2078!^ & v(:,:,:,nnew))
2079!^
2080 CALL exchange_v3d_tile (ng, tile, &
2081 & lbi, ubi, lbj, ubj, 1, n(ng), &
2082 & tl_v(:,:,:,nnew))
2083
2084 CALL exchange_u3d_tile (ng, tile, &
2085 & lbi, ubi, lbj, ubj, 1, n(ng), &
2086 & huon)
2087 CALL exchange_u3d_tile (ng, tile, &
2088 & lbi, ubi, lbj, ubj, 1, n(ng), &
2089 & tl_huon)
2090
2091 CALL exchange_v3d_tile (ng, tile, &
2092 & lbi, ubi, lbj, ubj, 1, n(ng), &
2093 & hvom)
2094 CALL exchange_v3d_tile (ng, tile, &
2095 & lbi, ubi, lbj, ubj, 1, n(ng), &
2096 & tl_hvom)
2097!
2098 DO k=1,2
2099!^ CALL exchange_u2d_tile (ng, tile, &
2100!^ & LBi, UBi, LBj, UBj, &
2101!^ & ubar(:,:,k))
2102!^
2103 CALL exchange_u2d_tile (ng, tile, &
2104 & lbi, ubi, lbj, ubj, &
2105 & tl_ubar(:,:,k))
2106!^ CALL exchange_v2d_tile (ng, tile, &
2107!^ & LBi, UBi, LBj, UBj, &
2108!^ & vbar(:,:,k))
2109!^
2110 CALL exchange_v2d_tile (ng, tile, &
2111 & lbi, ubi, lbj, ubj, &
2112 & tl_vbar(:,:,k))
2113 END DO
2114 END IF
2115
2116# ifdef DISTRIBUTE
2117!^ CALL mp_exchange3d (ng, tile, iNLM, 4, &
2118!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
2119!^ & NghostPoints, &
2120!^ & EWperiodic(ng), NSperiodic(ng), &
2121!^ & u(:,:,:,nnew), v(:,:,:,nnew), &
2122!^ & Huon, Hvom)
2123!^
2124 CALL mp_exchange3d (ng, tile, itlm, 4, &
2125 & lbi, ubi, lbj, ubj, 1, n(ng), &
2126 & nghostpoints, &
2127 & ewperiodic(ng), nsperiodic(ng), &
2128 & tl_u(:,:,:,nnew), tl_v(:,:,:,nnew), &
2129 & tl_huon, tl_hvom)
2130 CALL mp_exchange3d (ng, tile, itlm, 2, &
2131 & lbi, ubi, lbj, ubj, 1, n(ng), &
2132 & nghostpoints, &
2133 & ewperiodic(ng), nsperiodic(ng), &
2134 & huon, hvom)
2135!
2136!^ CALL mp_exchange2d (ng, tile, iNLM, 4, &
2137!^ & LBi, UBi, LBj, UBj, &
2138!^ & NghostPoints, &
2139!^ & EWperiodic(ng), NSperiodic(ng), &
2140!^ & ubar(:,:,1), vbar(:,:,1), &
2141!^ & ubar(:,:,2), vbar(:,:,2))
2142!^
2143 CALL mp_exchange2d (ng, tile, itlm, 4, &
2144 & lbi, ubi, lbj, ubj, &
2145 & nghostpoints, &
2146 & ewperiodic(ng), nsperiodic(ng), &
2147 & tl_ubar(:,:,1), tl_vbar(:,:,1), &
2148 & tl_ubar(:,:,2), tl_vbar(:,:,2))
2149# endif
2150# ifdef UV_DESTAGGERED
2151!
2152!-----------------------------------------------------------------------
2153! Compute tangent linear 3D momentum (tl_ua, tl_va) at RHO-points
2154! (A-Grid) for output purposes and data assimilation where the
2155! observations and state vector is located at the cell-center.
2156!-----------------------------------------------------------------------
2157!
2158!> CALL uv_C2A_grid (ng, tile, iNLM, nnew)
2159!>
2160 CALL tl_uv_c2a_grid (ng, tile, itlm, nnew)
2161# endif
2162!
2163 RETURN
2164 END SUBROUTINE tl_step3d_uv_tile
2165#endif
2166 END MODULE tl_step3d_uv_mod
exchange_2d_mod
Definition exchange_2d.F:2
exchange_2d_mod::exchange_u2d_tile
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:440
exchange_2d_mod::exchange_v2d_tile
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition exchange_2d.F:627
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_ncparam
Definition mod_ncparam.F:2
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_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::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::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
tl_step3d_uv_mod
Definition tl_step3d_uv.F:2
tl_step3d_uv_mod::tl_step3d_uv
subroutine, public tl_step3d_uv(ng, tile)
Definition tl_step3d_uv.F:56
tl_step3d_uv_mod::tl_step3d_uv_tile
subroutine tl_step3d_uv_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, nstp, nnew, umask, vmask, umask_wet, vmask_wet, om_v, on_u, pm, pn, hz, tl_hz, z_r, tl_z_r, z_w, tl_z_w, akv, tl_akv, du_avg1, tl_du_avg1, dv_avg1, tl_dv_avg1, du_avg2, tl_du_avg2, dv_avg2, tl_dv_avg2, tl_ru, tl_rv, u, tl_u, v, tl_v, tl_ubar, tl_vbar, ubar_stokes, tl_ubar_stokes, vbar_stokes, tl_vbar_stokes, u_stokes, tl_u_stokes, v_stokes, tl_v_stokes, huon, tl_huon, hvom, tl_hvom)
Definition tl_step3d_uv.F:183
tl_u3dbc_mod
Definition tl_u3dbc_im.F:2
tl_u3dbc_mod::tl_u3dbc_tile
subroutine, public tl_u3dbc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nstp, nout, tl_u)
Definition tl_u3dbc_im.F:57
tl_v3dbc_mod
Definition tl_v3dbc_im.F:2
tl_v3dbc_mod::tl_v3dbc_tile
subroutine, public tl_v3dbc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nstp, nout, tl_v)
Definition tl_v3dbc_im.F:57
uv_var_change_mod
Definition uv_var_change.F:2
uv_var_change_mod::tl_uv_c2a_grid
subroutine, public tl_uv_c2a_grid(ng, tile, model, ninp)
Definition uv_var_change.F:370
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