diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e167054cb6..ccfc5f4a7e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -274,6 +274,7 @@ target_sources(mom6shared PRIVATE tracer/MOM_tracer_diabatic.F90 tracer/MOM_tracer_flow_control.F90 tracer/MOM_tracer_hor_diff.F90 + tracer/MOM_tracer_numerical_mixing.F90 tracer/MOM_tracer_registry.F90 tracer/MOM_tracer_types.F90 tracer/MOM_tracer_Z_init.F90 diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 8bb85762bb..bd5d74d6d9 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -29,7 +29,7 @@ module MOM_diagnostics use MOM_interface_heights, only : find_eta, find_col_mass use MOM_spatial_means, only : global_area_mean, global_layer_mean use MOM_spatial_means, only : global_volume_mean, global_area_integral -use MOM_tracer_registry, only : tracer_registry_type, post_tracer_transport_diagnostics +use MOM_tracer_registry, only : tracer_type, tracer_registry_type, post_tracer_transport_diagnostics use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, ocean_internal_state, p3d use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, surface @@ -164,7 +164,6 @@ module MOM_diagnostics !>@} end type transport_diag_IDs - contains !> Diagnostics not more naturally calculated elsewhere are computed here. subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & @@ -1654,7 +1653,7 @@ subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dy real :: Idt ! The inverse of the time interval [T-1 ~> s-1] real :: H_to_RZ_dt ! A conversion factor from accumulated transports to fluxes ! [R Z H-1 T-1 ~> kg m-3 s-1 or s-1]. - integer :: i, j, k, is, ie, js, je, nz + integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Idt = 1. / dt_trans @@ -1705,7 +1704,7 @@ subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dy call post_data(IDs%id_dynamics_h_tendency, h_tend, diag, alt_h=diag_pre_dyn%h_state) endif - call post_tracer_transport_diagnostics(G, GV, Reg, diag_pre_dyn%h_state, diag) + call post_tracer_transport_diagnostics(G, GV, Reg, diag_pre_dyn%h_state, diag, uhtr, vhtr, h, dt_trans, Idt) call diag_restore_grids(diag) diff --git a/src/tracer/MOM_tracer_numerical_mixing.F90 b/src/tracer/MOM_tracer_numerical_mixing.F90 new file mode 100644 index 0000000000..90cbffba95 --- /dev/null +++ b/src/tracer/MOM_tracer_numerical_mixing.F90 @@ -0,0 +1,160 @@ +!> Functions and routines involved in calculating numerical mixing of tracers due to advection schemes. +module MOM_tracer_numerical_mixing + +use MOM_grid, only : ocean_grid_type +use MOM_tracer_types, only : tracer_type +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +#include + +public advection_scheme_variance_production + +contains + +!< Calculate the spurious variance production of tracer `Tr` due to the advection schemes. +subroutine advection_scheme_variance_production(G, GV, Tr, h_diag, h, dt_trans, Idt, uhtr, vhtr, asvp) + + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(tracer_type), intent(in) :: Tr !< Pointer to the tracer regsitry + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_diag !< Previous thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Updated layer thicknesses [H ~> m or kg m-2] + real, intent(in) :: dt_trans !< The transport time interval [T ~> s] + real, intent(in) :: Idt !< Inverse time interval [T-1 ~> s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: uhtr !< Accumulated zonal thickness fluxes + !! used to advect tracers [H L2 ~> m3 or kg] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vhtr !< Accumulated meridional thickness fluxes + !! used to advect tracers [H L2 ~> m3 or kg] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: asvp !< Advection scheme varianve production + !! diagnostic [CU2 H T-1 ~> conc2 m s-1] + + call thickness_weighted_variance_advection(G, GV, Tr, h_diag, h, dt_trans, Idt, asvp) + call thickness_weighted_variance_flux_divergence(G, Gv, Tr, uhtr, vhtr, Idt, asvp) + +end subroutine advection_scheme_variance_production + +!< Subroutine to calculate the thickness weighted variance advection over the transport timestep. +subroutine thickness_weighted_variance_advection(G, GV, Tr, h_diag, h, dt, Idt, asvp) + + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(tracer_type), intent(in) :: Tr !< Pointer to the tracer registry + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_diag !< Previous thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Updated thicknesses [H ~> m or kg m-2] + real, intent(in) :: dt !< Transport time interval [T ~> s] + real, intent(in) :: Idt !< Inverse time interval [T-1 ~> s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: asvp !< Array to store asvp in + !! [CU2 H T-1 ~> conc2 m s-1] + + !< Local variables + integer :: is, ie, js, je, nz !< Grid cell centre and layer indexes + integer :: i, j, k !< Counters + real :: Ih !< Inverse updated thickness [H-1 ~> m-1] + real :: ht_prev !< Thickness weighted tracer prior to dynamics [CU H ~> conc m] + real :: ht_adv !< Thickness weighted tracer after lateral advection [CU H ~> conc m] + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + do k=1,nz ; do j=js,je ; do i=is,ie + Ih = 1 / h(i,j,k) + ht_prev = h_diag(i,j,k) * Tr%t_prev(i,j,k) + ht_adv = ht_prev + dt * Tr%advection_xy(i,j,k) + asvp(i,j,k) = ( (Ih * (ht_adv*ht_adv)) - (ht_prev * Tr%t_prev(i,j,k)) ) * Idt + enddo ; enddo ; enddo + +end subroutine thickness_weighted_variance_advection + +!< Subroutine to calculate the divergence of the variance flux due to the advection scheme +subroutine thickness_weighted_variance_flux_divergence(G, Gv, Tr, uhtr, vhtr, Idt, asvp) + + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(tracer_type), intent(in) :: Tr !< Pointer to the tracer registry + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: uhtr !< Accumulated zonal thickness fluxes + !! used to advect tracers [H L2 ~> m3 or kg] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vhtr !< Accumulated meridional thickness fluxes + !! used to advect tracers [H L2 ~> m3 or kg] + real, intent(in) :: Idt !< Inverse time interval [T-1 ~> s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: asvp !< Array to store asvp in + !! [CU2 H T-1 ~> conc2 m s-1] + + !< Local variables + integer :: is, ie, js, je, nz !< Grid cell centre and layer indexes + integer :: i, j, k !< Counters + real :: east, west, north, south !< east, west, north and south faces of h point + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: x_upwind !< Zonal upwind values for tracer [CU ~> conc] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: y_upwind !< Meridional upwind values for tracer [CU ~> conc] + + x_upwind(:,:,:) = 0. + y_upwind(:,:,:) = 0. + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + call zonal_upwind_values(G, GV, Tr, uhtr, x_upwind) + call meridional_upwind_values(G, GV, Tr, vhtr, y_upwind) + + do k=1,nz ; do j=js,je ; do i=is,ie + east = (2 * (Tr%ad_x(I,j,k) *x_upwind(I,j,k))) - ((Idt*uhtr(I,j,k)) * (x_upwind(I,j,k) *x_upwind(I,j,k))) + west = (2 * (Tr%ad_x(I-1,j,k)*x_upwind(I-1,j,k))) - ((Idt*uhtr(I-1,j,k)) * (x_upwind(I-1,j,k)*x_upwind(I-1,j,k))) + north = (2 * (Tr%ad_y(i,J,k) *y_upwind(i,J,k))) - ((Idt*vhtr(i,J,k)) * (y_upwind(i,J,k) *y_upwind(i,J,k))) + south = (2 * (Tr%ad_y(i,J-1,k)*y_upwind(i,J-1,k))) - ((Idt*vhtr(i,J-1,k)) * (y_upwind(i,J-1,k)*y_upwind(i,J-1,k))) + asvp(i,j,k) = asvp(i,j,k) + (((east - west) + (north - south)) * G%IareaT(i,j)) + enddo ; enddo; enddo + +end subroutine thickness_weighted_variance_flux_divergence + +!< Subroutine to calculate upwind values in zonal direction +subroutine zonal_upwind_values(G, GV, Tr, uhtr, x_upwind) + + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(tracer_type), intent(in) :: Tr !< Pointer to the tracer registry + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: uhtr !< Accumulated zonal thickness fluxes + !! used to advect tracers [H L2 ~> m3 or kg] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: x_upwind !< zonal upwind tracer [CU ~> conc] + + !< Local variables + integer :: is, ie, js, je, nz !< Grid cell centre indexes + integer :: i, j, k !< Counters + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + do k=1,nz ; do j=js,je ; do I=is-1,ie + if (uhtr(I,j,k) >= 0) then + x_upwind(I,j,k) = Tr%t_prev(i,j,k) + elseif (uhtr(I,j,k) < 0) then + x_upwind(I,j,k) = Tr%t_prev(i+1,j,k) + endif + enddo ; enddo ; enddo + +end subroutine zonal_upwind_values + +!< Subroutine to calculate upwind values in the meridional direction +subroutine meridional_upwind_values(G, GV, Tr, vhtr, y_upwind) + + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure + type(tracer_type), intent(in) :: Tr !< Pointer to tracer registry + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: vhtr !< Accumulated meridional thickness fluxes used + !! to advect tracers [H L2 ~> m3 or kg] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)),intent(inout) :: y_upwind !< Meridional upwind tracer values [CU ~> conc] + + !< Local variables + integer :: is, ie, js, je, nz !< Grid cell centre indexes + integer :: i, j, k !< Counters + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + do k=1,nz ; do J=js-1,je ; do i=is,ie + if (vhtr(i,J,k) >= 0) then + y_upwind(i,J,k) = Tr%t_prev(i,j,k) + elseif (vhtr(i,J,k) < 0) then + y_upwind(i,J,k) = Tr%t_prev(i,j+1,k) + endif + enddo ; enddo ; enddo + +end subroutine meridional_upwind_values + +end module MOM_tracer_numerical_mixing diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 50220be343..c647c4afda 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -14,6 +14,7 @@ module MOM_tracer_registry use MOM_diag_mediator, only : diag_ctrl, register_diag_field, post_data, safe_alloc_ptr use MOM_diag_mediator, only : diag_grid_storage use MOM_diag_mediator, only : diag_copy_storage_to_diag, diag_save_grids, diag_restore_grids +use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_hor_index, only : hor_index_type @@ -25,6 +26,7 @@ module MOM_tracer_registry use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type use MOM_tracer_types, only : tracer_type, tracer_registry_type +use MOM_tracer_numerical_mixing, only : advection_scheme_variance_production implicit none ; private @@ -382,6 +384,10 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u "flux from the horizontal boundary diffusion scheme", trim(flux_units), & v_extensive=.true., & x_cell_method='sum', conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) + Tr%id_advection_scheme_variance_production = register_diag_field("ocean_model", & + trim(shortnm)//"_advection_scheme_variance_production", diag%axesTL, Time, & + "Spurious variance production of "//trim(shortnm)//" variance due to advection", & + trim(Tr%units)//"2 m s-1", conversion=(TR%conc_scale**2)*GV%H_to_MKS*US%s_to_T) else Tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", & diag%axesCuL, Time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), & @@ -407,6 +413,10 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u "Horizontal Boundary Diffusive Meridional Flux of "//trim(flux_longname), & flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T, & x_cell_method='sum') + Tr%id_advection_scheme_variance_production = register_diag_field("ocean_model", & + trim(shortnm)//"_advection_scheme_variance_production", diag%axesTL, Time, & + "Spurious variance production of "//trim(shortnm)//" variance due to advection", & + trim(Tr%units)//"2 m s-1", conversion=(TR%conc_scale**2)*GV%H_to_MKS*US%s_to_T) endif Tr%id_zint = register_diag_field("ocean_model", trim(shortnm)//"_zint", & diag%axesT1, Time, & @@ -418,8 +428,10 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u trim(units) // " m") Tr%id_surf = register_diag_field("ocean_model", trim(shortnm)//"_SURF", & diag%axesT1, Time, "Surface values of "// trim(longname), trim(units)) - if (Tr%id_adx > 0) call safe_alloc_ptr(Tr%ad_x,IsdB,IedB,jsd,jed,nz) - if (Tr%id_ady > 0) call safe_alloc_ptr(Tr%ad_y,isd,ied,JsdB,JedB,nz) + if ((Tr%id_adx > 0) .or. (Tr%id_advection_scheme_variance_production > 0)) & + call safe_alloc_ptr(Tr%ad_x,IsdB,IedB,jsd,jed,nz) + if ((Tr%id_ady > 0) .or. (Tr%id_advection_scheme_variance_production > 0)) & + call safe_alloc_ptr(Tr%ad_y,isd,ied,JsdB,JedB,nz) if (Tr%id_dfx > 0) call safe_alloc_ptr(Tr%df_x,IsdB,IedB,jsd,jed,nz) if (Tr%id_dfy > 0) call safe_alloc_ptr(Tr%df_y,isd,ied,JsdB,JedB,nz) if (Tr%id_hbd_dfx > 0) call safe_alloc_ptr(Tr%hbd_dfx,IsdB,IedB,jsd,jed,nz) @@ -468,7 +480,7 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u diag%axesT1, Time, & 'Vertical sum of horizontal convergence of residual mean advective fluxes of '//& trim(lowercase(flux_longname)), conv_units, conversion=Tr%conv_scale*US%s_to_T) - if ((Tr%id_adv_xy > 0) .or. (Tr%id_adv_xy_2d > 0)) & + if ((Tr%id_adv_xy > 0) .or. (Tr%id_adv_xy_2d > 0) .or. (Tr%id_advection_scheme_variance_production > 0)) & call safe_alloc_ptr(Tr%advection_xy,isd,ied,jsd,jed,nz) Tr%id_tendency = register_diag_field('ocean_model', trim(shortnm)//'_tendency', & @@ -476,9 +488,9 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u 'Net time tendency for '//trim(lowercase(longname)), & trim(units)//' s-1', conversion=Tr%conc_scale*US%s_to_T) - if (Tr%id_tendency > 0) then + if ((Tr%id_tendency > 0) .or. (Tr%id_advection_scheme_variance_production > 0)) then call safe_alloc_ptr(Tr%t_prev,isd,ied,jsd,jed,nz) - do k=1,nz ; do j=js,je ; do i=is,ie + do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2 Tr%t_prev(i,j,k) = Tr%t(i,j,k) enddo ; enddo ; enddo endif @@ -733,6 +745,12 @@ subroutine post_tracer_diagnostics_at_sync(Reg, h, diag_prev, diag, G, GV, dt) enddo ; enddo ; enddo call post_data(Tr%id_tendency, work3d, diag, alt_h=diag_prev%h_state) endif + if (Tr%id_advection_scheme_variance_production > 0) then + call pass_var(Tr%t, G%Domain, halo=2) + do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2 + tr%t_prev(i,j,k) = Tr%t(i,j,k) + enddo ; enddo ; enddo + endif if ((Tr%id_trxh_tendency > 0) .or. (Tr%id_trxh_tendency_2d > 0)) then do k=1,nz ; do j=js,je ; do i=is,ie work3d(i,j,k) = (Tr%t(i,j,k)*h(i,j,k) - Tr%Trxh_prev(i,j,k)) * Idt @@ -754,13 +772,23 @@ subroutine post_tracer_diagnostics_at_sync(Reg, h, diag_prev, diag, G, GV, dt) end subroutine post_tracer_diagnostics_at_sync !> Post the advective and diffusive tendencies -subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) +subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag, uhtr, vhtr, h, dt_trans, Idt) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(tracer_registry_type), pointer :: Reg !< pointer to the tracer registry real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(in) :: h_diag !< Layer thicknesses on which to post fields [H ~> m or kg m-2] type(diag_ctrl), intent(in) :: diag !< structure to regulate diagnostic output + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: uhtr !< Accumulated zonal thickness fluxes + !! used to advect tracers [H L2 ~> m3 or kg] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in) :: vhtr !< Accumulated meridional thickness fluxes + !! used to advect tracers [H L2 ~> m3 or kg] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< The updated layer thicknesses [H ~> m or kg m-2] + real, intent(in) :: dt_trans !< The transport time interval [T ~> s] + real, intent(in) :: Idt !< The inverse of the time interval [T-1 ~> s-1] integer :: i, j, k, is, ie, js, je, nz, m, khi real :: work2d(SZI_(G),SZJ_(G)) ! The vertically integrated convergence of lateral advective @@ -768,9 +796,14 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) real :: frac_under_100m(SZI_(G),SZJ_(G),SZK_(GV)) ! weights used to compute 100m vertical integrals [nondim] real :: ztop(SZI_(G),SZJ_(G)) ! position of the top interface [H ~> m or kg m-2] real :: zbot(SZI_(G),SZJ_(G)) ! position of the bottom interface [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: asvp ! Advection scheme varianve production of a + ! tracer [CU2 H T-1 ~> conc2 m s-1] + real :: H_to_RZ_dt ! A conversion factor from accumulated transports to fluxes + ! [R Z H-1 T-1 ~> kg m-3 s-1 or s-1]. type(tracer_type), pointer :: Tr=>NULL() is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + H_to_RZ_dt = GV%H_to_RZ * Idt ! If any tracers are posting 100m vertical integrals, compute weights frac_under_100m(:,:,:) = 0.0 @@ -821,6 +854,12 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) call post_data(Tr%id_adv_xy_2d, work2d, diag) endif + if (Tr%id_advection_scheme_variance_production > 0) then + asvp(:,:,:) = 0. + call advection_scheme_variance_production(G, GV, Tr, h_diag, h, dt_trans, Idt, uhtr, vhtr, asvp) + call post_data(Tr%id_advection_scheme_variance_production, asvp, diag, alt_h=h_diag) + endif + ! A few diagnostics introduce with MARBL driver ! Compute full-depth vertical integral if (Tr%id_zint > 0) then diff --git a/src/tracer/MOM_tracer_types.F90 b/src/tracer/MOM_tracer_types.F90 index 730a453695..2c43d4829a 100644 --- a/src/tracer/MOM_tracer_types.F90 +++ b/src/tracer/MOM_tracer_types.F90 @@ -120,6 +120,7 @@ module MOM_tracer_types integer :: id_tr_vardec = -1 integer :: id_zint = -1, id_zint_100m = -1, id_surf = -1 integer :: id_net_surfflux = -1, id_NLT_tendency = -1, id_NLT_budget = -1 + integer :: id_advection_scheme_variance_production = -1 !>@} end type tracer_type