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