ROMS
Loading...
Searching...
No Matches
ad_balance.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE ad_balance_mod
3
4#if defined TANGENT && defined BALANCE_OPERATOR
5!
6!git $Id$
7!================================================== Hernan G. Arango ===
8! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
9! Licensed under a MIT/X style license !
10! See License_ROMS.md !
11!=======================================================================
12! !
13! These routines impose a multivariate balance operator to constraint !
14! the background/model error covariance matrix, B, such that the !
15! unobserved variables information is extracted from observed data. !
16! It follows the approach proposed by Weaver et al. (2006). The state !
17! vector is split between balanced and unbalanced components, except !
18! for temperature, which is used to establish the balanced part of !
19! other state variables. !
20! !
21! The background/model error covariance is represented as: !
22! !
23! B = K Bu K^(T) !
24! !
25! where !
26! !
27! B : background/model error covariance matrix. !
28! Bu: unbalanced background/model error covariance matrix modeled !
29! with the generalized diffusion operator. !
30! K : balance matrix operator. !
31! !
32! Here, T denotes the transpose. !
33! !
34! The multivariate formulation is obtained by establishing balance !
35! relationships with the other state variables using T-S empirical !
36! formulas, hydrostatic balance, and geostrophic balance. !
37! !
38! Reference: !
39! !
40! Weaver, A.T., C. Deltel, E. Machu, S. Ricci, and N. Daget, 2005: !
41! A multivariate balance operator for variational data assimila- !
42! tion, Q. J. R. Meteorol. Soc., 131, 3605-3625. !
43! !
44! (See also, ECMWR Technical Memorandum # 491, April 2006) !
45! !
46!=======================================================================
47!
48 USE mod_kinds
49
50 implicit none
51
52 PRIVATE
53 PUBLIC :: ad_balance
54
55 CONTAINS
56!
57!***********************************************************************
58 SUBROUTINE ad_balance (ng, tile, Lbck, Linp)
59!***********************************************************************
60!
61 USE mod_param
62 USE mod_grid
63# ifdef SOLVE3D
64 USE mod_coupling
65 USE mod_mixing
66# endif
67# ifdef ZETA_ELLIPTIC
68 USE mod_fourdvar
69# endif
70 USE mod_ocean
71 USE mod_stepping
72# ifdef SOLVE3D
73!
74 USE rho_eos_mod
75 USE set_depth_mod
76# endif
77!
78! Imported variable declarations.
79!
80 integer, intent(in) :: ng, tile, lbck, linp
81!
82! Local variable declarations.
83!
84# include "tile.h"
85!
86# ifdef SOLVE3D
87!
88! Compute background state thickness, depth arrays, thermal expansion,
89! and saline contraction coefficients.
90!
91 coupling(ng) % Zt_avg1 =0.0_r8
92
93 CALL set_depth (ng, tile, iadm)
94 nrhs(ng)=lbck
95 CALL rho_eos (ng, tile, iadm)
96!
97# endif
98 CALL ad_balance_tile (ng, tile, &
99 & lbi, ubi, lbj, ubj, lbij, ubij, &
100 & imins, imaxs, jmins, jmaxs, &
101 & lbck, linp, &
102 & grid(ng) % f, &
103 & grid(ng) % pm, &
104 & grid(ng) % pn, &
105# ifdef ZETA_ELLIPTIC
106 & grid(ng) % pmon_u, &
107 & grid(ng) % pnom_v, &
108 & grid(ng) % h, &
109# endif
110# ifdef SOLVE3D
111 & grid(ng) % Hz, &
112 & grid(ng) % z_r, &
113 & grid(ng) % z_w, &
114# endif
115# ifdef MASKING
116 & grid(ng) % rmask, &
117 & grid(ng) % umask, &
118 & grid(ng) % vmask, &
119# endif
120# ifdef SOLVE3D
121 & mixing(ng) % alpha, &
122 & mixing(ng) % beta, &
123 & ocean(ng) % t, &
124# endif
125# ifdef ZETA_ELLIPTIC
126 & fourdvar(ng) % ad_zeta_ref, &
127 & fourdvar(ng) % ad_rhs_r2d, &
128 & fourdvar(ng) % pc_r2d, &
129 & fourdvar(ng) % r_r2d, &
130 & fourdvar(ng) % br_r2d, &
131 & fourdvar(ng) % p_r2d, &
132 & fourdvar(ng) % bp_r2d, &
133 & fourdvar(ng) % zdf1, &
134 & fourdvar(ng) % zdf2, &
135 & fourdvar(ng) % zdf3, &
136 & fourdvar(ng) % bc_ak, &
137 & fourdvar(ng) % bc_bk, &
138# endif
139# ifdef SOLVE3D
140 & ocean(ng) % ad_rho, &
141 & ocean(ng) % ad_t, &
142 & ocean(ng) % ad_u, &
143 & ocean(ng) % ad_v, &
144# endif
145 & ocean(ng) % ad_zeta)
146
147 RETURN
148 END SUBROUTINE ad_balance
149!
150!***********************************************************************
151 SUBROUTINE ad_balance_tile (ng, tile, &
152 & LBi, UBi, LBj, UBj, LBij, UBij, &
153 & IminS, ImaxS, JminS, JmaxS, &
154 & Lbck, Linp, &
155 & f, pm, pn, &
156# ifdef ZETA_ELLIPTIC
157 & pmon_u, pnom_v, h, &
158# endif
159# ifdef SOLVE3D
160 & Hz, z_r, z_w, &
161# endif
162# ifdef MASKING
163 & rmask, umask, vmask, &
164# endif
165# ifdef SOLVE3D
166 & alpha, beta, t, &
167# endif
168# ifdef ZETA_ELLIPTIC
169 & ad_zeta_ref, ad_rhs_r2d, &
170 & pc_r2d, r_r2d, br_r2d, &
171 & p_r2d, bp_r2d, &
172 & zdf1, zdf2, zdf3, bc_ak, bc_bk, &
173# endif
174# ifdef SOLVE3D
175 & ad_rho, ad_t, ad_u, ad_v, &
176# endif
177 & ad_zeta)
178!***********************************************************************
179!
180 USE mod_param
181 USE mod_ncparam
182 USE mod_scalars
183!
184 USE ad_bc_2d_mod
185 USE ad_exchange_2d_mod
186 USE ad_exchange_3d_mod
187# ifdef DISTRIBUTE
188 USE mp_exchange_mod, ONLY : ad_mp_exchange2d, ad_mp_exchange3d
189# endif
190# ifdef ZETA_ELLIPTIC
191 USE zeta_balance_mod, ONLY: ad_r2d_bc, u2d_bc, v2d_bc
192 USE zeta_balance_mod, ONLY: ad_biconj_tile
193# endif
194!
195! Imported variable declarations.
196!
197 integer, intent(in) :: ng, tile
198 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
199 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
200 integer, intent(in) :: Lbck, Linp
201!
202# ifdef ASSUMED_SHAPE
203 real(r8), intent(in) :: f(LBi:,LBj:)
204 real(r8), intent(in) :: pm(LBi:,LBj:)
205 real(r8), intent(in) :: pn(LBi:,LBj:)
206# ifdef ZETA_ELLIPTIC
207 real(r8), intent(in) :: h(LBi:,LBj:)
208 real(r8), intent(in) :: pmon_u(LBi:,LBj:)
209 real(r8), intent(in) :: pnom_v(LBi:,LBj:)
210# endif
211# ifdef SOLVE3D
212 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
213 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
214 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
215# endif
216# ifdef MASKING
217 real(r8), intent(in) :: rmask(LBi:,LBj:)
218 real(r8), intent(in) :: umask(LBi:,LBj:)
219 real(r8), intent(in) :: vmask(LBi:,LBj:)
220# endif
221# ifdef SOLVE3D
222 real(r8), intent(in) :: alpha(LBi:,LBj:)
223 real(r8), intent(in) :: beta(LBi:,LBj:)
224 real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
225 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
226 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
227 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
228# endif
229 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
230# ifdef SOLVE3D
231 real(r8), intent(out) :: ad_rho(LBi:,LBj:,:)
232# endif
233# ifdef ZETA_ELLIPTIC
234 real(r8), intent(in) :: bc_ak(:)
235 real(r8), intent(in) :: bc_bk(:)
236 real(r8), intent(in) :: zdf1(:)
237 real(r8), intent(in) :: zdf2(:)
238 real(r8), intent(in) :: zdf3(:)
239 real(r8), intent(inout) :: pc_r2d(LBi:,LBj:)
240 real(r8), intent(inout) :: r_r2d(LBi:,LBj:,:)
241 real(r8), intent(inout) :: br_r2d(LBi:,LBj:,:)
242 real(r8), intent(inout) :: p_r2d(LBi:,LBj:,:)
243 real(r8), intent(inout) :: bp_r2d(LBi:,LBj:,:)
244 real(r8), intent(inout) :: ad_rhs_r2d(LBi:,LBj:)
245 real(r8), intent(inout) :: ad_zeta_ref(LBi:,LBj:)
246# endif
247
248# else
249
250 real(r8), intent(in) :: f(LBi:UBi,LBj:UBj)
251 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
252 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
253# ifdef ZETA_ELLIPTIC
254 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
255 real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
256 real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
257# endif
258# ifdef SOLVE3D
259 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
260 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
261 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
262# endif
263# ifdef MASKING
264 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
265 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
266 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
267# endif
268# ifdef SOLVE3D
269 real(r8), intent(in) :: alpha(LBi:UBi,LBj:UBj)
270 real(r8), intent(in) :: beta(LBi:UBi,LBj:UBj)
271 real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
272 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
273 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,2,N(ng))
274 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,2,N(ng))
275# endif
276 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
277# ifdef SOLVE3D
278 real(r8), intent(out) :: ad_rho(LBi:UBi,LBj:UBj,N(ng))
279# endif
280# ifdef ZETA_ELLIPTIC
281 real(r8), intent(in) :: bc_ak(Nbico(ng))
282 real(r8), intent(in) :: bc_bk(Nbico(ng))
283 real(r8), intent(in) :: zdf1(Nbico(ng))
284 real(r8), intent(in) :: zdf2(Nbico(ng))
285 real(r8), intent(in) :: zdf3(Nbico(ng))
286 real(r8), intent(inout) :: pc_r2d(LBi:UBi,LBj:UBj)
287 real(r8), intent(inout) :: r_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
288 real(r8), intent(inout) :: br_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
289 real(r8), intent(inout) :: p_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
290 real(r8), intent(inout) :: bp_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
291 real(r8), intent(inout) :: ad_rhs_r2d(LBi:UBi,LBj:UBj)
292 real(r8), intent(inout) :: ad_zeta_ref(LBi:UBi,LBj:UBj)
293# endif
294
295# endif
296!
297! Local variable declarations.
298!
299 integer :: i, j, k, order
300
301 integer :: Norder = 2 ! Shapiro filter order
302
303 real(r8) :: fac, fac1, fac2, fac3, gamma
304 real(r8) :: cff, cff1, cff2, cff3, cff4
305 real(r8) :: ad_cff, ad_cff1, ad_cff2, adfac, adfac1, adfac2
306 real(r8) :: dzdT, zphi, zphi1, zwbot, zwtop
307
308 real(r8), dimension(20) :: filter_coef = &
309 & (/ 2.500000E-1_r8, 6.250000E-2_r8, 1.562500E-2_r8, &
310 & 3.906250E-3_r8, 9.765625E-4_r8, 2.44140625E-4_r8, &
311 & 6.103515625E-5_r8, 1.5258789063E-5_r8, 3.814697E-6_r8, &
312 & 9.536743E-7_r8, 2.384186E-7_r8, 5.960464E-8_r8, &
313 & 1.490116E-8_r8, 3.725290E-9_r8, 9.313226E-10_r8, &
314 & 2.328306E-10_r8, 5.820766E-11_r8, 1.455192E-11_r8, &
315 & 3.637979E-12_r8, 9.094947E-13_r8 /)
316
317 real(r8), dimension(N(ng)) :: dSdT, dSdT_filter
318
319 real(r8), dimension(IminS:ImaxS) :: ad_phie
320 real(r8), dimension(IminS:ImaxS) :: ad_phix
321
322# ifdef SALINITY
323 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
324
325 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: dTdz
326 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: dSdz
327# endif
328 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: ad_gradP
329
330# ifdef ZETA_ELLIPTIC
331 real(r8), dimension(IminS:ImaxS,N(ng)) :: ad_phi
332
333 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_gradPx
334 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_gradPy
335 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_r2d_rhs
336# endif
337
338# include "set_bounds.h"
339!
340!-----------------------------------------------------------------------
341! Initialize adjoint private variables.
342!-----------------------------------------------------------------------
343!
344 ad_cff=0.0_r8
345 ad_cff1=0.0_r8
346 ad_cff2=0.0_r8
347 DO i=imins,imaxs
348 ad_phie(i)=0.0_r8
349 ad_phix(i)=0.0_r8
350 END DO
351 DO k=1,n(ng)
352 DO j=lbj,ubj
353 DO i=lbi,ubi
354 ad_rho(i,j,k)=0.0_r8
355 END DO
356 END DO
357 DO j=jmins,jmaxs
358 DO i=imins,imaxs
359 ad_gradp(i,j,k)=0.0_r8
360 END DO
361 END DO
362 END DO
363# ifdef ZETA_ELLIPTIC
364 DO k=1,n(ng)
365 DO i=imins,imaxs
366 ad_phi(i,k)=0.0_r8
367 END DO
368 END DO
369 DO j=jmins,jmaxs
370 DO i=imins,imaxs
371 ad_gradpx(i,j)=0.0_r8
372 ad_gradpy(i,j)=0.0_r8
373 END DO
374 END DO
375# endif
376!
377!-----------------------------------------------------------------------
378! Add balanced free-surface contribution to unbalanced free-surface
379! increment.
380!-----------------------------------------------------------------------
381!
382 IF (balance(isfsur)) THEN
383
384# ifdef DISTRIBUTE
385!^ CALL mp_exchange2d (ng, tile, iTLM, 1, &
386!^ & LBi, UBi, LBj, UBj, &
387!^ & NghostPoints, &
388!^ & EWperiodic(ng), NSperiodic(ng), &
389!^ & tl_zeta(:,:,Linp))
390!^
391 CALL ad_mp_exchange2d (ng, tile, iadm, 1, &
392 & lbi, ubi, lbj, ubj, &
393 & nghostpoints, &
394 & ewperiodic(ng), nsperiodic(ng), &
395 & ad_zeta(:,:,linp))
396
397# endif
398 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
399!^ CALL exchange_r2d_tile (ng, tile, &
400!^ & LBi, UBi, LBj, UBj, &
401!^ & tl_zeta(:,:,Linp))
402!^
403 CALL ad_exchange_r2d_tile (ng, tile, &
404 & lbi, ubi, lbj, ubj, &
405 & ad_zeta(:,:,linp))
406 END IF
407# ifdef ZETA_ELLIPTIC
408!
409! Adjoint balanced free-surface contribution to unbalanced free-surface
410! increment: Solve elliptic equation (Fukumori et al., 1998).
411!
412 DO j=jstr,jend
413 DO i=istr,iend
414!^ tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)+tl_zeta_ref(i,j)
415!^
416 ad_zeta_ref(i,j)=ad_zeta(i,j,linp)
417 END DO
418 END DO
419!
420! Call the adjoint biconjugate gradient solver.
421!
422!^ CALL tl_biconj_tile (ng, tile, iTLM, &
423!^ & LBi, UBi, LBj, UBj, &
424!^ & IminS, ImaxS, JminS, JmaxS, &
425!^ & Lbck, &
426!^ & h, pmon_u, pnom_v, pm, pn, &
427# ifdef MASKING
428!^ & umask, vmask, rmask, &
429# endif
430!^ & bc_ak, bc_bk, zdf1, zdf2, zdf3, &
431!^ & pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, &
432!^ & tl_zeta_ref, tl_rhs_r2d)
433!^
434 CALL ad_biconj_tile (ng, tile, iadm, &
435 & lbi, ubi, lbj, ubj, &
436 & imins, imaxs, jmins, jmaxs, &
437 & lbck, &
438 & h, pmon_u, pnom_v, pm, pn, &
439# ifdef MASKING
440 & umask, vmask, rmask, &
441# endif
442 & bc_ak, bc_bk, zdf1, zdf2, zdf3, &
443 & pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, &
444 & ad_zeta_ref, ad_rhs_r2d)
445!
446 DO j=jstr,jend
447 DO i=istr,iend
448!^ tl_zeta_ref(i,j)=0.0_r8
449!^
450 ad_zeta_ref(i,j)=0.0_r8
451 END DO
452 END DO
453
454# ifdef DISTRIBUTE
455!^ CALL mp_exchange2d (ng, tile, iTLM, 1, &
456!^ & LBi, UBi, LBj, UBj, &
457!^ & NghostPoints, &
458!^ & EWperiodic(ng), NSperiodic(ng), &
459!^ & tl_rhs_r2d)
460!^
461 CALL ad_mp_exchange2d (ng, tile, iadm, 1, &
462 & lbi, ubi, lbj, ubj, &
463 & nghostpoints, &
464 & ewperiodic(ng), nsperiodic(ng), &
465 & ad_rhs_r2d)
466# endif
467!^ CALL r2d_bc (ng, tile, &
468!^ & LBi, UBi, LBj, UBj, &
469!^ & tl_rhs_r2d)
470!^
471 CALL ad_r2d_bc (ng, tile, &
472 & lbi, ubi, lbj, ubj, &
473 & ad_rhs_r2d)
474!
475! Adjoint of compute the RHS term for the biconjugate gradient solver.
476!
477 DO j=jstr,jend
478 DO i=istr,iend
479# ifdef MASKING
480!^ tl_rhs_r2d(i,j)=tl_rhs_r2d(i,j)*rmask(i,j)
481!^
482 ad_rhs_r2d(i,j)=ad_rhs_r2d(i,j)*rmask(i,j)
483# endif
484!^ tl_rhs_r2d(i,j)=-pm(i,j)*pn(i,j)* &
485!^ & (pmon_u(i+1,j)*tl_gradPx(i+1,j)- &
486!^ & pmon_u(i ,j)*tl_gradPx(i ,j)+ &
487!^ & pnom_v(i,j+1)*tl_gradPy(i,j+1)- &
488!^ & pnom_v(i,j )*tl_gradPy(i,j ))
489!^
490 adfac=-pm(i,j)*pn(i,j)*ad_rhs_r2d(i,j)
491 ad_gradpx(i ,j)=ad_gradpx(i ,j)-pmon_u(i ,j)*adfac
492 ad_gradpx(i+1,j)=ad_gradpx(i+1,j)+pmon_u(i+1,j)*adfac
493 ad_gradpy(i,j )=ad_gradpy(i,j )-pnom_v(i,j )*adfac
494 ad_gradpy(i,j+1)=ad_gradpy(i,j+1)+pnom_v(i,j+1)*adfac
495 ad_rhs_r2d(i,j)=0.0_r8
496 END DO
497 END DO
498!
499! Use forward boundary routines here since they are self-adjoint.
500!
501!^ CALL v2d_bc (ng, tile, &
502!^ & IminS, ImaxS, JminS, JmaxS, &
503!^ & tl_gradPy)
504!^
505 CALL v2d_bc (ng, tile, &
506 & imins, imaxs, jmins, jmaxs, &
507 & ad_gradpy)
508!^ CALL u2d_bc (ng, tile, &
509!^ & IminS, ImaxS, JminS, JmaxS, &
510!^ & tl_gradPx)
511!^
512 CALL u2d_bc (ng, tile, &
513 & imins, imaxs, jmins, jmaxs, &
514 & ad_gradpx)
515
516# else
517!
518! Integrate hydrostatic equation from bottom to surface.
519!
520 IF (lnm_flag.eq.0) THEN
521 cff1=1.0_r8/rho0
522 DO k=1,n(ng)
523 DO j=jstr,jend
524 DO i=istr,iend
525!^ tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)+tl_cff
526!^
527 ad_cff=ad_cff+ad_zeta(i,j,linp)
528# ifdef MASKING
529!^ tl_cff=tl_cff*rmask(i,j)
530!^
531 ad_cff=ad_cff*rmask(i,j)
532# endif
533!^ tl_cff=-cff1*tl_rho(i,j,k)*Hz(i,j,k)
534!^
535 ad_rho(i,j,k)=ad_rho(i,j,k)- &
536 & cff1*hz(i,j,k)*ad_cff
537 ad_cff=0.0_r8
538 END DO
539 END DO
540 END DO
541!
542! Integrate from level of no motion (LNM_depth) to surface or
543! integrate from local bottom if shallower than LNM_depth.
544! Notice that the level of motion depth is tested against the
545! bottom of the grid cell (at W-points) and bracketed between
546! W-points during the interpolation. Also positive depths are
547! used for clarity.
548!
549 ELSE IF (lnm_flag.eq.1) THEN
550 cff1=1.0_r8/rho0
551 DO j=jstr,jend
552 DO i=istr,iend
553 DO k=1,n(ng)
554 zwtop=abs(z_w(i,j,k ))
555 zwbot=abs(z_w(i,j,k-1))
556 IF (zwbot.lt.lnm_depth(ng)) THEN
557!^ tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)+tl_cff
558!^
559 ad_cff=ad_cff+ad_zeta(i,j,linp)
560# ifdef MASKING
561!^ tl_cff=tl_cff*rmask(i,j)
562!^
563 ad_cff=ad_cff*rmask(i,j)
564# endif
565!^ tl_cff=-cff1*tl_rho(i,j,k)*Hz(i,j,k)
566!^
567 ad_rho(i,j,k)=ad_rho(i,j,k)- &
568 & cff1*hz(i,j,k)*ad_cff
569 ad_cff=0.0_r8
570 ELSE IF ((zwbot.gt.lnm_depth(ng)).and. &
571 & (lnm_depth(ng).ge.zwtop)) THEN ! interpolate
572!^ tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)+tl_cff
573!^
574 ad_cff=ad_cff+ad_zeta(i,j,linp)
575# ifdef MASKING
576!^ tl_cff=tl_cff*rmask(i,j)
577!^
578 ad_cff=ad_cff*rmask(i,j)
579# endif
580 zphi=abs(z_r(i,j,k))
581 IF (zphi.ge.lnm_depth(ng)) THEN ! above cell center
582 IF (k.eq.n(ng)) THEN
583!^ tl_cff=-cff1*tl_rho(i,j,k)*Hz(i,j,k)
584!^
585 ad_rho(i,j,k)=ad_rho(i,j,k)- &
586 & cff1*hz(i,j,k)*ad_cff
587 ad_cff=0.0_r8
588 ELSE
589 zphi1=abs(z_r(i,j,k+1))
590 fac=(lnm_depth(ng)-zphi1)/(zphi-zphi1)
591!^ tl_cff=-cff1* &
592!^ & (tl_rho(i,j,k+1)+ &
593!^ & fac*(tl_rho(i,j,k)-tl_rho(i,j,k+1)))* &
594!^ & (LNM_depth(ng)-zwtop)
595!^
596 adfac=cff1*(lnm_depth(ng)-zwtop)*ad_cff
597 adfac1=fac*adfac
598 ad_rho(i,j,k )=ad_rho(i,j,k )-adfac1
599 ad_rho(i,j,k+1)=ad_rho(i,j,k+1)-adfac+adfac1
600 ad_cff=0.0_r8
601 END IF
602 ELSE ! below cell center
603 IF (k.eq.1) THEN
604!^ tl_cff=-cff1*tl_rho(i,j,k)*Hz(i,j,k)
605!^
606 ad_rho(i,j,k)=ad_rho(i,j,k)- &
607 & cff1*hz(i,j,k)*ad_cff
608 ad_cff=0.0_r8
609 ELSE
610 zphi1=abs(z_r(i,j,k-1))
611 fac=(lnm_depth(ng)-zphi)/(zphi1-zphi)
612!^ tl_cff=-cff1* &
613!^ & (tl_rho(i,j,k)+ &
614!^ & fac*(tl_rho(i,j,k-1)-tl_rho(i,j,k)))* &
615!^ & (zwbot-LNM_depth(ng))
616!^
617 adfac=cff1*(zwbot-lnm_depth(ng))*ad_cff
618 adfac1=fac*adfac
619 ad_rho(i,j,k-1)=ad_rho(i,j,k-1)-adfac1
620 ad_rho(i,j,k )=ad_rho(i,j,k )-adfac+adfac1
621 ad_cff=0.0_r8
622 END IF
623 END IF
624 END IF
625 END DO
626 END DO
627 END DO
628 END IF
629# endif
630 END IF
631!
632!-----------------------------------------------------------------------
633! Add balanced velocity contributions to unbalanced velocity
634! increments. Use linear pressure gradient formulation based
635! to that found in routine "prsgrd31.h".
636!-----------------------------------------------------------------------
637!
638 IF (balance(isvvel)) THEN
639!
640! Compute 3D V-velocity error covariance constraint.
641!
642# ifdef DISTRIBUTE
643!^ CALL mp_exchange3d (ng, tile, iTLM, 2, &
644!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
645!^ & NghostPoints, &
646!^ & EWperiodic(ng), NSperiodic(ng), &
647!^ & tl_u(:,:,:,Linp), tl_v(:,:,:,Linp))
648!^
649 CALL ad_mp_exchange3d (ng, tile, iadm, 2, &
650 & lbi, ubi, lbj, ubj, 1, n(ng), &
651 & nghostpoints, &
652 & ewperiodic(ng), nsperiodic(ng), &
653 & ad_u(:,:,:,linp), ad_v(:,:,:,linp))
654# endif
655 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
656!^ CALL exchange_v3d_tile (ng, tile, &
657!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
658!^ & tl_v(:,:,:,Linp))
659!^
660 CALL ad_exchange_v3d_tile (ng, tile, &
661 & lbi, ubi, lbj, ubj, 1, n(ng), &
662 & ad_v(:,:,:,linp))
663!^ CALL exchange_u3d_tile (ng, tile, &
664!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
665!^ & tl_u(:,:,:,Linp))
666!^
667 CALL ad_exchange_u3d_tile (ng, tile, &
668 & lbi, ubi, lbj, ubj, 1, n(ng), &
669 & ad_u(:,:,:,linp))
670 END IF
671
672 DO k=1,n(ng)
673 DO j=jstrv,jend
674 DO i=istr,iend
675# ifdef MASKING
676!^ tl_v(i,j,k,Linp)=tl_v(i,j,k,Linp)*vmask(i,j)
677!^
678 ad_v(i,j,k,linp)=ad_v(i,j,k,linp)*vmask(i,j)
679# endif
680!^ tl_v(i,j,k,Linp)=tl_v(i,j,k,Linp)+ &
681!^ & 0.25_r8*(tl_gradP(i ,j-1,k)+ &
682!^ & tl_gradP(i+1,j-1,k)+ &
683!^ & tl_gradP(i ,j ,k)+ &
684!^ & tl_gradP(i+1,j ,k))
685!^
686 adfac=0.25_r8*ad_v(i,j,k,linp)
687 ad_gradp(i ,j-1,k)=ad_gradp(i ,j-1,k)+adfac
688 ad_gradp(i+1,j-1,k)=ad_gradp(i+1,j-1,k)+adfac
689 ad_gradp(i ,j ,k)=ad_gradp(i ,j ,k)+adfac
690 ad_gradp(i+1,j ,k)=ad_gradp(i+1,j ,k)+adfac
691 END DO
692 END DO
693 END DO
694 END IF
695!
696! NOTE: fac2=0 because the balanced component should consist of the
697! baroclinic pressure gradient only.
698!
699 fac1=0.5_r8*g/rho0
700!! fac2=g
701 fac2=0.0_r8
702 fac3=0.25_r8*g/rho0
703!
704 IF (balance(isfsur).or.balance(isvvel)) THEN
705 DO j=jstr-1,jend
706
707# ifdef ZETA_ELLIPTIC
708!
709! Compute right-hand-side term used in the elliptic equation for
710! unbalanced free-surface. Integrate from bottom to surface.
711!
712 IF (balance(isfsur)) THEN
713 DO k=1,n(ng)
714 DO i=istr,iend+1
715!^ tl_gradPx(i,j)=tl_gradPx(i,j)+tl_cff
716!^
717 ad_cff=ad_cff+ad_gradpx(i,j)
718# ifdef MASKING
719!^ tl_cff=tl_cff*umask(i,j)
720!^
721 ad_cff=ad_cff*umask(i,j)
722# endif
723!^ tl_cff=0.5_r8*(Hz(i-1,j,k)+Hz(i,j,k))*tl_phi(i,k)
724!^
725 ad_phi(i,k)=ad_phi(i,k)+ &
726 & 0.5_r8*(hz(i-1,j,k)+hz(i,j,k))*ad_cff
727 ad_cff=0.0_r8
728 END DO
729 END DO
730 END IF
731# endif
732!
733! Compute balanced, interior V-momentum from baroclinic pressure
734! gradient (differentiate and then vertically integrate).
735!
736!^ DO k=N(ng)-1,1,-1
737!^
738 DO k=1,n(ng)-1
739 DO i=istr,iend+1
740 cff1=1.0_r8/((z_r(i ,j,k+1)-z_r(i ,j,k))* &
741 & (z_r(i-1,j,k+1)-z_r(i-1,j,k)))
742 cff2=z_r(i ,j,k )-z_r(i-1,j,k )+ &
743 & z_r(i ,j,k+1)-z_r(i-1,j,k+1)
744 cff3=z_r(i ,j,k+1)-z_r(i ,j,k )- &
745 & z_r(i-1,j,k+1)+z_r(i-1,j,k )
746 gamma=0.125_r8*cff1*cff2*cff3
747!
748 cff3=z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
749 & z_r(i,j,k )-z_r(i-1,j,k )
750 cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i-1,j,k+1))+ &
751 & (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i-1,j,k ))
752
753# ifdef ZETA_ELLIPTIC
754!^ tl_phi(i,k)=tl_phix(i)
755!^
756 ad_phix(i)=ad_phix(i)+ad_phi(i,k)
757 ad_phi(i,k)=0.0_r8
758# endif
759!^ tl_gradP(i,j,k)=0.5_r8*tl_phix(i)*(pm(i-1,j)+pm(i,j))/ &
760!^ & (f(i-1,j)+f(i,j))
761!^
762 ad_phix(i)=ad_phix(i)+ &
763 & ad_gradp(i,j,k)*0.5_r8*(pm(i-1,j)+pm(i,j))/ &
764 & (f(i-1,j)+f(i,j))
765 ad_gradp(i,j,k)=0.0_r8
766!^ tl_phix(i)=tl_phix(i)+ &
767!^ & fac3*(tl_cff1*cff3-tl_cff2*cff4)
768!^
769 ad_cff1=fac3*cff3*ad_phix(i)
770 ad_cff2=-fac3*cff4*ad_phix(i)
771!^ tl_cff2=tl_rho(i,j,k+1)+tl_rho(i-1,j,k+1)- &
772!^ & tl_rho(i,j,k )-tl_rho(i-1,j,k )
773!^ tl_cff1=(1.0_r8+gamma)*(tl_rho(i ,j,k+1)- &
774!^ & tl_rho(i-1,j,k+1))+ &
775!^ & (1.0_r8-gamma)*(tl_rho(i ,j,k )- &
776!^ & tl_rho(i-1,j,k ))
777!^
778 adfac1=(1.0_r8+gamma)*ad_cff1
779 adfac2=(1.0_r8-gamma)*ad_cff1
780 ad_rho(i-1,j,k )=ad_rho(i-1,j,k )-adfac2-ad_cff2
781 ad_rho(i ,j,k )=ad_rho(i ,j,k )+adfac2-ad_cff2
782 ad_rho(i-1,j,k+1)=ad_rho(i-1,j,k+1)-adfac1+ad_cff2
783 ad_rho(i ,j,k+1)=ad_rho(i ,j,k+1)+adfac1+ad_cff2
784 ad_cff1=0.0_r8
785 ad_cff2=0.0_r8
786 END DO
787 END DO
788!
789! Compute balanced, surface V-momentum from baroclinic and barotropic
790! surface pressure gradient.
791!
792 DO i=istr,iend+1
793 cff1=z_w(i ,j,n(ng))-z_r(i ,j,n(ng))+ &
794 & z_w(i-1,j,n(ng))-z_r(i-1,j,n(ng))
795# ifdef ZETA_ELLIPTIC
796!^ tl_phi(i,N(ng))=tl_phix(i)
797!^
798 ad_phix(i)=ad_phix(i)+ad_phi(i,n(ng))
799 ad_phi(i,n(ng))=0.0_r8
800# endif
801!^ tl_gradP(i,j,N(ng))=0.5_r8*tl_phix(i)*(pm(i-1,j)+pm(i,j))/ &
802!^ & (f(i-1,j)+f(i,j))
803!^
804 ad_phix(i)=ad_phix(i)+ &
805 & ad_gradp(i,j,n(ng))*0.5_r8*(pm(i-1,j)+pm(i,j))/ &
806 & (f(i-1,j)+f(i,j))
807 ad_gradp(i,j,n(ng))=0.0_r8
808!^ tl_phix(i)=fac1*(tl_rho(i ,j,N(ng))- &
809!^ tl_rho(i-1,j,N(ng)))*cff1+ &
810!^ & fac2*(tl_zeta(i ,j,Linp)- &
811!^ tl_zeta(i-1,j,Linp))
812!^
813 adfac1=fac1*cff1*ad_phix(i)
814 adfac2=fac2*ad_phix(i)
815 ad_rho(i-1,j,n(ng))=ad_rho(i-1,j,n(ng))-adfac1
816 ad_rho(i ,j,n(ng))=ad_rho(i ,j,n(ng))+adfac1
817 ad_zeta(i-1,j,linp)=ad_zeta(i-1,j,linp)-adfac2
818 ad_zeta(i ,j,linp)=ad_zeta(i ,j,linp)+adfac2
819 ad_phix(i)=0.0_r8
820 END DO
821 END DO
822 END IF
823!
824! Compute 3D U-velocity error covariance constraint.
825!
826 IF (balance(isuvel)) THEN
827 DO k=1,n(ng)
828 DO j=jstr,jend
829 DO i=istru,iend
830# ifdef MASKING
831!^ tl_u(i,j,k,Linp)=tl_u(i,j,k,Linp)*umask(i,j)
832!^
833 ad_u(i,j,k,linp)=ad_u(i,j,k,linp)*umask(i,j)
834# endif
835!^ tl_u(i,j,k,Linp)=tl_u(i,j,k,Linp)- &
836!^ & 0.25_r8*(tl_gradP(i-1,j ,k)+ &
837!^ & tl_gradP(i ,j ,k)+ &
838!^ & tl_gradP(i-1,j+1,k)+ &
839!^ & tl_gradP(i ,j+1,k))
840!^
841 adfac=0.25_r8*ad_u(i,j,k,linp)
842 ad_gradp(i-1,j ,k)=ad_gradp(i-1,j ,k)-adfac
843 ad_gradp(i ,j ,k)=ad_gradp(i ,j ,k)-adfac
844 ad_gradp(i-1,j+1,k)=ad_gradp(i-1,j+1,k)-adfac
845 ad_gradp(i ,j+1,k)=ad_gradp(i ,j+1,k)-adfac
846 END DO
847 END DO
848 END DO
849 END IF
850!
851 IF (balance(isfsur).or.balance(isuvel)) THEN
852 DO j=jstr,jend+1
853
854# ifdef ZETA_ELLIPTIC
855!
856! Compute right-hand-side term used in the elliptic equation for
857! unbalanced free-surface. Integrate from bottom to surface.
858!
859 IF (balance(isfsur)) THEN
860 DO k=1,n(ng)
861 DO i=istr-1,iend
862!^ tl_gradPy(i,j)=tl_gradPy(i,j)+tl_cff
863!^
864 ad_cff=ad_cff+ad_gradpy(i,j)
865# ifdef MASKING
866!^ tl_cff=tl_cff*vmask(i,j)
867!^
868 ad_cff=ad_cff*vmask(i,j)
869# endif
870!^ tl_cff=0.5_r8*(Hz(i,j-1,k)+Hz(i,j,k))*tl_phi(i,k)
871!^
872 ad_phi(i,k)=ad_phi(i,k)+ &
873 & 0.5_r8*(hz(i,j-1,k)+hz(i,j,k))*ad_cff
874 ad_cff=0.0_r8
875 END DO
876 END DO
877 END IF
878# endif
879!
880! Compute balanced, interior U-momentum from baroclinic pressure
881! gradient (differentiate and then vertically integrate).
882!
883!^ DO k=N(ng)-1,1,-1
884!^
885 DO k=1,n(ng)-1
886 DO i=istr-1,iend
887 cff1=1.0_r8/((z_r(i,j ,k+1)-z_r(i,j ,k))* &
888 & (z_r(i,j-1,k+1)-z_r(i,j-1,k)))
889 cff2=z_r(i,j ,k )-z_r(i,j-1,k )+ &
890 & z_r(i,j ,k+1)-z_r(i,j-1,k+1)
891 cff3=z_r(i,j ,k+1)-z_r(i,j ,k )- &
892 & z_r(i,j-1,k+1)+z_r(i,j-1,k )
893 gamma=0.125_r8*cff1*cff2*cff3
894!
895 cff3=z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
896 & z_r(i,j,k )-z_r(i,j-1,k )
897 cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i,j-1,k+1))+ &
898 & (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i,j-1,k ))
899# ifdef ZETA_ELLIPTIC
900!^ tl_phi(i,k)=tl_phie(i)
901!^
902 ad_phie(i)=ad_phie(i)+ad_phi(i,k)
903 ad_phi(i,k)=0.0_r8
904# endif
905!^ tl_gradP(i,j,k)=0.5_r8*tl_phie(i)*(pn(i,j-1)+pn(i,j))/ &
906!^ & (f(i,j-1)+f(i,j))
907!^
908 ad_phie(i)=ad_phie(i)+ &
909 & ad_gradp(i,j,k)*0.5_r8*(pn(i,j-1)+pn(i,j))/ &
910 & (f(i,j-1)+f(i,j))
911 ad_gradp(i,j,k)=0.0_r8
912!^ tl_phie(i)=tl_phie(i)+ &
913!^ & fac3*(tl_cff1*cff3-tl_cff2*cff4)
914!^
915 ad_cff1=fac3*cff3*ad_phie(i)
916 ad_cff2=-fac3*cff4*ad_phie(i)
917!^ tl_cff2=tl_rho(i,j,k+1)+tl_rho(i,j-1,k+1)- &
918!^ & tl_rho(i,j,k )-tl_rho(i,j-1,k )
919!^ tl_cff1=(1.0_r8+gamma)*(tl_rho(i,j ,k+1)- &
920!^ & tl_rho(i,j-1,k+1))+ &
921!^ & (1.0_r8-gamma)*(tl_rho(i,j ,k )- &
922!^ & tl_rho(i,j-1,k ))
923!^
924 adfac1=(1.0_r8+gamma)*ad_cff1
925 adfac2=(1.0_r8-gamma)*ad_cff1
926 ad_rho(i,j-1,k )=ad_rho(i,j-1,k )-adfac2-ad_cff2
927 ad_rho(i,j ,k )=ad_rho(i,j ,k )+adfac2-ad_cff2
928 ad_rho(i,j-1,k+1)=ad_rho(i,j-1,k+1)-adfac1+ad_cff2
929 ad_rho(i,j ,k+1)=ad_rho(i,j ,k+1)+adfac1+ad_cff2
930 ad_cff1=0.0_r8
931 ad_cff2=0.0_r8
932 END DO
933 END DO
934!
935! Compute balanced, surface U-momentum from baroclinic and barotropic
936! surface pressure gradient.
937!
938 DO i=istr-1,iend
939 cff1=z_w(i,j ,n(ng))-z_r(i,j ,n(ng))+ &
940 & z_w(i,j-1,n(ng))-z_r(i,j-1,n(ng))
941# ifdef ZETA_ELLIPTIC
942!^ tl_phi(i,N(ng))=tl_phie(i)
943!^
944 ad_phie(i)=ad_phie(i)+ad_phi(i,n(ng))
945 ad_phi(i,n(ng))=0.0_r8
946# endif
947!^ tl_gradP(i,j,N(ng))=0.5_r8*tl_phie(i)*(pn(i,j-1)+pn(i,j))/ &
948!^ & (f(i,j-1)+f(i,j))
949!^
950 ad_phie(i)=ad_phie(i)+ &
951 & ad_gradp(i,j,n(ng))*0.5_r8*(pn(i,j-1)+pn(i,j))/ &
952 & (f(i,j-1)+f(i,j))
953 ad_gradp(i,j,n(ng))=0.0_r8
954!^ tl_phie(i)=fac1*(tl_rho(i,j ,N(ng))- &
955!^ & tl_rho(i,j-1,N(ng)))*cff1+ &
956!^ & fac2*(tl_zeta(i,j ,Linp)- &
957!^ & tl_zeta(i,j-1,Linp))
958!^
959 adfac1=fac1*cff1*ad_phie(i)
960 adfac2=fac2*ad_phie(i)
961 ad_rho(i,j-1,n(ng))=ad_rho(i,j-1,n(ng))-adfac1
962 ad_rho(i,j ,n(ng))=ad_rho(i,j ,n(ng))+adfac1
963 ad_zeta(i,j-1,linp)=ad_zeta(i,j-1,linp)-adfac2
964 ad_zeta(i,j ,linp)=ad_zeta(i,j ,linp)+adfac2
965 ad_phie(i)=0.0_r8
966 END DO
967 END DO
968
969# ifdef ZETA_ELLIPTIC
970!
971! Initialize vertical intergal local variables.
972!
973 DO j=jmins,jmaxs
974 DO i=imins,imaxs
975!^ tl_gradPy(i,j)=0.0_r8
976!^
977 ad_gradpy(i,j)=0.0_r8
978!^ tl_gradPx(i,j)=0.0_r8
979!^
980 ad_gradpx(i,j)=0.0_r8
981 END DO
982 END DO
983# endif
984 END IF
985!
986!-----------------------------------------------------------------------
987! Compute balanced density anomaly increment using linearized equation
988! of state. The thermal expansion and saline contraction coefficients
989! are computed from the background state.
990!-----------------------------------------------------------------------
991!
992 IF (balance(isfsur).or.balance(isvvel)) THEN
993
994# ifdef DISTRIBUTE
995!^ CALL mp_exchange3d (ng, tile, iTLM, 1, &
996!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
997!^ & NghostPoints, &
998!^ & EWperiodic(ng), NSperiodic(ng), &
999!^ & tl_rho)
1000!^
1001 CALL ad_mp_exchange3d (ng, tile, iadm, 1, &
1002 & lbi, ubi, lbj, ubj, 1, n(ng), &
1003 & nghostpoints, &
1004 & ewperiodic(ng), nsperiodic(ng), &
1005 & ad_rho)
1006# endif
1007 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1008!^ CALL exchange_r3d_tile (ng, tile, &
1009!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
1010!^ & tl_rho)
1011!^
1012 CALL ad_exchange_r3d_tile (ng, tile, &
1013 & lbi, ubi, lbj, ubj, 1, n(ng), &
1014 & ad_rho)
1015 END IF
1016
1017 DO j=jstrr,jendr
1018 DO k=1,n(ng)
1019 DO i=istrr,iendr
1020# ifdef MASKING
1021!^ tl_rho(i,j,k)=tl_rho(i,j,k)*rmask(i,j)
1022!^
1023 ad_rho(i,j,k)=ad_rho(i,j,k)*rmask(i,j)
1024# endif
1025# ifdef SALINITY
1026!^ tl_rho(i,j,k)=tl_rho(i,j,k)+ &
1027!^ & rho0*beta(i,j)*tl_t(i,j,k,Linp,isalt)
1028!^
1029 ad_t(i,j,k,linp,isalt)=ad_t(i,j,k,linp,isalt)+ &
1030 & rho0*beta(i,j)*ad_rho(i,j,k)
1031# endif
1032!^ tl_rho(i,j,k)=-rho0*alpha(i,j)*tl_t(i,j,k,Linp,itemp)
1033!^
1034 ad_t(i,j,k,linp,itemp)=ad_t(i,j,k,linp,itemp)- &
1035 & rho0*alpha(i,j)*ad_rho(i,j,k)
1036 ad_rho(i,j,k)=0.0
1037 END DO
1038 END DO
1039 END DO
1040 END IF
1041
1042# ifdef SALINITY
1043!
1044!-----------------------------------------------------------------------
1045! Compute balance salinity contribution.
1046!-----------------------------------------------------------------------
1047!
1048 IF (balance(istvar(isalt))) THEN
1049
1050# ifdef DISTRIBUTE
1051!^ CALL mp_exchange3d (ng, tile, iTLM, 1, &
1052!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
1053!^ & NghostPoints, &
1054!^ & EWperiodic(ng), NSperiodic(ng), &
1055!^ & tl_t(:,:,:,Linp,isalt))
1056!^
1057 CALL ad_mp_exchange3d (ng, tile, iadm, 1, &
1058 & lbi, ubi, lbj, ubj, 1, n(ng), &
1059 & nghostpoints, &
1060 & ewperiodic(ng), nsperiodic(ng), &
1061 & ad_t(:,:,:,linp,isalt))
1062# endif
1063 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1064!^ CALL exchange_r3d_tile (ng, tile, &
1065!^ & LBi, UBi, LBj, UBj, 1, N(ng), &
1066!^ & tl_t(:,:,:,Linp,isalt))
1067!^
1068 CALL ad_exchange_r3d_tile (ng, tile, &
1069 & lbi, ubi, lbj, ubj, 1, n(ng), &
1070 & ad_t(:,:,:,linp,isalt))
1071 END IF
1072!
1073! Compute temperature (dTdz) and salinity (dSdz) shears.
1074!
1075 DO j=jstrr,jendr
1076 DO i=istrr,iendr
1077 fc(i,0)=0.0_r8
1078 dtdz(i,j,0)=0.0_r8
1079 dsdz(i,j,0)=0.0_r8
1080 END DO
1081 DO k=1,n(ng)-1
1082 DO i=istrr,iendr
1083 cff=1.0_r8/(2.0_r8*hz(i,j,k+1)+ &
1084 & hz(i,j,k)*(2.0_r8-fc(i,k-1)))
1085 fc(i,k)=cff*hz(i,j,k+1)
1086 dtdz(i,j,k)=cff*(6.0_r8*(t(i,j,k+1,lbck,itemp)- &
1087 & t(i,j,k ,lbck,itemp))- &
1088 & hz(i,j,k)*dtdz(i,j,k-1))
1089 dsdz(i,j,k)=cff*(6.0_r8*(t(i,j,k+1,lbck,isalt)- &
1090 & t(i,j,k ,lbck,isalt))- &
1091 & hz(i,j,k)*dsdz(i,j,k-1))
1092 END DO
1093 END DO
1094 DO i=istrr,iendr
1095 dtdz(i,j,n(ng))=0.0_r8
1096 dsdz(i,j,n(ng))=0.0_r8
1097 END DO
1098 DO k=n(ng)-1,1,-1
1099 DO i=istrr,iendr
1100 dtdz(i,j,k)=dtdz(i,j,k)-fc(i,k)*dtdz(i,j,k+1)
1101 dsdz(i,j,k)=dsdz(i,j,k)-fc(i,k)*dsdz(i,j,k+1)
1102 END DO
1103 END DO
1104 END DO
1105!
1106! Add balanced salinity (deltaS_b) contribution to unbalanced salinity
1107! increment. The unbalanced salinity increment is related related to
1108! temperature increment:
1109!
1110! deltaS_b = cff * dS/dz * dz/dT * deltaT
1111!
1112! Here, cff is a coefficient that depends on the mixed layer depth:
1113!
1114! cff = 1.0 - EXP (z_r / ml_depth)
1115!
1116! the coefficient is smoothly reduced to zero at the surface and below
1117! the mixed layer.
1118!
1119 DO j=jstrr,jendr
1120 DO i=istrr,iendr
1121 DO k=1,n(ng)
1122 cff=0.5_r8*(dtdz(i,j,k-1)+dtdz(i,j,k))
1123 IF (abs(cff).lt.dtdz_min(ng)) THEN
1124 dzdt=0.0_r8
1125 ELSE
1126 dzdt=1.0_r8/cff
1127 END IF
1128 dsdt(k)=(0.5_r8*(dsdz(i,j,k-1)+ &
1129 & dsdz(i,j,k )))*dzdt
1130 END DO
1131!
1132! Shapiro filter.
1133!
1134 DO order=1,norder/2
1135 IF (order.ne.norder/2) THEN
1136 dsdt_filter(1)=2.0_r8*(dsdt(1)-dsdt(2))
1137 dsdt_filter(n(ng))=2.0_r8*(dsdt(n(ng))-dsdt(n(ng)-1))
1138 ELSE
1139 dsdt_filter(1)=0.0_r8
1140 dsdt_filter(n(ng))=0.0_r8
1141 END IF
1142 DO k=2,n(ng)-1
1143 dsdt_filter(k)=2.0_r8*dsdt(k)-dsdt(k-1)-dsdt(k+1)
1144 END DO
1145 DO k=1,n(ng)
1146 dsdt(k)=dsdt(k)-filter_coef(norder/2)*dsdt_filter(k)
1147 END DO
1148 END DO
1149
1150 DO k=1,n(ng)
1151 cff=(1.0_r8-exp(z_r(i,j,k)/ml_depth(ng)))*dsdt(k)
1152# ifdef MASKING
1153!^ tl_t(i,j,k,Linp,isalt)=tl_t(i,j,k,Linp,isalt)*rmask(i,j)
1154!^
1155 ad_t(i,j,k,linp,isalt)=ad_t(i,j,k,linp,isalt)*rmask(i,j)
1156# endif
1157!^ tl_t(i,j,k,Linp,isalt)=tl_t(i,j,k,Linp,isalt)+ &
1158!^ & cff*tl_t(i,j,k,Linp,itemp)
1159!^
1160 ad_t(i,j,k,linp,itemp)=ad_t(i,j,k,linp,itemp)+ &
1161 & cff*ad_t(i,j,k,linp,isalt)
1162 END DO
1163 END DO
1164 END DO
1165 END IF
1166# endif
1167
1168 RETURN
1169 END SUBROUTINE ad_balance_tile
1170
1171#endif
1172 END MODULE ad_balance_mod
ad_balance_mod
Definition ad_balance.F:2
ad_balance_mod::ad_balance
subroutine, public ad_balance(ng, tile, lbck, linp)
Definition ad_balance.F:59
ad_balance_mod::ad_balance_tile
subroutine ad_balance_tile(ng, tile, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, lbck, linp, f, pm, pn, pmon_u, pnom_v, h, hz, z_r, z_w, rmask, umask, vmask, alpha, beta, t, ad_zeta_ref, ad_rhs_r2d, pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, zdf1, zdf2, zdf3, bc_ak, bc_bk, ad_rho, ad_t, ad_u, ad_v, ad_zeta)
Definition ad_balance.F:178
ad_bc_2d_mod
Definition ad_bc_2d.F:2
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
ad_exchange_3d_mod
Definition ad_exchange_3d.F:2
ad_exchange_3d_mod::ad_exchange_v3d_tile
subroutine ad_exchange_v3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, ad_a)
Definition ad_exchange_3d.F:1283
ad_exchange_3d_mod::ad_exchange_r3d_tile
subroutine ad_exchange_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, ad_a)
Definition ad_exchange_3d.F:461
ad_exchange_3d_mod::ad_exchange_u3d_tile
subroutine ad_exchange_u3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, ad_a)
Definition ad_exchange_3d.F:872
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_kinds
Definition mod_kinds.F:2
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ncparam
Definition mod_ncparam.F:2
mod_ncparam::isvvel
integer isvvel
Definition mod_ncparam.F:507
mod_ncparam::istvar
integer, dimension(:), allocatable istvar
Definition mod_ncparam.F:520
mod_ncparam::isuvel
integer isuvel
Definition mod_ncparam.F:506
mod_ncparam::isfsur
integer isfsur
Definition mod_ncparam.F:503
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::iadm
integer, parameter iadm
Definition mod_param.F:665
mod_scalars
Definition mod_scalars.F:2
mod_scalars::ml_depth
real(r8), dimension(:), allocatable ml_depth
Definition mod_scalars.F:1295
mod_scalars::dtdz_min
real(r8), dimension(:), allocatable dtdz_min
Definition mod_scalars.F:1294
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::balance
logical, dimension(:), allocatable balance
Definition mod_scalars.F:1314
mod_scalars::lnm_depth
real(r8), dimension(:), allocatable lnm_depth
Definition mod_scalars.F:1300
mod_scalars::isalt
integer isalt
Definition mod_scalars.F:127
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::lnm_flag
integer lnm_flag
Definition mod_scalars.F:1309
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_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
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
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::v2d_bc
subroutine, public v2d_bc(ng, tile, lbi, ubi, lbj, ubj, a)
Definition zeta_balance.F:1253
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::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