ROMS
Loading...
Searching...
No Matches
ad_t3dmix4_iso.h
Go to the documentation of this file.
1 MODULE ad_t3dmix4_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! This subroutine computes adjoint horizontal biharmonic mixing of !
11! tracers along isopycnic surfaces. !
12! !
13! BASIC STATE variables needed: diff4, Hz, t, rho, z_r !
14! !
15!=======================================================================
16!
17 implicit none
18!
19 PRIVATE
20 PUBLIC ad_t3dmix4
21!
22 CONTAINS
23!
24!***********************************************************************
25 SUBROUTINE ad_t3dmix4 (ng, tile)
26!***********************************************************************
27!
28 USE mod_param
29#ifdef TS_MIX_CLIMA
30 USE mod_clima
31#endif
32#ifdef DIAGNOSTICS_TS
33!! USE mod_diags
34#endif
35 USE mod_grid
36 USE mod_mixing
37 USE mod_ocean
38 USE mod_stepping
39!
40! Imported variable declarations.
41!
42 integer, intent(in) :: ng, tile
43!
44! Local variable declarations.
45!
46 character (len=*), parameter :: MyFile = &
47 & __FILE__
48!
49#include "tile.h"
50!
51#ifdef PROFILE
52 CALL wclock_on (ng, iadm, 29, __line__, myfile)
53#endif
54 CALL ad_t3dmix4_iso_tile (ng, tile, &
55 & lbi, ubi, lbj, ubj, &
56 & imins, imaxs, jmins, jmaxs, &
57 & nrhs(ng), nstp(ng), nnew(ng), &
58#ifdef MASKING
59 & grid(ng) % umask, &
60 & grid(ng) % vmask, &
61#endif
62#ifdef WET_DRY_NOT_YET
63 & grid(ng) % umask_wet, &
64 & grid(ng) % vmask_wet, &
65#endif
66 & grid(ng) % om_v, &
67 & grid(ng) % on_u, &
68 & grid(ng) % pm, &
69 & grid(ng) % pn, &
70 & grid(ng) % Hz, &
71 & grid(ng) % ad_Hz, &
72 & grid(ng) % z_r, &
73 & grid(ng) % ad_z_r, &
74#ifdef DIFF_3DCOEF
75# ifdef TS_U3ADV_SPLIT
76 & mixing(ng) % diff3d_u, &
77 & mixing(ng) % diff3d_v, &
78# else
79 & mixing(ng) % diff3d_r, &
80# endif
81#else
82 & mixing(ng) % diff4, &
83#endif
84 & ocean(ng) % pden, &
85 & ocean(ng) % ad_pden, &
86#ifdef TS_MIX_CLIMA
87 & clima(ng) % tclm, &
88#endif
89#ifdef DIAGNOSTICS_TS
90!! & DIAGS(ng) % DiaTwrk, &
91#endif
92 & ocean(ng) % t, &
93 & ocean(ng) % ad_t)
94#ifdef PROFILE
95 CALL wclock_off (ng, iadm, 29, __line__, myfile)
96#endif
97!
98 RETURN
99 END SUBROUTINE ad_t3dmix4
100!
101!***********************************************************************
102 SUBROUTINE ad_t3dmix4_iso_tile (ng, tile, &
103 & LBi, UBi, LBj, UBj, &
104 & IminS, ImaxS, JminS, JmaxS, &
105 & nrhs, nstp, nnew, &
106#ifdef MASKING
107 & umask, vmask, &
108#endif
109#ifdef WET_DRY_NOT_YET
110 & umask_wet, vmask_wet, &
111#endif
112 & om_v, on_u, pm, pn, &
113 & Hz, ad_Hz, &
114 & z_r, ad_z_r, &
115#ifdef DIFF_3DCOEF
116# ifdef TS_U3ADV_SPLIT
117 & diff3d_u, diff3d_v, &
118# else
119 & diff3d_r, &
120# endif
121#else
122 & diff4, &
123#endif
124 & pden, ad_pden, &
125#ifdef TS_MIX_CLIMA
126 & tclm, &
127#endif
128#ifdef DIAGNOSTICS_TS
129!! & DiaTwrk, &
130#endif
131 & t, ad_t)
132!***********************************************************************
133!
134 USE mod_param
135 USE mod_ncparam
136 USE mod_scalars
137!
138! Imported variable declarations.
139!
140 integer, intent(in) :: ng, tile
141 integer, intent(in) :: LBi, UBi, LBj, UBj
142 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
143 integer, intent(in) :: nrhs, nstp, nnew
144
145#ifdef ASSUMED_SHAPE
146# ifdef MASKING
147 real(r8), intent(in) :: umask(LBi:,LBj:)
148 real(r8), intent(in) :: vmask(LBi:,LBj:)
149# endif
150# ifdef WET_DRY_NOT_YET
151 real(r8), intent(in) :: umask_wet(LBi:,LBj:)
152 real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
153# endif
154# ifdef DIFF_3DCOEF
155# ifdef TS_U3ADV_SPLIT
156 real(r8), intent(in) :: diff3d_u(LBi:,LBj:,:)
157 real(r8), intent(in) :: diff3d_v(LBi:,LBj:,:)
158# else
159 real(r8), intent(in) :: diff3d_r(LBi:,LBj:,:)
160# endif
161# else
162 real(r8), intent(in) :: diff4(LBi:,LBj:,:)
163# endif
164 real(r8), intent(in) :: om_v(LBi:,LBj:)
165 real(r8), intent(in) :: on_u(LBi:,LBj:)
166 real(r8), intent(in) :: pm(LBi:,LBj:)
167 real(r8), intent(in) :: pn(LBi:,LBj:)
168 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
169 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
170 real(r8), intent(in) :: pden(LBi:,LBj:,:)
171 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
172# ifdef TS_MIX_CLIMA
173 real(r8), intent(in) :: tclm(LBi:,LBj:,:,:)
174# endif
175# ifdef DIAGNOSTICS_TS
176 real(r8), intent(inout) :: DiaTwrk(LBi:,LBj:,:,:,:)
177# endif
178 real(r8), intent(inout) :: ad_Hz(LBi:,LBj:,:)
179 real(r8), intent(inout) :: ad_z_r(LBi:,LBj:,:)
180 real(r8), intent(inout) :: ad_pden(LBi:,LBj:,:)
181 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
182#else
183# ifdef MASKING
184 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
185 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
186# endif
187# ifdef WET_DRY_NOT_YET
188 real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
189 real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
190# endif
191# ifdef DIFF_3DCOEF
192# ifdef TS_U3ADV_SPLIT
193 real(r8), intent(in) :: diff3d_u(LBi:UBi,LBj:UBj,N(ng))
194 real(r8), intent(in) :: diff3d_v(LBi:UBi,LBj:UBj,N(ng))
195# else
196 real(r8), intent(in) :: diff3d_r(LBi:UBi,LBj:UBj,N(ng))
197# endif
198# else
199 real(r8), intent(in) :: diff4(LBi:UBi,LBj:UBj,NT(ng))
200# endif
201 real(r8), intent(in) :: om_v(LBi:UBi,LBj:UBj)
202 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
203 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
204 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
205 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
206 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
207 real(r8), intent(in) :: pden(LBi:UBi,LBj:UBj,N(ng))
208 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
209# ifdef TS_MIX_CLIMA
210 real(r8), intent(in) :: tclm(LBi:UBi,LBj:UBj,N(ng),NT(ng))
211# endif
212# ifdef DIAGNOSTICS_TS
213!! real(r8), intent(inout) :: DiaTwrk(LBi:UBi,LBj:UBj,N(ng),NT(ng), &
214!! & NDT)
215# endif
216 real(r8), intent(inout) :: ad_Hz(LBi:UBi,LBj:UBj,N(ng))
217 real(r8), intent(inout) :: ad_z_r(LBi:UBi,LBj:UBj,N(ng))
218 real(r8), intent(inout) :: ad_pden(LBi:UBi,LBj:UBj,N(ng))
219 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
220#endif
221!
222! Local variable declarations.
223!
224 integer :: Imin, Imax, Jmin, Jmax
225 integer :: i, itrc, j, k, kk, kt, k1, k1b, k2, k2b
226
227 real(r8), parameter :: eps = 0.5_r8
228 real(r8), parameter :: small = 1.0e-14_r8
229 real(r8), parameter :: slope_max = 0.0001_r8
230 real(r8), parameter :: strat_min = 0.1_r8
231
232 real(r8) :: cff, cff1, cff2, cff3, cff4, dife, difx
233 real(r8) :: ad_cff, ad_cff1, ad_cff2, ad_cff3, ad_cff4
234 real(r8) :: adfac, adfac1, adfac2, adfac3, adfac4
235
236 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: LapT
237
238 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: ad_LapT
239
240 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
241 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
242
243 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_FE
244 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_FX
245
246 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FS
247 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: FS1
248 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dRde
249 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dRdx
250 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTde
251 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdr
252 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: dTdx
253
254 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_FS
255 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dRde
256 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dRdx
257 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dTde
258 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dTdr
259 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,2) :: ad_dTdx
260
261#include "set_bounds.h"
262!
263!-----------------------------------------------------------------------
264! Initialize adjoint private variables.
265!-----------------------------------------------------------------------
266!
267 ad_cff=0.0_r8
268 ad_cff1=0.0_r8
269 ad_cff2=0.0_r8
270 ad_cff3=0.0_r8
271 ad_cff4=0.0_r8
272
273 ad_fe(imins:imaxs,jmins:jmaxs)=0.0_r8
274 ad_fx(imins:imaxs,jmins:jmaxs)=0.0_r8
275
276 ad_fs(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
277
278 ad_drde(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
279 ad_drdx(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
280 ad_dtde(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
281 ad_dtdr(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
282 ad_dtdx(imins:imaxs,jmins:jmaxs,1:2)=0.0_r8
283
284 ad_lapt(imins:imaxs,jmins:jmaxs,1:n(ng))=0.0_r8
285!
286!----------------------------------------------------------------------
287! Compute horizontal biharmonic diffusion along isopycnic surfaces.
288! The biharmonic operator is computed by applying the harmonic
289! operator twice.
290!----------------------------------------------------------------------
291!
292! Set local I- and J-ranges.
293!
294 IF (ewperiodic(ng)) THEN
295 imin=istr-1
296 imax=iend+1
297 ELSE
298 imin=max(istr-1,1)
299 imax=min(iend+1,lm(ng))
300 END IF
301 IF (nsperiodic(ng)) THEN
302 jmin=jstr-1
303 jmax=jend+1
304 ELSE
305 jmin=max(jstr-1,1)
306 jmax=min(jend+1,mm(ng))
307 END IF
308!
309! Compute horizontal and density gradients for the BASIC STATE. Notice
310! the recursive blocking sequence. The vertical placement of the
311! gradients is:
312!
313! dTdx,dTde(:,:,k1) k rho-points
314! dTdx,dTde(:,:,k2) k+1 rho-points
315! FS,dTdr(:,:,k1) k-1/2 W-points
316! FS,dTdr(:,:,k2) k+1/2 W-points
317!
318#ifdef TS_MIX_STABILITY
319! In order to increase stability, the biharmonic operator is applied
320! as: 3/4 t(:,:,:,nrhs,:) + 1/4 t(:,:,:,nstp,:).
321!
322#endif
323
324 t_loop : DO itrc=1,nt(ng)
325 k2=1
326 k_loop1 : DO k=0,n(ng)
327 k1=k2
328 k2=3-k1
329 IF (k.lt.n(ng)) THEN
330 DO j=jmin,jmax
331 DO i=imin,imax+1
332 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
333#ifdef MASKING
334 cff=cff*umask(i,j)
335#endif
336#ifdef WET_DRY_NOT_YET
337 cff=cff*umask_wet(i,j)
338#endif
339 drdx(i,j,k2)=cff*(pden(i ,j,k+1)- &
340 & pden(i-1,j,k+1))
341#if defined TS_MIX_STABILITY
342 dtdx(i,j,k2)=cff*(0.75_r8*(t(i ,j,k+1,nrhs,itrc)- &
343 & t(i-1,j,k+1,nrhs,itrc))+ &
344 & 0.25_r8*(t(i ,j,k+1,nstp,itrc)- &
345 & t(i-1,j,k+1,nstp,itrc)))
346#elif defined TS_MIX_CLIMA
347 IF (ltracerclm(itrc,ng)) THEN
348 dtdx(i,j,k2)=cff*((t(i ,j,k+1,nrhs,itrc)- &
349 & tclm(i ,j,k+1,itrc))- &
350 & (t(i-1,j,k+1,nrhs,itrc)- &
351 & tclm(i-1,j,k+1,itrc)))
352 ELSE
353 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
354 & t(i-1,j,k+1,nrhs,itrc))
355 END IF
356#else
357 dtdx(i,j,k2)=cff*(t(i ,j,k+1,nrhs,itrc)- &
358 & t(i-1,j,k+1,nrhs,itrc))
359#endif
360 END DO
361 END DO
362 DO j=jmin,jmax+1
363 DO i=imin,imax
364 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
365#ifdef MASKING
366 cff=cff*vmask(i,j)
367#endif
368#ifdef WET_DRY_NOT_YET
369 cff=cff*vmask_wet(i,j)
370#endif
371 drde(i,j,k2)=cff*(pden(i,j ,k+1)- &
372 & pden(i,j-1,k+1))
373#if defined TS_MIX_STABILITY
374 dtde(i,j,k2)=cff*(0.75_r8*(t(i,j ,k+1,nrhs,itrc)- &
375 & t(i,j-1,k+1,nrhs,itrc))+ &
376 & 0.25_r8*(t(i,j ,k+1,nstp,itrc)- &
377 & t(i,j-1,k+1,nstp,itrc)))
378#elif defined TS_MIX_CLIMA
379 IF (ltracerclm(itrc,ng)) THEN
380 dtde(i,j,k2)=cff*((t(i,j ,k+1,nrhs,itrc)- &
381 & tclm(i,j ,k+1,itrc))- &
382 & (t(i,j-1,k+1,nrhs,itrc)- &
383 & tclm(i,j-1,k+1,itrc)))
384 ELSE
385 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
386 & t(i,j-1,k+1,nrhs,itrc))
387 END IF
388#else
389 dtde(i,j,k2)=cff*(t(i,j ,k+1,nrhs,itrc)- &
390 & t(i,j-1,k+1,nrhs,itrc))
391#endif
392 END DO
393 END DO
394 END IF
395 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
396 DO j=-1+jmin,jmax+1
397 DO i=-1+imin,imax+1
398 dtdr(i,j,k2)=0.0_r8
399 fs(i,j,k2)=0.0_r8
400 END DO
401 END DO
402 ELSE
403 DO j=-1+jmin,jmax+1
404 DO i=-1+imin,imax+1
405#if defined TS_MIX_MAX_SLOPE
406 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
407 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
408 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
409 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
410 cff2=0.25_r8*slope_max* &
411 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
412 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
413 cff4=max(cff2,cff3)
414 cff=-1.0_r8/cff4
415#elif defined TS_MIX_MIN_STRAT
416 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
417 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
418 cff=-1.0_r8/cff1
419#else
420 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
421 cff=-1.0_r8/cff1
422#endif
423#if defined TS_MIX_STABILITY
424 dtdr(i,j,k2)=cff*(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
425 & t(i,j,k ,nrhs,itrc))+ &
426 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
427 & t(i,j,k ,nstp,itrc)))
428#elif defined TS_MIX_CLIMA
429 IF (ltracerclm(itrc,ng)) THEN
430 dtdr(i,j,k2)=cff*((t(i,j,k+1,nrhs,itrc)- &
431 & tclm(i,j,k+1,itrc))- &
432 & (t(i,j,k ,nrhs,itrc)- &
433 & tclm(i,j,k ,itrc)))
434 ELSE
435 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
436 & t(i,j,k ,nrhs,itrc))
437 END IF
438#else
439 dtdr(i,j,k2)=cff*(t(i,j,k+1,nrhs,itrc)- &
440 & t(i,j,k ,nrhs,itrc))
441#endif
442 fs(i,j,k2)=cff*(z_r(i,j,k+1)- &
443 & z_r(i,j,k ))
444 END DO
445 END DO
446 END IF
447 IF (k.gt.0) THEN
448 DO j=jmin,jmax
449 DO i=imin,imax+1
450#ifdef DIFF_3DCOEF
451# ifdef TS_U3ADV_SPLIT
452 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
453# else
454 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
455 & on_u(i,j)
456# endif
457#else
458 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
459 & on_u(i,j)
460#endif
461 fx(i,j)=cff* &
462 & (hz(i,j,k)+hz(i-1,j,k))* &
463 & (dtdx(i,j,k1)- &
464 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
465 & (dtdr(i-1,j,k1)+ &
466 & dtdr(i ,j,k2))+ &
467 & min(drdx(i,j,k1),0.0_r8)* &
468 & (dtdr(i-1,j,k2)+ &
469 & dtdr(i ,j,k1))))
470 END DO
471 END DO
472 DO j=jmin,jmax+1
473 DO i=imin,imax
474#ifdef DIFF_3DCOEF
475# ifdef TS_U3ADV_SPLIT
476 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
477# else
478 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
479 & om_v(i,j)
480# endif
481#else
482 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
483 & om_v(i,j)
484#endif
485 fe(i,j)=cff* &
486 & (hz(i,j,k)+hz(i,j-1,k))* &
487 & (dtde(i,j,k1)- &
488 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
489 & (dtdr(i,j-1,k1)+ &
490 & dtdr(i,j ,k2))+ &
491 & min(drde(i,j,k1),0.0_r8)* &
492 & (dtdr(i,j-1,k2)+ &
493 & dtdr(i,j ,k1))))
494 END DO
495 END DO
496 IF (k.lt.n(ng)) THEN
497 DO j=jmin,jmax
498 DO i=imin,imax
499#ifdef DIFF_3DCOEF
500# ifdef TS_U3ADV_SPLIT
501 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
502 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
503 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
504 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
505# else
506 difx=0.5_r8*diff3d_r(i,j,k)
507 dife=difx
508# endif
509#else
510 difx=0.5_r8*diff4(i,j,itrc)
511 dife=difx
512#endif
513 cff1=max(drdx(i ,j,k1),0.0_r8)
514 cff2=max(drdx(i+1,j,k2),0.0_r8)
515 cff3=min(drdx(i ,j,k2),0.0_r8)
516 cff4=min(drdx(i+1,j,k1),0.0_r8)
517 cff=difx* &
518 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
519 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
520 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
521 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
522 cff1=max(drde(i,j ,k1),0.0_r8)
523 cff2=max(drde(i,j+1,k2),0.0_r8)
524 cff3=min(drde(i,j ,k2),0.0_r8)
525 cff4=min(drde(i,j+1,k1),0.0_r8)
526 cff=cff+ &
527 & dife* &
528 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
529 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
530 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
531 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
532 fs(i,j,k2)=cff*fs(i,j,k2)
533 END DO
534 END DO
535 END IF
536!
537! Compute first BASIC STATE harmonic operator, without mixing
538! coefficient. Multiply by the metrics of the second harmonic
539! operator. Save into work array "LapT".
540!
541 DO j=jmin,jmax
542 DO i=imin,imax
543 cff=pm(i,j)*pn(i,j)
544 cff1=1.0_r8/hz(i,j,k)
545 lapt(i,j,k)=cff1*(cff* &
546 & (fx(i+1,j)-fx(i,j)+ &
547 & fe(i,j+1)-fe(i,j))+ &
548 & (fs(i,j,k2)-fs(i,j,k1)))
549 END DO
550 END DO
551 END IF
552 END DO k_loop1
553!
554! Apply boundary conditions (except periodic; closed or gradient)
555! to the first BASIC STATE harmonic operator.
556!
557 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
558 IF (domain(ng)%Western_Edge(tile)) THEN
559 IF (ad_lbc(iwest,istvar(itrc),ng)%closed) THEN
560 DO k=1,n(ng)
561 DO j=jmin,jmax
562 lapt(istr-1,j,k)=0.0_r8
563 END DO
564 END DO
565 ELSE
566 DO k=1,n(ng)
567 DO j=jmin,jmax
568 lapt(istr-1,j,k)=lapt(istr,j,k)
569 END DO
570 END DO
571 END IF
572 END IF
573 END IF
574!
575 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
576 IF (domain(ng)%Eastern_Edge(tile)) THEN
577 IF (ad_lbc(ieast,istvar(itrc),ng)%closed) THEN
578 DO k=1,n(ng)
579 DO j=jmin,jmax
580 lapt(iend+1,j,k)=0.0_r8
581 END DO
582 END DO
583 ELSE
584 DO k=1,n(ng)
585 DO j=jmin,jmax
586 lapt(iend+1,j,k)=lapt(iend,j,k)
587 END DO
588 END DO
589 END IF
590 END IF
591 END IF
592!
593 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
594 IF (domain(ng)%Southern_Edge(tile)) THEN
595 IF (ad_lbc(isouth,istvar(itrc),ng)%closed) THEN
596 DO k=1,n(ng)
597 DO i=imin,imax
598 lapt(i,jstr-1,k)=0.0_r8
599 END DO
600 END DO
601 ELSE
602 DO k=1,n(ng)
603 DO i=imin,imax
604 lapt(i,jstr-1,k)=lapt(i,jstr,k)
605 END DO
606 END DO
607 END IF
608 END IF
609 END IF
610!
611 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
612 IF (domain(ng)%Northern_Edge(tile)) THEN
613 IF (ad_lbc(inorth,istvar(itrc),ng)%closed) THEN
614 DO k=1,n(ng)
615 DO i=imin,imax
616 lapt(i,jend+1,k)=0.0_r8
617 END DO
618 END DO
619 ELSE
620 DO k=1,n(ng)
621 DO i=imin,imax
622 lapt(i,jend+1,k)=lapt(i,jend,k)
623 END DO
624 END DO
625 END IF
626 END IF
627 END IF
628!
629 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
630 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
631 IF (domain(ng)%SouthWest_Corner(tile)) THEN
632 DO k=1,n(ng)
633 lapt(istr-1,jstr-1,k)=0.5_r8* &
634 & (lapt(istr ,jstr-1,k)+ &
635 & lapt(istr-1,jstr ,k))
636 END DO
637 END IF
638 END IF
639
640 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
641 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
642 IF (domain(ng)%SouthEast_Corner(tile)) THEN
643 DO k=1,n(ng)
644 lapt(iend+1,jstr-1,k)=0.5_r8* &
645 & (lapt(iend ,jstr-1,k)+ &
646 & lapt(iend+1,jstr ,k))
647 END DO
648 END IF
649 END IF
650
651 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
652 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
653 IF (domain(ng)%NorthWest_Corner(tile)) THEN
654 DO k=1,n(ng)
655 lapt(istr-1,jend+1,k)=0.5_r8* &
656 & (lapt(istr ,jend+1,k)+ &
657 & lapt(istr-1,jend ,k))
658 END DO
659 END IF
660 END IF
661
662 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
663 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
664 IF (domain(ng)%NorthEast_Corner(tile)) THEN
665 DO k=1,n(ng)
666 lapt(iend+1,jend+1,k)=0.5_r8* &
667 & (lapt(iend ,jend+1,k)+ &
668 & lapt(iend+1,jend ,k))
669 END DO
670 END IF
671 END IF
672!
673! Compute adjoint of starting storage recursive indices k1 and k2.
674!
675 k1=2
676 k2=1
677 DO k=0,n(ng)
678!!
679!! Note: The following code is equivalent to
680!!
681!! kt=k1
682!! k1=k2
683!! k2=kt
684!!
685!! We use the adjoint of the above code.
686!!
687 k1=k2
688 k2=3-k1
689 END DO
690!
691! Compute required basic state fields. Need to look forward in
692! recursive kk index.
693!
694 k_loop2: DO k=n(ng),0,-1
695 k2b=1
696 DO kk=0,k
697 k1b=k2b
698 k2b=3-k1b
699!
700! Compute components of the rotated tracer flux (T m3/s) along
701! isopycnic surfaces (required BASIC STATE fields).
702!
703 IF (kk.lt.n(ng)) THEN
704 DO j=jstr,jend
705 DO i=istr,iend+1
706 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
707#ifdef MASKING
708 cff=cff*umask(i,j)
709#endif
710#ifdef WET_DRY_NOT_YET
711 cff=cff*umask_wet(i,j)
712#endif
713 drdx(i,j,k2b)=cff*(pden(i ,j,kk+1)- &
714 & pden(i-1,j,kk+1))
715 dtdx(i,j,k2b)=cff*(lapt(i ,j,kk+1)- &
716 & lapt(i-1,j,kk+1))
717 END DO
718 END DO
719 IF (kk.eq.0) THEN
720 DO j=jstr,jend
721 DO i=istr,iend+1
722 drdx(i,j,k1b)=0.0_r8
723 dtdx(i,j,k1b)=0.0_r8
724 END DO
725 END DO
726 END IF
727 DO j=jstr,jend+1
728 DO i=istr,iend
729 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
730#ifdef MASKING
731 cff=cff*vmask(i,j)
732#endif
733#ifdef WET_DRY_NOT_YET
734 cff=cff*vmask_wet(i,j)
735#endif
736 drde(i,j,k2b)=cff*(pden(i,j ,kk+1)- &
737 & pden(i,j-1,kk+1))
738 dtde(i,j,k2b)=cff*(lapt(i,j ,kk+1)- &
739 & lapt(i,j-1,kk+1))
740 END DO
741 END DO
742 IF (kk.eq.0) THEN
743 DO j=jstr,jend+1
744 DO i=istr,iend
745 drde(i,j,k1b)=0.0_r8
746 dtde(i,j,k1b)=0.0_r8
747 END DO
748 END DO
749 END IF
750 END IF
751 IF ((kk.eq.0).or.(kk.eq.n(ng))) THEN
752 DO j=jstr-1,jend+1
753 DO i=istr-1,iend+1
754 dtdr(i,j,k2b)=0.0_r8
755 fs(i,j,k2b)=0.0_r8
756 END DO
757 END DO
758 IF (kk.eq.0) THEN
759 DO j=jstr-1,jend+1
760 DO i=istr-1,iend+1
761 dtdr(i,j,k1b)=0.0_r8
762 fs(i,j,k1b)=0.0_r8
763 END DO
764 END DO
765 END IF
766 ELSE
767 DO j=jstr-1,jend+1
768 DO i=istr-1,iend+1
769#if defined TS_MIX_MAX_SLOPE
770 cff1=sqrt(drdx(i,j,k2b)**2+drdx(i+1,j,k2b)**2+ &
771 & drdx(i,j,k1b)**2+drdx(i+1,j,k1b)**2+ &
772 & drde(i,j,k2b)**2+drde(i,j+1,k2b)**2+ &
773 & drde(i,j,k1b)**2+drde(i,j+1,k1b)**2)
774 cff2=0.25_r8*slope_max* &
775 & (z_r(i,j,kk+1)-z_r(i,j,kk))*cff1
776 cff3=max(pden(i,j,kk)-pden(i,j,kk+1),small)
777 cff4=max(cff2,cff3)
778 cff=-1.0_r8/cff4
779#elif defined TS_MIX_MIN_STRAT
780 cff1=max(pden(i,j,kk)-pden(i,j,kk+1), &
781 & strat_min*(z_r(i,j,kk+1)-z_r(i,j,kk)))
782 cff=-1.0_r8/cff1
783#else
784 cff1=max(pden(i,j,kk)-pden(i,j,kk+1),eps)
785 cff=-1.0_r8/cff1
786#endif
787 dtdr(i,j,k2b)=cff*(lapt(i,j,kk+1)- &
788 & lapt(i,j,kk ))
789 fs(i,j,k2b)=cff*(z_r(i,j,kk+1)- &
790 & z_r(i,j,kk ))
791 END DO
792 END DO
793 END IF
794 END DO
795!
796 IF (k.gt.0) THEN
797!
798! Time-step biharmonic, isopycnal diffusion term (m Tunits).
799!
800 DO j=jstr,jend
801 DO i=istr,iend
802#ifdef DIAGNOSTICS_TS
803!! DiaTwrk(i,j,k,itrc,iThdif)=-cff
804#endif
805!^ tl_t(i,j,k,nnew,itrc)=tl_t(i,j,k,nnew,itrc)-tl_cff
806!^
807 ad_cff=ad_cff-ad_t(i,j,k,nnew,itrc)
808!^ tl_cff=dt(ng)*pm(i,j)*pn(i,j)* &
809!^ & (tl_FX(i+1,j)-tl_FX(i,j)+ &
810!^ & tl_FE(i,j+1)-tl_FE(i,j))+ &
811!^ & dt(ng)*(tl_FS(i,j,k2)-tl_FS(i,j,k1))
812!^
813 adfac=dt(ng)*ad_cff
814 adfac1=adfac*pm(i,j)*pn(i,j)
815 ad_fs(i,j,k1)=ad_fs(i,j,k1)-adfac
816 ad_fs(i,j,k2)=ad_fs(i,j,k2)+adfac
817 ad_fx(i ,j)=ad_fx(i ,j)-adfac1
818 ad_fx(i+1,j)=ad_fx(i+1,j)+adfac1
819 ad_fe(i,j )=ad_fe(i,j )-adfac1
820 ad_fe(i,j+1)=ad_fe(i,j+1)+adfac1
821 ad_cff=0.0_r8
822 END DO
823 END DO
824 IF (k.lt.n(ng)) THEN
825 DO j=jstr,jend
826 DO i=istr,iend
827#ifdef DIFF_3DCOEF
828# ifdef TS_U3ADV_SPLIT
829 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
830 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
831 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
832 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
833# else
834 difx=0.5_r8*diff3d_r(i,j,k)
835 dife=difx
836# endif
837#else
838 difx=0.5_r8*diff4(i,j,itrc)
839 dife=difx
840#endif
841 cff1=max(drdx(i ,j,k1),0.0_r8)
842 cff2=max(drdx(i+1,j,k2),0.0_r8)
843 cff3=min(drdx(i ,j,k2),0.0_r8)
844 cff4=min(drdx(i+1,j,k1),0.0_r8)
845 cff=difx* &
846 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
847 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
848 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
849 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
850 cff1=max(drde(i,j ,k1),0.0_r8)
851 cff2=max(drde(i,j+1,k2),0.0_r8)
852 cff3=min(drde(i,j ,k2),0.0_r8)
853 cff4=min(drde(i,j+1,k1),0.0_r8)
854 cff=cff+ &
855 & dife* &
856 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
857 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
858 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
859 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
860!^ tl_FS(i,j,k2)=tl_cff*FS(i,j,k2)+ &
861!^ & cff*tl_FS(i,j,k2
862!^
863 ad_cff=ad_cff+fs(i,j,k2)*ad_fs(i,j,k2)
864 ad_fs(i,j,k2)=cff*ad_fs(i,j,k2)
865!^ tl_cff=tl_cff+ &
866!^ & dife* &
867!^ & (tl_cff1*(cff1*dTdr(i,j,k2)- &
868!^ & dTde(i,j ,k1))+ &
869!^ & tl_cff2*(cff2*dTdr(i,j,k2)- &
870!^ & dTde(i,j+1,k2))+ &
871!^ & tl_cff3*(cff3*dTdr(i,j,k2)- &
872!^ & dTde(i,j ,k2))+ &
873!^ & tl_cff4*(cff4*dTdr(i,j,k2)- &
874!^ & dTde(i,j+1,k1))+ &
875!^ & cff1*(tl_cff1*dTdr(i,j,k2)+ &
876!^ & cff1*tl_dTdr(i,j,k2)- &
877!^ & tl_dTde(i,j ,k1))+ &
878!^ & cff2*(tl_cff2*dTdr(i,j,k2)+ &
879!^ & cff2*tl_dTdr(i,j,k2)- &
880!^ & tl_dTde(i,j+1,k2))+ &
881!^ & cff3*(tl_cff3*dTdr(i,j,k2)+ &
882!^ & cff3*tl_dTdr(i,j,k2)- &
883!^ & tl_dTde(i,j ,k2))+ &
884!^ & cff4*(tl_cff4*dTdr(i,j,k2)+ &
885!^ & cff4*tl_dTdr(i,j,k2)- &
886!^ & tl_dTde(i,j+1,k1)))
887!^
888 adfac=dife*ad_cff
889 ad_cff1=ad_cff1+ &
890 & (2.0_r8*cff1*dtdr(i,j,k2)-dtde(i,j ,k1))* &
891 & adfac
892 ad_cff2=ad_cff2+ &
893 & (2.0_r8*cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))* &
894 & adfac
895 ad_cff3=ad_cff3+ &
896 & (2.0_r8*cff3*dtdr(i,j,k2)-dtde(i,j ,k2))* &
897 & adfac
898 ad_cff4=ad_cff4+ &
899 & (2.0_r8*cff4*dtdr(i,j,k2)-dtde(i,j+1,k1))* &
900 & adfac
901 ad_dtdr(i,j,k2)=ad_dtdr(i,j,k2)+ &
902 & (cff1*cff1+ &
903 & cff2*cff2+ &
904 & cff3*cff3+ &
905 & cff4*cff4)*adfac
906 ad_dtde(i,j ,k1)=ad_dtde(i,j ,k1)-cff1*adfac
907 ad_dtde(i,j+1,k2)=ad_dtde(i,j+1,k2)-cff2*adfac
908 ad_dtde(i,j ,k2)=ad_dtde(i,j ,k2)-cff3*adfac
909 ad_dtde(i,j+1,k1)=ad_dtde(i,j+1,k1)-cff4*adfac
910!^ tl_cff4=(0.5_r8+SIGN(0.5_r8,-dRde(i,j+1,k1)))* &
911!^ & tl_dRde(i,j+1,k1)
912!^
913 ad_drde(i,j+1,k1)=ad_drde(i,j+1,k1)+ &
914 & (0.5_r8+sign(0.5_r8, &
915 & -drde(i,j+1,k1)))* &
916 & ad_cff4
917 ad_cff4=0.0_r8
918!^ tl_cff3=(0.5_r8+SIGN(0.5_r8,-dRde(i,j ,k2)))* &
919!^ & tl_dRde(i,j ,k2)
920!^
921 ad_drde(i,j ,k2)=ad_drde(i,j ,k2)+ &
922 & (0.5_r8+sign(0.5_r8, &
923 & -drde(i,j ,k2)))* &
924 & ad_cff3
925 ad_cff3=0.0_r8
926!^ tl_cff2=(0.5_r8+SIGN(0.5_r8, dRde(i,j+1,k2)))* &
927!^ & tl_dRde(i,j+1,k2)
928!^
929 ad_drde(i,j+1,k2)=ad_drde(i,j+1,k2)+ &
930 & (0.5_r8+sign(0.5_r8, &
931 & drde(i,j+1,k2)))* &
932 & ad_cff2
933 ad_cff2=0.0_r8
934!^ tl_cff1=(0.5_r8+SIGN(0.5_r8, dRde(i,j ,k1)))* &
935!^ & tl_dRde(i,j ,k1)
936!^
937 ad_drde(i,j ,k1)=ad_drde(i,j ,k1)+ &
938 & (0.5_r8+sign(0.5_r8, &
939 & drde(i,j ,k1)))* &
940 & ad_cff1
941 ad_cff1=0.0_r8
942
943 cff1=max(drdx(i ,j,k1),0.0_r8)
944 cff2=max(drdx(i+1,j,k2),0.0_r8)
945 cff3=min(drdx(i ,j,k2),0.0_r8)
946 cff4=min(drdx(i+1,j,k1),0.0_r8)
947!^ tl_cff=difx* &
948!^ & (tl_cff1*(cff1*dTdr(i,j,k2)- &
949!^ & dTdx(i ,j,k1))+ &
950!^ & tl_cff2*(cff2*dTdr(i,j,k2)- &
951!^ & dTdx(i+1,j,k2))+ &
952!^ & tl_cff3*(cff3*dTdr(i,j,k2)- &
953!^ & dTdx(i ,j,k2))+ &
954!^ & tl_cff4*(cff4*dTdr(i,j,k2)- &
955!^ & dTdx(i+1,j,k1))+ &
956!^ & cff1*(tl_cff1*dTdr(i,j,k2)+ &
957!^ & cff1*tl_dTdr(i,j,k2)- &
958!^ & tl_dTdx(i ,j,k1))+ &
959!^ & cff2*(tl_cff2*dTdr(i,j,k2)+ &
960!^ & cff2*tl_dTdr(i,j,k2)- &
961!^ & tl_dTdx(i+1,j,k2))+ &
962!^ & cff3*(tl_cff3*dTdr(i,j,k2)+ &
963!^ & cff3*tl_dTdr(i,j,k2)- &
964!^ & tl_dTdx(i ,j,k2))+ &
965!^ & cff4*(tl_cff4*dTdr(i,j,k2)+ &
966!^ & cff4*tl_dTdr(i,j,k2)- &
967!^ & tl_dTdx(i+1,j,k1)))
968!^
969 adfac=difx*ad_cff
970 ad_cff1=ad_cff1+ &
971 & (2.0_r8*cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))* &
972 & adfac
973 ad_cff2=ad_cff2+ &
974 & (2.0_r8*cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))* &
975 & adfac
976 ad_cff3=ad_cff3+ &
977 & (2.0_r8*cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))* &
978 & adfac
979 ad_cff4=ad_cff4+ &
980 & (2.0_r8*cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1))* &
981 & adfac
982 ad_dtdr(i,j,k2)=ad_dtdr(i,j,k2)+ &
983 & (cff1*cff1+ &
984 & cff2*cff2+ &
985 & cff3*cff3+ &
986 & cff4*cff4)*adfac
987 ad_dtdx(i ,j,k1)=ad_dtdx(i ,j,k1)-cff1*adfac
988 ad_dtdx(i+1,j,k2)=ad_dtdx(i+1,j,k2)-cff2*adfac
989 ad_dtdx(i ,j,k2)=ad_dtdx(i ,j,k2)-cff3*adfac
990 ad_dtdx(i+1,j,k1)=ad_dtdx(i+1,j,k1)-cff4*adfac
991 ad_cff=0.0_r8
992!^ tl_cff4=(0.5_r8+SIGN(0.5_r8,-dRdx(i+1,j,k1)))* &
993!^ & tl_dRdx(i+1,j,k1)
994!^
995 ad_drdx(i+1,j,k1)=ad_drdx(i+1,j,k1)+ &
996 & (0.5_r8+sign(0.5_r8, &
997 & -drdx(i+1,j,k1)))* &
998 & ad_cff4
999 ad_cff4=0.0_r8
1000!^ tl_cff3=(0.5_r8+SIGN(0.5_r8,-dRdx(i ,j,k2)))* &
1001!^ & tl_dRdx(i ,j,k2)
1002!^
1003 ad_drdx(i ,j,k2)=ad_drdx(i ,j,k2)+ &
1004 & (0.5_r8+sign(0.5_r8, &
1005 & -drdx(i ,j,k2)))* &
1006 & ad_cff3
1007 ad_cff3=0.0_r8
1008!^ tl_cff2=(0.5_r8+SIGN(0.5_r8, dRdx(i+1,j,k2)))* &
1009!^ & tl_dRdx(i+1,j,k2)
1010!^
1011 ad_drdx(i+1,j,k2)=ad_drdx(i+1,j,k2)+ &
1012 & (0.5_r8+sign(0.5_r8, &
1013 & drdx(i+1,j,k2)))* &
1014 & ad_cff2
1015 ad_cff2=0.0_r8
1016!^ tl_cff1=(0.5_r8+SIGN(0.5_r8, dRdx(i ,j,k1)))* &
1017!^ & tl_dRdx(i ,j,k1)
1018!^
1019 ad_drdx(i ,j,k1)=ad_drdx(i ,j,k1)+ &
1020 & (0.5_r8+sign(0.5_r8, &
1021 & drdx(i ,j,k1)))* &
1022 & ad_cff1
1023 ad_cff1=0.0_r8
1024 END DO
1025 END DO
1026 END IF
1027 DO j=jstr,jend+1
1028 DO i=istr,iend
1029#ifdef DIFF_3DCOEF
1030# ifdef TS_U3ADV_SPLIT
1031 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1032# else
1033 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1034 & om_v(i,j)
1035# endif
1036#else
1037 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1038 & om_v(i,j)
1039#endif
1040!^ tl_FE(i,j)=cff* &
1041!^ & ((tl_Hz(i,j,k)+tl_Hz(i,j-1,k))* &
1042!^ & (dTde(i,j,k1)- &
1043!^ & 0.5_r8*(MAX(dRde(i,j,k1),0.0_r8)* &
1044!^ & (dTdr(i,j-1,k1)+ &
1045!^ & dTdr(i,j ,k2))+ &
1046!^ & MIN(dRde(i,j,k1),0.0_r8)* &
1047!^ & (dTdr(i,j-1,k2)+ &
1048!^ & dTdr(i,j ,k1))))+ &
1049!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
1050!^ & (tl_dTde(i,j,k1)- &
1051!^ & 0.5_r8*(MAX(dRde(i,j,k1),0.0_r8)* &
1052!^ & (tl_dTdr(i,j-1,k1)+ &
1053!^ & tl_dTdr(i ,j,k2))+ &
1054!^ & MIN(dRde(i,j,k1),0.0_r8)* &
1055!^ & (tl_dTdr(i,j-1,k2)+ &
1056!^ & tl_dTdr(i,j ,k1)))- &
1057!^ & 0.5_r8*((0.5_r8+ &
1058!^ & SIGN(0.5_r8, dRde(i,j,k1)))* &
1059!^ & tl_dRde(i,j,k1)* &
1060!^ & (dTdr(i,j-1,k1)+dTdr(i,j,k2))+ &
1061!^ & (0.5_r8+ &
1062!^ & SIGN(0.5_r8,-dRde(i,j,k1)))* &
1063!^ & tl_dRde(i,j,k1)* &
1064!^ & (dTdr(i,j-1,k2)+dTdr(i,j,k1)))))
1065!^
1066 adfac=cff*ad_fe(i,j)
1067 adfac1=adfac*(dtde(i,j,k1)- &
1068 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
1069 & (dtdr(i,j-1,k1)+ &
1070 & dtdr(i,j ,k2))+ &
1071 & min(drde(i,j,k1),0.0_r8)* &
1072 & (dtdr(i,j-1,k2)+ &
1073 & dtdr(i,j ,k1))))
1074 adfac2=adfac*(hz(i,j,k)+hz(i,j-1,k))
1075 adfac3=adfac2*0.5_r8*max(drde(i,j,k1),0.0_r8)
1076 adfac4=adfac2*0.5_r8*min(drde(i,j,k1),0.0_r8)
1077 ad_hz(i,j-1,k)=ad_hz(i,j-1,k)+adfac1
1078 ad_hz(i,j ,k)=ad_hz(i,j ,k)+adfac1
1079 ad_dtde(i,j,k1)=ad_dtde(i,j,k1)+adfac2
1080 ad_dtdr(i,j-1,k1)=ad_dtdr(i,j-1,k1)-adfac3
1081 ad_dtdr(i,j ,k2)=ad_dtdr(i,j ,k2)-adfac3
1082 ad_dtdr(i,j-1,k2)=ad_dtdr(i,j-1,k2)-adfac4
1083 ad_dtdr(i,j ,k1)=ad_dtdr(i,j ,k1)-adfac4
1084 ad_drde(i,j,k1)=ad_drde(i,j,k1)- &
1085 & adfac2*0.5_r8* &
1086 & ((0.5_r8+sign(0.5_r8, drde(i,j,k1)))* &
1087 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
1088 & (0.5_r8+sign(0.5_r8,-drde(i,j,k1)))* &
1089 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))
1090 ad_fe(i,j)=0.0_r8
1091 END DO
1092 END DO
1093 DO j=jstr,jend
1094 DO i=istr,iend+1
1095#ifdef DIFF_3DCOEF
1096# ifdef TS_U3ADV_SPLIT
1097 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
1098# else
1099 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
1100 & on_u(i,j)
1101# endif
1102#else
1103 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
1104 & on_u(i,j)
1105#endif
1106!^ tl_FX(i,j)=cff* &
1107!^ & ((tl_Hz(i,j,k)+tl_Hz(i-1,j,k))* &
1108!^ & (dTdx(i,j,k1)- &
1109!^ & 0.5_r8*(MAX(dRdx(i,j,k1),0.0_r8)* &
1110!^ & (dTdr(i-1,j,k1)+ &
1111!^ & dTdr(i ,j,k2))+ &
1112!^ & MIN(dRdx(i,j,k1),0.0_r8)* &
1113!^ & (dTdr(i-1,j,k2)+ &
1114!^ & dTdr(i ,j,k1))))+ &
1115!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
1116!^ & (tl_dTdx(i,j,k1)- &
1117!^ & 0.5_r8*(MAX(dRdx(i,j,k1),0.0_r8)* &
1118!^ & (tl_dTdr(i-1,j,k1)+ &
1119!^ & tl_dTdr(i ,j,k2))+ &
1120!^ & MIN(dRdx(i,j,k1),0.0_r8)* &
1121!^ & (tl_dTdr(i-1,j,k2)+ &
1122!^ & tl_dTdr(i ,j,k1)))- &
1123!^ & 0.5_r8*((0.5_r8+ &
1124!^ & SIGN(0.5_r8, dRdx(i,j,k1)))* &
1125!^ & tl_dRdx(i,j,k1)* &
1126!^ & (dTdr(i-1,j,k1)+dTdr(i,j,k2))+ &
1127!^ & (0.5_r8+ &
1128!^ & SIGN(0.5_r8,-dRdx(i,j,k1)))* &
1129!^ & tl_dRdx(i,j,k1)* &
1130!^ & (dTdr(i-1,j,k2)+dTdr(i,j,k1)))))
1131!^
1132 adfac=cff*ad_fx(i,j)
1133 adfac1=adfac*(dtdx(i ,j,k1)- &
1134 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
1135 & (dtdr(i-1,j,k1)+ &
1136 & dtdr(i ,j,k2))+ &
1137 & min(drdx(i,j,k1),0.0_r8)* &
1138 & (dtdr(i-1,j,k2)+ &
1139 & dtdr(i ,j,k1))))
1140 adfac2=adfac*(hz(i,j,k)+hz(i-1,j,k))
1141 adfac3=adfac2*0.5_r8*max(drdx(i,j,k1),0.0_r8)
1142 adfac4=adfac2*0.5_r8*min(drdx(i,j,k1),0.0_r8)
1143 ad_hz(i-1,j,k)=ad_hz(i-1,j,k)+adfac1
1144 ad_hz(i ,j,k)=ad_hz(i ,j,k)+adfac1
1145 ad_dtdx(i ,j,k1)=ad_dtdx(i ,j,k1)+adfac2
1146 ad_dtdr(i-1,j,k1)=ad_dtdr(i-1,j,k1)-adfac3
1147 ad_dtdr(i ,j,k2)=ad_dtdr(i ,j,k2)-adfac3
1148 ad_dtdr(i-1,j,k2)=ad_dtdr(i-1,j,k2)-adfac4
1149 ad_dtdr(i ,j,k1)=ad_dtdr(i ,j,k1)-adfac4
1150 ad_drdx(i,j,k1)=ad_drdx(i,j,k1)- &
1151 & adfac2*0.5_r8* &
1152 & ((0.5_r8+sign(0.5_r8, drdx(i,j,k1)))* &
1153 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
1154 & (0.5_r8+sign(0.5_r8,-drdx(i,j,k1)))* &
1155 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))
1156 ad_fx(i,j)=0.0_r8
1157 END DO
1158 END DO
1159 END IF
1160 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
1161 DO j=jstr-1,jend+1
1162 DO i=istr-1,iend+1
1163!^ tl_FS(i,j,k2)=0.0_r8
1164!^
1165 ad_fs(i,j,k2)=0.0_r8
1166!^ tl_dTdr(i,j,k2)=0.0_r8
1167!^
1168 ad_dtdr(i,j,k2)=0.0_r8
1169 END DO
1170 END DO
1171 ELSE
1172 DO j=jstr-1,jend+1
1173 DO i=istr-1,iend+1
1174#if defined TS_MIX_MAX_SLOPE
1175 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
1176 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
1177 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
1178 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
1179 cff2=0.25_r8*slope_max* &
1180 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
1181 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
1182 cff4=max(cff2,cff3)
1183 cff=-1.0_r8/cff4
1184#elif defined TS_MIX_MIN_STRAT
1185 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
1186 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
1187 cff=-1.0_r8/cff1
1188#else
1189 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
1190 cff=-1.0_r8/cff1
1191#endif
1192!^ tl_FS(i,j,k2)=tl_cff*(z_r(i,j,k+1)- &
1193!^ & z_r(i,j,k ))+ &
1194!^ & cff*(tl_z_r(i,j,k+1)- &
1195!^ & tl_z_r(i,j,k ))
1196!^
1197 adfac=cff*ad_fs(i,j,k2)
1198 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac
1199 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac
1200 ad_cff=ad_cff+(z_r(i,j,k+1)- &
1201 & z_r(i,j,k ))*ad_fs(i,j,k2)
1202 ad_fs(i,j,k2)=0.0_r8
1203!^ tl_dTdr(i,j,k2)=tl_cff*(LapT(i,j,k+1)- &
1204!^ & LapT(i,j,k ))+ &
1205!^ & cff*(tl_LapT(i,j,k+1)- &
1206!^ & tl_LapT(i,j,k ))
1207!^
1208 adfac=cff*ad_dtdr(i,j,k2)
1209 ad_lapt(i,j,k )=ad_lapt(i,j,k )-adfac
1210 ad_lapt(i,j,k+1)=ad_lapt(i,j,k+1)+adfac
1211 ad_cff=ad_cff+(lapt(i,j,k+1)- &
1212 & lapt(i,j,k ))*ad_dtdr(i,j,k2)
1213 ad_dtdr(i,j,k2)=0.0_r8
1214#if defined TS_MIX_MAX_SLOPE
1215!^ tl_cff=cff*cff*tl_cff4
1216!^
1217 ad_cff4=ad_cff4+cff*cff*ad_cff
1218 ad_cff=0.0_r8
1219!^ tl_cff4=(0.5_r8+SIGN(0.5_r8,cff2-cff3))*tl_cff2+ &
1220!^ & (0.5_r8-SIGN(0.5_r8,cff2-cff3))*tl_cff3
1221!^
1222 ad_cff3=ad_cff3+ &
1223 & (0.5_r8-sign(0.5_r8,cff2-cff3))*ad_cff4
1224 ad_cff2=ad_cff2+ &
1225 & (0.5_r8+sign(0.5_r8,cff2-cff3))*ad_cff4
1226 ad_cff4=0.0_r8
1227!^ tl_cff3=(0.5_r8+SIGN(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
1228!^ & small))* &
1229!^ & (tl_pden(i,j,k)-tl_pden(i,j,k+1))
1230!^
1231 adfac=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
1232 & small))*ad_cff3
1233 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac
1234 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac
1235 ad_cff3=0.0_r8
1236!^ tl_cff2=0.25_r8*slope_max* &
1237!^ & ((tl_z_r(i,j,k+1)-tl_z_r(i,j,k))*cff1+ &
1238!^ & (z_r(i,j,k+1)-z_r(i,j,k))*tl_cff1)
1239!^
1240 adfac=0.25_r8*slope_max*ad_cff2
1241 adfac1=adfac*cff1
1242 ad_cff1=ad_cff1+(z_r(i,j,k+1)-z_r(i,j,k))*adfac
1243 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac1
1244 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac1
1245 ad_cff2=0.0_r8
1246 IF (cff1.ne.0.0_r8) THEN
1247!^ tl_cff1=(dRdx(i ,j,k2)*tl_dRdx(i ,j,k2)+ &
1248!^ & dRdx(i+1,j,k2)*tl_dRdx(i+1,j,k2)+ &
1249!^ & dRdx(i ,j,k1)*tl_dRdx(i ,j,k1)+ &
1250!^ & dRdx(i+1,j,k1)*tl_dRdx(i+1,j,k1)+ &
1251!^ & dRde(i,j ,k2)*tl_dRde(i,j ,k2)+ &
1252!^ & dRde(i,j+1,k2)*tl_dRde(i,j+1,k2)+ &
1253!^ & dRde(i,j ,k1)*tl_dRde(i,j ,k1)+ &
1254!^ & dRde(i,j+1,k1)*tl_dRde(i,j+1,k1))/cff1
1255!^
1256 adfac=ad_cff1/cff1
1257 ad_drdx(i ,j,k1)=ad_drdx(i ,j,k1)+ &
1258 & drdx(i ,j,k1)*adfac
1259 ad_drdx(i+1,j,k1)=ad_drdx(i+1,j,k1)+ &
1260 & drdx(i+1,j,k1)*adfac
1261 ad_drdx(i ,j,k2)=ad_drdx(i ,j,k2)+ &
1262 & drdx(i ,j,k2)*adfac
1263 ad_drdx(i+1,j,k2)=ad_drdx(i+1,j,k2)+ &
1264 & drdx(i+1,j,k2)*adfac
1265 ad_drde(i,j ,k2)=ad_drde(i,j ,k2)+ &
1266 & drde(i,j ,k2)*adfac
1267 ad_drde(i,j+1,k2)=ad_drde(i,j+1,k2)+ &
1268 & drde(i,j+1,k2)*adfac
1269 ad_drde(i,j ,k1)=ad_drde(i,j ,k1)+ &
1270 & drde(i,j ,k1)*adfac
1271 ad_drde(i,j+1,k1)=ad_drde(i,j+1,k1)+ &
1272 & drde(i,j+1,k1)*adfac
1273 ad_cff1=0.0_r8
1274 ELSE
1275!^ tl_cff1=0.0_r8
1276!^
1277 ad_cff1=0.0_r8
1278 END IF
1279#elif defined TS_MIX_MIN_STRAT
1280!^ tl_cff=cff*cff*tl_cff1
1281!^
1282 ad_cff1=ad_cff1+cff*cff*ad_cff
1283 ad_cff=0.0_r8
1284!^ tl_cff1=(0.5_r8+SIGN(0.5_r8, &
1285!^ & pden(i,j,k)-pden(i,j,k+1)- &
1286!^ & strat_min*(z_r(i,j,k+1)- &
1287!^ & z_r(i,j,k ))))* &
1288!^ & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
1289!^ & (0.5_r8-SIGN(0.5_r8, &
1290!^ & pden(i,j,k)-pden(i,j,k+1)- &
1291!^ & strat_min*(z_r(i,j,k+1)- &
1292!^ & z_r(i,j,k ))))* &
1293!^ & (strat_min*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k )))
1294!^
1295 adfac1=(0.5_r8+sign(0.5_r8, &
1296 & pden(i,j,k)-pden(i,j,k+1)- &
1297 & strat_min*(z_r(i,j,k+1)- &
1298 & z_r(i,j,k ))))* &
1299 & ad_cff1
1300 adfac2=(0.5_r8-sign(0.5_r8, &
1301 & pden(i,j,k)-pden(i,j,k+1)- &
1302 & strat_min*(z_r(i,j,k+1)- &
1303 & z_r(i,j,k ))))* &
1304 & strat_min*ad_cff1
1305 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac1
1306 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac1
1307 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac2
1308 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac2
1309 ad_cff1=0.0_r8
1310#else
1311!^ tl_cff=cff*cff*tl_cff1
1312!^
1313 ad_cff1=ad_cff1+cff*cff*ad_cff
1314 ad_cff=0.0_r8
1315!^ tl_cff1=(0.5_r8+SIGN(0.5_r8, &
1316!^ & pden(i,j,k)-pden(i,j,k+1)-eps))* &
1317!^ & (tl_pden(i,j,k)-tl_pden(i,j,k+1))
1318!^
1319 adfac=(0.5_r8+sign(0.5_r8, &
1320 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
1321 & ad_cff1
1322 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac
1323 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac
1324 ad_cff1=0.0_r8
1325#endif
1326 END DO
1327 END DO
1328 END IF
1329 IF (k.lt.n(ng)) THEN
1330 DO j=jstr,jend+1
1331 DO i=istr,iend
1332 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
1333#ifdef MASKING
1334 cff=cff*vmask(i,j)
1335#endif
1336#ifdef WET_DRY_NOT_YET
1337 cff=cff*vmask_wet(i,j)
1338#endif
1339!^ tl_dTde(i,j,k2)=cff*(tl_LapT(i,j ,k+1)- &
1340!^ & tl_LapT(i,j-1,k+1))
1341!^
1342 adfac=cff*ad_dtde(i,j,k2)
1343 ad_lapt(i,j-1,k+1)=ad_lapt(i,j-1,k+1)-adfac
1344 ad_lapt(i,j ,k+1)=ad_lapt(i,j ,k+1)+adfac
1345 ad_dtde(i,j,k2)=0.0_r8
1346!^ tl_dRde(i,j,k2)=cff*(tl_pden(i,j ,k+1)- &
1347!^ & tl_pden(i,j-1,k+1))
1348!^
1349 adfac=cff*ad_drde(i,j,k2)
1350 ad_pden(i,j-1,k+1)=ad_pden(i,j-1,k+1)-adfac
1351 ad_pden(i,j ,k+1)=ad_pden(i,j ,k+1)+adfac
1352 ad_drde(i,j,k2)=0.0_r8
1353 END DO
1354 END DO
1355 DO j=jstr,jend
1356 DO i=istr,iend+1
1357 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
1358#ifdef MASKING
1359 cff=cff*umask(i,j)
1360#endif
1361#ifdef WET_DRY_NOT_YET
1362 cff=cff*umask_wet(i,j)
1363#endif
1364!^ tl_dTdx(i,j,k2)=cff*(tl_LapT(i ,j,k+1)- &
1365!^ & tl_LapT(i-1,j,k+1))
1366!^
1367 adfac=cff*ad_dtdx(i,j,k2)
1368 ad_lapt(i-1,j,k+1)=ad_lapt(i-1,j,k+1)-adfac
1369 ad_lapt(i ,j,k+1)=ad_lapt(i ,j,k+1)+adfac
1370 ad_dtdx(i,j,k2)=0.0_r8
1371!^ tl_dRdx(i,j,k2)=cff*(tl_pden(i ,j,k+1)- &
1372!^ & tl_pden(i-1,j,k+1))
1373!^
1374 adfac=cff*ad_drdx(i,j,k2)
1375 ad_pden(i-1,j,k+1)=ad_pden(i-1,j,k+1)-adfac
1376 ad_pden(i ,j,k+1)=ad_pden(i ,j,k+1)+adfac
1377 ad_drdx(i,j,k2)=0.0_r8
1378 END DO
1379 END DO
1380 END IF
1381!
1382! Compute new storage recursive indices.
1383!
1384 kt=k2
1385 k2=k1
1386 k1=kt
1387 END DO k_loop2
1388!
1389! Apply adjoint boundary conditions (except periodic; closed or
1390! gradient) to the first harmonic operator.
1391!
1392 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
1393 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
1394 IF (domain(ng)%NorthEast_Corner(tile)) THEN
1395 DO k=1,n(ng)
1396!^ tl_LapT(Iend+1,Jend+1,k)=0.5_r8* &
1397!^ & (tl_LapT(Iend ,Jend+1,k)+ &
1398!^ & tl_LapT(Iend+1,Jend ,k))
1399!^
1400 adfac=0.5_r8*ad_lapt(iend+1,jend+1,k)
1401 ad_lapt(iend+1,jend ,k)=ad_lapt(iend+1,jend ,k)+adfac
1402 ad_lapt(iend ,jend+1,k)=ad_lapt(iend ,jend+1,k)+adfac
1403 ad_lapt(iend+1,jend+1,k)=0.0_r8
1404 END DO
1405 END IF
1406 END IF
1407
1408 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng).or. &
1409 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
1410 IF (domain(ng)%NorthWest_Corner(tile)) THEN
1411 DO k=1,n(ng)
1412!^ tl_LapT(Istr-1,Jend+1,k)=0.5_r8* &
1413!^ & (tl_LapT(Istr ,Jend+1,k)+ &
1414!^ & tl_LapT(Istr-1,Jend ,k))
1415!^
1416 adfac=0.5_r8*ad_lapt(istr-1,jend+1,k)
1417 ad_lapt(istr-1,jend ,k)=ad_lapt(istr-1,jend ,k)+adfac
1418 ad_lapt(istr ,jend+1,k)=ad_lapt(istr ,jend+1,k)+adfac
1419 ad_lapt(istr-1,jend+1,k)=0.0_r8
1420 END DO
1421 END IF
1422 END IF
1423
1424 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
1425 & compositegrid(ieast ,ng).or.ewperiodic(ng))) THEN
1426 IF (domain(ng)%SouthEast_Corner(tile)) THEN
1427 DO k=1,n(ng)
1428!^ tl_LapT(Iend+1,Jstr-1,k)=0.5_r8* &
1429!^ & (tl_LapT(Iend ,Jstr-1,k)+ &
1430!^ & tl_LapT(Iend+1,Jstr ,k))
1431!^
1432 adfac=0.5_r8*ad_lapt(iend+1,jstr-1,k)
1433 ad_lapt(iend ,jstr-1,k)=ad_lapt(iend ,jstr-1,k)+adfac
1434 ad_lapt(iend+1,jstr ,k)=ad_lapt(iend+1,jstr ,k)+adfac
1435 ad_lapt(iend+1,jstr-1,k)=0.0_r8
1436 END DO
1437 END IF
1438 END IF
1439
1440 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng).or. &
1441 & compositegrid(iwest ,ng).or.ewperiodic(ng))) THEN
1442 IF (domain(ng)%SouthWest_Corner(tile)) THEN
1443 DO k=1,n(ng)
1444!^ tl_LapT(Istr-1,Jstr-1,k)=0.5_r8* &
1445!^ & (tl_LapT(Istr ,Jstr-1,k)+ &
1446!^ tl_LapT(Istr-1,Jstr ,k))
1447!^
1448 adfac=0.5_r8*ad_lapt(istr-1,jstr-1,k)
1449 ad_lapt(istr ,jstr-1,k)=ad_lapt(istr ,jstr-1,k)+adfac
1450 ad_lapt(istr-1,jstr ,k)=ad_lapt(istr-1,jstr ,k)+adfac
1451 ad_lapt(istr-1,jstr-1,k)=0.0_r8
1452 END DO
1453 END IF
1454 END IF
1455!
1456 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
1457 IF (domain(ng)%Northern_Edge(tile)) THEN
1458 IF (ad_lbc(inorth,istvar(itrc),ng)%closed) THEN
1459 DO k=1,n(ng)
1460 DO i=imin,imax
1461!^ tl_LapT(i,Jend+1,k)=0.0_r8
1462!^
1463 ad_lapt(i,jend+1,k)=0.0_r8
1464 END DO
1465 END DO
1466 ELSE
1467 DO k=1,n(ng)
1468 DO i=imin,imax
1469!^ tl_LapT(i,Jend+1,k)=tl_LapT(i,Jend,k)
1470!^
1471 ad_lapt(i,jend,k)=ad_lapt(i,jend,k)+ &
1472 & ad_lapt(i,jend+1,k)
1473 ad_lapt(i,jend+1,k)=0.0_r8
1474 END DO
1475 END DO
1476 END IF
1477 END IF
1478 END IF
1479!
1480 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
1481 IF (domain(ng)%Southern_Edge(tile)) THEN
1482 IF (ad_lbc(isouth,istvar(itrc),ng)%closed) THEN
1483 DO k=1,n(ng)
1484 DO i=imin,imax
1485!^ tl_LapT(i,Jstr-1,k)=0.0_r8
1486!^
1487 ad_lapt(i,jstr-1,k)=0.0_r8
1488 END DO
1489 END DO
1490 ELSE
1491 DO k=1,n(ng)
1492 DO i=imin,imax
1493!^ tl_LapT(i,Jstr-1,k)=tl_LapT(i,Jstr,k)
1494!^
1495 ad_lapt(i,jstr,k)=ad_lapt(i,jstr,k)+ &
1496 & ad_lapt(i,jstr-1,k)
1497 ad_lapt(i,jstr-1,k)=0.0_r8
1498 END DO
1499 END DO
1500 END IF
1501 END IF
1502 END IF
1503!
1504 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
1505 IF (domain(ng)%Eastern_Edge(tile)) THEN
1506 IF (ad_lbc(ieast,istvar(itrc),ng)%closed) THEN
1507 DO k=1,n(ng)
1508 DO j=jmin,jmax
1509!^ tl_LapT(Iend+1,j,k)=0.0_r8
1510!^
1511 ad_lapt(iend+1,j,k)=0.0_r8
1512 END DO
1513 END DO
1514 ELSE
1515 DO k=1,n(ng)
1516 DO j=jmin,jmax
1517!^ tl_LapT(Iend+1,j,k)=tl_LapT(Iend,j,k)
1518!^
1519 ad_lapt(iend,j,k)=ad_lapt(iend,j,k)+ &
1520 & ad_lapt(iend+1,j,k)
1521 ad_lapt(iend+1,j,k)=0.0_r8
1522 END DO
1523 END DO
1524 END IF
1525 END IF
1526 END IF
1527!
1528 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
1529 IF (domain(ng)%Western_Edge(tile)) THEN
1530 IF (ad_lbc(iwest,istvar(itrc),ng)%closed) THEN
1531 DO k=1,n(ng)
1532 DO j=jmin,jmax
1533!^ tl_LapT(Istr-1,j,k)=0.0_r8
1534!^
1535 ad_lapt(istr-1,j,k)=0.0_r8
1536 END DO
1537 END DO
1538 ELSE
1539 DO k=1,n(ng)
1540 DO j=jmin,jmax
1541!^ tl_LapT(Istr-1,j,k)=tl_LapT(Istr,j,k)
1542!^
1543 ad_lapt(istr,j,k)=ad_lapt(istr,j,k)+ &
1544 & ad_lapt(istr-1,j,k)
1545 ad_lapt(istr-1,j,k)=0.0_r8
1546 END DO
1547 END DO
1548 END IF
1549 END IF
1550 END IF
1551!
1552!-----------------------------------------------------------------------
1553! Compute first adjoint harmonic operator, without mixing coefficient.
1554! Multiply by the metrics of the second harmonic operator.
1555!-----------------------------------------------------------------------
1556!
1557! Compute adjoint of starting recursive indices k1 and k2.
1558!
1559 k1=2
1560 k2=1
1561 DO k=0,n(ng)
1562!!
1563!! Note: The following code is equivalent to
1564!!
1565!! kt=k1
1566!! k1=k2
1567!! k2=kt
1568!!
1569!! We use the adjoint of above code.
1570!!
1571 k1=k2
1572 k2=3-k1
1573 END DO
1574!
1575! Compute required basic state fields. Need to look forward in "kk"
1576! index.
1577!
1578 k_loop3: DO k=n(ng),0,-1
1579 k2b=1
1580 DO kk=0,k
1581 k1b=k2b
1582 k2b=3-k1b
1583 IF (kk.lt.n(ng)) THEN
1584 DO j=jmin,jmax
1585 DO i=imin,imax+1
1586 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
1587#ifdef MASKING
1588 cff=cff*umask(i,j)
1589#endif
1590#ifdef WET_DRY_NOT_YET
1591 cff=cff*umask_wet(i,j)
1592#endif
1593 drdx(i,j,k2b)=cff*(pden(i ,j,kk+1)- &
1594 & pden(i-1,j,kk+1))
1595#if defined TS_MIX_STABILITY
1596 dtdx(i,j,k2b)=cff*(0.75_r8*(t(i ,j,kk+1,nrhs,itrc)- &
1597 & t(i-1,j,kk+1,nrhs,itrc))+ &
1598 & 0.25_r8*(t(i ,j,kk+1,nstp,itrc)- &
1599 & t(i-1,j,kk+1,nstp,itrc)))
1600#elif defined TS_MIX_CLIMA
1601 IF (ltracerclm(itrc,ng)) THEN
1602 dtdx(i,j,k2b)=cff*((t(i ,j,kk+1,nrhs,itrc)- &
1603 & tclm(i ,j,kk+1,itrc))- &
1604 & (t(i-1,j,kk+1,nrhs,itrc)- &
1605 & tclm(i-1,j,kk+1,itrc)))
1606 ELSE
1607 dtdx(i,j,k2b)=cff*(t(i ,j,kk+1,nrhs,itrc)- &
1608 & t(i-1,j,kk+1,nrhs,itrc))
1609 END IF
1610#else
1611 dtdx(i,j,k2b)=cff*(t(i ,j,kk+1,nrhs,itrc)- &
1612 & t(i-1,j,kk+1,nrhs,itrc))
1613#endif
1614 END DO
1615 END DO
1616 IF (kk.eq.0) THEN
1617 DO j=jmin,jmax
1618 DO i=imin,imax+1
1619 drdx(i,j,k1b)=0.0_r8
1620 dtdx(i,j,k1b)=0.0_r8
1621 END DO
1622 END DO
1623 END IF
1624 DO j=jmin,jmax+1
1625 DO i=imin,imax
1626 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
1627#ifdef MASKING
1628 cff=cff*vmask(i,j)
1629#endif
1630#ifdef WET_DRY_NOT_YET
1631 cff=cff*vmask_wet(i,j)
1632#endif
1633 drde(i,j,k2b)=cff*(pden(i,j ,kk+1)- &
1634 & pden(i,j-1,kk+1))
1635#if defined TS_MIX_STABILITY
1636 dtde(i,j,k2b)=cff*(0.75_r8*(t(i,j ,kk+1,nrhs,itrc)- &
1637 & t(i,j-1,kk+1,nrhs,itrc))+ &
1638 & 0.25_r8*(t(i,j ,kk+1,nstp,itrc)- &
1639 & t(i,j-1,kk+1,nstp,itrc)))
1640#elif defined TS_MIX_CLIMA
1641 IF (ltracerclm(itrc,ng)) THEN
1642 dtde(i,j,k2b)=cff*((t(i,j ,kk+1,nrhs,itrc)- &
1643 & tclm(i,j ,kk+1,itrc))- &
1644 & (t(i,j-1,kk+1,nrhs,itrc)- &
1645 & tclm(i,j-1,kk+1,itrc)))
1646 ELSE
1647 dtde(i,j,k2b)=cff*(t(i,j ,kk+1,nrhs,itrc)- &
1648 & t(i,j-1,kk+1,nrhs,itrc))
1649 END IF
1650#else
1651 dtde(i,j,k2b)=cff*(t(i,j ,kk+1,nrhs,itrc)- &
1652 & t(i,j-1,kk+1,nrhs,itrc))
1653#endif
1654 END DO
1655 END DO
1656 IF (kk.eq.0) THEN
1657 DO j=jmin,jmax+1
1658 DO i=imin,imax
1659 drde(i,j,k1b)=0.0_r8
1660 dtde(i,j,k1b)=0.0_r8
1661 END DO
1662 END DO
1663 END IF
1664 END IF
1665 IF ((kk.eq.0).or.(kk.eq.n(ng))) THEN
1666 DO j=-1+jmin,jmax+1
1667 DO i=-1+imin,imax+1
1668 dtdr(i,j,k2b)=0.0_r8
1669 fs(i,j,k2b)=0.0_r8
1670 END DO
1671 END DO
1672 IF (kk.eq.0) THEN
1673 DO j=-1+jmin,jmax+1
1674 DO i=-1+imin,imax+1
1675 dtdr(i,j,k1b)=0.0_r8
1676 fs(i,j,k1b)=0.0_r8
1677 END DO
1678 END DO
1679 END IF
1680 ELSE
1681 DO j=-1+jmin,jmax+1
1682 DO i=-1+imin,imax+1
1683#if defined TS_MIX_MAX_SLOPE
1684 cff1=sqrt(drdx(i,j,k2b)**2+drdx(i+1,j,k2b)**2+ &
1685 & drdx(i,j,k1b)**2+drdx(i+1,j,k1b)**2+ &
1686 & drde(i,j,k2b)**2+drde(i,j+1,k2b)**2+ &
1687 & drde(i,j,k1b)**2+drde(i,j+1,k1b)**2)
1688 cff2=0.25_r8*slope_max* &
1689 & (z_r(i,j,kk+1)-z_r(i,j,kk))*cff1
1690 cff3=max(pden(i,j,kk)-pden(i,j,kk+1),small)
1691 cff4=max(cff2,cff3)
1692 cff=-1.0_r8/cff4
1693#elif defined TS_MIX_MIN_STRAT
1694 cff1=max(pden(i,j,kk)-pden(i,j,kk+1), &
1695 & strat_min*(z_r(i,j,kk+1)-z_r(i,j,kk)))
1696 cff=-1.0_r8/cff1
1697#else
1698 cff1=max(pden(i,j,kk)-pden(i,j,kk+1),eps)
1699 cff=-1.0_r8/cff1
1700#endif
1701#if defined TS_MIX_STABILITY
1702 dtdr(i,j,k2b)=cff*(0.75_r8*(t(i,j,kk+1,nrhs,itrc)- &
1703 & t(i,j,kk ,nrhs,itrc))+ &
1704 & 0.25_r8*(t(i,j,kk+1,nstp,itrc)- &
1705 & t(i,j,kk ,nstp,itrc)))
1706#elif defined TS_MIX_CLIMA
1707 IF (ltracerclm(itrc,ng)) THEN
1708 dtdr(i,j,k2b)=cff*((t(i,j,kk+1,nrhs,itrc)- &
1709 & tclm(i,j,kk+1,itrc))- &
1710 & (t(i,j,kk ,nrhs,itrc)- &
1711 & tclm(i,j,kk ,itrc)))
1712 ELSE
1713 dtdr(i,j,k2b)=cff*(t(i,j,kk+1,nrhs,itrc)- &
1714 & t(i,j,kk ,nrhs,itrc))
1715 END IF
1716#else
1717 dtdr(i,j,k2b)=cff*(t(i,j,kk+1,nrhs,itrc)- &
1718 & t(i,j,kk ,nrhs,itrc))
1719#endif
1720 fs(i,j,k2b)=cff*(z_r(i,j,kk+1)- &
1721 & z_r(i,j,kk ))
1722 END DO
1723 END DO
1724 END IF
1725 IF (kk.gt.0) THEN
1726 DO j=jmin,jmax
1727 DO i=imin,imax+1
1728#ifdef DIFF_3DCOEF
1729# ifdef TS_U3ADV_SPLIT
1730 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
1731# else
1732 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
1733 & on_u(i,j)
1734# endif
1735#else
1736 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
1737 & on_u(i,j)
1738#endif
1739 fx(i,j)=cff* &
1740 & (hz(i,j,kk)+hz(i-1,j,kk))* &
1741 & (dtdx(i ,j,k1b)- &
1742 & 0.5_r8*(max(drdx(i,j,k1b),0.0_r8)* &
1743 & (dtdr(i-1,j,k1b)+ &
1744 & dtdr(i ,j,k2b))+ &
1745 & min(drdx(i,j,k1b),0.0_r8)* &
1746 & (dtdr(i-1,j,k2b)+ &
1747 & dtdr(i ,j,k1b))))
1748 END DO
1749 END DO
1750 DO j=jmin,jmax+1
1751 DO i=imin,imax
1752#ifdef DIFF_3DCOEF
1753# ifdef TS_U3ADV_SPLIT
1754 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
1755# else
1756 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
1757 & om_v(i,j)
1758# endif
1759#else
1760 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
1761 & om_v(i,j)
1762#endif
1763 fe(i,j)=cff* &
1764 & (hz(i,j,kk)+hz(i,j-1,kk))* &
1765 & (dtde(i,j,k1b)- &
1766 & 0.5_r8*(max(drde(i,j,k1b),0.0_r8)* &
1767 & (dtdr(i,j-1,k1b)+ &
1768 & dtdr(i,j ,k2b))+ &
1769 & min(drde(i,j,k1b),0.0_r8)* &
1770 & (dtdr(i,j-1,k2b)+ &
1771 & dtdr(i,j ,k1b))))
1772 END DO
1773 END DO
1774 IF (kk.lt.n(ng)) THEN
1775 DO j=jmin,jmax
1776 DO i=imin,imax
1777#ifdef DIFF_3DCOEF
1778# ifdef TS_U3ADV_SPLIT
1779 difx=0.125_r8* &
1780 & (diff3d_u(i,j,kk )+diff3d_u(i+1,j,kk )+ &
1781 & diff3d_u(i,j,kk+1)+diff3d_u(i+1,j,kk+1))
1782 dife=0.125_r8* &
1783 & (diff3d_v(i,j,kk )+diff3d_v(i,j+1,kk )+ &
1784 & diff3d_v(i,j,kk+1)+diff3d_v(i,j+1,kk+1))
1785# else
1786 difx=0.5_r8*diff3d_r(i,j,kk)
1787 dife=difx
1788# endif
1789#else
1790 difx=0.5_r8*diff4(i,j,itrc)
1791 dife=difx
1792#endif
1793 cff1=max(drdx(i ,j,k1b),0.0_r8)
1794 cff2=max(drdx(i+1,j,k2b),0.0_r8)
1795 cff3=min(drdx(i ,j,k2b),0.0_r8)
1796 cff4=min(drdx(i+1,j,k1b),0.0_r8)
1797 cff=difx* &
1798 & (cff1*(cff1*dtdr(i,j,k2b)-dtdx(i ,j,k1b))+ &
1799 & cff2*(cff2*dtdr(i,j,k2b)-dtdx(i+1,j,k2b))+ &
1800 & cff3*(cff3*dtdr(i,j,k2b)-dtdx(i ,j,k2b))+ &
1801 & cff4*(cff4*dtdr(i,j,k2b)-dtdx(i+1,j,k1b)))
1802 cff1=max(drde(i,j ,k1b),0.0_r8)
1803 cff2=max(drde(i,j+1,k2b),0.0_r8)
1804 cff3=min(drde(i,j ,k2b),0.0_r8)
1805 cff4=min(drde(i,j+1,k1b),0.0_r8)
1806 cff=cff+ &
1807 & dife* &
1808 & (cff1*(cff1*dtdr(i,j,k2b)-dtde(i,j ,k1b))+ &
1809 & cff2*(cff2*dtdr(i,j,k2b)-dtde(i,j+1,k2b))+ &
1810 & cff3*(cff3*dtdr(i,j,k2b)-dtde(i,j ,k2b))+ &
1811 & cff4*(cff4*dtdr(i,j,k2b)-dtde(i,j+1,k1b)))
1812 fs1(i,j,k2b)=fs(i,j,k2b) ! recursive
1813 fs(i,j,k2b)=cff*fs(i,j,k2b)
1814 END DO
1815 END DO
1816 IF (kk.eq.0) THEN
1817 DO j=jmin,jmax
1818 DO i=imin,imax
1819 fs1(i,j,k1b)=0.0_r8
1820 fs(i,j,k1b)=0.0_r8
1821 END DO
1822 END DO
1823 END IF
1824 END IF
1825 END IF
1826 END DO
1827!
1828! Compute adjoint first harmonic operator, without mixing coefficient.
1829! Multiply by the metrics of the second harmonic operator. Save
1830! into work array "LapT".
1831!
1832 IF (k.gt.0) THEN
1833 DO j=jmin,jmax
1834 DO i=imin,imax
1835 cff=pm(i,j)*pn(i,j)
1836 cff1=1.0_r8/hz(i,j,k)
1837!^ tl_LapT(i,j,k)=tl_cff1*(cff* &
1838!^ & (FX(i+1,j)-FX(i,j)+ &
1839!^ & FE(i,j+1)-FE(i,j))+ &
1840!^ & (FS(i,j,k2)-FS(i,j,k1)))+ &
1841!^ & cff1*(cff* &
1842!^ & (tl_FX(i+1,j)-tl_FX(i,j)+ &
1843!^ & tl_FE(i,j+1)-tl_FE(i,j))+ &
1844!^ & (tl_FS(i,j,k2)-tl_FS(i,j,k1)))
1845!^
1846 adfac=cff1*ad_lapt(i,j,k)
1847 adfac1=adfac*cff
1848 ad_fs(i,j,k1)=ad_fs(i,j,k1)-adfac
1849 ad_fs(i,j,k2)=ad_fs(i,j,k2)+adfac
1850 ad_fe(i,j )=ad_fe(i,j )-adfac1
1851 ad_fe(i,j+1)=ad_fe(i,j+1)+adfac1
1852 ad_fx(i ,j)=ad_fx(i ,j)-adfac1
1853 ad_fx(i+1,j)=ad_fx(i+1,j)+adfac1
1854 ad_cff1=ad_cff1+(cff* &
1855 & (fx(i+1,j)-fx(i,j)+ &
1856 & fe(i,j+1)-fe(i,j))+ &
1857 & (fs(i,j,k2)-fs(i,j,k1)))* &
1858 & ad_lapt(i,j,k)
1859 ad_lapt(i,j,k)=0.0_r8
1860!^ tl_cff1=-cff1*cff1*tl_Hz(i,j,k)
1861!^
1862 ad_hz(i,j,k)=ad_hz(i,j,k)-cff1*cff1*ad_cff1
1863 ad_cff1=0.0_r8
1864 END DO
1865 END DO
1866 IF (k.lt.n(ng)) THEN
1867 DO j=jmin,jmax
1868 DO i=imin,imax
1869#ifdef DIFF_3DCOEF
1870# ifdef TS_U3ADV_SPLIT
1871 difx=0.125_r8*(diff3d_u(i,j,k )+diff3d_u(i+1,j,k )+ &
1872 & diff3d_u(i,j,k+1)+diff3d_u(i+1,j,k+1))
1873 dife=0.125_r8*(diff3d_v(i,j,k )+diff3d_v(i,j+1,k )+ &
1874 & diff3d_v(i,j,k+1)+diff3d_v(i,j+1,k+1))
1875# else
1876 difx=0.5_r8*diff3d_r(i,j,k)
1877 dife=difx
1878# endif
1879#else
1880 difx=0.5_r8*diff4(i,j,itrc)
1881 dife=difx
1882#endif
1883 cff1=max(drdx(i ,j,k1),0.0_r8)
1884 cff2=max(drdx(i+1,j,k2),0.0_r8)
1885 cff3=min(drdx(i ,j,k2),0.0_r8)
1886 cff4=min(drdx(i+1,j,k1),0.0_r8)
1887 cff=difx* &
1888 & (cff1*(cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))+ &
1889 & cff2*(cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))+ &
1890 & cff3*(cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))+ &
1891 & cff4*(cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1)))
1892 cff1=max(drde(i,j ,k1),0.0_r8)
1893 cff2=max(drde(i,j+1,k2),0.0_r8)
1894 cff3=min(drde(i,j ,k2),0.0_r8)
1895 cff4=min(drde(i,j+1,k1),0.0_r8)
1896 cff=cff+ &
1897 & dife* &
1898 & (cff1*(cff1*dtdr(i,j,k2)-dtde(i,j ,k1))+ &
1899 & cff2*(cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))+ &
1900 & cff3*(cff3*dtdr(i,j,k2)-dtde(i,j ,k2))+ &
1901 & cff4*(cff4*dtdr(i,j,k2)-dtde(i,j+1,k1)))
1902!^ tl_FS(i,j,k2)=tl_cff*FS(i,j,k2)+ &
1903!^ & cff*tl_FS(i,j,k2) ! recursive
1904!^ ! use FS1
1905!^
1906 ad_fs(i,j,k2)=cff*ad_fs(i,j,k2)
1907 ad_cff=ad_cff+fs1(i,j,k2)*ad_fs(i,j,k2)
1908!^ tl_cff=tl_cff+ &
1909!^ & dife* &
1910!^ & (tl_cff1*(cff1*dTdr(i,j,k2)- &
1911!^ & dTde(i,j ,k1))+ &
1912!^ & tl_cff2*(cff2*dTdr(i,j,k2)- &
1913!^ & dTde(i,j+1,k2))+ &
1914!^ & tl_cff3*(cff3*dTdr(i,j,k2)- &
1915!^ & dTde(i,j ,k2))+ &
1916!^ & tl_cff4*(cff4*dTdr(i,j,k2)- &
1917!^ & dTde(i,j+1,k1))+ &
1918!^ & cff1*(tl_cff1*dTdr(i,j,k2)+ &
1919!^ & cff1*tl_dTdr(i,j,k2)- &
1920!^ & tl_dTde(i,j ,k1))+ &
1921!^ & cff2*(tl_cff2*dTdr(i,j,k2)+ &
1922!^ & cff2*tl_dTdr(i,j,k2)- &
1923!^ & tl_dTde(i,j+1,k2))+ &
1924!^ & cff3*(tl_cff3*dTdr(i,j,k2)+ &
1925!^ & cff3*tl_dTdr(i,j,k2)- &
1926!^ & tl_dTde(i,j ,k2))+ &
1927!^ & cff4*(tl_cff4*dTdr(i,j,k2)+ &
1928!^ & cff4*tl_dTdr(i,j,k2)- &
1929!^ & tl_dTde(i,j+1,k1)))
1930!^
1931 adfac*dife*ad_cff
1932 ad_cff1=ad_cff1+ &
1933 & (2.0_r8*cff1*dtdr(i,j,k2)-dtde(i,j ,k1))* &
1934 & adfac
1935 ad_cff2=ad_cff2+ &
1936 & (2.0_r8*cff2*dtdr(i,j,k2)-dtde(i,j+1,k2))* &
1937 & adfac
1938 ad_cff3=ad_cff3+ &
1939 & (2.0_r8*cff3*dtdr(i,j,k2)-dtde(i,j ,k2))* &
1940 & adfac
1941 ad_cff4=ad_cff4+ &
1942 & (2.0_r8*cff4*dtdr(i,j,k2)-dtde(i,j+1,k1))* &
1943 & adfac
1944 ad_dtdr(i,j,k2)=ad_dtdr(i,j,k2)+ &
1945 & (cff1*cff1+ &
1946 & cff2*cff2+ &
1947 & cff3*cff3+ &
1948 & cff4*cff4)*adfac
1949 ad_dtde(i,j ,k1)=ad_dtde(i,j ,k1)-cff1*adfac
1950 ad_dtde(i,j+1,k2)=ad_dtde(i,j+1,k2)-cff2*adfac
1951 ad_dtde(i,j ,k2)=ad_dtde(i,j ,k2)-cff3*adfac
1952 ad_dtde(i,j+1,k1)=ad_dtde(i,j+1,k1)-cff4*adfac
1953!^ tl_cff4=(0.5_r8+SIGN(0.5_r8,-dRde(i,j+1,k1)))* &
1954!^ & tl_dRde(i,j+1,k1)
1955!^
1956 ad_drde(i,j+1,k1)=ad_drde(i,j+1,k1)+ &
1957 & (0.5_r8+sign(0.5_r8, &
1958 & -drde(i,j+1,k1)))* &
1959 & ad_cff4
1960 ad_cff4=0.0_r8
1961!^ tl_cff3=(0.5_r8+SIGN(0.5_r8,-dRde(i,j ,k2)))* &
1962!^ & tl_dRde(i,j ,k2)
1963!^
1964 ad_drde(i,j ,k2)=ad_drde(i,j ,k2)+ &
1965 & (0.5_r8+sign(0.5_r8, &
1966 & -drde(i,j ,k2)))* &
1967 & ad_cff3
1968 ad_cff3=0.0_r8
1969!^ tl_cff2=(0.5_r8+SIGN(0.5_r8, dRde(i,j+1,k2)))* &
1970!^ & tl_dRde(i,j+1,k2)
1971!^
1972 ad_drde(i,j+1,k2)=ad_drde(i,j+1,k2)+ &
1973 & (0.5_r8+sign(0.5_r8, &
1974 & drde(i,j+1,k2)))* &
1975 & ad_cff2
1976 ad_cff2=0.0_r8
1977!^ tl_cff1=(0.5_r8+SIGN(0.5_r8, dRde(i,j ,k1)))* &
1978!^ & tl_dRde(i,j ,k1)
1979!^
1980 ad_drde(i ,j,k1)=ad_drde(i ,j,k1)+ &
1981 & (0.5_r8+sign(0.5_r8, &
1982 & drde(i ,j,k1)))* &
1983 & ad_cff1
1984 ad_cff1=0.0_r8
1985 cff1=max(drdx(i ,j,k1),0.0_r8)
1986 cff2=max(drdx(i+1,j,k2),0.0_r8)
1987 cff3=min(drdx(i ,j,k2),0.0_r8)
1988 cff4=min(drdx(i+1,j,k1),0.0_r8)
1989!^ tl_cff=difx* &
1990!^ & (tl_cff1*(cff1*dTdr(i ,j,k2)- &
1991!^ & dTdx(i ,j,k1))+ &
1992!^ & tl_cff2*(cff2*dTdr(i,j,k2)- &
1993!^ & dTdx(i+1,j,k2))+ &
1994!^ & tl_cff3*(cff3*dTdr(i,j,k2)- &
1995!^ & dTdx(i ,j,k2))+ &
1996!^ & tl_cff4*(cff4*dTdr(i,j,k2)- &
1997!^ & dTdx(i+1,j,k1))+ &
1998!^ & cff1*(tl_cff1*dTdr(i,j,k2)+ &
1999!^ & cff1*tl_dTdr(i,j,k2)- &
2000!^ & tl_dTdx(i ,j,k1))+ &
2001!^ & cff2*(tl_cff2*dTdr(i,j,k2)+ &
2002!^ & cff2*tl_dTdr(i,j,k2)- &
2003!^ & tl_dTdx(i+1,j,k2))+ &
2004!^ & cff3*(tl_cff3*dTdr(i,j,k2)+ &
2005!^ & cff3*tl_dTdr(i,j,k2)- &
2006!^ & tl_dTdx(i ,j,k2))+ &
2007!^ & cff4*(tl_cff4*dTdr(i,j,k2)+ &
2008!^ & cff4*tl_dTdr(i,j,k2)- &
2009!^ & tl_dTdx(i+1,j,k1)))
2010!^
2011 adfac=difx*ad_cff
2012 ad_cff1=ad_cff1+ &
2013 & (2.0_r8*cff1*dtdr(i,j,k2)-dtdx(i ,j,k1))* &
2014 & adfac
2015 ad_cff2=ad_cff2+ &
2016 & (2.0_r8*cff2*dtdr(i,j,k2)-dtdx(i+1,j,k2))* &
2017 & adfac
2018 ad_cff3=ad_cff3+ &
2019 & (2.0_r8*cff3*dtdr(i,j,k2)-dtdx(i ,j,k2))* &
2020 & adfac
2021 ad_cff4=ad_cff4+ &
2022 & (2.0_r8*cff4*dtdr(i,j,k2)-dtdx(i+1,j,k1))* &
2023 & adfac
2024 ad_dtdr(i,j,k2)=ad_dtdr(i,j,k2)+ &
2025 & (cff1*cff1+ &
2026 & cff2*cff2+ &
2027 & cff3*cff3+ &
2028 & cff4*cff4)*adfac
2029 ad_dtdx(i ,j,k1)=ad_dtdx(i ,j,k1)-cff1*adfac
2030 ad_dtdx(i+1,j,k2)=ad_dtdx(i+1,j,k2)-cff2*adfac
2031 ad_dtdx(i ,j,k2)=ad_dtdx(i ,j,k2)-cff3*adfac
2032 ad_dtdx(i+1,j,k1)=ad_dtdx(i+1,j,k1)-cff4*adfac
2033 ad_cff=0.0_r8
2034!^ tl_cff4=(0.5_r8+SIGN(0.5_r8,-dRdx(i+1,j,k1)))* &
2035!^ & tl_dRdx(i+1,j,k1)
2036!^
2037 ad_drdx(i+1,j,k1)=ad_drdx(i+1,j,k1)+ &
2038 & (0.5_r8+sign(0.5_r8, &
2039 & -drdx(i+1,j,k1)))* &
2040 & ad_cff4
2041 ad_cff4=0.0_r8
2042!^ tl_cff3=(0.5_r8+SIGN(0.5_r8,-dRdx(i ,j,k2)))* &
2043!^ & tl_dRdx(i ,j,k2)
2044!^
2045 ad_drdx(i ,j,k2)=ad_drdx(i ,j,k2)+ &
2046 & (0.5_r8+sign(0.5_r8, &
2047 & -drdx(i ,j,k2)))* &
2048 & ad_cff3
2049 ad_cff3=0.0_r8
2050!^ tl_cff2=(0.5_r8+SIGN(0.5_r8, dRdx(i+1,j,k2)))* &
2051!^ & tl_dRdx(i+1,j,k2)
2052!^
2053 ad_drdx(i+1,j,k2)=ad_drdx(i+1,j,k2)+ &
2054 & (0.5_r8+sign(0.5_r8, &
2055 & drdx(i+1,j,k2)))* &
2056 & ad_cff2
2057 ad_cff2=0.0_r8
2058!^ tl_cff1=(0.5_r8+SIGN(0.5_r8, dRdx(i ,j,k1)))* &
2059!^ & tl_dRdx(i ,j,k1)
2060!^
2061 ad_drdx(i ,j,k1)=ad_drdx(i ,j,k1)+ &
2062 & (0.5_r8+sign(0.5_r8, &
2063 & drdx(i ,j,k1)))* &
2064 & ad_cff1
2065 ad_cff1=0.0_r8
2066 END DO
2067 END DO
2068 END IF
2069 DO j=jmin,jmax+1
2070 DO i=imin,imax
2071#ifdef DIFF_3DCOEF
2072# ifdef TS_U3ADV_SPLIT
2073 cff=0.5_r8*diff3d_v(i,j,k)*om_v(i,j)
2074# else
2075 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i,j-1,k))* &
2076 & om_v(i,j)
2077# endif
2078#else
2079 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i,j-1,itrc))* &
2080 & om_v(i,j)
2081#endif
2082!^ tl_FE(i,j)=cff* &
2083!^ & ((tl_Hz(i,j,k)+tl_Hz(i,j-1,k))* &
2084!^ & (dTde(i,j,k1)- &
2085!^ & 0.5_r8*(MAX(dRde(i,j,k1),0.0_r8)* &
2086!^ & (dTdr(i,j-1,k1)+ &
2087!^ & dTdr(i,j ,k2))+ &
2088!^ & MIN(dRde(i,j,k1),0.0_r8)* &
2089!^ & (dTdr(i,j-1,k2)+ &
2090!^ & dTdr(i,j ,k1))))+ &
2091!^ & (Hz(i,j,k)+Hz(i,j-1,k))* &
2092!^ & (tl_dTde(i,j,k1)- &
2093!^ & 0.5_r8*(MAX(dRde(i,j,k1),0.0_r8)* &
2094!^ & (tl_dTdr(i,j-1,k1)+ &
2095!^ & tl_dTdr(i,j ,k2))+ &
2096!^ & MIN(dRde(i,j,k1),0.0_r8)* &
2097!^ & (tl_dTdr(i,j-1,k2)+ &
2098!^ & tl_dTdr(i,j ,k1)))- &
2099!^ & 0.5_r8*((0.5_r8+ &
2100!^ & SIGN(0.5_r8, dRde(i,j,k1)))* &
2101!^ & tl_dRde(i,j,k1)* &
2102!^ & (dTdr(i,j-1,k1)+dTdr(i,j,k2))+ &
2103!^ & (0.5_r8+ &
2104!^ & SIGN(0.5_r8,-dRde(i,j,k1)))* &
2105!^ & tl_dRde(i,j,k1)* &
2106!^ & (dTdr(i,j-1,k2)+dTdr(i,j,k1)))))
2107!^
2108 adfac=cff*ad_fe(i,j)
2109 adfac1=adfac*(dtde(i,j,k1)- &
2110 & 0.5_r8*(max(drde(i,j,k1),0.0_r8)* &
2111 & (dtdr(i,j-1,k1)+ &
2112 & dtdr(i ,j,k2))+ &
2113 & min(drde(i,j,k1),0.0_r8)* &
2114 & (dtdr(i,j-1,k2)+ &
2115 & dtdr(i ,j,k1))))
2116 adfac2=adfac*(hz(i,j,k)+hz(i,j-1,k))
2117 adfac3=adfac2*0.5_r8*max(drde(i,j,k1),0.0_r8)
2118 adfac4=adfac2*0.5_r8*min(drde(i,j,k1),0.0_r8)
2119 ad_hz(i,j-1,k)=ad_hz(i,j-1,k)+adfac1
2120 ad_hz(i,j ,k)=ad_hz(i,j ,k)+adfac1
2121 ad_dtde(i,j,k1)=ad_dtde(i,j,k1)+adfac2
2122 ad_dtdr(i,j-1,k1)=ad_dtdr(i,j-1,k1)-adfac3
2123 ad_dtdr(i,j, k2)=ad_dtdr(i,j ,k2)-adfac3
2124 ad_dtdr(i,j-1,k2)=ad_dtdr(i,j-1,k2)-adfac4
2125 ad_dtdr(i,j ,k1)=ad_dtdr(i,j ,k1)-adfac4
2126 ad_drde(i,j,k1)=ad_drde(i,j,k1)- &
2127 & adfac2*0.5_r8* &
2128 & ((0.5_r8+sign(0.5_r8, drde(i,j,k1)))* &
2129 & (dtdr(i,j-1,k1)+dtdr(i,j,k2))+ &
2130 & (0.5_r8+sign(0.5_r8,-drde(i,j,k1)))* &
2131 & (dtdr(i,j-1,k2)+dtdr(i,j,k1)))
2132 ad_fe(i,j)=0.0_r8
2133 END DO
2134 END DO
2135 DO j=jmin,jmax
2136 DO i=imin,imax+1
2137#ifdef DIFF_3DCOEF
2138# ifdef TS_U3ADV_SPLIT
2139 cff=0.5_r8*diff3d_u(i,j,k)*on_u(i,j)
2140# else
2141 cff=0.25_r8*(diff3d_r(i,j,k)+diff3d_r(i-1,j,k))* &
2142 & on_u(i,j)
2143# endif
2144#else
2145 cff=0.25_r8*(diff4(i,j,itrc)+diff4(i-1,j,itrc))* &
2146 & on_u(i,j)
2147#endif
2148!^ tl_FX(i,j)=cff* &
2149!^ & ((tl_Hz(i,j,k)+tl_Hz(i-1,j,k))* &
2150!^ & (dTdx(i,j,k1)- &
2151!^ & 0.5_r8*(MAX(dRdx(i,j,k1),0.0_r8)* &
2152!^ & (dTdr(i-1,j,k1)+ &
2153!^ & dTdr(i ,j,k2))+ &
2154!^ & MIN(dRdx(i,j,k1),0.0_r8)* &
2155!^ & (dTdr(i-1,j,k2)+ &
2156!^ & dTdr(i ,j,k1))))+ &
2157!^ & (Hz(i,j,k)+Hz(i-1,j,k))* &
2158!^ & (tl_dTdx(i,j,k1)- &
2159!^ & 0.5_r8*(MAX(dRdx(i,j,k1),0.0_r8)* &
2160!^ & (tl_dTdr(i-1,j,k1)+ &
2161!^ & tl_dTdr(i ,j,k2))+ &
2162!^ & MIN(dRdx(i,j,k1),0.0_r8)* &
2163!^ & (tl_dTdr(i-1,j,k2)+ &
2164!^ & tl_dTdr(i ,j,k1)))- &
2165!^ & 0.5_r8*((0.5_r8+ &
2166!^ & SIGN(0.5_r8, dRdx(i,j,k1)))* &
2167!^ & tl_dRdx(i,j,k1)* &
2168!^ & (dTdr(i-1,j,k1)+dTdr(i,j,k2))+ &
2169!^ & (0.5_r8+ &
2170!^ & SIGN(0.5_r8,-dRdx(i,j,k1)))* &
2171!^ & tl_dRdx(i,j,k1)* &
2172!^ & (dTdr(i-1,j,k2)+dTdr(i,j,k1)))))
2173!^
2174 adfac=cff*ad_fx(i,j)
2175 adfac1=adfac*(dtdx(i,j,k1)- &
2176 & 0.5_r8*(max(drdx(i,j,k1),0.0_r8)* &
2177 & (dtdr(i-1,j,k1)+ &
2178 & dtdr(i ,j,k2))+ &
2179 & min(drdx(i,j,k1),0.0_r8)* &
2180 & (dtdr(i-1,j,k2)+ &
2181 & dtdr(i ,j,k1))))
2182 adfac2=adfac*(hz(i,j,k)+hz(i-1,j,k))
2183 adfac3=adfac2*0.5_r8*max(drdx(i,j,k1),0.0_r8)
2184 adfac4=adfac2*0.5_r8*min(drdx(i,j,k1),0.0_r8)
2185 ad_hz(i-1,j,k)=ad_hz(i-1,j,k)+adfac1
2186 ad_hz(i ,j,k)=ad_hz(i ,j,k)+adfac1
2187 ad_dtdx(i,j,k1)=ad_dtdx(i,j,k1)+adfac2
2188 ad_dtdr(i-1,j,k1)=ad_dtdr(i-1,j,k1)-adfac3
2189 ad_dtdr(i ,j,k2)=ad_dtdr(i ,j,k2)-adfac3
2190 ad_dtdr(i-1,j,k2)=ad_dtdr(i-1,j,k2)-adfac4
2191 ad_dtdr(i ,j,k1)=ad_dtdr(i ,j,k1)-adfac4
2192 ad_drdx(i,j,k1)=ad_drdx(i,j,k1)- &
2193 & adfac2*0.5_r8* &
2194 & ((0.5_r8+sign(0.5_r8, drdx(i,j,k1)))* &
2195 & (dtdr(i-1,j,k1)+dtdr(i,j,k2))+ &
2196 & (0.5_r8+sign(0.5_r8,-drdx(i,j,k1)))* &
2197 & (dtdr(i-1,j,k2)+dtdr(i,j,k1)))
2198 ad_fx(i,j)=0.0_r8
2199 END DO
2200 END DO
2201 END IF
2202 IF ((k.eq.0).or.(k.eq.n(ng))) THEN
2203 DO j=-1+jmin,jmax+1
2204 DO i=-1+imin,imax+1
2205!^ tl_dTdr(i,j,k2)=0.0_r8
2206!^
2207 ad_dtdr(i,j,k2)=0.0_r8
2208!^ tl_FS(i,j,k2)=0.0_r8
2209!^
2210 ad_fs(i,j,k2)=0.0_r8
2211 END DO
2212 END DO
2213 ELSE
2214 DO j=-1+jmin,jmax+1
2215 DO i=-1+imin,imax+1
2216#if defined TS_MIX_MAX_SLOPE
2217 cff1=sqrt(drdx(i,j,k2)**2+drdx(i+1,j,k2)**2+ &
2218 & drdx(i,j,k1)**2+drdx(i+1,j,k1)**2+ &
2219 & drde(i,j,k2)**2+drde(i,j+1,k2)**2+ &
2220 & drde(i,j,k1)**2+drde(i,j+1,k1)**2)
2221 cff2=0.25_r8*slope_max* &
2222 & (z_r(i,j,k+1)-z_r(i,j,k))*cff1
2223 cff3=max(pden(i,j,k)-pden(i,j,k+1),small)
2224 cff4=max(cff2,cff3)
2225 cff=-1.0_r8/cff4
2226#elif defined TS_MIX_MIN_STRAT
2227 cff1=max(pden(i,j,k)-pden(i,j,k+1), &
2228 & strat_min*(z_r(i,j,k+1)-z_r(i,j,k)))
2229 cff=-1.0_r8/cff1
2230#else
2231 cff1=max(pden(i,j,k)-pden(i,j,k+1),eps)
2232 cff=-1.0_r8/cff1
2233#endif
2234!^ tl_FS(i,j,k2)=tl_cff*(z_r(i,j,k+1)- &
2235!^ & z_r(i,j,k ))+ &
2236!^ & cff*(tl_z_r(i,j,k+1)- &
2237!^ & tl_z_r(i,j,k ))
2238!^
2239 adfac=cff*ad_fs(i,j,k2)
2240 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac
2241 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac
2242 ad_cff=ad_cff+(z_r(i,j,k+1)- &
2243 & z_r(i,j,k ))*ad_fs(i,j,k2)
2244 ad_fs(i,j,k2)=0.0_r8
2245#if defined TS_MIX_STABILITY
2246!^ tl_dTdr(i,j,k2)=tl_cff* &
2247!^ & (0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
2248!^ & t(i,j,k ,nrhs,itrc))+ &
2249!^ & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
2250!^ & t(i,j,k ,nstp,itrc)))+ &
2251!^ & cff*
2252!^ & (0.75_r8*(tl_t(i,j,k+1,nrhs,itrc)- &
2253!^ & tl_t(i,j,k ,nrhs,itrc))+ &
2254!^ & 0.25_r8*(tl_t(i,j,k+1,nstp,itrc)- &
2255!^ & tl_t(i,j,k ,nstp,itrc)))
2256!^
2257 adfac=cff*ad_dtdr(i,j,k2)
2258 adfac1=adfac*0.75_r8
2259 adfac2=adfac*0.25_r8
2260 ad_t(i,j,k ,nrhs,itrc)=ad_t(i,j,k ,nrhs,itrc)-adfac1
2261 ad_t(i,j,k+1,nrhs,itrc)=ad_t(i,j,k+1,nrhs,itrc)+adfac1
2262 ad_t(i,j,k ,nstp,itrc)=ad_t(i,j,k ,nstp,itrc)-adfac2
2263 ad_t(i,j,k+1,nstp,itrc)=ad_t(i,j,k+1,nstp,itrc)+adfac2
2264 ad_cff=ad_cff+(0.75_r8*(t(i,j,k+1,nrhs,itrc)- &
2265 & t(i,j,k ,nrhs,itrc))+ &
2266 & 0.25_r8*(t(i,j,k+1,nstp,itrc)- &
2267 & t(i,j,k ,nstp,itrc)))* &
2268 & ad_dtdz(i,j,k2)
2269 ad_dtdz(i,j,k2)=0.0_r8
2270#elif defined TS_MIX_CLIMA
2271 IF (ltracerclm(itrc,ng)) THEN
2272!^ tl_dTdr(i,j,k2)=tl_cff*((t(i,j,k+1,nrhs,itrc)- &
2273!^ & tclm(i,j,k+1,itrc))- &
2274!^ & (t(i,j,k ,nrhs,itrc)- &
2275!^ & tclm(i,j,k ,itrc)))+ &
2276!^ & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
2277!^ & tl_t(i,j,k ,nrhs,itrc))
2278!^
2279 adfac=cff*ad_dtdr(i,j,k2)
2280 ad_t(i,j,k ,nrhs,itrc)=ad_t(i,j,k ,nrhs,itrc)-adfac
2281 ad_t(i,j,k+1,nrhs,itrc)=ad_t(i,j,k+1,nrhs,itrc)+adfac
2282 ad_cff=ad_cff+((t(i,j,k+1,nrhs,itrc)- &
2283 & tclm(i,j,k+1,itrc))- &
2284 & (t(i,j,k ,nrhs,itrc)- &
2285 & tclm(i,j,k ,itrc)))*ad_dtdr(i,j,k2)
2286 ad_dtdr(i,j,k2)=0.0_r8
2287 ELSE
2288!^ tl_dTdr(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
2289!^ & t(i,j,k ,nrhs,itrc))+ &
2290!^ & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
2291!^ & tl_t(i,j,k ,nrhs,itrc))
2292!^
2293 END IF
2294#else
2295!^ tl_dTdr(i,j,k2)=tl_cff*(t(i,j,k+1,nrhs,itrc)- &
2296!^ & t(i,j,k ,nrhs,itrc))+ &
2297!^ & cff*(tl_t(i,j,k+1,nrhs,itrc)- &
2298!^ & tl_t(i,j,k ,nrhs,itrc))
2299!^
2300 adfac=cff*ad_dtdr(i,j,k2)
2301 ad_t(i,j,k ,nrhs,itrc)=ad_t(i,j,k ,nrhs,itrc)-adfac
2302 ad_t(i,j,k+1,nrhs,itrc)=ad_t(i,j,k+1,nrhs,itrc)+adfac
2303 ad_cff=ad_cff+(t(i,j,k+1,nrhs,itrc)- &
2304 & t(i,j,k ,nrhs,itrc))*ad_dtdr(i,j,k2)
2305 ad_dtdr(i,j,k2)=0.0_r8
2306#endif
2307#if defined TS_MIX_MAX_SLOPE
2308!^ tl_cff=cff*cff*tl_cff4
2309!^
2310 ad_cff4=ad_cff4+cff*cff*ad_cff
2311 ad_cff=0.0_r8
2312!^ tl_cff4=(0.5_r8+SIGN(0.5_r8,cff2-cff3))*tl_cff2+ &
2313!^ & (0.5_r8-SIGN(0.5_r8,cff2-cff3))*tl_cff3
2314!^
2315 ad_cff3=ad_cff3+ &
2316 & (0.5_r8-sign(0.5_r8,cff2-cff3))*ad_cff4
2317 ad_cff2=ad_cff2+ &
2318 & (0.5_r8+sign(0.5_r8,cff2-cff3))*ad_cff4
2319 ad_cff4=0.0_r8
2320!^ tl_cff3=(0.5_r8+SIGN(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
2321!^ & small))* &
2322!^ & (tl_pden(i,j,k)-tl_pden(i,j,k+1))
2323!^
2324 adfac=(0.5_r8+sign(0.5_r8,pden(i,j,k)-pden(i,j,k+1)- &
2325 & small))*ad_cff3
2326 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac
2327 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac
2328 ad_cff3=0.0_r8
2329!^ tl_cff2=0.25_r8*slope_max* &
2330!^ & ((tl_z_r(i,j,k+1)-tl_z_r(i,j,k))*cff1+ &
2331!^ & (z_r(i,j,k+1)-z_r(i,j,k))*tl_cff1)
2332!^
2333 adfac=0.25_r8*slope_max*ad_cff2
2334 adfac1=adfac*cff1
2335 ad_cff1=ad_cff1+(z_r(i,j,k+1)-z_r(i,j,k))*adfac
2336 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac1
2337 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac1
2338 ad_cff2=0.0_r8
2339 IF (cff1.ne.0.0_r8) THEN
2340!^ tl_cff1=(dRdx(i ,j,k2)*tl_dRdx(i ,j,k2)+ &
2341!^ & dRdx(i+1,j,k2)*tl_dRdx(i+1,j,k2)+ &
2342!^ & dRdx(i ,j,k1)*tl_dRdx(i ,j,k1)+ &
2343!^ & dRdx(i+1,j,k1)*tl_dRdx(i+1,j,k1)+ &
2344!^ & dRde(i,j ,k2)*tl_dRde(i,j ,k2)+ &
2345!^ & dRde(i,j+1,k2)*tl_dRde(i,j+1,k2)+ &
2346!^ & dRde(i,j ,k1)*tl_dRde(i,j ,k1)+ &
2347!^ & dRde(i,j+1,k1)*tl_dRde(i,j+1,k1))/cff1
2348!^
2349 adfac=ad_cff1/cff1
2350 ad_drdx(i ,j,k1)=ad_drdx(i ,j,k1)+ &
2351 & drdx(i ,j,k1)*adfac
2352 ad_drdx(i+1,j,k1)=ad_drdx(i+1,j,k1)+ &
2353 & drdx(i+1,j,k1)*adfac
2354 ad_drdx(i ,j,k2)=ad_drdx(i ,j,k2)+ &
2355 & drdx(i ,j,k2)*adfac
2356 ad_drdx(i+1,j,k2)=ad_drdx(i+1,j,k2)+ &
2357 & drdx(i+1,j,k2)*adfac
2358 ad_drde(i,j ,k2)=ad_drde(i,j ,k2)+ &
2359 & drde(i,j ,k2)*adfac
2360 ad_drde(i,j+1,k2)=ad_drde(i,j+1,k2)+ &
2361 & drde(i,j+1,k2)*adfac
2362 ad_drde(i,j ,k1)=ad_drde(i,j ,k1)+ &
2363 & drde(i,j ,k1)*adfac
2364 ad_drde(i,j+1,k1)=ad_drde(i,j+1,k1)+ &
2365 & drde(i,j+1,k1)*adfac
2366 ad_cff1=0.0_r8
2367 ELSE
2368!^ tl_cff1=0.0_r8
2369!^
2370 ad_cff1=0.0_r8
2371 END IF
2372#elif defined TS_MIX_MIN_STRAT
2373!^ tl_cff=cff*cff*tl_cff1
2374!^
2375 ad_cff1=ad_cff1+cff*cff*ad_cff
2376 ad_cff=0.0_r8
2377!^ tl_cff1=(0.5_r8+SIGN(0.5_r8, &
2378!^ & pden(i,j,k)-pden(i,j,k+1)- &
2379!^ & strat_min*(z_r(i,j,k+1)- &
2380!^ & z_r(i,j,k ))))* &
2381!^ & (tl_pden(i,j,k)-tl_pden(i,j,k+1))+ &
2382!^ & (0.5_r8-SIGN(0.5_r8, &
2383!^ & pden(i,j,k)-pden(i,j,k+1)- &
2384!^ & strat_min*(z_r(i,j,k+1)- &
2385!^ & z_r(i,j,k ))))* &
2386!^ & (strat_min*(tl_z_r(i,j,k+1)-tl_z_r(i,j,k )))
2387!^
2388 adfac1=(0.5_r8+sign(0.5_r8, &
2389 & pden(i,j,k)-pden(i,j,k+1)- &
2390 & strat_min*(z_r(i,j,k+1)- &
2391 & z_r(i,j,k ))))* &
2392 & ad_cff1
2393 adfac2=(0.5_r8-sign(0.5_r8, &
2394 & pden(i,j,k)-pden(i,j,k+1)- &
2395 & strat_min*(z_r(i,j,k+1)- &
2396 & z_r(i,j,k ))))* &
2397 & strat_min*ad_cff1
2398 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac1
2399 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac1
2400 ad_z_r(i,j,k )=ad_z_r(i,j,k )-adfac2
2401 ad_z_r(i,j,k+1)=ad_z_r(i,j,k+1)+adfac2
2402 ad_cff1=0.0_r8
2403#else
2404!^ tl_cff=cff*cff*tl_cff1
2405!^
2406 ad_cff1=ad_cff1+cff*cff*ad_cff
2407 ad_cff=0.0_r8
2408!^ tl_cff1=(0.5_r8+SIGN(0.5_r8, &
2409!^ & pden(i,j,k)-pden(i,j,k+1)-eps))* &
2410!^ & (tl_pden(i,j,k)-tl_pden(i,j,k+1))
2411!^
2412 adfac=(0.5_r8+sign(0.5_r8, &
2413 & pden(i,j,k)-pden(i,j,k+1)-eps))* &
2414 & ad_cff1
2415 ad_pden(i,j,k )=ad_pden(i,j,k )+adfac
2416 ad_pden(i,j,k+1)=ad_pden(i,j,k+1)-adfac
2417 ad_cff1=0.0_r8
2418#endif
2419 END DO
2420 END DO
2421 END IF
2422 IF (k.lt.n(ng)) THEN
2423 DO j=jmin,jmax+1
2424 DO i=imin,imax
2425 cff=0.5_r8*(pn(i,j)+pn(i,j-1))
2426#ifdef MASKING
2427 cff=cff*vmask(i,j)
2428#endif
2429#ifdef WET_DRY_NOT_YET
2430 cff=cff*vmask_wet(i,j)
2431#endif
2432#if defined TS_MIX_STABILITY
2433!^ tl_dTde(i,j,k2)=cff* &
2434!^ & (0.75_r8*(tl_t(i,j ,k+1,nrhs,itrc)- &
2435!^ & tl_t(i,j-1,k+1,nrhs,itrc))+ &
2436!^ & 0.25_r8*(tl_t(i,j ,k+1,nstp,itrc)- &
2437!^ & tl_t(i,j-1,k+1,nstp,itrc)))
2438!^
2439 adfac=cff*ad_dtde(i,j,k2)
2440 adfac1=adfac*0.75_r8
2441 adfac2=adfac*0.25_r8
2442 ad_t(i,j-1,k+1,nrhs,itrc)=ad_t(i,j-1,k+1,nrhs,itrc)- &
2443 & adfac1
2444 ad_t(i,j ,k+1,nrhs,itrc)=ad_t(i,j ,k+1,nrhs,itrc)+ &
2445 & adfac1
2446 ad_t(i,j-1,k+1,nstp,itrc)=ad_t(i,j-1,k+1,nstp,itrc)- &
2447 & adfac2
2448 ad_t(i,j ,k+1,nstp,itrc)=ad_t(i,j ,k+1,nstp,itrc)+ &
2449 & adfac2
2450 ad_dtde(i,j,k2)=0.0_r8
2451#elif defined TS_MIX_CLIMA
2452!^ tl_dTde(i,j,k2)=cff*(tl_t(i,j ,k+1,nrhs,itrc)- &
2453!^ & tl_t(i,j-1,k+1,nrhs,itrc))
2454!^
2455 adfac=cff*ad_dtde(i,j,k2)
2456 ad_t(i,j-1,k+1,nrhs,itrc)=ad_t(i,j-1,k+1,nrhs,itrc)- &
2457 & adfac
2458 ad_t(i,j ,k+1,nrhs,itrc)=ad_t(i,j ,k+1,nrhs,itrc)+ &
2459 & adfac
2460 ad_dtde(i,j,k2)=0.0_r8
2461#else
2462!^ tl_dTde(i,j,k2)=cff*(tl_t(i,j ,k+1,nrhs,itrc)- &
2463!^ & tl_t(i,j-1,k+1,nrhs,itrc))
2464!^
2465 adfac=cff*ad_dtde(i,j,k2)
2466 ad_t(i,j-1,k+1,nrhs,itrc)=ad_t(i,j-1,k+1,nrhs,itrc)- &
2467 & adfac
2468 ad_t(i,j ,k+1,nrhs,itrc)=ad_t(i,j ,k+1,nrhs,itrc)+ &
2469 & adfac
2470 ad_dtde(i,j,k2)=0.0_r8
2471#endif
2472!^ tl_dRde(i,j,k2)=cff*(tl_pden(i,j ,k+1)- &
2473!^ & tl_pden(i,j-1,k+1))
2474!^
2475 adfac=cff*ad_drde(i,j,k2)
2476 ad_pden(i,j-1,k+1)=ad_pden(i,j-1,k+1)-adfac
2477 ad_pden(i,j ,k+1)=ad_pden(i,j ,k+1)+adfac
2478 ad_drde(i,j,k2)=0.0_r8
2479 END DO
2480 END DO
2481 DO j=jmin,jmax
2482 DO i=imin,imax+1
2483 cff=0.5_r8*(pm(i,j)+pm(i-1,j))
2484#ifdef MASKING
2485 cff=cff*umask(i,j)
2486#endif
2487#ifdef WET_DRY_NOT_YET
2488 cff=cff*umask_wet(i,j)
2489#endif
2490#if defined TS_MIX_STABILITY
2491!^ tl_dTdx(i,j,k2)=cff* &
2492!^ & (0.75_r8*(tl_t(i ,j,k+1,nrhs,itrc)- &
2493!^ & tl_t(i-1,j,k+1,nrhs,itrc))+ &
2494!^ & 0.25_r8*(tl_t(i ,j,k+1,nstp,itrc)- &
2495!^ & tl_t(i-1,j,k+1,nstp,itrc)))
2496!^
2497 adfac=cff*ad_dtdx(i,j,k2)
2498 adfac1=adfac*0.75_r8
2499 adfac2=adfac*0.25_r8
2500 ad_t(i-1,j,k+1,nrhs,itrc)=ad_t(i-1,j,k+1,nrhs,itrc)- &
2501 & adfac1
2502 ad_t(i ,j,k+1,nrhs,itrc)=ad_t(i ,j,k+1,nrhs,itrc)+ &
2503 & adfac1
2504 ad_t(i-1,j,k+1,nstp,itrc)=ad_t(i-1,j,k+1,nstp,itrc)- &
2505 & adfac2
2506 ad_t(i ,j,k+1,nstp,itrc)=ad_t(i ,j,k+1,nstp,itrc)+ &
2507 & adfac2
2508 ad_dtdx(i,j,k2)=0.0_r8
2509#elif defined TS_MIX_CLIMA
2510!^ tl_dTdx(i,j,k2)=cff*(tl_t(i ,j,k+1,nrhs,itrc)- &
2511!^ & tl_t(i-1,j,k+1,nrhs,itrc))
2512!^
2513 adfac=cff*ad_dtdx(i,j,k2)
2514 ad_t(i-1,j,k+1,nrhs,itrc)=ad_t(i-1,j,k+1,nrhs,itrc)- &
2515 & adfac
2516 ad_t(i ,j,k+1,nrhs,itrc)=ad_t(i ,j,k+1,nrhs,itrc)+ &
2517 & adfac
2518 ad_dtdx(i,j,k2)=0.0_r8
2519#else
2520!^ tl_dTdx(i,j,k2)=cff*(tl_t(i ,j,k+1,nrhs,itrc)- &
2521!^ & tl_t(i-1,j,k+1,nrhs,itrc))
2522!^
2523 adfac=cff*ad_dtdx(i,j,k2)
2524 ad_t(i-1,j,k+1,nrhs,itrc)=ad_t(i-1,j,k+1,nrhs,itrc)- &
2525 & adfac
2526 ad_t(i ,j,k+1,nrhs,itrc)=ad_t(i ,j,k+1,nrhs,itrc)+ &
2527 & adfac
2528 ad_dtdx(i,j,k2)=0.0_r8
2529#endif
2530!^ tl_dRdx(i,j,k2)=cff*(tl_pden(i ,j,k+1)- &
2531!^ & tl_pden(i-1,j,k+1))
2532!^
2533 adfac=cff*ad_drdx(i,j,k2)
2534 ad_pden(i-1,j,k+1)=ad_pden(i-1,j,k+1)-adfac
2535 ad_pden(i ,j,k+1)=ad_pden(i ,j,k+1)+adfac
2536 ad_drdx(i,j,k2)=0.0_r8
2537 END DO
2538 END DO
2539 END IF
2540!
2541! Compute new storage recursive indices.
2542!
2543 kt=k2
2544 k2=k1
2545 k1=kt
2546 END DO k_loop3
2547 END DO t_loop
2548!
2549 RETURN
2550 END SUBROUTINE ad_t3dmix4_iso_tile
2551
2552 END MODULE ad_t3dmix4_mod
ad_t3dmix4_mod
Definition ad_t3dmix4_geo.h:1
ad_t3dmix4_mod::ad_t3dmix4
subroutine, public ad_t3dmix4(ng, tile)
Definition ad_t3dmix4_geo.h:26
ad_t3dmix4_mod::ad_t3dmix4_iso_tile
subroutine ad_t3dmix4_iso_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, ad_hz, z_r, ad_z_r, diff3d_u, diff3d_v, diff3d_r, diff4, pden, ad_pden, tclm, t, ad_t)
Definition ad_t3dmix4_iso.h:132
mod_clima
Definition mod_clima.F:2
mod_clima::clima
type(t_clima), dimension(:), allocatable clima
Definition mod_clima.F:153
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_ncparam::istvar
integer, dimension(:), allocatable istvar
Definition mod_ncparam.F:520
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::ad_lbc
type(t_lbc), dimension(:,:,:), allocatable ad_lbc
Definition mod_param.F:378
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::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
mod_scalars
Definition mod_scalars.F:2
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
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::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::ltracerclm
logical, dimension(:,:), allocatable ltracerclm
Definition mod_scalars.F:528
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
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
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