From c1d5fc5b7e6b6e568d11a7e77f14b72483fbec1f Mon Sep 17 00:00:00 2001 From: Nora Loose Date: Tue, 9 Apr 2024 17:22:15 -0400 Subject: [PATCH] Fix bug and add option to include memory in eddy forcing --- .../lateral/MOM_Zanna_Bolton.F90 | 304 +++++++++++++++--- 1 file changed, 267 insertions(+), 37 deletions(-) diff --git a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 index 60a587fe64..b2d0c359bf 100644 --- a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 +++ b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 @@ -55,9 +55,10 @@ module MOM_Zanna_Bolton !! 2 - Eulerian memory of eddy forcing integer :: ZB_memory_location !< Select where memory is incorporated within ZB model: !! 0 - no memory - !! 1 - memory for eddy forcing + !! 1 - memory for eddy fluxes !! 2 - memory for large scale flow - !! 3 - memory for both eddy forcing and large scale flow + !! 3 - memory for both eddy fluxes and large scale flow + !! 4 - memory for eddy forcing integer :: ZB_memory_dynamical_scale !< Select which dynamical length scale to infer memory time !! scale from: !! 0 - no memory @@ -98,7 +99,9 @@ module MOM_Zanna_Bolton sh_xx_with_memory, & !< Memory-enhanced horizontal tension incl metric terms [T-1 ~> s-1] sh_xy_with_memory, & !< Memory-enhanced shearing strain incl metric terms [T-1 ~> s-1] vort_xy_with_memory, & !< Memory-enhanced vorticity incl metric terms [T-1 ~> s-1] - ZBtau_large_scale !< Memory timescale for velocity gradients [T ~> s] + ZBtau_large_scale, & !< Memory timescale for velocity gradients [T ~> s] + ZB2020u_with_memory, & !< Memory-enhanced ZB subgrid forcing in zonal direction [L T-2 ~> m s-2] + ZB2020v_with_memory !< Memory-enhanced ZB subgrid forcing in meridional direction [L T-2 ~> m s-2] real, dimension(:,:), allocatable :: & kappa_h, & !< Scaling coefficient in h points [L2 ~> m2] @@ -255,9 +258,10 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) call get_param(param_file, mdl, "ZB_MEMORY_LOCATION", CS%ZB_memory_location, & "Select where memory is incorporated within ZB model:\n" //& "\t 0 - no memory\n" //& - "\t 1 - memory for eddy forcing\n" //& + "\t 1 - memory for eddy fluxes\n" //& "\t 2 - memory for large scale flow\n" //& - "\t 3 - memory for both eddy forcing and large scale flow", default=0) + "\t 3 - memory for both eddy fluxes and large scale flow\n" //& + "\t 4 - memory for eddy forcing", default=0) call get_param(param_file, mdl, "ZB_MEMORY_DYNAMICAL_SCALE", CS%ZB_memory_dynamical_scale, & "Select which dynamical length scale to infer memory time scale from:\n" //& "\t 0 - no memory\n" //& @@ -428,6 +432,11 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) allocate(CS%vort_xy_with_memory(SZIB_(G),SZJB_(G),SZK_(GV))); CS%vort_xy_with_memory(:,:,:) = 0. allocate(CS%ZBtau_large_scale(SZI_(G),SZJ_(G),SZK_(GV))); CS%ZBtau_large_scale(:,:,:) = 0. endif + if (CS%ZB_memory_location == 4) then + allocate(CS%ZB2020u_with_memory(SZIB_(G),SZJ_(G),SZK_(GV))); CS%ZB2020u_with_memory(:,:,:) = 0. + allocate(CS%ZB2020v_with_memory(SZI_(G),SZJB_(G),SZK_(GV))); CS%ZB2020v_with_memory(:,:,:) = 0. + allocate(CS%ZBtau_eddy(SZI_(G),SZJ_(G),SZK_(GV))); CS%ZBtau_eddy(:,:,:) = 0. + endif endif @@ -528,6 +537,11 @@ subroutine ZB_2020_end(CS) deallocate(CS%vort_xy_with_memory) deallocate(CS%ZBtau_large_scale) endif + if (CS%ZB_memory_location == 4) then + deallocate(CS%ZB2020u_with_memory) + deallocate(CS%ZB2020v_with_memory) + deallocate(CS%ZBtau_eddy) + endif endif if (CS%use_ann==1) then @@ -676,6 +690,9 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & dx2h, dy2h, dx2q, dy2q, & G, GV, CS) + if (CS%ZB_memory_location == 4) then + if (CS%ZB_memory_type > 0) call step_forward_ZB_memory_forcing(ZB2020u, ZB2020v, u, v, h, G, GV, CS) + endif ! Update the eddy viscosity acceleration with ZB model call cpu_clock_begin(CS%id_clock_upd) @@ -730,7 +747,7 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) intent(inout) :: v !< The meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. - logical, intent(in) :: eddy !< If true, step forward memory on eddy forcing; + logical, intent(in) :: eddy !< If true, step forward memory on eddies fluxes; ! if false, step forward memory on velocity gradients ! Local variables @@ -775,16 +792,16 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) call pass_var(CS%Txx_with_memory, G%Domain) call pass_var(CS%Tyy_with_memory, G%Domain) call pass_var(CS%Txy_with_memory, G%Domain, position=CORNER) - call advection(u, v, G, GV, CS, h=CS%Txx_with_memory) - call advection(u, v, G, GV, CS, h=CS%Tyy_with_memory) - call advection(u, v, G, GV, CS, q=CS%Txy_with_memory) + call advection(u, v, G, GV, CS, hfld=CS%Txx_with_memory) + call advection(u, v, G, GV, CS, hfld=CS%Tyy_with_memory) + call advection(u, v, G, GV, CS, qfld=CS%Txy_with_memory) else call pass_var(CS%sh_xx_with_memory, G%Domain) call pass_var(CS%sh_xy_with_memory, G%Domain, position=CORNER) call pass_var(CS%vort_xy_with_memory, G%Domain, position=CORNER) - call advection(u, v, G, GV, CS, h=CS%sh_xx_with_memory) - call advection(u, v, G, GV, CS, q=CS%sh_xy_with_memory) - call advection(u, v, G, GV, CS, q=CS%vort_xy_with_memory) + call advection(u, v, G, GV, CS, hfld=CS%sh_xx_with_memory) + call advection(u, v, G, GV, CS, qfld=CS%sh_xy_with_memory) + call advection(u, v, G, GV, CS, qfld=CS%vort_xy_with_memory) endif endif @@ -823,7 +840,6 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) tau = max(tau, CS%ZB_min_memory) CS%ZBtau_eddy(i,j,k) = tau - ! Eulerian memory with implicit time stepping CS%Txx_with_memory(i,j,k) = tau / (tau + CS%dt) * CS%Txx_with_memory(i,J,k) & + CS%dt / (tau + CS%dt) * CS%Txx(i,j,k) CS%Tyy_with_memory(i,j,k) = tau / (tau + CS%dt) * CS%Tyy_with_memory(i,J,k) & @@ -831,7 +847,6 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) CS%Txx(i,j,k) = CS%Txx_with_memory(i,j,k) CS%Tyy(i,j,k) = CS%Tyy_with_memory(i,j,k) - else if (CS%ZB_large_scale_memory_const > 0.0) then tau = CS%ZB_large_scale_memory_const * tau0 @@ -850,23 +865,16 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) enddo ; enddo - do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 - ! compute eddy memory time scale at cell corner - sh_xx_q = 0.25 * ( (CS%sh_xx(i+1,j+1,k) + CS%sh_xx(i,j,k)) & - + (CS%sh_xx(i+1,j,k) + CS%sh_xx(i,j+1,k))) - D = sqrt(CS%sh_xy(I,J,k) * CS%sh_xy(I,J,k) + sh_xx_q * sh_xx_q) + do J=Jsq,Jeq ; do I=Isq,Ieq if (eddy) then - tau = CS%ZB_eddy_memory_const / (D + 1.0/CS%ZB_max_memory) - !tau = 7776000.0 ! 90 days + tau = 0.25 * (CS%ZBtau_eddy(i,j,k) + CS%ZBtau_eddy(i+1,j,k) + CS%ZBtau_eddy(i,j+1,k) + CS%ZBtau_eddy(i+1,j+1,k)) CS%Txy_with_memory(I,J,k) = tau / (tau + CS%dt) * CS%Txy_with_memory(i,J,k) & + CS%dt / (tau + CS%dt) * CS%Txy(I,J,k) CS%Txy(I,J,k) = CS%Txy_with_memory(I,J,k) - else - tau = CS%ZB_large_scale_memory_const / (D + 1.0/CS%ZB_max_memory) - !tau = 7776000.0 ! 90 days + tau = 0.25 * (CS%ZBtau_large_scale(i,j,k) + CS%ZBtau_large_scale(i+1,j,k) + CS%ZBtau_large_scale(i,j+1,k) + CS%ZBtau_large_scale(i+1,j+1,k)) CS%sh_xy_with_memory(I,J,k) = tau / (tau + CS%dt) * CS%sh_xy_with_memory(i,J,k) & + CS%dt / (tau + CS%dt) * CS%sh_xy(I,J,k) @@ -894,8 +902,110 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) end subroutine step_forward_ZB_memory +subroutine step_forward_ZB_memory_forcing(ZB2020u, ZB2020v, u, v, h, G, GV, CS) + type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: ZB2020u !< Zonal acceleration due to convergence of along- + !! coordinate stress tensor in fixed layer k [L T-2 ~> m s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(inout) :: ZB2020v !< Meridional acceleration due to convergence of along- + !! coordinate stress tensor in fixed layer k [L T-2 ~> m s-2] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: u !< The zonal velocity [L T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(inout) :: v !< The meridional velocity [L T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + + ! Local variables -subroutine advection(u, v, G, GV, CS, h, q) + real :: sh_xy_h !< Shearing strain interpolated to h point [T-1 ~> s-1] + real :: vort_xy_h !< Relative vorticity interpolated to h point [T-1 ~> s-1] + real :: sh_xx_q !< Horizontal tension interpolated to q point [T-1 ~> s-1] + real :: D !< deformation rate [T-1 ~> s-1] + real :: tau, tau0 !< eddy memory time scale tau [T ~> s] + real :: Rd !< first Rossby radius of deformation in baroclinic DG [L ~> m] + integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq + integer :: i, j, k + character(len=160) :: mesg ! The text of an error message + real :: subroundoff_Cor ! A negligible parameter which avoids division by zero + ! but small compared to Coriolis parameter [T-1 ~> s-1] + real :: ICoriolis ! inverse of Coriolis parameter [T ~> s] + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + + subroundoff_Cor = 1e-30 + + call pass_vector(ZB2020u, ZB2020v, G%Domain) + + if (CS%ZB_memory_type == 1) then + call pass_vector(CS%ZB2020u_with_memory, CS%ZB2020v_with_memory, G%Domain) + call advection(u, v, G, GV, CS, ufld=CS%ZB2020u_with_memory) + call advection(u, v, G, GV, CS, vfld=CS%ZB2020v_with_memory) + endif + + do k=1,nz + do j=js-1,je+1 ; do i=is-1,ie+1 + + ! compute dynamical time scale tau0 + if (CS%ZB_memory_dynamical_scale > 0) then + ! compute eddy memory time scale at cell center + sh_xy_h = 0.25 * ( (CS%sh_xy(I-1,J-1,k) + CS%sh_xy(I,J,k)) & + + (CS%sh_xy(I-1,J,k) + CS%sh_xy(I,J-1,k)) ) + vort_xy_h = 0.25 * ( (CS%vort_xy(I-1,J-1,k) + CS%vort_xy(I,J,k)) & + + (CS%vort_xy(I-1,J,k) + CS%vort_xy(I,J-1,k)) ) + if (CS%ZB_memory_dynamical_scale == 1) then + D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h) + tau0 = 1 / (D + 1.0/CS%ZB_max_memory) + elseif (CS%ZB_memory_dynamical_scale == 2) then + D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h + vort_xy_h * vort_xy_h) + tau0 = 1 / (D + 1.0/CS%ZB_max_memory) + elseif (CS%ZB_memory_dynamical_scale == 3) then + ICoriolis = 1. / (0.25 * ((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) & + + (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) + subroundoff_Cor) + Rd = sqrt(GV%g_prime(2) * h(i,j,1) * h(i,j,2) / (h(i,j,1)+h(i,j,2))) * ICoriolis + D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h + vort_xy_h * vort_xy_h) + tau0 = Rd / sqrt(G%areaT(i,j)) / (D + Rd / sqrt(G%areaT(i,j)) /CS%ZB_max_memory) + endif + + endif + + if (CS%ZB_eddy_memory_const > 0.0) then + tau = CS%ZB_eddy_memory_const * tau0 + else + tau = CS%ZB_eddy_memory + endif + tau = max(tau, CS%ZB_min_memory) + CS%ZBtau_eddy(i,j,k) = tau + enddo; enddo + + + ! introduce memory for zonal ZB subgrid forcing + do j=js,je ; do I=Isq,Ieq + tau = 0.5 * (CS%ZBtau_eddy(i,j,k) + CS%ZBtau_eddy(i+1,j,k)) + + ! apply Eulerian memory with implicit time stepping to ZB2020u + CS%ZB2020u_with_memory(I,j,k) = tau / (tau + CS%dt) * CS%ZB2020u_with_memory(I,j,k) & + + CS%dt / (tau + CS%dt) * G%mask2dCu(I,j) * ZB2020u(I,j,k) + ZB2020u(I,j,k) = CS%ZB2020u_with_memory(I,j,k) + enddo ; enddo + + ! introduce memory for meridional ZB subgrid forcing + do J=Jsq,Jeq ; do i=is,ie + tau = 0.5 * (CS%ZBtau_eddy(i,j,k) + CS%ZBtau_eddy(i,j+1,k)) + + ! apply Eulerian memory with implicit time stepping to h*ZB2020v + CS%ZB2020v_with_memory(i,J,k) = tau / (tau + CS%dt) * CS%ZB2020v_with_memory(i,J,k) & + + CS%dt / (tau + CS%dt) * G%mask2dCv(i,J) * ZB2020v(i,J,k) + ZB2020v(i,J,k) = CS%ZB2020v_with_memory(i,J,k) + enddo ; enddo + enddo ! end of k loop + +end subroutine step_forward_ZB_memory_forcing + +subroutine advection(u, v, G, GV, CS, hfld, qfld, ufld, vfld) type(ZB2020_CS), intent(inout) :: CS !< ZB2020 control structure. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -904,18 +1014,28 @@ subroutine advection(u, v, G, GV, CS, h, q) real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & intent(inout) :: v !< The meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, & - intent(inout) :: h !< Input/output array at h points [optional] + intent(inout) :: hfld !< Input/output array at h points [optional] real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)), optional, & - intent(inout) :: q !< Input/output array at q points [optional] + intent(inout) :: qfld !< Input/output array at q points [optional] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), optional, & + intent(inout) :: ufld !< Input/output array at u points [optional] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), optional, & + intent(inout) :: vfld !< Input/output array at v points [optional] ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & - h_advected !< advected quantity at h-points + hfld_advected !< advected quantity at h-points real, dimension(SZIB_(G),SZJB_(G)) :: & - q_advected !< advected quantity at q-points + qfld_advected !< advected quantity at q-points + real, dimension(SZIB_(G),SZJ_(G)) :: & + ufld_advected !< advected quantity at u-points + real, dimension(SZI_(G),SZJB_(G)) :: & + vfld_advected !< advected quantity at v-points real :: uT, vT ! Interpolation of advective velocity to T points real :: uQ, vQ ! Interpolation of advective velocity to Q points + real :: uU, vU ! Interpolation of advective velocity to U points + real :: uV, vV ! Interpolation of advective velocity to V points real :: dx, dy ! Displacement of the fluid particle backward in time in metres real :: x, y ! Nondimension coordinate on unit square of the interpolation point integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq @@ -929,13 +1049,15 @@ subroutine advection(u, v, G, GV, CS, h, q) call pass_vector(u,v,G%Domain) call pass_var(G%dxCv,G%Domain) + call pass_var(G%dyCv,G%Domain) + call pass_var(G%dxCu,G%Domain) call pass_var(G%dyCu,G%Domain) call pass_var(G%dyBu,G%Domain) call pass_var(G%dxBu,G%Domain) call pass_var(G%dyT,G%Domain) call pass_var(G%dxT,G%Domain) - if (present(h)) then + if (present(hfld)) then do k=1,nz do j=js-1,je+1 ; do i=is-1,ie+1 @@ -977,19 +1099,19 @@ subroutine advection(u, v, G, GV, CS, h, q) y = 1. + dy / G%dyBu(i0,j0) endif - h_advected(i,j) = bilin_interp(h(i0,j0,k), h(i0,j0+1,k), & - h(i0+1,j0,k), h(i0+1,j0+1,k), x, y) + hfld_advected(i,j) = bilin_interp(hfld(i0,j0,k), hfld(i0,j0+1,k), & + hfld(i0+1,j0,k), hfld(i0+1,j0+1,k), x, y) enddo; enddo do j=js-1,je+1 ; do i=is-1,ie+1 ! Save computations to storage array and enforcing zero B.C. - h(i,j,k) = h_advected(i,j) * G%mask2dT(i,j) + hfld(i,j,k) = hfld_advected(i,j) * G%mask2dT(i,j) enddo ; enddo enddo ! k loop endif - if (present(q)) then + if (present(qfld)) then do k=1,nz do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 @@ -1031,18 +1153,126 @@ subroutine advection(u, v, G, GV, CS, h, q) y = 1. + dy / G%dyT(I0+1,j0+1) endif - q_advected(I,J) = bilin_interp(q(I0,J0,k), q(I0,J0+1,k), & - q(I0+1,J0,k), q(I0+1,J0+1,k), x, y) + qfld_advected(I,J) = bilin_interp(qfld(I0,J0,k), qfld(I0,J0+1,k), & + qfld(I0+1,J0,k), qfld(I0+1,J0+1,k), x, y) enddo; enddo do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1 ! Save computations to storage array and enforcing zero B.C. - q(I,J,k) = q_advected(I,J) * G%mask2dBu(I,J) + qfld(I,J,k) = qfld_advected(I,J) * G%mask2dBu(I,J) enddo ; enddo enddo ! k loop endif + + if (present(ufld)) then + do k=1,nz + + do j=js-1,je+1 ; do I=Isq-1,Ieq+1 + uU = u(I,j,k) + vU = (v(i,J,k) + v(i,J-1,k)+ v(i+1,J,k) + v(i+1,J-1,k)) * 0.25 + + ! Displacement backward in time, so the sign is minus + dx = - uU * CS%dt + dy = - vU * CS%dt + + ! Now we determine the left bottom corner of the unit square + ! used for interpolation. + if (dx > 0) then + I0 = I + else + I0 = I-1 + endif + + if (dy>0) then + j0 = j + else + j0 = j-1 + endif + + ! Here we determine the relative position of the interpolation point + ! Note in every case the corner point is the center of the square + ! on which bilinear interpolation is performed. + ! If dx<0, we need to convert the non-dimensional coordinate to positive one, + ! which is w.r.t. leftmost corner of the unit square + if (dx > 0.) then + x = dx / G%dxCu(I0,j0) + else + x = 1. + dx / G%dxCu(I0,j0) ! this is smaller than 1 + endif + + if (dy > 0.) then + y = dy / G%dyCu(I0,j0) + else + y = 1. + dy / G%dyCu(I0,j0) + endif + + ufld_advected(I,j) = bilin_interp(ufld(I0,j0,k), ufld(I0,j0+1,k), & + ufld(I0+1,j0,k), ufld(I0+1,j0+1,k), x, y) + enddo; enddo + + do j=js-1,je+1 ; do i=Isq-1,Ieq+1 + ! Save computations to storage array and enforcing zero B.C. + ufld(I,j,k) = ufld_advected(I,j) * G%mask2dCu(I,j) + enddo ; enddo + + enddo ! k loop + endif + + if (present(vfld)) then + do k=1,nz + + do J=Jsq-1,Jeq+1 ; do i=is-1,ie+1 + uV = (u(I,j,k) + u(I-1,j,k) + u(I,j+1,k) + u(I-1,j+1,k)) * 0.25 + vV = v(i,J,k) + + ! Displacement backward in time, so the sign is minus + dx = - uV * CS%dt + dy = - vV * CS%dt + + ! Now we determine the left bottom corner of the unit square + ! used for interpolation. + if (dx > 0) then + i0 = i + else + i0 = i-1 + endif + + if (dy>0) then + J0 = J + else + J0 = J-1 + endif + + ! Here we determine the relative position of the interpolation point + ! Note in every case the corner point is the center of the square + ! on which bilinear interpolation is performed. + ! If dx<0, we need to convert the non-dimensional coordinate to positive one, + ! which is w.r.t. leftmost corner of the unit square + if (dx > 0.) then + x = dx / G%dxCv(i0,J0) + else + x = 1. + dx / G%dxCv(i0,J0) ! this is smaller than 1 + endif + + if (dy > 0.) then + y = dy / G%dyCv(i0,J0) + else + y = 1. + dy / G%dyCv(i0,J0) + endif + + vfld_advected(i,J) = bilin_interp(vfld(i0,J0,k), vfld(i0,J0+1,k), & + vfld(i0+1,J0,k), vfld(i0+1,J0+1,k), x, y) + enddo; enddo + + do j=Jsq-1,Jeq+1 ; do i=is-1,ie+1 + ! Save computations to storage array and enforcing zero B.C. + vfld(i,J,k) = vfld_advected(i,J) * G%mask2dCv(i,J) + enddo ; enddo + + enddo ! k loop + endif end subroutine advection