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)