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