From b9f73c5f07367d2d313f23bf73265aa979b47836 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 23 Apr 2025 17:14:00 -0400 Subject: [PATCH] +Add post_tracer_integral_diagnostics Added the new routine post_tracer_integral_diagnostics to calculate vertically integrated tracer diagnostics. This includes properly setting the vertical layer extents from layer thicknesses via a call to thickness_to_dz. This new routine is being called from within step_MOM(), immediately after another closely related call to write out other tracer diagnostics. In so doing, it changes the thicknesses used for these diagnostics to be consistent with the state of the tracers. This commit also sets the appropriate unit conversion factors in the register_diag_field calls for three recently diagnostics of surface tracer concentrations and the integrated tracers amounts in the full water column and in the topmost 100 m. All solutions are bitwise identical, but this will correct the thicknesses used in the calculation of two groups of tracer diagnostics and properly implements the dimensional rescaling of these diagnostics when tracers (like temperature and salinity) when they are being rescaled. --- src/core/MOM.F90 | 3 + src/tracer/MOM_tracer_registry.F90 | 119 +++++++++++++++++++---------- 2 files changed, 83 insertions(+), 39 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 156a397ff6..33e08a4a2a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -139,6 +139,7 @@ module MOM use MOM_tracer_registry, only : tracer_registry_type, register_tracer, tracer_registry_init use MOM_tracer_registry, only : register_tracer_diagnostics, post_tracer_diagnostics_at_sync use MOM_tracer_registry, only : post_tracer_transport_diagnostics, MOM_tracer_chksum +use MOM_tracer_registry, only : post_tracer_integral_diagnostics use MOM_tracer_registry, only : preALE_tracer_diagnostics, postALE_tracer_diagnostics use MOM_tracer_registry, only : lock_tracer_registry, tracer_registry_end use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_CS @@ -1053,6 +1054,8 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS CS%CDp, p_surf, CS%t_dyn_rel_diag, CS%diag_pre_sync,& G, GV, US, CS%diagnostics_CSp) call post_tracer_diagnostics_at_sync(CS%Tracer_reg, h, CS%diag_pre_sync, CS%diag, G, GV, CS%t_dyn_rel_diag) + call post_tracer_integral_diagnostics(G, GV, US, CS%Tracer_reg, h, CS%tv, CS%diag) + call diag_copy_diag_to_storage(CS%diag_pre_sync, h, CS%diag) if (showCallTree) call callTree_waypoint("finished calculate_diagnostic_fields (step_MOM)") call disable_averaging(CS%diag) diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 9c69a06c7c..a60542355c 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -16,11 +16,13 @@ module MOM_tracer_registry use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_hor_index, only : hor_index_type use MOM_grid, only : ocean_grid_type +use MOM_interface_heights, only : thickness_to_dz use MOM_io, only : vardesc, query_vardesc, cmor_long_std use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_tracer_types, only : tracer_type, tracer_registry_type @@ -32,6 +34,7 @@ module MOM_tracer_registry public MOM_tracer_chksum, MOM_tracer_chkinv public register_tracer_diagnostics public post_tracer_diagnostics_at_sync, post_tracer_transport_diagnostics +public post_tracer_integral_diagnostics public preALE_tracer_diagnostics, postALE_tracer_diagnostics public tracer_registry_init, lock_tracer_registry, tracer_registry_end public tracer_name_lookup @@ -409,13 +412,13 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u Tr%id_zint = register_diag_field("ocean_model", trim(shortnm)//"_zint", & diag%axesT1, Time, & "Thickness-weighted integral of " // trim(longname), & - trim(units) // " m") + trim(units) // " m", conversion=Tr%conc_scale*US%Z_to_m) Tr%id_zint_100m = register_diag_field("ocean_model", trim(shortnm)//"_zint_100m", & diag%axesT1, Time, & "Thickness-weighted integral of "// trim(longname) // " over top 100m", & - trim(units) // " m") + trim(units) // " m", conversion=Tr%conc_scale*US%Z_to_m) Tr%id_surf = register_diag_field("ocean_model", trim(shortnm)//"_SURF", & - diag%axesT1, Time, "Surface values of "// trim(longname), trim(units)) + diag%axesT1, Time, "Surface values of "// trim(longname), trim(units), conversion=Tr%conc_scale) 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_dfx > 0) call safe_alloc_ptr(Tr%df_x,IsdB,IedB,jsd,jed,nz) @@ -760,45 +763,13 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) 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 - integer :: i, j, k, is, ie, js, je, nz, m, khi + integer :: i, j, k, is, ie, js, je, nz, m real :: work2d(SZI_(G),SZJ_(G)) ! The vertically integrated convergence of lateral advective ! tracer fluxes [CU H T-1 ~> conc m s-1 or conc kg m-2 s-1] - 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] type(tracer_type), pointer :: Tr=>NULL() is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - ! If any tracers are posting 100m vertical integrals, compute weights - frac_under_100m(:,:,:) = 0.0 - ! khi will be the largest layer index corresponding where ztop < 100m and ztop >= 100m - ! in any column (we can reduce computation of 100m integrals by only looping through khi - ! rather than GV%ke) - khi = 0 - do m=1,Reg%ntr ; if (Reg%Tr(m)%registry_diags) then - Tr => Reg%Tr(m) - if (Tr%id_zint_100m > 0) then - zbot(:,:) = 0.0 - do k=1, nz - do j=js,je ; do i=is,ie - ztop(i,j) = zbot(i,j) - zbot(i,j) = ztop(i,j) + h_diag(i,j,k)*GV%H_to_m - if (zbot(i,j) <= 100.0) then - frac_under_100m(i,j,k) = 1.0 - elseif (ztop(i,j) < 100.0) then - frac_under_100m(i,j,k) = (100.0 - ztop(i,j)) / (zbot(i,j) - ztop(i,j)) - else - frac_under_100m(i,j,k) = 0.0 - endif - ! frac_under_100m(i,j,k) = max(0, min(1.0, (100.0 - ztop(i,j)) / (zbot(i,j) - ztop(i,j)))) - enddo ; enddo - if (any(frac_under_100m(:,:,k) > 0)) khi = k - enddo - exit - endif - endif; enddo - do m=1,Reg%ntr ; if (Reg%Tr(m)%registry_diags) then Tr => Reg%Tr(m) if (Tr%id_tr_post_horzn> 0) call post_data(Tr%id_tr_post_horzn, Tr%t, diag) @@ -818,13 +789,83 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) enddo ; enddo ; enddo call post_data(Tr%id_adv_xy_2d, work2d, diag) endif + endif ; enddo +end subroutine post_tracer_transport_diagnostics + +!> Post diagnostics of vertically integrated tracer amouints +subroutine post_tracer_integral_diagnostics(G, GV, US, Reg, h_diag, tv, diag) + 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(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + 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(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables. + type(diag_ctrl), intent(in) :: diag !< structure to regulate diagnostic output + + integer :: i, j, k, is, ie, js, je, nz, m, khi + real :: work2d(SZI_(G),SZJ_(G)) ! The vertically integrated tracer amounts [CU Z T-1 ~> conc m] + real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) !< Geometric layer thicknesses in height units [Z ~> m] + 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 [Z ~> m] + real :: zbot(SZI_(G),SZJ_(G)) ! position of the bottom interface [Z ~> m] + real :: Z_100 ! 100 m in depth units [Z ~> m] + logical :: dz_needed, dz100_used + type(tracer_type), pointer :: Tr=>NULL() + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + dz_needed = .false. + dz100_used = .false. + do m=1,Reg%ntr ; if (Reg%Tr(m)%registry_diags) then + if (Reg%Tr(m)%id_zint_100m > 0) dz100_used = .true. + if (Reg%Tr(m)%id_zint > 0) dz_needed = .true. + endif ; enddo + if (dz100_used) dz_needed = .true. + + if (dz_needed) then + ! Convert the layer thicknesses into geometric depths, using the pre-stored layer-mean specific + ! volumes when in non-Boussinesq mode. + call thickness_to_dz(h_diag, tv, dz, G, GV, US) + endif + + if (dz100_used) then + ! If any tracers are posting 100m vertical integrals, compute weights + frac_under_100m(:,:,:) = 0.0 + ! khi will be the largest layer index corresponding where ztop < 100m and ztop >= 100m + ! in any column (we can reduce computation of 100m integrals by only looping through khi + ! rather than GV%ke) + khi = 0 + + Z_100 = 100.0*US%m_to_Z + zbot(:,:) = 0.0 + do k=1,nz + do j=js,je ; do i=is,ie + ztop(i,j) = zbot(i,j) + zbot(i,j) = ztop(i,j) + dz(i,j,k) + if (zbot(i,j) <= Z_100) then + frac_under_100m(i,j,k) = 1.0 + elseif (ztop(i,j) < Z_100) then + frac_under_100m(i,j,k) = (Z_100 - ztop(i,j)) / (zbot(i,j) - ztop(i,j)) + else + frac_under_100m(i,j,k) = 0.0 + endif + ! frac_under_100m(i,j,k) = max(0, min(1.0, (Z_100 - ztop(i,j)) / (zbot(i,j) - ztop(i,j)))) + enddo ; enddo + if (any(frac_under_100m(:,:,k) > 0)) khi = k + enddo + endif + + do m=1,Reg%ntr ; if (Reg%Tr(m)%registry_diags) then + Tr => Reg%Tr(m) ! A few diagnostics introduce with MARBL driver ! Compute full-depth vertical integral if (Tr%id_zint > 0) then work2d(:,:) = 0.0 do k=1,nz ; do j=js,je ; do i=is,ie - work2d(i,j) = work2d(i,j) + (h_diag(i,j,k)*GV%H_to_m)*tr%t(i,j,k) + work2d(i,j) = work2d(i,j) + dz(i,j,k)*tr%t(i,j,k) enddo ; enddo ; enddo call post_data(Tr%id_zint, work2d, diag) endif @@ -833,7 +874,7 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) if (Tr%id_zint_100m > 0) then work2d(:,:) = 0.0 do k=1,khi ; do j=js,je ; do i=is,ie - work2d(i,j) = work2d(i,j) + frac_under_100m(i,j,k)*((h_diag(i,j,k)*GV%H_to_m)*tr%t(i,j,k)) + work2d(i,j) = work2d(i,j) + frac_under_100m(i,j,k) * dz(i,j,k)*Tr%t(i,j,k) enddo ; enddo ; enddo call post_data(Tr%id_zint_100m, work2d, diag) endif @@ -842,7 +883,7 @@ subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag) if (Tr%id_SURF > 0) call post_data(Tr%id_SURF, Tr%t(:,:,1), diag) endif ; enddo -end subroutine post_tracer_transport_diagnostics +end subroutine post_tracer_integral_diagnostics !> This subroutine writes out chksums for the first ntr registered tracers. subroutine tracer_array_chksum(mesg, Tr, ntr, G)