From ddefcc71993e1b3cadad58612363bf26fdbea081 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 10 Sep 2019 11:31:15 -0400 Subject: [PATCH 1/5] (+) MKS thickness scaling of horizontal avg diag The horizontal_average_diag_field function uses thickness when computing layer averages in some cases. Reproducibility requires that reproducible sums be calculated in MKS units. While the grids have been converted to MKS lengths, the thickeness have not. This patch rescales the thicknesses to volumes, which restores the reproducibility of the dimentionally scaled quantities involving H. --- src/framework/MOM_diag_mediator.F90 | 5 +++-- src/framework/MOM_diag_remap.F90 | 25 +++++++++++++++---------- 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 54f1934abd..674046a750 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -1765,7 +1765,7 @@ subroutine post_xy_average(diag_cs, diag, field) staggered_in_y = diag%axes%is_v_point .or. diag%axes%is_q_point if (diag%axes%is_native) then - call horizontally_average_diag_field(diag_cs%G, diag_cs%h, & + call horizontally_average_diag_field(diag_cs%G, diag_cs%GV, diag_cs%h, & staggered_in_x, staggered_in_y, & diag%axes%is_layer, diag%v_extensive, & diag_cs%missing_value, field, & @@ -1783,7 +1783,8 @@ subroutine post_xy_average(diag_cs, diag, field) call assert(IMPLIES(.not. diag%axes%is_layer, nz == remap_nz+1), & 'post_xy_average: interface field dimension mismatch.') - call horizontally_average_diag_field(diag_cs%G, diag_cs%diag_remap_cs(coord)%h, & + call horizontally_average_diag_field(diag_cs%G, diag_cs%GV, & + diag_cs%diag_remap_cs(coord)%h, & staggered_in_x, staggered_in_y, & diag%axes%is_layer, diag%v_extensive, & diag_cs%missing_value, field, & diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 8f1d309b06..372d6d65cc 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -637,11 +637,12 @@ subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, sta end subroutine vertically_interpolate_diag_field !> Horizontally average field -subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, & +subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_in_y, & is_layer, is_extensive, & missing_value, field, averaged_field, & averaged_mask) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean vertical grid structure real, dimension(:,:,:), intent(in) :: h !< The current thicknesses logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points @@ -663,6 +664,7 @@ subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, ! TODO: These averages could potentially be modified to use the function in ! the MOM_spatial_means module. + ! NOTE: Reproducible sums must be computed in the original MKS units if (staggered_in_x .and. .not. staggered_in_y) then if (is_layer) then @@ -673,14 +675,15 @@ subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, if (is_extensive) then do j=G%jsc, G%jec ; do I=G%isc, G%iec I1 = I - G%isdB + 1 - volume(I,j,k) = G%US%L_to_m**2*G%areaCu(I,j) * G%mask2dCu(I,j) + volume(I,j,k) = (G%US%L_to_m**2 * G%areaCu(I,j)) * G%mask2dCu(I,j) stuff(I,j,k) = volume(I,j,k) * field(I1,j,k) enddo ; enddo else ! Intensive do j=G%jsc, G%jec ; do I=G%isc, G%iec I1 = i - G%isdB + 1 height = 0.5 * (h(i,j,k) + h(i+1,j,k)) - volume(I,j,k) = G%US%L_to_m**2*G%areaCu(I,j) * height * G%mask2dCu(I,j) + volume(I,j,k) = (G%US%L_to_m**2 * G%areaCu(I,j)) & + * (GV%H_to_m * height) * G%mask2dCu(I,j) stuff(I,j,k) = volume(I,j,k) * field(I1,j,k) enddo ; enddo endif @@ -689,7 +692,7 @@ subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, do k=1,nz do j=G%jsc, G%jec ; do I=G%isc, G%iec I1 = I - G%isdB + 1 - volume(I,j,k) = G%US%L_to_m**2*G%areaCu(I,j) * G%mask2dCu(I,j) + volume(I,j,k) = (G%US%L_to_m**2 * G%areaCu(I,j)) * G%mask2dCu(I,j) stuff(I,j,k) = volume(I,j,k) * field(I1,j,k) enddo ; enddo enddo @@ -701,14 +704,15 @@ subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, if (is_extensive) then do J=G%jsc, G%jec ; do i=G%isc, G%iec J1 = J - G%jsdB + 1 - volume(i,J,k) = G%US%L_to_m**2*G%areaCv(i,J) * G%mask2dCv(i,J) + volume(i,J,k) = (G%US%L_to_m**2 * G%areaCv(i,J)) * G%mask2dCv(i,J) stuff(i,J,k) = volume(i,J,k) * field(i,J1,k) enddo ; enddo else ! Intensive do J=G%jsc, G%jec ; do i=G%isc, G%iec J1 = J - G%jsdB + 1 height = 0.5 * (h(i,j,k) + h(i,j+1,k)) - volume(i,J,k) = G%US%L_to_m**2*G%areaCv(i,J) * height * G%mask2dCv(i,J) + volume(i,J,k) = (G%US%L_to_m**2 * G%areaCv(i,J)) & + * (GV%H_to_m * height) * G%mask2dCv(i,J) stuff(i,J,k) = volume(i,J,k) * field(i,J1,k) enddo ; enddo endif @@ -717,7 +721,7 @@ subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, do k=1,nz do J=G%jsc, G%jec ; do i=G%isc, G%iec J1 = J - G%jsdB + 1 - volume(i,J,k) = G%US%L_to_m**2*G%areaCv(i,J) * G%mask2dCv(i,J) + volume(i,J,k) = (G%US%L_to_m**2 * G%areaCv(i,J)) * G%mask2dCv(i,J) stuff(i,J,k) = volume(i,J,k) * field(i,J1,k) enddo ; enddo enddo @@ -729,7 +733,7 @@ subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, if (is_extensive) then do j=G%jsc, G%jec ; do i=G%isc, G%iec if (h(i,j,k) > 0.) then - volume(i,j,k) = G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j) + volume(i,j,k) = (G%US%L_to_m**2 * G%areaT(i,j)) * G%mask2dT(i,j) stuff(i,j,k) = volume(i,j,k) * field(i,j,k) else volume(i,j,k) = 0. @@ -738,7 +742,8 @@ subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, enddo ; enddo else ! Intensive do j=G%jsc, G%jec ; do i=G%isc, G%iec - volume(i,j,k) = G%US%L_to_m**2*G%areaT(i,j) * h(i,j,k) * G%mask2dT(i,j) + volume(i,j,k) = (G%US%L_to_m**2 * G%areaT(i,j)) & + * (GV%H_to_m * h(i,j,k)) * G%mask2dT(i,j) stuff(i,j,k) = volume(i,j,k) * field(i,j,k) enddo ; enddo endif @@ -746,7 +751,7 @@ subroutine horizontally_average_diag_field(G, h, staggered_in_x, staggered_in_y, else ! Interface do k=1,nz do j=G%jsc, G%jec ; do i=G%isc, G%iec - volume(i,j,k) = G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j) + volume(i,j,k) = (G%US%L_to_m**2 * G%areaT(i,j)) * G%mask2dT(i,j) stuff(i,j,k) = volume(i,j,k) * field(i,j,k) enddo ; enddo enddo From af62e2e678df28357a7cc9726347c9b3ddcf9de5 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 10 Sep 2019 15:31:20 -0400 Subject: [PATCH 2/5] Diagnostic dimensional scaling fixes This fixes up the dimensional scaling of several diagnostics: * dynamics_h_tendency * MLD_Restrat * ea, eb * h_predia * frazil_h We also include a change to the flux_to_kg_per_s scaling factor. Previously, this was dynamically set the flux units to either m3 s-1 or kg s-1 depending on the Boussinesq mode. However, the units of the diagnostics using this factor (uhml, vhml) were always explicitly set to kg -1. Assuming that the scaled fluxes were H L2 T-1, this would have given the wrong scaling in non-Boussinesq mode. We resolve this by removing the conditional and always using GV%H_to_kg_m2 when defining this scaling factor. --- src/diagnostics/MOM_diagnostics.F90 | 2 +- .../lateral/MOM_mixed_layer_restrat.F90 | 20 ++++++++--------- .../vertical/MOM_diabatic_driver.F90 | 22 ++++++++++++++----- 3 files changed, 27 insertions(+), 17 deletions(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 4aa6fad710..8fa106c4e0 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1845,7 +1845,7 @@ subroutine register_transport_diags(Time, G, GV, US, IDs, diag) 'm s-1', v_extensive=.true., conversion=GV%H_to_m) IDs%id_dynamics_h_tendency = register_diag_field('ocean_model','dynamics_h_tendency', & diag%axesTl, Time, 'Change in layer thicknesses due to horizontal dynamics', & - 'm s-1', v_extensive=.true.) + 'm s-1', v_extensive=.true., conversion=GV%H_to_m) end subroutine register_transport_diags diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 6543ab7af6..ba241ea4b1 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -884,30 +884,30 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, CS%diag => diag - if (GV%Boussinesq) then ; flux_to_kg_per_s = GV%H_to_kg_m2 * US%L_to_m**2 * US%s_to_T - else ; flux_to_kg_per_s = GV%H_to_m * US%L_to_m**2 * US%s_to_T ; endif + flux_to_kg_per_s = GV%H_to_kg_m2 * US%L_to_m**2 * US%s_to_T CS%id_uhml = register_diag_field('ocean_model', 'uhml', diag%axesCuL, Time, & - 'Zonal Thickness Flux to Restratify Mixed Layer', 'kg s-1', conversion=flux_to_kg_per_s, & - y_cell_method='sum', v_extensive=.true.) + 'Zonal Thickness Flux to Restratify Mixed Layer', 'kg s-1', & + conversion=flux_to_kg_per_s, y_cell_method='sum', v_extensive=.true.) CS%id_vhml = register_diag_field('ocean_model', 'vhml', diag%axesCvL, Time, & - 'Meridional Thickness Flux to Restratify Mixed Layer', 'kg s-1', conversion=flux_to_kg_per_s, & - x_cell_method='sum', v_extensive=.true.) + 'Meridional Thickness Flux to Restratify Mixed Layer', 'kg s-1', & + conversion=flux_to_kg_per_s, x_cell_method='sum', v_extensive=.true.) CS%id_urestrat_time = register_diag_field('ocean_model', 'MLu_restrat_time', diag%axesCu1, Time, & 'Mixed Layer Zonal Restratification Timescale', 's', conversion=US%T_to_s) CS%id_vrestrat_time = register_diag_field('ocean_model', 'MLv_restrat_time', diag%axesCv1, Time, & 'Mixed Layer Meridional Restratification Timescale', 's', conversion=US%T_to_s) CS%id_MLD = register_diag_field('ocean_model', 'MLD_restrat', diag%axesT1, Time, & - 'Mixed Layer Depth as used in the mixed-layer restratification parameterization', 'm') + 'Mixed Layer Depth as used in the mixed-layer restratification parameterization', 'm', & + conversion=GV%H_to_m) CS%id_Rml = register_diag_field('ocean_model', 'ML_buoy_restrat', diag%axesT1, Time, & 'Mixed Layer Buoyancy as used in the mixed-layer restratification parameterization', & - 'm s2', conversion=US%m_to_Z*US%L_to_m**2*US%s_to_T**2) + 'm s2', conversion=US%m_to_Z*(US%L_to_m**2)*(US%s_to_T**2)) CS%id_uDml = register_diag_field('ocean_model', 'udml_restrat', diag%axesCu1, Time, & 'Transport stream function amplitude for zonal restratification of mixed layer', & - 'm3 s-1', conversion=GV%H_to_m*US%L_to_m**2*US%s_to_T) + 'm3 s-1', conversion=GV%H_to_m*(US%L_to_m**2)*US%s_to_T) CS%id_vDml = register_diag_field('ocean_model', 'vdml_restrat', diag%axesCv1, Time, & 'Transport stream function amplitude for meridional restratification of mixed layer', & - 'm3 s-1', conversion=GV%H_to_m*US%L_to_m**2*US%s_to_T) + 'm3 s-1', conversion=GV%H_to_m*(US%L_to_m**2)*US%s_to_T) CS%id_uml = register_diag_field('ocean_model', 'uml_restrat', diag%axesCu1, Time, & 'Surface zonal velocity component of mixed layer restratification', & 'm s-1', conversion=US%L_T_to_m_s) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e6f644d210..2354fcf933 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3203,6 +3203,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di real :: Kd integer :: num_mode logical :: use_temperature, differentialDiffusion + real :: H_to_MKS !< Conversion factor from H to equivalent MKS unit ! This "include" declares and sets the variable "version". #include "version_variable.h" @@ -3361,8 +3362,13 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! Register all available diagnostics for this module. - if (GV%Boussinesq) then ; thickness_units = "m" - else ; thickness_units = "kg m-2" ; endif + if (GV%Boussinesq) then + thickness_units = "m" + H_to_MKS = GV%H_to_m + else + thickness_units = "kg m-2" + H_to_MKS = GV%H_to_kg_m2 + endif CS%id_ea_t = register_diag_field('ocean_model','ea_t',diag%axesTL,Time, & 'Layer (heat) entrainment from above per timestep','m') @@ -3374,9 +3380,11 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Layer (salt) entrainment from below per timestep', 'm') ! used by layer diabatic CS%id_ea = register_diag_field('ocean_model','ea',diag%axesTL,Time, & - 'Layer entrainment from above per timestep','m') + 'Layer entrainment from above per timestep','m', & + conversion=GV%H_to_m) CS%id_eb = register_diag_field('ocean_model','eb',diag%axesTL,Time, & - 'Layer entrainment from below per timestep', 'm') + 'Layer entrainment from below per timestep', 'm', & + conversion=GV%H_to_m) CS%id_wd = register_diag_field('ocean_model','wd',diag%axesTi,Time, & 'Diapycnal velocity', 'm s-1') if (CS%id_wd > 0) call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) @@ -3446,7 +3454,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di CS%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, Time, & 'Meridional velocity before diabatic forcing', 'm s-1', conversion=US%L_T_to_m_s) CS%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, Time, & - 'Layer Thickness before diabatic forcing', thickness_units, v_extensive=.true.) + 'Layer Thickness before diabatic forcing', thickness_units, v_extensive=.true., & + conversion=H_to_MKS) CS%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, Time, & 'Interface Heights before diabatic forcing', 'm') if (use_temperature) then @@ -3640,7 +3649,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! diagnostics for tendencies of temp and heat due to frazil CS%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, Time, & - long_name = 'Cell Thickness', standard_name='cell_thickness', units='m', v_extensive=.true.) + long_name = 'Cell Thickness', standard_name='cell_thickness', units='m', & + v_extensive=.true., conversion=GV%H_to_m) ! diagnostic for tendency of temp due to frazil CS%id_frazil_temp_tend = register_diag_field('ocean_model',& From ff12685c3836a3fc2a5f8ad33a56fbc25d3e2274 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 11 Sep 2019 17:31:00 -0400 Subject: [PATCH 3/5] (+) Dimentional scaling fixes; H_to_MKS scaling This patch fixes the H dimensional scaling of several diagnostics. * h_preale * h_predia * Tflx_dia_diff * Tflx_dia_adv * Sflx_dia_diff * Sflx_dia_adv * diabatic_diff_h * boundary_forcing_h * frazil_h * internal_heat_h_tendency Two interface changes were introduced to complement these fixes. - The vertical grid struct (GV) was added to geothermal_init - We have introduced the H_to_MKS scaling factor, which is set to either H_to_m or H_to_kg_m2 based on whether the model is in Boussinesq mode. This to accommodate diagnostics whose units are based on the Boussinesq mode. The vert_remap_h diagnostic has also been partially rescaled and has been modified by this patch, but is not yet invariant to H scaling. --- src/ALE/MOM_ALE.F90 | 8 ++- src/core/MOM_verticalGrid.F90 | 4 ++ .../vertical/MOM_diabatic_driver.F90 | 63 +++++++++---------- .../vertical/MOM_geothermal.F90 | 10 ++- 4 files changed, 47 insertions(+), 38 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index def4148f25..63c0f027b0 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -250,7 +250,8 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS) CS%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, Time, & 'Meridional velocity before remapping', 'm s-1', conversion=US%L_T_to_m_s) CS%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, Time, & - 'Layer Thickness before remapping', get_thickness_units(GV), v_extensive=.true.) + 'Layer Thickness before remapping', get_thickness_units(GV), & + conversion=GV%H_to_MKS, v_extensive=.true.) CS%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, Time, & 'Temperature before remapping', 'degC') CS%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, Time, & @@ -260,8 +261,9 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS) CS%id_dzRegrid = register_diag_field('ocean_model','dzRegrid',diag%axesTi,Time, & 'Change in interface height due to ALE regridding', 'm') - cs%id_vert_remap_h = register_diag_field('ocean_model','vert_remap_h',diag%axestl,time, & - 'layer thicknesses after ALE regridding and remapping', 'm', v_extensive = .true.) + cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', & + diag%axestl, time, 'layer thicknesses after ALE regridding and remapping', 'm', & + conversion=GV%H_to_m, v_extensive=.true.) cs%id_vert_remap_h_tendency = register_diag_field('ocean_model','vert_remap_h_tendency',diag%axestl,time, & 'Layer thicknesses tendency due to ALE regridding and remapping', 'm', v_extensive = .true.) diff --git a/src/core/MOM_verticalGrid.F90 b/src/core/MOM_verticalGrid.F90 index c11de0d0dd..43c673a592 100644 --- a/src/core/MOM_verticalGrid.F90 +++ b/src/core/MOM_verticalGrid.F90 @@ -61,6 +61,8 @@ module MOM_verticalGrid real :: H_to_Pa !< A constant that translates the units of thickness to pressure [Pa]. real :: H_to_Z !< A constant that translates thickness units to the units of depth. real :: Z_to_H !< A constant that translates depth units to thickness units. + real :: H_to_MKS !< A constant that translates thickness units to its + !! MKS unit (m or kg m-2) based on GV%Boussinesq real :: m_to_H_restart = 0.0 !< A copy of the m_to_H that is used in restart files. end type verticalGrid_type @@ -143,11 +145,13 @@ subroutine verticalGridInit( param_file, GV, US ) GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2 GV%m_to_H = 1.0 / GV%H_to_m GV%Angstrom_H = GV%m_to_H * GV%Angstrom_m + GV%H_to_MKS = GV%H_to_m else GV%kg_m2_to_H = 1.0 / GV%H_to_kg_m2 GV%m_to_H = GV%Rho0 * GV%kg_m2_to_H GV%H_to_m = GV%H_to_kg_m2 / GV%Rho0 GV%Angstrom_H = GV%Angstrom_m*1000.0*GV%kg_m2_to_H + GV%H_to_MKS = GV%H_to_kg_m2 endif GV%H_subroundoff = 1e-20 * max(GV%Angstrom_H,GV%m_to_H*1e-17) GV%H_to_Pa = GV%mks_g_Earth * GV%H_to_kg_m2 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 2354fcf933..d6c38fa93e 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -66,7 +66,7 @@ module MOM_diabatic_driver use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d -use MOM_verticalGrid, only : verticalGrid_type +use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_wave_speed, only : wave_speeds use MOM_wave_interface, only : wave_parameters_CS @@ -508,10 +508,10 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim Kd_heat, & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1] Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1] Kd_ePBL, & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1] - Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC m s-1] - Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC m s-1] - Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt m s-1] - Sadv_flx ! advective diapycnal salt flux across interfaces [ppt m s-1] + Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] + Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] + Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] + Sadv_flx ! advective diapycnal salt flux across interfaces [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below, ! where massive is defined as sufficiently thick that @@ -1291,10 +1291,10 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, Kd_heat, & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1] Kd_salt, & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1] Kd_ePBL, & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1] - Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC m s-1] - Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC m s-1] - Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt m s-1] - Sadv_flx ! advective diapycnal salt flux across interfaces [ppt m s-1] + Tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] + Tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H s-1 ~> degC m s-1 or degC kg m-2 s-1] + Sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] + Sadv_flx ! advective diapycnal salt flux across interfaces [ppt H s-1 ~> ppt m s-1 or ppt kg m-2 s-1] logical :: in_boundary(SZI_(G)) ! True if there are no massive layers below, ! where massive is defined as sufficiently thick that @@ -2705,7 +2705,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e !$OMP parallel do default(shared) do j=js,je do K=2,nz ; do i=is,ie - CDp%diapyc_vel(i,j,K) = Idt * (GV%H_to_m * (ea(i,j,k) - eb(i,j,k-1))) + CDp%diapyc_vel(i,j,K) = Idt * (ea(i,j,k) - eb(i,j,k-1)) enddo ; enddo do i=is,ie CDp%diapyc_vel(i,j,1) = 0.0 @@ -3203,7 +3203,6 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di real :: Kd integer :: num_mode logical :: use_temperature, differentialDiffusion - real :: H_to_MKS !< Conversion factor from H to equivalent MKS unit ! This "include" declares and sets the variable "version". #include "version_variable.h" @@ -3362,22 +3361,20 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! Register all available diagnostics for this module. - if (GV%Boussinesq) then - thickness_units = "m" - H_to_MKS = GV%H_to_m - else - thickness_units = "kg m-2" - H_to_MKS = GV%H_to_kg_m2 - endif + thickness_units = get_thickness_units(GV) CS%id_ea_t = register_diag_field('ocean_model','ea_t',diag%axesTL,Time, & - 'Layer (heat) entrainment from above per timestep','m') + 'Layer (heat) entrainment from above per timestep','m', & + conversion=GV%H_to_m) CS%id_eb_t = register_diag_field('ocean_model','eb_t',diag%axesTL,Time, & - 'Layer (heat) entrainment from below per timestep', 'm') + 'Layer (heat) entrainment from below per timestep', 'm', & + conversion=GV%H_to_m) CS%id_ea_s = register_diag_field('ocean_model','ea_s',diag%axesTL,Time, & - 'Layer (salt) entrainment from above per timestep','m') + 'Layer (salt) entrainment from above per timestep','m', & + conversion=GV%H_to_m) CS%id_eb_s = register_diag_field('ocean_model','eb_s',diag%axesTL,Time, & - 'Layer (salt) entrainment from below per timestep', 'm') + 'Layer (salt) entrainment from below per timestep', 'm', & + conversion=GV%H_to_m) ! used by layer diabatic CS%id_ea = register_diag_field('ocean_model','ea',diag%axesTL,Time, & 'Layer entrainment from above per timestep','m', & @@ -3410,16 +3407,16 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di if (use_temperature) then CS%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff",diag%axesTi, & Time, "Diffusive diapycnal temperature flux across interfaces", & - "degC m s-1") + "degC m s-1", conversion=GV%H_to_m) CS%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv",diag%axesTi, & Time, "Advective diapycnal temperature flux across interfaces", & - "degC m s-1") + "degC m s-1", conversion=GV%H_to_m) CS%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff",diag%axesTi, & Time, "Diffusive diapycnal salnity flux across interfaces", & - "psu m s-1") + "psu m s-1", conversion=GV%H_to_m) CS%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv",diag%axesTi, & Time, "Advective diapycnal salnity flux across interfaces", & - "psu m s-1") + "psu m s-1", conversion=GV%H_to_m) CS%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, Time, & 'Mixed layer depth (delta rho = 0.03)', 'm', conversion=US%Z_to_m, & cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', & @@ -3454,8 +3451,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di CS%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, Time, & 'Meridional velocity before diabatic forcing', 'm s-1', conversion=US%L_T_to_m_s) CS%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, Time, & - 'Layer Thickness before diabatic forcing', thickness_units, v_extensive=.true., & - conversion=H_to_MKS) + 'Layer Thickness before diabatic forcing', trim(thickness_units), & + conversion=GV%H_to_MKS, v_extensive=.true.) CS%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, Time, & 'Interface Heights before diabatic forcing', 'm') if (use_temperature) then @@ -3519,7 +3516,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! available only for ALE algorithm. ! diagnostics for tendencies of temp and heat due to frazil CS%id_diabatic_diff_h = register_diag_field('ocean_model', 'diabatic_diff_h', diag%axesTL, Time, & - long_name = 'Cell thickness used during diabatic diffusion', units='m', v_extensive=.true.) + long_name = 'Cell thickness used during diabatic diffusion', units='m', & + conversion=GV%H_to_m, v_extensive=.true.) if (CS%useALEalgorithm) then CS%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', & 'diabatic_diff_temp_tendency', diag%axesTL, Time, & @@ -3591,7 +3589,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! available only for ALE algorithm. ! diagnostics for tendencies of temp and heat due to frazil CS%id_boundary_forcing_h = register_diag_field('ocean_model', 'boundary_forcing_h', diag%axesTL, Time, & - long_name = 'Cell thickness after applying boundary forcing', units='m', v_extensive=.true.) + long_name = 'Cell thickness after applying boundary forcing', units='m', & + conversion=GV%H_to_m, v_extensive=.true.) CS%id_boundary_forcing_h_tendency = register_diag_field('ocean_model', & 'boundary_forcing_h_tendency', diag%axesTL, Time, & 'Cell thickness tendency due to boundary forcing', 'm s-1', & @@ -3650,7 +3649,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! diagnostics for tendencies of temp and heat due to frazil CS%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, Time, & long_name = 'Cell Thickness', standard_name='cell_thickness', units='m', & - v_extensive=.true., conversion=GV%H_to_m) + conversion=GV%H_to_m, v_extensive=.true.) ! diagnostic for tendency of temp due to frazil CS%id_frazil_temp_tend = register_diag_field('ocean_model',& @@ -3693,7 +3692,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! initialize the geothermal heating module if (CS%use_geothermal) & - call geothermal_init(Time, G, param_file, diag, CS%geothermal_CSp) + call geothermal_init(Time, G, GV, param_file, diag, CS%geothermal_CSp) ! initialize module for internal tide induced mixing if (CS%use_int_tides) then diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90 index d3d0f2bc6e..5fefbf199e 100644 --- a/src/parameterizations/vertical/MOM_geothermal.F90 +++ b/src/parameterizations/vertical/MOM_geothermal.F90 @@ -11,7 +11,7 @@ module MOM_geothermal use MOM_io, only : MOM_read_data, slasher use MOM_grid, only : ocean_grid_type use MOM_variables, only : thermo_var_ptrs -use MOM_verticalGrid, only : verticalGrid_type +use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units use MOM_EOS, only : calculate_density, calculate_density_derivs implicit none ; private @@ -371,9 +371,10 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo) end subroutine geothermal !> Initialize parameters and allocate memory associated with the geothermal heating module. -subroutine geothermal_init(Time, G, param_file, diag, CS) +subroutine geothermal_init(Time, G, GV, param_file, diag, CS) type(time_type), target, intent(in) :: Time !< Current model time. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time !! parameters. type(diag_ctrl), target, intent(inout) :: diag !< Structure used to regulate diagnostic output. @@ -382,6 +383,7 @@ subroutine geothermal_init(Time, G, param_file, diag, CS) ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_geothermal" ! module name + character(len=48) :: thickness_units ! Local variables character(len=200) :: inputdir, geo_file, filename, geotherm_var real :: scale @@ -442,6 +444,8 @@ subroutine geothermal_init(Time, G, param_file, diag, CS) endif call pass_var(CS%geo_heat, G%domain) + thickness_units = get_thickness_units(GV) + ! post the static geothermal heating field id = register_static_field('ocean_model', 'geo_heat', diag%axesT1, & 'Geothermal heat flux into ocean', 'W m-2', & @@ -463,7 +467,7 @@ subroutine geothermal_init(Time, G, param_file, diag, CS) CS%id_internal_heat_h_tendency=register_diag_field('ocean_model', & 'internal_heat_h_tendency', diag%axesTL, Time, & 'Thickness tendency (in 3D) due to internal (geothermal) sources', & - 'm OR kg m-2', v_extensive=.true.) + trim(thickness_units), conversion=GV%H_to_MKS, v_extensive=.true.) end subroutine geothermal_init From f1b866da3b12457a25baf4ef776c43536637ab7c Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 12 Sep 2019 11:50:57 -0400 Subject: [PATCH 4/5] Diagnostic scaling; global_spatial_mean in MKS The dimensional scaling of three diagnostics has been amended: * dzRegrid * vert_remap_h_tendency * wd The spatial averaging function global_spatial_mean has also been modified to run entirely in MKS units, due to the reproducing_sum function requirement. A previous change fixed a bug in this function, where only part of the solution was being scaled. This new fix restores the original scaling, but de-scales the final result in order to retain reproducibility. --- src/ALE/MOM_ALE.F90 | 6 ++++-- src/framework/MOM_spatial_means.F90 | 8 ++++---- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 2 +- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 63c0f027b0..8eed4aa925 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -260,12 +260,14 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS) 'Interface Heights before remapping', 'm', conversion=US%Z_to_m) CS%id_dzRegrid = register_diag_field('ocean_model','dzRegrid',diag%axesTi,Time, & - 'Change in interface height due to ALE regridding', 'm') + 'Change in interface height due to ALE regridding', 'm', & + conversion=GV%H_to_m) cs%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', & diag%axestl, time, 'layer thicknesses after ALE regridding and remapping', 'm', & conversion=GV%H_to_m, v_extensive=.true.) cs%id_vert_remap_h_tendency = register_diag_field('ocean_model','vert_remap_h_tendency',diag%axestl,time, & - 'Layer thicknesses tendency due to ALE regridding and remapping', 'm', v_extensive = .true.) + 'Layer thicknesses tendency due to ALE regridding and remapping', 'm', & + conversion=GV%H_to_m, v_extensive = .true.) end subroutine ALE_register_diags diff --git a/src/framework/MOM_spatial_means.F90 b/src/framework/MOM_spatial_means.F90 index f6282faa52..5a84ca0001 100644 --- a/src/framework/MOM_spatial_means.F90 +++ b/src/framework/MOM_spatial_means.F90 @@ -36,9 +36,9 @@ function global_area_mean(var,G) tmpForSumming(:,:) = 0. do j=js,je ; do i=is,ie - tmpForSumming(i,j) = var(i,j) * (G%areaT(i,j) * G%mask2dT(i,j)) + tmpForSumming(i,j) = var(i,j) * (G%US%L_to_m**2 * G%areaT(i,j) * G%mask2dT(i,j)) enddo ; enddo - global_area_mean = reproducing_sum( tmpForSumming ) * G%IareaT_global + global_area_mean = reproducing_sum(tmpForSumming) * (G%US%m_to_L**2 * G%IareaT_global) end function global_area_mean @@ -54,9 +54,9 @@ function global_area_integral(var,G) tmpForSumming(:,:) = 0. do j=js,je ; do i=is, ie - tmpForSumming(i,j) = ( var(i,j) * (G%US%L_to_m**2*G%areaT(i,j) * G%mask2dT(i,j)) ) + tmpForSumming(i,j) = var(i,j) * (G%US%L_to_m**2 * G%areaT(i,j) * G%mask2dT(i,j)) enddo ; enddo - global_area_integral = reproducing_sum( tmpForSumming ) + global_area_integral = reproducing_sum(tmpForSumming) end function global_area_integral diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index d6c38fa93e..369ee5da40 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3383,7 +3383,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Layer entrainment from below per timestep', 'm', & conversion=GV%H_to_m) CS%id_wd = register_diag_field('ocean_model','wd',diag%axesTi,Time, & - 'Diapycnal velocity', 'm s-1') + 'Diapycnal velocity', 'm s-1', conversion=GV%H_to_m) if (CS%id_wd > 0) call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) CS%id_dudt_dia = register_diag_field('ocean_model','dudt_dia',diag%axesCuL,Time, & From a045278a3b08c1529e4837591e61a358d8eeadda Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 19 Sep 2019 18:03:57 +0000 Subject: [PATCH 5/5] Dimension scaling of massout_flux, tracers This patch fixes the dimensional scaling of the following diagnostics: - massout_flux - ideal_age - various DOME derived tracers - Derived tracer diagnostics: - *_adx - *_ady - *_dfx - *_dfx - *h_tendency_2d Many tracer fixes rely on the "flux_" and "conv_" arguments for rescaling, which is probably not the original intended use, and should probably be revised in the future. --- src/core/MOM_forcing_type.F90 | 2 +- src/tracer/DOME_tracer.F90 | 4 ++-- src/tracer/MOM_tracer_registry.F90 | 20 ++++++++++++-------- src/tracer/ideal_age_example.F90 | 3 ++- 4 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 7df4213a2f..a8c6f7bf1a 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1323,7 +1323,7 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_massout_flux = register_diag_field('ocean_model', 'massout_flux', diag%axesT1, Time, & 'Net mass flux of freshwater out of the ocean (used in the boundary flux calculation)', & - 'kg m-2') + 'kg m-2', conversion=diag%GV%H_to_kg_m2) handles%id_massin_flux = register_diag_field('ocean_model', 'massin_flux', diag%axesT1, Time, & 'Net mass flux of freshwater into the ocean (used in boundary flux calculation)', 'kg m-2') diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90 index d820ecf36a..7589f04ed0 100644 --- a/src/tracer/DOME_tracer.F90 +++ b/src/tracer/DOME_tracer.F90 @@ -123,8 +123,8 @@ function register_DOME_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Register the tracer for horizontal advection, diffusion, and restarts. call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & name=name, longname=longname, units="kg kg-1", & - registry_diags=.true., flux_units=flux_units, & - restart_CS=restart_CS) + registry_diags=.true., restart_CS=restart_CS, & + flux_units=trim(flux_units), flux_scale=GV%H_to_MKS) ! Set coupled_tracers to be true (hard-coded above) to provide the surface ! values to the coupler (if any). This is meta-code and its arguments will diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 41299be3e8..6a2dd79b5b 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -392,18 +392,20 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE) if (Tr%diag_form == 1) then Tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", & diag%axesCuL, Time, trim(flux_longname)//" advective zonal flux" , & - trim(flux_units), v_extensive = .true., y_cell_method = 'sum') + trim(flux_units), v_extensive = .true., y_cell_method = 'sum', & + conversion=Tr%flux_scale) Tr%id_ady = register_diag_field("ocean_model", trim(shortnm)//"_ady", & diag%axesCvL, Time, trim(flux_longname)//" advective meridional flux" , & - trim(flux_units), v_extensive = .true., x_cell_method = 'sum') + trim(flux_units), v_extensive = .true., x_cell_method = 'sum', & + conversion=Tr%flux_scale) Tr%id_dfx = register_diag_field("ocean_model", trim(shortnm)//"_dfx", & diag%axesCuL, Time, trim(flux_longname)//" diffusive zonal flux" , & - trim(flux_units), v_extensive = .true., conversion=US%L_to_m**2, & - y_cell_method = 'sum') + trim(flux_units), v_extensive = .true., y_cell_method = 'sum', & + conversion=(US%L_to_m**2)*Tr%flux_scale) Tr%id_dfy = register_diag_field("ocean_model", trim(shortnm)//"_dfy", & diag%axesCvL, Time, trim(flux_longname)//" diffusive zonal flux" , & - trim(flux_units), v_extensive = .true., conversion=US%L_to_m**2, & - x_cell_method = 'sum') + trim(flux_units), v_extensive = .true., x_cell_method = 'sum', & + conversion=(US%L_to_m**2)*Tr%flux_scale) 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), & @@ -508,9 +510,11 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE) var_lname = "Net time tendency for "//lowercase(flux_longname) if (len_trim(Tr%cmor_tendprefix) == 0) then Tr%id_trxh_tendency = register_diag_field('ocean_model', trim(shortnm)//'h_tendency', & - diag%axesTL, Time, var_lname, conv_units, v_extensive=.true.) + diag%axesTL, Time, var_lname, conv_units, v_extensive=.true., & + conversion=Tr%conv_scale) Tr%id_trxh_tendency_2d = register_diag_field('ocean_model', trim(shortnm)//'h_tendency_2d', & - diag%axesT1, Time, "Vertical sum of "//trim(lowercase(var_lname)), conv_units) + diag%axesT1, Time, "Vertical sum of "//trim(lowercase(var_lname)), conv_units, & + conversion=Tr%conv_scale) else cmor_var_lname = "Tendency of "//trim(cmor_longname)//" Expressed as "//& trim(flux_longname)//" Content" diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90 index 35975bccb0..a46e42f415 100644 --- a/src/tracer/ideal_age_example.F90 +++ b/src/tracer/ideal_age_example.F90 @@ -176,7 +176,8 @@ function register_ideal_age_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Register the tracer for horizontal advection, diffusion, and restarts. call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, tr_desc=CS%tr_desc(m), & registry_diags=.true., restart_CS=restart_CS, & - mandatory=.not.CS%tracers_may_reinit) + mandatory=.not.CS%tracers_may_reinit, & + flux_scale=GV%H_to_m) ! Set coupled_tracers to be true (hard-coded above) to provide the surface ! values to the coupler (if any). This is meta-code and its arguments will