ROMS
Loading...
Searching...
No Matches
tl_balance.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE tl_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, 2006: !
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 :: tl_balance
54!
55 CONTAINS
56!
57!***********************************************************************
58 SUBROUTINE tl_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, itlm)
94 nrhs(ng)=lbck
95 CALL rho_eos (ng, tile, itlm)
96!
97# endif
98 CALL tl_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) % tl_zeta_ref, &
127 & fourdvar(ng) % tl_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) % tl_rho, &
141 & ocean(ng) % tl_t, &
142 & ocean(ng) % tl_u, &
143 & ocean(ng) % tl_v, &
144# endif
145 & ocean(ng) % tl_zeta)
146
147 RETURN
148 END SUBROUTINE tl_balance
149!
150!***********************************************************************
151 SUBROUTINE tl_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 & tl_zeta_ref, tl_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 & tl_rho, tl_t, tl_u, tl_v, &
176# endif
177 & tl_zeta)
178!***********************************************************************
179!
180 USE mod_param
181 USE mod_ncparam
182 USE mod_scalars
183!
184 USE bc_2d_mod
185 USE exchange_2d_mod
186 USE exchange_3d_mod
187# ifdef DISTRIBUTE
188 USE mp_exchange_mod, ONLY : mp_exchange2d, mp_exchange3d
189# endif
190# ifdef ZETA_ELLIPTIC
191 USE zeta_balance_mod, ONLY: r2d_bc, u2d_bc, v2d_bc
192 USE zeta_balance_mod, ONLY: tl_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) :: tl_t(LBi:,LBj:,:,:,:)
226 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
227 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
228# endif
229 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
230# ifdef SOLVE3D
231 real(r8), intent(out) :: tl_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) :: tl_rhs_r2d(LBi:,LBj:)
245 real(r8), intent(inout) :: tl_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) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
273 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,2,N(ng))
274 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,2,N(ng))
275# endif
276 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
277# ifdef SOLVE3D
278 real(r8), intent(out) :: tl_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) :: tl_rhs_r2d(LBi:UBi,LBj:UBj)
292 real(r8), intent(inout) :: tl_zeta_ref(LBi:UBi,LBj:UBj)
293# endif
294
295# endif
296!
297! Local variable declarations.
298!
299 integer :: i, j, k, order
300 integer :: Norder = 2 ! Shapiro filter order
301
302 real(r8) :: fac, fac1, fac2, fac3, gamma
303 real(r8) :: cff, cff1, cff2, cff3, cff4
304 real(r8) :: tl_cff, tl_cff1, tl_cff2
305 real(r8) :: dzdT, zphi, zphi1, zwbot, zwtop
306
307 real(r8), dimension(20) :: filter_coef = &
308 & (/ 2.500000E-1_r8, 6.250000E-2_r8, 1.562500E-2_r8, &
309 & 3.906250E-3_r8, 9.765625E-4_r8, 2.44140625E-4_r8, &
310 & 6.103515625E-5_r8, 1.5258789063E-5_r8, 3.814697E-6_r8, &
311 & 9.536743E-7_r8, 2.384186E-7_r8, 5.960464E-8_r8, &
312 & 1.490116E-8_r8, 3.725290E-9_r8, 9.313226E-10_r8, &
313 & 2.328306E-10_r8, 5.820766E-11_r8, 1.455192E-11_r8, &
314 & 3.637979E-12_r8, 9.094947E-13_r8 /)
315
316 real(r8), dimension(N(ng)) :: dSdT, dSdT_filter
317
318 real(r8), dimension(IminS:ImaxS) :: tl_phie
319 real(r8), dimension(IminS:ImaxS) :: tl_phix
320
321# ifdef SALINITY
322 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
323
324 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: dTdz
325 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: dSdz
326# endif
327 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: tl_gradP
328
329# ifdef ZETA_ELLIPTIC
330 real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_phi
331
332 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_gradPx
333 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_gradPy
334 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_r2d_rhs
335# endif
336
337# include "set_bounds.h"
338
339# ifdef SALINITY
340!
341!-----------------------------------------------------------------------
342! Compute balance salinity contribution.
343!-----------------------------------------------------------------------
344!
345 IF (balance(istvar(isalt))) THEN
346!
347! Compute temperature (dTdz) and salinity (dSdz) shears.
348!
349 DO j=jstrr,jendr
350 DO i=istrr,iendr
351 fc(i,0)=0.0_r8
352 dtdz(i,j,0)=0.0_r8
353 dsdz(i,j,0)=0.0_r8
354 END DO
355 DO k=1,n(ng)-1
356 DO i=istrr,iendr
357 cff=1.0_r8/(2.0_r8*hz(i,j,k+1)+ &
358 & hz(i,j,k)*(2.0_r8-fc(i,k-1)))
359 fc(i,k)=cff*hz(i,j,k+1)
360 dtdz(i,j,k)=cff*(6.0_r8*(t(i,j,k+1,lbck,itemp)- &
361 & t(i,j,k ,lbck,itemp))- &
362 & hz(i,j,k)*dtdz(i,j,k-1))
363 dsdz(i,j,k)=cff*(6.0_r8*(t(i,j,k+1,lbck,isalt)- &
364 & t(i,j,k ,lbck,isalt))- &
365 & hz(i,j,k)*dsdz(i,j,k-1))
366 END DO
367 END DO
368 DO i=istrr,iendr
369 dtdz(i,j,n(ng))=0.0_r8
370 dsdz(i,j,n(ng))=0.0_r8
371 END DO
372 DO k=n(ng)-1,1,-1
373 DO i=istrr,iendr
374 dtdz(i,j,k)=dtdz(i,j,k)-fc(i,k)*dtdz(i,j,k+1)
375 dsdz(i,j,k)=dsdz(i,j,k)-fc(i,k)*dsdz(i,j,k+1)
376 END DO
377 END DO
378 END DO
379!
380! Add balanced salinity (deltaS_b) contribution to unbalanced salinity
381! increment. The unbalanced salinity increment is related related to
382! temperature increment:
383!
384! deltaS_b = cff * dS/dz * dz/dT * deltaT
385!
386! Here, cff is a coefficient that depends on the mixed layer depth:
387!
388! cff = 1.0 - EXP (z_r / ml_depth)
389!
390! the coefficient is smoothly reduced to zero at the surface and below
391! the mixed layer.
392!
393 DO j=jstrr,jendr
394 DO i=istrr,iendr
395 DO k=1,n(ng)
396 cff=0.5_r8*(dtdz(i,j,k-1)+dtdz(i,j,k))
397 IF (abs(cff).lt.dtdz_min(ng)) THEN
398 dzdt=0.0_r8
399 ELSE
400 dzdt=1.0_r8/cff
401 END IF
402 dsdt(k)=(0.5_r8*(dsdz(i,j,k-1)+ &
403 & dsdz(i,j,k )))*dzdt
404 END DO
405!
406! Shapiro filter.
407!
408 DO order=1,norder/2
409 IF (order.ne.norder/2) THEN
410 dsdt_filter(1)=2.0_r8*(dsdt(1)-dsdt(2))
411 dsdt_filter(n(ng))=2.0_r8*(dsdt(n(ng))-dsdt(n(ng)-1))
412 ELSE
413 dsdt_filter(1)=0.0_r8
414 dsdt_filter(n(ng))=0.0_r8
415 END IF
416 DO k=2,n(ng)-1
417 dsdt_filter(k)=2.0_r8*dsdt(k)-dsdt(k-1)-dsdt(k+1)
418 END DO
419 DO k=1,n(ng)
420 dsdt(k)=dsdt(k)-filter_coef(norder/2)*dsdt_filter(k)
421 END DO
422 END DO
423
424 DO k=1,n(ng)
425 cff=(1.0_r8-exp(z_r(i,j,k)/ml_depth(ng)))*dsdt(k)
426 tl_t(i,j,k,linp,isalt)=tl_t(i,j,k,linp,isalt)+ &
427 & cff*tl_t(i,j,k,linp,itemp)
428# ifdef MASKING
429 tl_t(i,j,k,linp,isalt)=tl_t(i,j,k,linp,isalt)*rmask(i,j)
430# endif
431 END DO
432 END DO
433 END DO
434
435 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
436 CALL exchange_r3d_tile (ng, tile, &
437 & lbi, ubi, lbj, ubj, 1, n(ng), &
438 & tl_t(:,:,:,linp,isalt))
439 END IF
440# ifdef DISTRIBUTE
441 CALL mp_exchange3d (ng, tile, itlm, 1, &
442 & lbi, ubi, lbj, ubj, 1, n(ng), &
443 & nghostpoints, &
444 & ewperiodic(ng), nsperiodic(ng), &
445 & tl_t(:,:,:,linp,isalt))
446# endif
447 END IF
448# endif
449!
450!-----------------------------------------------------------------------
451! Compute balanced density anomaly increment using linearized equation
452! of state. The thermal expansion and saline contraction coefficients
453! are computed from the background state.
454!-----------------------------------------------------------------------
455!
456 IF (balance(isfsur).or.balance(isvvel)) THEN
457 DO j=jstrr,jendr
458 DO k=1,n(ng)
459 DO i=istrr,iendr
460 tl_rho(i,j,k)=-rho0*alpha(i,j)*tl_t(i,j,k,linp,itemp)
461# ifdef SALINITY
462 tl_rho(i,j,k)=tl_rho(i,j,k)+ &
463 & rho0*beta(i,j)*tl_t(i,j,k,linp,isalt)
464# endif
465# ifdef MASKING
466 tl_rho(i,j,k)=tl_rho(i,j,k)*rmask(i,j)
467# endif
468 END DO
469 END DO
470 END DO
471
472 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
473 CALL exchange_r3d_tile (ng, tile, &
474 & lbi, ubi, lbj, ubj, 1, n(ng), &
475 & tl_rho)
476 END IF
477# ifdef DISTRIBUTE
478 CALL mp_exchange3d (ng, tile, itlm, 1, &
479 & lbi, ubi, lbj, ubj, 1, n(ng), &
480 & nghostpoints, &
481 & ewperiodic(ng), nsperiodic(ng), &
482 & tl_rho)
483# endif
484 END IF
485!
486!-----------------------------------------------------------------------
487! Add balanced velocity contributions to unbalanced velocity
488! increments. Use linear pressure gradient formulation based
489! to that found in routine "prsgrd31.h".
490!-----------------------------------------------------------------------
491!
492! NOTE: fac2=0 because the balanced component should consist of the
493! baroclinic pressure gradient only.
494!
495 fac1=0.5_r8*g/rho0
496!! fac2=g
497 fac2=0.0_r8
498 fac3=0.25_r8*g/rho0
499
500# ifdef ZETA_ELLIPTIC
501!
502! Initialize vertical intergal local variables.
503!
504 DO j=jmins,jmaxs
505 DO i=imins,imaxs
506 tl_gradpx(i,j)=0.0_r8
507 tl_gradpy(i,j)=0.0_r8
508 END DO
509 END DO
510# endif
511!
512! Compute balanced, surface U-momentum from baroclinic and barotropic
513! surface pressure gradient.
514!
515 IF (balance(isfsur).or.balance(isuvel)) THEN
516 DO j=jstr,jend+1
517 DO i=istr-1,iend
518 cff1=z_w(i,j ,n(ng))-z_r(i,j ,n(ng))+ &
519 & z_w(i,j-1,n(ng))-z_r(i,j-1,n(ng))
520 tl_phie(i)=fac1*(tl_rho(i,j ,n(ng))- &
521 & tl_rho(i,j-1,n(ng)))*cff1+ &
522 & fac2*(tl_zeta(i,j ,linp)- &
523 & tl_zeta(i,j-1,linp))
524 tl_gradp(i,j,n(ng))=0.5_r8*tl_phie(i)* &
525 & (pn(i,j-1)+pn(i,j))/(f(i,j-1)+f(i,j))
526# ifdef ZETA_ELLIPTIC
527 tl_phi(i,n(ng))=tl_phie(i)
528# endif
529 END DO
530!
531! Compute balance, interior U-momentum from baroclinic pressure
532! gradient (differentiate and then vertically integrate).
533!
534 DO k=n(ng)-1,1,-1
535 DO i=istr-1,iend
536 cff1=1.0_r8/((z_r(i,j ,k+1)-z_r(i,j ,k))* &
537 & (z_r(i,j-1,k+1)-z_r(i,j-1,k)))
538 cff2=z_r(i,j ,k )-z_r(i,j-1,k )+ &
539 & z_r(i,j ,k+1)-z_r(i,j-1,k+1)
540 cff3=z_r(i,j ,k+1)-z_r(i,j ,k )- &
541 & z_r(i,j-1,k+1)+z_r(i,j-1,k )
542 gamma=0.125_r8*cff1*cff2*cff3
543!
544 tl_cff1=(1.0_r8+gamma)*(tl_rho(i,j ,k+1)- &
545 & tl_rho(i,j-1,k+1))+ &
546 & (1.0_r8-gamma)*(tl_rho(i,j ,k )- &
547 & tl_rho(i,j-1,k ))
548 tl_cff2=tl_rho(i,j,k+1)+tl_rho(i,j-1,k+1)- &
549 & tl_rho(i,j,k )-tl_rho(i,j-1,k )
550 cff3=z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
551 & z_r(i,j,k )-z_r(i,j-1,k )
552 cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i,j-1,k+1))+ &
553 & (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i,j-1,k ))
554 tl_phie(i)=tl_phie(i)+ &
555 & fac3*(tl_cff1*cff3-tl_cff2*cff4)
556 tl_gradp(i,j,k)=0.5_r8*tl_phie(i)* &
557 & (pn(i,j-1)+pn(i,j))/(f(i,j-1)+f(i,j))
558# ifdef ZETA_ELLIPTIC
559 tl_phi(i,k)=tl_phie(i)
560# endif
561 END DO
562 END DO
563
564# ifdef ZETA_ELLIPTIC
565!
566! Compute right-hand-side term used in the elliptic equation for
567! unbalanced free-surface. Integrate from bottom to surface.
568!
569 IF (balance(isfsur)) THEN
570 DO k=1,n(ng)
571 DO i=istr-1,iend
572 tl_cff=0.5_r8*(hz(i,j-1,k)+hz(i,j,k))*tl_phi(i,k)
573# ifdef MASKING
574 tl_cff=tl_cff*vmask(i,j)
575# endif
576 tl_gradpy(i,j)=tl_gradpy(i,j)+tl_cff
577 END DO
578 END DO
579 END IF
580# endif
581 END DO
582 END IF
583!
584! Compute 3D U-velocity error covariance constraint.
585!
586 IF (balance(isuvel)) THEN
587 DO k=1,n(ng)
588 DO j=jstr,jend
589 DO i=istru,iend
590 tl_u(i,j,k,linp)=tl_u(i,j,k,linp)- &
591 & 0.25_r8*(tl_gradp(i-1,j ,k)+ &
592 & tl_gradp(i ,j ,k)+ &
593 & tl_gradp(i-1,j+1,k)+ &
594 & tl_gradp(i ,j+1,k))
595# ifdef MASKING
596 tl_u(i,j,k,linp)=tl_u(i,j,k,linp)*umask(i,j)
597# endif
598 END DO
599 END DO
600 END DO
601 END IF
602!
603! Compute balanced, surface V-momentum from baroclinic and barotropic
604! surface pressure gradient.
605!
606 IF (balance(isfsur).or.balance(isvvel)) THEN
607 DO j=jstr-1,jend
608 DO i=istr,iend+1
609 cff1=z_w(i ,j,n(ng))-z_r(i ,j,n(ng))+ &
610 & z_w(i-1,j,n(ng))-z_r(i-1,j,n(ng))
611 tl_phix(i)=fac1*(tl_rho(i ,j,n(ng))- &
612 & tl_rho(i-1,j,n(ng)))*cff1+ &
613 & fac2*(tl_zeta(i ,j,linp)- &
614 & tl_zeta(i-1,j,linp))
615 tl_gradp(i,j,n(ng))=0.5_r8*tl_phix(i)* &
616 & (pm(i-1,j)+pm(i,j))/(f(i-1,j)+f(i,j))
617# ifdef ZETA_ELLIPTIC
618 tl_phi(i,n(ng))=tl_phix(i)
619# endif
620 END DO
621!
622! Compute balance, interior V-momentum from baroclinic pressure
623! gradient (differentiate and then vertically integrate).
624!
625 DO k=n(ng)-1,1,-1
626 DO i=istr,iend+1
627 cff1=1.0_r8/((z_r(i ,j,k+1)-z_r(i ,j,k))* &
628 & (z_r(i-1,j,k+1)-z_r(i-1,j,k)))
629 cff2=z_r(i ,j,k )-z_r(i-1,j,k )+ &
630 & z_r(i ,j,k+1)-z_r(i-1,j,k+1)
631 cff3=z_r(i ,j,k+1)-z_r(i ,j,k )- &
632 & z_r(i-1,j,k+1)+z_r(i-1,j,k )
633 gamma=0.125_r8*cff1*cff2*cff3
634!
635 tl_cff1=(1.0_r8+gamma)*(tl_rho(i ,j,k+1)- &
636 & tl_rho(i-1,j,k+1))+ &
637 & (1.0_r8-gamma)*(tl_rho(i ,j,k )- &
638 & tl_rho(i-1,j,k ))
639 tl_cff2=tl_rho(i,j,k+1)+tl_rho(i-1,j,k+1)- &
640 & tl_rho(i,j,k )-tl_rho(i-1,j,k )
641 cff3=z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
642 & z_r(i,j,k )-z_r(i-1,j,k )
643 cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i-1,j,k+1))+ &
644 & (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i-1,j,k ))
645 tl_phix(i)=tl_phix(i)+ &
646 & fac3*(tl_cff1*cff3-tl_cff2*cff4)
647 tl_gradp(i,j,k)=0.5_r8*tl_phix(i)* &
648 & (pm(i-1,j)+pm(i,j))/(f(i-1,j)+f(i,j))
649# ifdef ZETA_ELLIPTIC
650 tl_phi(i,k)=tl_phix(i)
651# endif
652 END DO
653 END DO
654
655# ifdef ZETA_ELLIPTIC
656!
657! Compute right-hand-side term used in the elliptic equation for
658! unbalanced free-surface. Integrate from bottom to surface.
659!
660 IF (balance(isfsur)) THEN
661 DO k=1,n(ng)
662 DO i=istr,iend+1
663 tl_cff=0.5_r8*(hz(i-1,j,k)+hz(i,j,k))*tl_phi(i,k)
664# ifdef MASKING
665 tl_cff=tl_cff*umask(i,j)
666# endif
667 tl_gradpx(i,j)=tl_gradpx(i,j)+tl_cff
668 END DO
669 END DO
670 END IF
671# endif
672 END DO
673 END IF
674!
675! Compute 3D V-velocity error covariance constraint.
676!
677 IF (balance(isvvel)) THEN
678 DO k=1,n(ng)
679 DO j=jstrv,jend
680 DO i=istr,iend
681 tl_v(i,j,k,linp)=tl_v(i,j,k,linp)+ &
682 & 0.25_r8*(tl_gradp(i ,j-1,k)+ &
683 & tl_gradp(i+1,j-1,k)+ &
684 & tl_gradp(i ,j ,k)+ &
685 & tl_gradp(i+1,j ,k))
686# ifdef MASKING
687 tl_v(i,j,k,linp)=tl_v(i,j,k,linp)*vmask(i,j)
688# endif
689 END DO
690 END DO
691 END DO
692
693 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
694 CALL exchange_u3d_tile (ng, tile, &
695 & lbi, ubi, lbj, ubj, 1, n(ng), &
696 & tl_u(:,:,:,linp))
697 CALL exchange_v3d_tile (ng, tile, &
698 & lbi, ubi, lbj, ubj, 1, n(ng), &
699 & tl_v(:,:,:,linp))
700 END IF
701# ifdef DISTRIBUTE
702 CALL mp_exchange3d (ng, tile, itlm, 2, &
703 & lbi, ubi, lbj, ubj, 1, n(ng), &
704 & nghostpoints, &
705 & ewperiodic(ng), nsperiodic(ng), &
706 & tl_u(:,:,:,linp), tl_v(:,:,:,linp))
707# endif
708 END IF
709!
710!-----------------------------------------------------------------------
711! Add balanced free-surface contribution to unbalanced free-surface
712! increment.
713!-----------------------------------------------------------------------
714!
715 IF (balance(isfsur)) THEN
716
717# ifdef ZETA_ELLIPTIC
718!
719! Solve elliptic equation (Fukumori et al., 1998).
720!
721 CALL u2d_bc (ng, tile, &
722 & imins, imaxs, jmins, jmaxs, &
723 & tl_gradpx)
724 CALL v2d_bc (ng, tile, &
725 & imins, imaxs, jmins, jmaxs, &
726 & tl_gradpy)
727!
728! Compute the RHS term for the biconjugate gradient solver.
729!
730 DO j=jstr,jend
731 DO i=istr,iend
732 tl_rhs_r2d(i,j)=-pm(i,j)*pn(i,j)* &
733 & (pmon_u(i+1,j)*tl_gradpx(i+1,j)- &
734 & pmon_u(i ,j)*tl_gradpx(i ,j)+ &
735 & pnom_v(i,j+1)*tl_gradpy(i,j+1)- &
736 & pnom_v(i,j )*tl_gradpy(i,j ))
737# ifdef MASKING
738 tl_rhs_r2d(i,j)=tl_rhs_r2d(i,j)*rmask(i,j)
739# endif
740 END DO
741 END DO
742
743 CALL r2d_bc (ng, tile, &
744 & lbi, ubi, lbj, ubj, &
745 & tl_rhs_r2d)
746# ifdef DISTRIBUTE
747 CALL mp_exchange2d (ng, tile, itlm, 1, &
748 & lbi, ubi, lbj, ubj, &
749 & nghostpoints, &
750 & ewperiodic(ng), nsperiodic(ng), &
751 & tl_rhs_r2d)
752# endif
753!
754! Choose the starting value of zeta_ref.
755!
756 DO j=jstr,jend
757 DO i=istr,iend
758 tl_zeta_ref(i,j)=0.0_r8
759 END DO
760 END DO
761!
762! Call the tangent linear biconjugate gradient solver.
763!
764 CALL tl_biconj_tile (ng, tile, itlm, &
765 & lbi, ubi, lbj, ubj, &
766 & imins, imaxs, jmins, jmaxs, &
767 & lbck, &
768 & h, pmon_u, pnom_v, pm, pn, &
769# ifdef MASKING
770 & umask, vmask, rmask, &
771# endif
772 & bc_ak, bc_bk, zdf1, zdf2, zdf3, &
773 & pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, &
774 & tl_zeta_ref, tl_rhs_r2d)
775!
776! Add balanced free-surface contribution to unbalanced free-surface
777! increment.
778!
779 DO j=jstr,jend
780 DO i=istr,iend
781 tl_zeta(i,j,linp)=tl_zeta(i,j,linp)+tl_zeta_ref(i,j)
782 END DO
783 END DO
784
785# else
786!
787! Integrate hydrostatic equation from bottom to surface.
788!
789 IF (lnm_flag.eq.0) THEN
790 cff1=1.0_r8/rho0
791 DO k=1,n(ng)
792 DO j=jstr,jend
793 DO i=istr,iend
794 tl_cff=-cff1*tl_rho(i,j,k)*hz(i,j,k)
795# ifdef MASKING
796 tl_cff=tl_cff*rmask(i,j)
797# endif
798 tl_zeta(i,j,linp)=tl_zeta(i,j,linp)+tl_cff
799 END DO
800 END DO
801 END DO
802!
803! Integrate from level of no motion (LNM_depth) to surface or
804! integrate from local bottom if shallower than LNM_depth.
805! Notice that the level of motion depth is tested against the
806! bottom of the grid cell (at W-points) and bracketed between
807! W-points during the interpolation. Also positive depths are
808! used for clarity.
809!
810 ELSE IF (lnm_flag.eq.1) THEN
811 cff1=1.0_r8/rho0
812 DO j=jstr,jend
813 DO i=istr,iend
814 DO k=1,n(ng)
815 zwtop=abs(z_w(i,j,k ))
816 zwbot=abs(z_w(i,j,k-1))
817 IF (zwbot.lt.lnm_depth(ng)) THEN
818 tl_cff=-cff1*tl_rho(i,j,k)*hz(i,j,k)
819# ifdef MASKING
820 tl_cff=tl_cff*rmask(i,j)
821# endif
822 tl_zeta(i,j,linp)=tl_zeta(i,j,linp)+tl_cff
823 ELSE IF ((zwbot.gt.lnm_depth(ng)).and. &
824 & (lnm_depth(ng).ge.zwtop)) THEN ! interpolate
825 zphi=abs(z_r(i,j,k))
826 IF (zphi.ge.lnm_depth(ng)) THEN ! above cell center
827 IF (k.eq.n(ng)) THEN
828 tl_cff=-cff1*tl_rho(i,j,k)*hz(i,j,k)
829 ELSE
830 zphi1=abs(z_r(i,j,k+1))
831 fac=(lnm_depth(ng)-zphi1)/(zphi-zphi1)
832 tl_cff=-cff1* &
833 & (tl_rho(i,j,k+1)+ &
834 & fac*(tl_rho(i,j,k)-tl_rho(i,j,k+1)))* &
835 & (lnm_depth(ng)-zwtop)
836 END IF
837 ELSE ! below cell center
838 IF (k.eq.1) THEN
839 tl_cff=-cff1*tl_rho(i,j,k)*hz(i,j,k)
840 ELSE
841 zphi1=abs(z_r(i,j,k-1))
842 fac=(lnm_depth(ng)-zphi)/(zphi1-zphi)
843 tl_cff=-cff1* &
844 & (tl_rho(i,j,k)+ &
845 & fac*(tl_rho(i,j,k-1)-tl_rho(i,j,k)))* &
846 & (zwbot-lnm_depth(ng))
847 END IF
848 END IF
849# ifdef MASKING
850 tl_cff=tl_cff*rmask(i,j)
851# endif
852 tl_zeta(i,j,linp)=tl_zeta(i,j,linp)+tl_cff
853 END IF
854 END DO
855 END DO
856 END DO
857 END IF
858# endif
859
860 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
861 CALL exchange_r2d_tile (ng, tile, &
862 & lbi, ubi, lbj, ubj, &
863 & tl_zeta(:,:,linp))
864 END IF
865# ifdef DISTRIBUTE
866 CALL mp_exchange2d (ng, tile, itlm, 1, &
867 & lbi, ubi, lbj, ubj, &
868 & nghostpoints, &
869 & ewperiodic(ng), nsperiodic(ng), &
870 & tl_zeta(:,:,linp))
871# endif
872 END IF
873
874 RETURN
875 END SUBROUTINE tl_balance_tile
876
877#endif
878 END MODULE tl_balance_mod
bc_2d_mod
Definition bc_2d.F:2
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_3d_mod
Definition exchange_3d.F:2
exchange_3d_mod::exchange_u3d_tile
subroutine exchange_u3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:495
exchange_3d_mod::exchange_r3d_tile
subroutine exchange_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:283
exchange_3d_mod::exchange_v3d_tile
subroutine exchange_v3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:707
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::itlm
integer, parameter itlm
Definition mod_param.F:663
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::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
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
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
tl_balance_mod
Definition tl_balance.F:2
tl_balance_mod::tl_balance_tile
subroutine tl_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, tl_zeta_ref, tl_rhs_r2d, pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, zdf1, zdf2, zdf3, bc_ak, bc_bk, tl_rho, tl_t, tl_u, tl_v, tl_zeta)
Definition tl_balance.F:178
tl_balance_mod::tl_balance
subroutine, public tl_balance(ng, tile, lbck, linp)
Definition tl_balance.F:59
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::v2d_bc
subroutine, public v2d_bc(ng, tile, lbi, ubi, lbj, ubj, a)
Definition zeta_balance.F:1253
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