diff --git a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 index 3dd595a1a1..2255c6c6b1 100644 --- a/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 +++ b/src/parameterizations/lateral/MOM_Zanna_Bolton.F90 @@ -53,9 +53,17 @@ module MOM_Zanna_Bolton !! 0 - no memory !! 1 - Lagrangian memory of eddy forcing !! 2 - Eulerian memory of eddy forcing - real :: ZB_memory_const !< The non-dimensional scaling factor for the eddy memory time scale + integer :: ZB_memory_location !< Select where memory is incorporated within ZB model: + !! 0 - no memory + !! 1 - memory for eddy forcing + !! 2 - memory for large scale flow + !! 3 - memory for both eddy forcing and large scale flow + real :: ZB_eddy_memory_const !< The non-dimensional scaling factor for the eddy memory time scale !! [nondim] - !! Typical range: 1-2 + !! Typical range: 1-10 + real :: ZB_large_scale_memory_const !< The non-dimensional scaling factor for the memory time scale + !! for the velocity gradients [nondim] + !! Typical range: 1-10 real :: ZB_max_memory !< The maximum eddy memory time scale [T ~> s] real :: DT !< The baroclinic model time step [T ~> s] @@ -77,7 +85,11 @@ module MOM_Zanna_Bolton Txx_with_memory, & !< Memory-enhanced subgrid stress xx component [L3 T-2 ~> m3 s-2] Tyy_with_memory, & !< Memory-enhanced subgrid stress yy component [L3 T-2 ~> m3 s-2] Txy_with_memory, & !< Memory-enhanced subgrid stress xy component [L3 T-2 ~> m3 s-2] - ZBtau !< Memory timescale for ZB forcing [T ~> s] + ZBtau_eddy, & !< Memory timescale for ZB forcing [T ~> s] + 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] real, dimension(:,:), allocatable :: & kappa_h, & !< Scaling coefficient in h points [L2 ~> m2] @@ -106,7 +118,7 @@ module MOM_Zanna_Bolton integer :: id_Txx = -1 integer :: id_Tyy = -1 integer :: id_Txy = -1 - integer :: id_tau = -1 + integer :: id_tau_eddy = -1, id_tau_large_scale = -1 integer :: id_cdiss = -1 ! Debug fields integer :: id_sh_xx = -1, id_sh_xy = -1, id_vort_xy = -1 @@ -231,10 +243,21 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) "\t 1 - Lagrangian memory of eddy forcing\n" //& "\t 2 - Eulerian memory of eddy forcing", default=0) - call get_param(param_file, mdl, "ZB_MEMORY_CONST", CS%ZB_memory_const, & + 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 2 - memory for large scale flow\n" //& + "\t 3 - memory for both eddy forcing and large scale flow", default=0) + + call get_param(param_file, mdl, "ZB_EDDY_MEMORY_CONST", CS%ZB_eddy_memory_const, & "The non-dimensional scaling factor for the eddy memory time "//& "scale. A value of zero means no memory.", units="nondim", default=0.0) + call get_param(param_file, mdl, "ZB_LARGE_SCALE_MEMORY_CONST", CS%ZB_large_scale_memory_const, & + "The non-dimensional scaling factor for the large scale memory time "//& + "scale. A value of zero means no memory.", units="nondim", default=0.0) + call get_param(param_file, mdl, "ZB_MAX_MEMORY", CS%ZB_max_memory, & "The maximal eddy memory time. Default is 1 year", units="s", scale=US%s_to_T, & default=31536000.0) @@ -243,12 +266,24 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) "The (baroclinic) dynamics time step.", units="s", scale=US%s_to_T, & fail_if_missing=.true.) - if (CS%ZB_memory_type == 0) then - if (CS%ZB_memory_const > 0) call MOM_error(FATAL, & - "MOM_Zanna_Bolton: ZB_MEMORY_TYPE must be 1 or 2 if ZB_MEMORY_CONST > 0") + if (CS%ZB_memory_location > 0) then + if (CS%ZB_memory_type == 0) call MOM_error(FATAL, & + "MOM_Zanna_Bolton: ZB_MEMORY_TYPE must be 1 or 2 if ZB_MEMORY_LOCATION > 0") + if (CS%ZB_memory_location == 1 .OR. CS%ZB_memory_location == 3) then + if (CS%ZB_eddy_memory_const <= 0) call MOM_error(FATAL, & + "MOM_Zanna_Bolton: ZB_EDDY_MEMORY_CONST must be larger than 0 if ZB_MEMORY_LOCATION == 1 or 3") + endif + if (CS%ZB_memory_location == 2 .OR. CS%ZB_memory_location == 3) then + if (CS%ZB_large_scale_memory_const <= 0) call MOM_error(FATAL, & + "MOM_Zanna_Bolton: ZB_LARGE_SCALE_MEMORY_CONST must be larger than 0 if ZB_MEMORY_LOCATION == 1 or 3") + endif else - if (CS%ZB_memory_const <= 0.0) call MOM_error(FATAL, & - "MOM_Zanna_Bolton: ZB_MEMORY_CONST must be > 0") + if (CS%ZB_eddy_memory_const > 0.0) call MOM_error(FATAL, & + "MOM_Zanna_Bolton: ZB_MEMORY_LOCATION must be 1 or 3 if ZB_EDDY_MEMORY_CONST > 0") + if (CS%ZB_large_scale_memory_const > 0.0) call MOM_error(FATAL, & + "MOM_Zanna_Bolton: ZB_MEMORY_LOCATION must be 2 or 3 if ZB_LARGE_SCALE_MEMORY_CONST > 0") + if (CS%ZB_memory_type > 0) call MOM_error(FATAL, & + "MOM_Zanna_Bolton: ZB_MEMORY_LOCATION must be 1, 2 or 3 if ZB_MEMORY_TYPE > 0") endif ! Register fields for output from this module. @@ -273,8 +308,14 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) if (CS%ZB_memory_type == 1) then - CS%id_tau = register_diag_field('ocean_model', 'ZBtau', diag%axesTL, Time, & - 'Memory timescale on Zanna-Bolton forcing', 's', conversion=US%T_to_s) + if (CS%ZB_memory_location == 1 .OR. CS%ZB_memory_location == 3) then + CS%id_tau_eddy = register_diag_field('ocean_model', 'ZBtau_eddy', diag%axesTL, Time, & + 'Memory timescale on Zanna-Bolton forcing', 's', conversion=US%T_to_s) + endif + if (CS%ZB_memory_location == 2 .OR. CS%ZB_memory_location == 3) then + CS%id_tau_large_scale = register_diag_field('ocean_model', 'ZBtau_large_scale', diag%axesTL, Time, & + 'Memory timescale on large-scale flow in Zanna-Bolton closure', 's', conversion=US%T_to_s) + endif endif if (CS%Klower_R_diss > 0) then @@ -334,10 +375,18 @@ subroutine ZB_2020_init(Time, G, GV, US, param_file, diag, CS, use_ZB2020) allocate(CS%kappa_q(SZIB_(G),SZJB_(G))); CS%kappa_q(:,:) = 0. if (CS%ZB_memory_type > 0) then + if (CS%ZB_memory_location == 1 .OR. CS%ZB_memory_location == 3) then allocate(CS%Txx_with_memory(SZI_(G),SZJ_(G),SZK_(GV))); CS%Txx_with_memory(:,:,:) = 0. allocate(CS%Tyy_with_memory(SZI_(G),SZJ_(G),SZK_(GV))); CS%Tyy_with_memory(:,:,:) = 0. allocate(CS%Txy_with_memory(SZIB_(G),SZJB_(G),SZK_(GV))); CS%Txy_with_memory(:,:,:) = 0. - allocate(CS%ZBtau(SZI_(G),SZJ_(G),SZK_(GV))); CS%ZBtau(:,:,:) = 0. + allocate(CS%ZBtau_eddy(SZI_(G),SZJ_(G),SZK_(GV))); CS%ZBtau_eddy(:,:,:) = 0. + endif + if (CS%ZB_memory_location == 2 .OR. CS%ZB_memory_location == 3) then + allocate(CS%sh_xx_with_memory(SZI_(G),SZJ_(G),SZK_(GV))); CS%sh_xx_with_memory(:,:,:) = 0. + allocate(CS%sh_xy_with_memory(SZIB_(G),SZJB_(G),SZK_(GV))); CS%sh_xy_with_memory(:,:,:) = 0. + 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 endif @@ -426,10 +475,18 @@ subroutine ZB_2020_end(CS) endif if (CS%ZB_memory_type > 0) then - deallocate(CS%Txx_with_memory) - deallocate(CS%Tyy_with_memory) - deallocate(CS%Txy_with_memory) - deallocate(CS%ZBtau) + if (CS%ZB_memory_location == 1 .OR. CS%ZB_memory_location == 3) then + deallocate(CS%Txx_with_memory) + deallocate(CS%Tyy_with_memory) + deallocate(CS%Txy_with_memory) + deallocate(CS%ZBtau_eddy) + endif + if (CS%ZB_memory_location == 2 .OR. CS%ZB_memory_location == 3) then + deallocate(CS%sh_xx_with_memory) + deallocate(CS%sh_xy_with_memory) + deallocate(CS%vort_xy_with_memory) + deallocate(CS%ZBtau_large_scale) + endif endif if (CS%use_ann==1) then @@ -552,6 +609,11 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & ! Sharpen velocity gradients if specified call filter_velocity_gradients(G, GV, CS) + + ! Introduce memory into velocity gradients if specified + if (CS%ZB_memory_location == 2 .OR. CS%ZB_memory_location == 3) then + if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, eddy=.false.) + endif ! Compute the stress tensor given the ! (optionally sharpened) velocity gradients @@ -564,8 +626,8 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & ! Smooth the stress tensor if specified call filter_stress(G, GV, CS) - if (CS%ZB_memory_type > 0) then - call step_forward_ZB_memory(u, v, h, G, GV, CS) + if (CS%ZB_memory_location == 1 .OR. CS%ZB_memory_location == 3) then + if (CS%ZB_memory_type > 0) call step_forward_ZB_memory(u, v, h, G, GV, CS, eddy=.true.) endif ! Compute acceleration from a given stress tensor @@ -595,7 +657,8 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & if (CS%id_cdiss>0) call post_data(CS%id_cdiss, CS%c_diss, CS%diag) - if (CS%id_tau>0) call post_data(CS%id_tau, CS%ZBtau, CS%diag) + if (CS%id_tau_eddy>0) call post_data(CS%id_tau_eddy, CS%ZBtau_eddy, CS%diag) + if (CS%id_tau_large_scale>0) call post_data(CS%id_tau_large_scale, CS%ZBtau_large_scale, CS%diag) ! Debug fields if (CS%id_sh_xx>0) call post_data(CS%id_sh_xx, CS%sh_xx, CS%diag) if (CS%id_sh_xy>0) call post_data(CS%id_sh_xy, CS%sh_xy, CS%diag) @@ -614,7 +677,9 @@ subroutine Zanna_Bolton_2020(u, v, h, diffu, diffv, G, GV, CS, & end subroutine Zanna_Bolton_2020 -subroutine step_forward_ZB_memory(u, v, h, G, GV, CS) + + +subroutine step_forward_ZB_memory(u, v, h, G, GV, CS, eddy) 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. @@ -624,8 +689,10 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS) 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 + logical, intent(in) :: eddy !< If true, step forward memory on eddy forcing; + ! if false, step forward memory on velocity gradients + ! Local variables real :: sh_xy_h !< Shearing strain interpolated to h point [T-1 ~> s-1] real :: sh_xx_q !< Horizontal tension interpolated to q point [T-1 ~> s-1] @@ -638,10 +705,15 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS) 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 - - call pass_var(CS%Txx, G%Domain) - call pass_var(CS%Tyy, G%Domain) - call pass_var(CS%Txy, G%Domain, position=CORNER) + if (eddy) then + call pass_var(CS%Txx, G%Domain) + call pass_var(CS%Tyy, G%Domain) + call pass_var(CS%Txy, G%Domain, position=CORNER) + else + call pass_var(CS%sh_xx, G%Domain) + call pass_var(CS%sh_xy, G%Domain, position=CORNER) + call pass_var(CS%vort_xy, G%Domain, position=CORNER) + endif !if (CS%Txx_with_memory(20,20,1) == 0.) then ! ! This is lazy initial condition @@ -652,29 +724,57 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS) !endif if (CS%ZB_memory_type == 1) then - call tensor_advection(u, v, G, GV, CS) + if (eddy) then + 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) + 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) + endif endif do k=1,nz - ! relax Txx, Tyy to ZB: this is the Eulerian part of the averaging + ! relaxation: this is the Eulerian part of the averaging do j=js-1,je+1 ; do i=is-1,ie+1 ! 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)) ) D = sqrt(CS%sh_xx(i,j,k) * CS%sh_xx(i,j,k) + sh_xy_h * sh_xy_h) - tau = CS%ZB_memory_const / (D + 1.0/CS%ZB_max_memory) - !tau = 7776000.0 ! 90 days - CS%ZBtau(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) & - + CS%dt / (tau + CS%dt) * CS%Tyy(i,j,k) - - CS%Txx(i,j,k) = CS%Txx_with_memory(i,j,k) - CS%Tyy(i,j,k) = CS%Tyy_with_memory(i,j,k) + + if (eddy) then + tau = CS%ZB_eddy_memory_const / (D + 1.0/CS%ZB_max_memory) + !tau = 7776000.0 ! 90 days + 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) & + + CS%dt / (tau + CS%dt) * CS%Tyy(i,j,k) + + 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 + tau = CS%ZB_large_scale_memory_const / (D + 1.0/CS%ZB_max_memory) + !tau = 7776000.0 ! 90 days + CS%ZBtau_large_scale(i,j,k) = tau + + ! Eulerian memory with implicit time stepping + CS%sh_xx_with_memory(i,j,k) = tau / (tau + CS%dt) * CS%sh_xx_with_memory(i,J,k) & + + CS%dt / (tau + CS%dt) * CS%sh_xx(i,j,k) + + CS%sh_xx(i,j,k) = CS%sh_xx_with_memory(i,j,k) + endif enddo ; enddo @@ -683,27 +783,47 @@ subroutine step_forward_ZB_memory(u, v, h, G, GV, CS) 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) - tau = CS%ZB_memory_const / (D + 1.0/CS%ZB_max_memory) - !tau = 7776000.0 ! 90 days - 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) + if (eddy) then + tau = CS%ZB_eddy_memory_const / (D + 1.0/CS%ZB_max_memory) + !tau = 7776000.0 ! 90 days - enddo ; enddo + 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 + + 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) + CS%vort_xy_with_memory(I,J,k) = tau / (tau + CS%dt) * CS%vort_xy_with_memory(i,J,k) & + + CS%dt / (tau + CS%dt) * CS%vort_xy(I,J,k) + + CS%sh_xy(I,J,k) = CS%sh_xy_with_memory(I,J,k) + CS%vort_xy(I,J,k) = CS%vort_xy_with_memory(I,J,k) + endif + enddo ; enddo enddo ! end of k loop - call pass_var(CS%Txx, G%Domain) - call pass_var(CS%Tyy, G%Domain) - call pass_var(CS%Txy, G%Domain, position=CORNER) + if (eddy) then + call pass_var(CS%Txx, G%Domain) + call pass_var(CS%Tyy, G%Domain) + call pass_var(CS%Txy, G%Domain, position=CORNER) + else + call pass_var(CS%sh_xx, G%Domain) + call pass_var(CS%sh_xy, G%Domain, position=CORNER) + call pass_var(CS%vort_xy, G%Domain, position=CORNER) + endif end subroutine step_forward_ZB_memory -subroutine tensor_advection(u, v, G, GV, CS) + +subroutine advection(u, v, G, GV, CS, h, q) 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. @@ -711,13 +831,16 @@ subroutine tensor_advection(u, v, G, GV, CS) 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)), optional, & + intent(inout) :: h !< 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] ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & - Txx_with_memory, & !< advected tensor-component Txx [L3 T-2 ~> m3 s-2] - Tyy_with_memory !< advected tensor-component Tyy [L3 T-2 ~> m3 s-2] + h_advected !< advected quantity at h-points real, dimension(SZIB_(G),SZJB_(G)) :: & - Txy_with_memory !< advected tensor-component Txy [L3 T-2 ~> m3 s-2] + q_advected !< advected quantity at q-points real :: uT, vT ! Interpolation of advective velocity to T points real :: uQ, vQ ! Interpolation of advective velocity to Q points @@ -732,11 +855,6 @@ subroutine tensor_advection(u, v, G, GV, CS) 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 - - - 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 pass_vector(u,v,G%Domain) call pass_var(G%dxCv,G%Domain) call pass_var(G%dyCu,G%Domain) @@ -745,58 +863,62 @@ subroutine tensor_advection(u, v, G, GV, CS) call pass_var(G%dyT,G%Domain) call pass_var(G%dxT,G%Domain) - do k=1,nz - - do j=js-1,je+1 ; do i=is-1,ie+1 - uT = (u(I,j,k) + u(I-1,j,k)) * 0.5 - vT = (v(i,J,k) + v(i,J-1,k)) * 0.5 - - ! Displacement backward in time, so the sign is minus - dx = - uT * CS%dt - dy = - vT * 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 (present(h)) then + do k=1,nz - if (dy>0) then - j0 = j - else - j0 = j-1 - endif + do j=js-1,je+1 ; do i=is-1,ie+1 + uT = (u(I,j,k) + u(I-1,j,k)) * 0.5 + vT = (v(i,J,k) + v(i,J-1,k)) * 0.5 + + ! Displacement backward in time, so the sign is minus + dx = - uT * CS%dt + dy = - vT * 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%dxBu(i0,j0) + else + x = 1. + dx / G%dxBu(i0,j0) ! this is smaller than 1 + endif + + if (dy > 0.) then + y = dy / G%dyBu(i0,j0) + else + 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) + 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) + enddo ; enddo - ! 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%dxBu(i0,j0) - else - x = 1. + dx / G%dxBu(i0,j0) ! this is smaller than 1 - endif - - if (dy > 0.) then - y = dy / G%dyBu(i0,j0) - else - y = 1. + dy / G%dyBu(i0,j0) - endif + enddo ! k loop + endif - Txx_with_memory(i,j) = bilin_interp(CS%Txx_with_memory(i0,j0,k), CS%Txx_with_memory(i0,j0+1,k), & - CS%Txx_with_memory(i0+1,j0,k), CS%Txx_with_memory(i0+1,j0+1,k), x, y) - Tyy_with_memory(i,j) = bilin_interp(CS%Tyy_with_memory(i0,j0,k), CS%Tyy_with_memory(i0,j0+1,k), & - CS%Tyy_with_memory(i0+1,j0,k), CS%Tyy_with_memory(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. - CS%Txx_with_memory(i,j,k) = Txx_with_memory(i,j) * G%mask2dT(i,j) - CS%Tyy_with_memory(i,j,k) = Tyy_with_memory(i,j) * G%mask2dT(i,j) - enddo ; enddo + if (present(q)) then + do k=1,nz do J=Jsq-1,Jeq+1 ; do I=Isq-1,Ieq+1 uQ = (u(I,j,k) + u(I,j+1,k)) * 0.5 @@ -837,19 +959,20 @@ subroutine tensor_advection(u, v, G, GV, CS) y = 1. + dy / G%dyT(I0+1,j0+1) endif - Txy_with_memory(I,J) = bilin_interp(CS%Txy_with_memory(I0,J0,k), CS%Txy_with_memory(I0,J0+1,k), & - CS%Txy_with_memory(I0+1,J0,k), CS%Txy_with_memory(I0+1,J0+1,k), x, y) + 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) enddo; enddo do J=Jsq-1,Jeq+1 ; do i=Isq-1,Ieq+1 ! Save computations to storage array and enforcing zero B.C. - CS%Txy_with_memory(I,J,k) = Txy_with_memory(I,J) * G%mask2dBu(I,J) + q(I,J,k) = q_advected(I,J) * G%mask2dBu(I,J) enddo ; enddo - enddo ! k loop + enddo ! k loop + endif -end subroutine tensor_advection +end subroutine advection ! Interpolation on a unit square is explained here: