#include "cppdefs.h" MODULE tl_balance_mod #if defined TANGENT && defined BALANCE_OPERATOR ! !git $Id$ !svn $Id: tl_balance.F 1180 2023-07-13 02:42:10Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2023 The ROMS/TOMS Group Andrew M. Moore ! ! Licensed under a MIT/X style license ! ! See License_ROMS.md ! !======================================================================= ! ! ! These routines impose a multivariate balance operator to constraint ! ! the background/model error covariance matrix, B, such that the ! ! unobserved variables information is extracted from observed data. ! ! It follows the approach proposed by Weaver et al. (2006). The state ! ! vector is split between balanced and unbalanced components, except ! ! for temperature, which is used to establish the balanced part of ! ! other state variables. ! ! ! ! The background/model error covariance is represented as: ! ! ! ! B = K Bu K^(T) ! ! ! ! where ! ! ! ! B : background/model error covariance matrix. ! ! Bu: unbalanced background/model error covariance matrix modeled ! ! with the generalized diffusion operator. ! ! K : balance matrix operator. ! ! ! ! Here, T denotes the transpose. ! ! ! ! The multivariate formulation is obtained by establishing balance ! ! relationships with the other state variables using T-S empirical ! ! formulas, hydrostatic balance, and geostrophic balance. ! ! ! ! Reference: ! ! ! ! Weaver, A.T., C. Deltel, E. Machu, S. Ricci, and N. Daget, 2006: ! ! A multivariate balance operator for variational data assimila- ! ! tion, Q. J. R. Meteorol. Soc., 131, 3605-3625. ! ! ! ! (See also, ECMWR Technical Memorandum # 491, April 2006) ! ! ! !======================================================================= ! USE mod_kinds ! implicit none ! PRIVATE PUBLIC :: tl_balance ! CONTAINS ! !*********************************************************************** SUBROUTINE tl_balance (ng, tile, Lbck, Linp) !*********************************************************************** ! USE mod_param USE mod_grid # ifdef SOLVE3D USE mod_coupling USE mod_mixing # endif # ifdef ZETA_ELLIPTIC USE mod_fourdvar # endif USE mod_ocean USE mod_stepping # ifdef SOLVE3D ! USE rho_eos_mod USE set_depth_mod # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, Lbck, Linp ! ! Local variable declarations. ! # include "tile.h" ! # ifdef SOLVE3D ! ! Compute background state thickness, depth arrays, thermal expansion, ! and saline contraction coefficients. ! COUPLING(ng) % Zt_avg1 = 0.0_r8 CALL set_depth (ng, tile, iTLM) nrhs(ng)=Lbck CALL rho_eos (ng, tile, iTLM) ! # endif CALL tl_balance_tile (ng, tile, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & Lbck, Linp, & & GRID(ng) % f, & & GRID(ng) % pm, & & GRID(ng) % pn, & # ifdef ZETA_ELLIPTIC & GRID(ng) % pmon_u, & & GRID(ng) % pnom_v, & & GRID(ng) % h, & # endif # ifdef SOLVE3D & GRID(ng) % Hz, & & GRID(ng) % z_r, & & GRID(ng) % z_w, & # endif # ifdef MASKING & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif # ifdef SOLVE3D & MIXING(ng) % alpha, & & MIXING(ng) % beta, & & OCEAN(ng) % t, & # endif # ifdef ZETA_ELLIPTIC & FOURDVAR(ng) % tl_zeta_ref, & & FOURDVAR(ng) % tl_rhs_r2d, & & FOURDVAR(ng) % pc_r2d, & & FOURDVAR(ng) % r_r2d, & & FOURDVAR(ng) % br_r2d, & & FOURDVAR(ng) % p_r2d, & & FOURDVAR(ng) % bp_r2d, & & FOURDVAR(ng) % zdf1, & & FOURDVAR(ng) % zdf2, & & FOURDVAR(ng) % zdf3, & & FOURDVAR(ng) % bc_ak, & & FOURDVAR(ng) % bc_bk, & # endif # ifdef SOLVE3D & OCEAN(ng) % tl_rho, & & OCEAN(ng) % tl_t, & & OCEAN(ng) % tl_u, & & OCEAN(ng) % tl_v, & # endif & OCEAN(ng) % tl_zeta) RETURN END SUBROUTINE tl_balance ! !*********************************************************************** SUBROUTINE tl_balance_tile (ng, tile, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & Lbck, Linp, & & f, pm, pn, & # ifdef ZETA_ELLIPTIC & pmon_u, pnom_v, h, & # endif # ifdef SOLVE3D & Hz, z_r, z_w, & # endif # ifdef MASKING & rmask, umask, vmask, & # endif # ifdef SOLVE3D & alpha, beta, t, & # endif # ifdef ZETA_ELLIPTIC & tl_zeta_ref, tl_rhs_r2d, & & pc_r2d, r_r2d, br_r2d, & & p_r2d, bp_r2d, & & zdf1, zdf2, zdf3, bc_ak, bc_bk, & # endif # ifdef SOLVE3D & tl_rho, tl_t, tl_u, tl_v, & # endif & tl_zeta) !*********************************************************************** ! USE mod_param USE mod_ncparam USE mod_scalars ! USE bc_2d_mod USE exchange_2d_mod USE exchange_3d_mod # ifdef DISTRIBUTE USE mp_exchange_mod, ONLY : mp_exchange2d, mp_exchange3d # endif # ifdef ZETA_ELLIPTIC USE zeta_balance_mod, ONLY: r2d_bc, u2d_bc, v2d_bc USE zeta_balance_mod, ONLY: tl_biconj_tile # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Lbck, Linp ! # ifdef ASSUMED_SHAPE real(r8), intent(in) :: f(LBi:,LBj:) real(r8), intent(in) :: pm(LBi:,LBj:) real(r8), intent(in) :: pn(LBi:,LBj:) # ifdef ZETA_ELLIPTIC real(r8), intent(in) :: h(LBi:,LBj:) real(r8), intent(in) :: pmon_u(LBi:,LBj:) real(r8), intent(in) :: pnom_v(LBi:,LBj:) # endif # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(in) :: z_r(LBi:,LBj:,:) real(r8), intent(in) :: z_w(LBi:,LBj:,0:) # endif # ifdef MASKING real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif # ifdef SOLVE3D real(r8), intent(in) :: alpha(LBi:,LBj:) real(r8), intent(in) :: beta(LBi:,LBj:) real(r8), intent(in) :: t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:) real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:) # endif real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:) # ifdef SOLVE3D real(r8), intent(out) :: tl_rho(LBi:,LBj:,:) # endif # ifdef ZETA_ELLIPTIC real(r8), intent(in) :: bc_ak(:) real(r8), intent(in) :: bc_bk(:) real(r8), intent(in) :: zdf1(:) real(r8), intent(in) :: zdf2(:) real(r8), intent(in) :: zdf3(:) real(r8), intent(inout) :: pc_r2d(LBi:,LBj:) real(r8), intent(inout) :: r_r2d(LBi:,LBj:,:) real(r8), intent(inout) :: br_r2d(LBi:,LBj:,:) real(r8), intent(inout) :: p_r2d(LBi:,LBj:,:) real(r8), intent(inout) :: bp_r2d(LBi:,LBj:,:) real(r8), intent(inout) :: tl_rhs_r2d(LBi:,LBj:) real(r8), intent(inout) :: tl_zeta_ref(LBi:,LBj:) # endif # else real(r8), intent(in) :: f(LBi:UBi,LBj:UBj) real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj) real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj) # ifdef ZETA_ELLIPTIC real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj) real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj) # endif # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng)) # endif # ifdef MASKING real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif # ifdef SOLVE3D real(r8), intent(in) :: alpha(LBi:UBi,LBj:UBj) real(r8), intent(in) :: beta(LBi:UBi,LBj:UBj) real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,2,N(ng)) real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,2,N(ng)) # endif real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:) # ifdef SOLVE3D real(r8), intent(out) :: tl_rho(LBi:UBi,LBj:UBj,N(ng)) # endif # ifdef ZETA_ELLIPTIC real(r8), intent(in) :: bc_ak(Nbico(ng)) real(r8), intent(in) :: bc_bk(Nbico(ng)) real(r8), intent(in) :: zdf1(Nbico(ng)) real(r8), intent(in) :: zdf2(Nbico(ng)) real(r8), intent(in) :: zdf3(Nbico(ng)) real(r8), intent(inout) :: pc_r2d(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: r_r2d(LBi:UBi,LBj:UBj,Nbico(ng)) real(r8), intent(inout) :: br_r2d(LBi:UBi,LBj:UBj,Nbico(ng)) real(r8), intent(inout) :: p_r2d(LBi:UBi,LBj:UBj,Nbico(ng)) real(r8), intent(inout) :: bp_r2d(LBi:UBi,LBj:UBj,Nbico(ng)) real(r8), intent(inout) :: tl_rhs_r2d(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: tl_zeta_ref(LBi:UBi,LBj:UBj) # endif # endif ! ! Local variable declarations. ! integer :: i, j, k, order integer :: Norder = 2 ! Shapiro filter order real(r8) :: fac, fac1, fac2, fac3, gamma real(r8) :: cff, cff1, cff2, cff3, cff4 real(r8) :: tl_cff, tl_cff1, tl_cff2 real(r8) :: dzdT, zphi, zphi1, zwbot, zwtop real(r8), dimension(20) :: filter_coef = & & (/ 2.500000E-1_r8, 6.250000E-2_r8, 1.562500E-2_r8, & & 3.906250E-3_r8, 9.765625E-4_r8, 2.44140625E-4_r8, & & 6.103515625E-5_r8, 1.5258789063E-5_r8, 3.814697E-6_r8, & & 9.536743E-7_r8, 2.384186E-7_r8, 5.960464E-8_r8, & & 1.490116E-8_r8, 3.725290E-9_r8, 9.313226E-10_r8, & & 2.328306E-10_r8, 5.820766E-11_r8, 1.455192E-11_r8, & & 3.637979E-12_r8, 9.094947E-13_r8 /) real(r8), dimension(N(ng)) :: dSdT, dSdT_filter real(r8), dimension(IminS:ImaxS) :: tl_phie real(r8), dimension(IminS:ImaxS) :: tl_phix # ifdef SALINITY real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: dTdz real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: dSdz # endif real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: tl_gradP # ifdef ZETA_ELLIPTIC real(r8), dimension(IminS:ImaxS,N(ng)) :: tl_phi real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_gradPx real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_gradPy real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_r2d_rhs # endif # include "set_bounds.h" # ifdef SALINITY ! !----------------------------------------------------------------------- ! Compute balance salinity contribution. !----------------------------------------------------------------------- ! IF (balance(isTvar(isalt))) THEN ! ! Compute temperature (dTdz) and salinity (dSdz) shears. ! DO j=JstrR,JendR DO i=IstrR,IendR FC(i,0)=0.0_r8 dTdz(i,j,0)=0.0_r8 dSdz(i,j,0)=0.0_r8 END DO DO k=1,N(ng)-1 DO i=IstrR,IendR cff=1.0_r8/(2.0_r8*Hz(i,j,k+1)+ & & Hz(i,j,k)*(2.0_r8-FC(i,k-1))) FC(i,k)=cff*Hz(i,j,k+1) dTdz(i,j,k)=cff*(6.0_r8*(t(i,j,k+1,Lbck,itemp)- & & t(i,j,k ,Lbck,itemp))- & & Hz(i,j,k)*dTdz(i,j,k-1)) dSdz(i,j,k)=cff*(6.0_r8*(t(i,j,k+1,Lbck,isalt)- & & t(i,j,k ,Lbck,isalt))- & & Hz(i,j,k)*dSdz(i,j,k-1)) END DO END DO DO i=IstrR,IendR dTdz(i,j,N(ng))=0.0_r8 dSdz(i,j,N(ng))=0.0_r8 END DO DO k=N(ng)-1,1,-1 DO i=IstrR,IendR dTdz(i,j,k)=dTdz(i,j,k)-FC(i,k)*dTdz(i,j,k+1) dSdz(i,j,k)=dSdz(i,j,k)-FC(i,k)*dSdz(i,j,k+1) END DO END DO END DO ! ! Add balanced salinity (deltaS_b) contribution to unbalanced salinity ! increment. The unbalanced salinity increment is related related to ! temperature increment: ! ! deltaS_b = cff * dS/dz * dz/dT * deltaT ! ! Here, cff is a coefficient that depends on the mixed layer depth: ! ! cff = 1.0 - EXP (z_r / ml_depth) ! ! the coefficient is smoothly reduced to zero at the surface and below ! the mixed layer. ! DO j=JstrR,JendR DO i=IstrR,IendR DO k=1,N(ng) cff=0.5_r8*(dTdz(i,j,k-1)+dTdz(i,j,k)) IF (ABS(cff).lt.dTdz_min(ng)) THEN dzdT=0.0_r8 ELSE dzdT=1.0_r8/cff END IF dSdT(k)=(0.5_r8*(dSdz(i,j,k-1)+ & & dSdz(i,j,k )))*dzdT END DO ! ! Shapiro filter. ! DO order=1,Norder/2 IF (order.ne.Norder/2) THEN dSdT_filter(1)=2.0_r8*(dSdT(1)-dSdT(2)) dSdT_filter(N(ng))=2.0_r8*(dSdT(N(ng))-dSdT(N(ng)-1)) ELSE dSdT_filter(1)=0.0_r8 dSdT_filter(N(ng))=0.0_r8 END IF DO k=2,N(ng)-1 dSdT_filter(k)=2.0_r8*dSdT(k)-dSdT(k-1)-dSdT(k+1) END DO DO k=1,N(ng) dSdT(k)=dSdT(k)-filter_coef(Norder/2)*dSdT_filter(k) END DO END DO DO k=1,N(ng) cff=(1.0_r8-EXP(z_r(i,j,k)/ml_depth(ng)))*dSdT(k) tl_t(i,j,k,Linp,isalt)=tl_t(i,j,k,Linp,isalt)+ & & cff*tl_t(i,j,k,Linp,itemp) # ifdef MASKING tl_t(i,j,k,Linp,isalt)=tl_t(i,j,k,Linp,isalt)*rmask(i,j) # endif END DO END DO END DO IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & tl_t(:,:,:,Linp,isalt)) END IF # ifdef DISTRIBUTE CALL mp_exchange3d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_t(:,:,:,Linp,isalt)) # endif END IF # endif ! !----------------------------------------------------------------------- ! Compute balanced density anomaly increment using linearized equation ! of state. The thermal expansion and saline contraction coefficients ! are computed from the background state. !----------------------------------------------------------------------- ! IF (balance(isFsur).or.balance(isVvel)) THEN DO j=JstrR,JendR DO k=1,N(ng) DO i=IstrR,IendR tl_rho(i,j,k)=-rho0*alpha(i,j)*tl_t(i,j,k,Linp,itemp) # ifdef SALINITY tl_rho(i,j,k)=tl_rho(i,j,k)+ & & rho0*beta(i,j)*tl_t(i,j,k,Linp,isalt) # endif # ifdef MASKING tl_rho(i,j,k)=tl_rho(i,j,k)*rmask(i,j) # endif END DO END DO END DO IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & tl_rho) END IF # ifdef DISTRIBUTE CALL mp_exchange3d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_rho) # endif END IF ! !----------------------------------------------------------------------- ! Add balanced velocity contributions to unbalanced velocity ! increments. Use linear pressure gradient formulation based ! to that found in routine "prsgrd31.h". !----------------------------------------------------------------------- ! ! NOTE: fac2=0 because the balanced component should consist of the ! baroclinic pressure gradient only. ! fac1=0.5_r8*g/rho0 !! fac2=g fac2=0.0_r8 fac3=0.25_r8*g/rho0 # ifdef ZETA_ELLIPTIC ! ! Initialize vertical intergal local variables. ! DO j=JminS,JmaxS DO i=IminS,ImaxS tl_gradPx(i,j)=0.0_r8 tl_gradPy(i,j)=0.0_r8 END DO END DO # endif ! ! Compute balanced, surface U-momentum from baroclinic and barotropic ! surface pressure gradient. ! IF (balance(isFsur).or.balance(isUvel)) THEN DO j=Jstr,Jend+1 DO i=Istr-1,Iend cff1=z_w(i,j ,N(ng))-z_r(i,j ,N(ng))+ & & z_w(i,j-1,N(ng))-z_r(i,j-1,N(ng)) tl_phie(i)=fac1*(tl_rho(i,j ,N(ng))- & & tl_rho(i,j-1,N(ng)))*cff1+ & & fac2*(tl_zeta(i,j ,Linp)- & & tl_zeta(i,j-1,Linp)) tl_gradP(i,j,N(ng))=0.5_r8*tl_phie(i)* & & (pn(i,j-1)+pn(i,j))/(f(i,j-1)+f(i,j)) # ifdef ZETA_ELLIPTIC tl_phi(i,N(ng))=tl_phie(i) # endif END DO ! ! Compute balance, interior U-momentum from baroclinic pressure ! gradient (differentiate and then vertically integrate). ! DO k=N(ng)-1,1,-1 DO i=Istr-1,Iend cff1=1.0_r8/((z_r(i,j ,k+1)-z_r(i,j ,k))* & & (z_r(i,j-1,k+1)-z_r(i,j-1,k))) cff2=z_r(i,j ,k )-z_r(i,j-1,k )+ & & z_r(i,j ,k+1)-z_r(i,j-1,k+1) cff3=z_r(i,j ,k+1)-z_r(i,j ,k )- & & z_r(i,j-1,k+1)+z_r(i,j-1,k ) gamma=0.125_r8*cff1*cff2*cff3 ! tl_cff1=(1.0_r8+gamma)*(tl_rho(i,j ,k+1)- & & tl_rho(i,j-1,k+1))+ & & (1.0_r8-gamma)*(tl_rho(i,j ,k )- & & tl_rho(i,j-1,k )) tl_cff2=tl_rho(i,j,k+1)+tl_rho(i,j-1,k+1)- & & tl_rho(i,j,k )-tl_rho(i,j-1,k ) cff3=z_r(i,j,k+1)+z_r(i,j-1,k+1)- & & z_r(i,j,k )-z_r(i,j-1,k ) cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i,j-1,k+1))+ & & (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i,j-1,k )) tl_phie(i)=tl_phie(i)+ & & fac3*(tl_cff1*cff3-tl_cff2*cff4) tl_gradP(i,j,k)=0.5_r8*tl_phie(i)* & & (pn(i,j-1)+pn(i,j))/(f(i,j-1)+f(i,j)) # ifdef ZETA_ELLIPTIC tl_phi(i,k)=tl_phie(i) # endif END DO END DO # ifdef ZETA_ELLIPTIC ! ! Compute right-hand-side term used in the elliptic equation for ! unbalanced free-surface. Integrate from bottom to surface. ! IF (balance(isFsur)) THEN DO k=1,N(ng) DO i=Istr-1,Iend tl_cff=0.5_r8*(Hz(i,j-1,k)+Hz(i,j,k))*tl_phi(i,k) # ifdef MASKING tl_cff=tl_cff*vmask(i,j) # endif tl_gradPy(i,j)=tl_gradPy(i,j)+tl_cff END DO END DO END IF # endif END DO END IF ! ! Compute 3D U-velocity error covariance constraint. ! IF (balance(isUvel)) THEN DO k=1,N(ng) DO j=Jstr,Jend DO i=IstrU,Iend tl_u(i,j,k,Linp)=tl_u(i,j,k,Linp)- & & 0.25_r8*(tl_gradP(i-1,j ,k)+ & & tl_gradP(i ,j ,k)+ & & tl_gradP(i-1,j+1,k)+ & & tl_gradP(i ,j+1,k)) # ifdef MASKING tl_u(i,j,k,Linp)=tl_u(i,j,k,Linp)*umask(i,j) # endif END DO END DO END DO END IF ! ! Compute balanced, surface V-momentum from baroclinic and barotropic ! surface pressure gradient. ! IF (balance(isFsur).or.balance(isVvel)) THEN DO j=Jstr-1,Jend DO i=Istr,Iend+1 cff1=z_w(i ,j,N(ng))-z_r(i ,j,N(ng))+ & & z_w(i-1,j,N(ng))-z_r(i-1,j,N(ng)) tl_phix(i)=fac1*(tl_rho(i ,j,N(ng))- & & tl_rho(i-1,j,N(ng)))*cff1+ & & fac2*(tl_zeta(i ,j,Linp)- & & tl_zeta(i-1,j,Linp)) tl_gradP(i,j,N(ng))=0.5_r8*tl_phix(i)* & & (pm(i-1,j)+pm(i,j))/(f(i-1,j)+f(i,j)) # ifdef ZETA_ELLIPTIC tl_phi(i,N(ng))=tl_phix(i) # endif END DO ! ! Compute balance, interior V-momentum from baroclinic pressure ! gradient (differentiate and then vertically integrate). ! DO k=N(ng)-1,1,-1 DO i=Istr,Iend+1 cff1=1.0_r8/((z_r(i ,j,k+1)-z_r(i ,j,k))* & & (z_r(i-1,j,k+1)-z_r(i-1,j,k))) cff2=z_r(i ,j,k )-z_r(i-1,j,k )+ & & z_r(i ,j,k+1)-z_r(i-1,j,k+1) cff3=z_r(i ,j,k+1)-z_r(i ,j,k )- & & z_r(i-1,j,k+1)+z_r(i-1,j,k ) gamma=0.125_r8*cff1*cff2*cff3 ! tl_cff1=(1.0_r8+gamma)*(tl_rho(i ,j,k+1)- & & tl_rho(i-1,j,k+1))+ & & (1.0_r8-gamma)*(tl_rho(i ,j,k )- & & tl_rho(i-1,j,k )) tl_cff2=tl_rho(i,j,k+1)+tl_rho(i-1,j,k+1)- & & tl_rho(i,j,k )-tl_rho(i-1,j,k ) cff3=z_r(i,j,k+1)+z_r(i-1,j,k+1)- & & z_r(i,j,k )-z_r(i-1,j,k ) cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i-1,j,k+1))+ & & (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i-1,j,k )) tl_phix(i)=tl_phix(i)+ & & fac3*(tl_cff1*cff3-tl_cff2*cff4) tl_gradP(i,j,k)=0.5_r8*tl_phix(i)* & & (pm(i-1,j)+pm(i,j))/(f(i-1,j)+f(i,j)) # ifdef ZETA_ELLIPTIC tl_phi(i,k)=tl_phix(i) # endif END DO END DO # ifdef ZETA_ELLIPTIC ! ! Compute right-hand-side term used in the elliptic equation for ! unbalanced free-surface. Integrate from bottom to surface. ! IF (balance(isFsur)) THEN DO k=1,N(ng) DO i=Istr,Iend+1 tl_cff=0.5_r8*(Hz(i-1,j,k)+Hz(i,j,k))*tl_phi(i,k) # ifdef MASKING tl_cff=tl_cff*umask(i,j) # endif tl_gradPx(i,j)=tl_gradPx(i,j)+tl_cff END DO END DO END IF # endif END DO END IF ! ! Compute 3D V-velocity error covariance constraint. ! IF (balance(isVvel)) THEN DO k=1,N(ng) DO j=JstrV,Jend DO i=Istr,Iend tl_v(i,j,k,Linp)=tl_v(i,j,k,Linp)+ & & 0.25_r8*(tl_gradP(i ,j-1,k)+ & & tl_gradP(i+1,j-1,k)+ & & tl_gradP(i ,j ,k)+ & & tl_gradP(i+1,j ,k)) # ifdef MASKING tl_v(i,j,k,Linp)=tl_v(i,j,k,Linp)*vmask(i,j) # endif END DO END DO END DO IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_u3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & tl_u(:,:,:,Linp)) CALL exchange_v3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & tl_v(:,:,:,Linp)) END IF # ifdef DISTRIBUTE CALL mp_exchange3d (ng, tile, iTLM, 2, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_u(:,:,:,Linp), tl_v(:,:,:,Linp)) # endif END IF ! !----------------------------------------------------------------------- ! Add balanced free-surface contribution to unbalanced free-surface ! increment. !----------------------------------------------------------------------- ! IF (balance(isFsur)) THEN # ifdef ZETA_ELLIPTIC ! ! Solve elliptic equation (Fukumori et al., 1998). ! CALL u2d_bc (ng, tile, & & IminS, ImaxS, JminS, JmaxS, & & tl_gradPx) CALL v2d_bc (ng, tile, & & IminS, ImaxS, JminS, JmaxS, & & tl_gradPy) ! ! Compute the RHS term for the biconjugate gradient solver. ! DO j=Jstr,Jend DO i=Istr,Iend tl_rhs_r2d(i,j)=-pm(i,j)*pn(i,j)* & & (pmon_u(i+1,j)*tl_gradPx(i+1,j)- & & pmon_u(i ,j)*tl_gradPx(i ,j)+ & & pnom_v(i,j+1)*tl_gradPy(i,j+1)- & & pnom_v(i,j )*tl_gradPy(i,j )) # ifdef MASKING tl_rhs_r2d(i,j)=tl_rhs_r2d(i,j)*rmask(i,j) # endif END DO END DO CALL r2d_bc (ng, tile, & & LBi, UBi, LBj, UBj, & & tl_rhs_r2d) # ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_rhs_r2d) # endif ! ! Choose the starting value of zeta_ref. ! DO j=Jstr,Jend DO i=Istr,Iend tl_zeta_ref(i,j)=0.0_r8 END DO END DO ! ! Call the tangent linear biconjugate gradient solver. ! CALL tl_biconj_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & Lbck, & & h, pmon_u, pnom_v, pm, pn, & # ifdef MASKING & umask, vmask, rmask, & # endif & bc_ak, bc_bk, zdf1, zdf2, zdf3, & & pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, & & tl_zeta_ref, tl_rhs_r2d) ! ! Add balanced free-surface contribution to unbalanced free-surface ! increment. ! DO j=Jstr,Jend DO i=Istr,Iend tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)+tl_zeta_ref(i,j) END DO END DO # else ! ! Integrate hydrostatic equation from bottom to surface. ! IF (LNM_flag.eq.0) THEN cff1=1.0_r8/rho0 DO k=1,N(ng) DO j=Jstr,Jend DO i=Istr,Iend tl_cff=-cff1*tl_rho(i,j,k)*Hz(i,j,k) # ifdef MASKING tl_cff=tl_cff*rmask(i,j) # endif tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)+tl_cff END DO END DO END DO ! ! Integrate from level of no motion (LNM_depth) to surface or ! integrate from local bottom if shallower than LNM_depth. ! Notice that the level of motion depth is tested against the ! bottom of the grid cell (at W-points) and bracketed between ! W-points during the interpolation. Also positive depths are ! used for clarity. ! ELSE IF (LNM_flag.eq.1) THEN cff1=1.0_r8/rho0 DO j=Jstr,Jend DO i=Istr,Iend DO k=1,N(ng) zwtop=ABS(z_w(i,j,k )) zwbot=ABS(z_w(i,j,k-1)) IF (zwbot.lt.LNM_depth(ng)) THEN tl_cff=-cff1*tl_rho(i,j,k)*Hz(i,j,k) # ifdef MASKING tl_cff=tl_cff*rmask(i,j) # endif tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)+tl_cff ELSE IF ((Zwbot.gt.LNM_depth(ng)).and. & & (LNM_depth(ng).ge.Zwtop)) THEN ! interpolate zphi=ABS(z_r(i,j,k)) IF (zphi.ge.LNM_depth(ng)) THEN ! above cell center IF (k.eq.N(ng)) THEN tl_cff=-cff1*tl_rho(i,j,k)*Hz(i,j,k) ELSE zphi1=ABS(z_r(i,j,k+1)) fac=(LNM_depth(ng)-zphi1)/(zphi-zphi1) tl_cff=-cff1* & & (tl_rho(i,j,k+1)+ & & fac*(tl_rho(i,j,k)-tl_rho(i,j,k+1)))* & & (LNM_depth(ng)-zwtop) END IF ELSE ! below cell center IF (k.eq.1) THEN tl_cff=-cff1*tl_rho(i,j,k)*Hz(i,j,k) ELSE zphi1=ABS(z_r(i,j,k-1)) fac=(LNM_depth(ng)-zphi)/(zphi1-zphi) tl_cff=-cff1* & & (tl_rho(i,j,k)+ & & fac*(tl_rho(i,j,k-1)-tl_rho(i,j,k)))* & & (zwbot-LNM_depth(ng)) END IF END IF # ifdef MASKING tl_cff=tl_cff*rmask(i,j) # endif tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)+tl_cff END IF END DO END DO END DO END IF # endif IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & tl_zeta(:,:,Linp)) END IF # ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_zeta(:,:,Linp)) # endif END IF RETURN END SUBROUTINE tl_balance_tile #endif END MODULE tl_balance_mod