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