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