From dfc2d0d80321d396f2a2e40baddc398d32350df9 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 30 Aug 2019 17:20:47 -0400 Subject: [PATCH 1/9] Diagnostic dimensional scaling fixes This patch adds and fixes the scaling factors for several diagnostics. --- src/core/MOM_barotropic.F90 | 35 +++++++++++-------- src/diagnostics/MOM_diagnostics.F90 | 31 ++++++++-------- .../vertical/MOM_bulk_mixed_layer.F90 | 21 +++++------ .../vertical/MOM_vert_friction.F90 | 12 ++++--- 4 files changed, 55 insertions(+), 44 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 8d48ebbb0b..347b5a1936 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -4205,13 +4205,13 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, 'Barotropic meridional acceleration from baroclinic terms', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_eta_bt = register_diag_field('ocean_model', 'eta_bt', diag%axesT1, Time, & - 'Barotropic end SSH', thickness_units) + 'Barotropic end SSH', thickness_units, conversion=GV%H_to_m) CS%id_ubt = register_diag_field('ocean_model', 'ubt', diag%axesCu1, Time, & 'Barotropic end zonal velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_vbt = register_diag_field('ocean_model', 'vbt', diag%axesCv1, Time, & 'Barotropic end meridional velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_eta_st = register_diag_field('ocean_model', 'eta_st', diag%axesT1, Time, & - 'Barotropic start SSH', thickness_units) + 'Barotropic start SSH', thickness_units, conversion=GV%H_to_m) CS%id_ubt_st = register_diag_field('ocean_model', 'ubt_st', diag%axesCu1, Time, & 'Barotropic start zonal velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_vbt_st = register_diag_field('ocean_model', 'vbt_st', diag%axesCv1, Time, & @@ -4221,31 +4221,34 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, CS%id_vbtav = register_diag_field('ocean_model', 'vbtav', diag%axesCv1, Time, & 'Barotropic time-average meridional velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_eta_cor = register_diag_field('ocean_model', 'eta_cor', diag%axesT1, Time, & - 'Corrective mass flux', 'm s-1') + 'Corrective mass flux', 'm s-1', conversion=GV%H_to_m) CS%id_visc_rem_u = register_diag_field('ocean_model', 'visc_rem_u', diag%axesCuL, Time, & 'Viscous remnant at u', 'nondim') CS%id_visc_rem_v = register_diag_field('ocean_model', 'visc_rem_v', diag%axesCvL, Time, & 'Viscous remnant at v', 'nondim') CS%id_gtotn = register_diag_field('ocean_model', 'gtot_n', diag%axesT1, Time, & - 'gtot to North', 'm s-2', conversion=US%L_T_to_m_s**2) + 'gtot to North', 'm s-2', conversion=GV%m_to_H*(US%L_T_to_m_s**2)) CS%id_gtots = register_diag_field('ocean_model', 'gtot_s', diag%axesT1, Time, & - 'gtot to South', 'm s-2', conversion=US%L_T_to_m_s**2) + 'gtot to South', 'm s-2', conversion=GV%m_to_H*(US%L_T_to_m_s**2)) CS%id_gtote = register_diag_field('ocean_model', 'gtot_e', diag%axesT1, Time, & - 'gtot to East', 'm s-2', conversion=US%L_T_to_m_s**2) + 'gtot to East', 'm s-2', conversion=GV%m_to_H*(US%L_T_to_m_s**2)) CS%id_gtotw = register_diag_field('ocean_model', 'gtot_w', diag%axesT1, Time, & - 'gtot to West', 'm s-2', conversion=US%L_T_to_m_s**2) + 'gtot to West', 'm s-2', conversion=GV%m_to_H*(US%L_T_to_m_s**2)) CS%id_eta_hifreq = register_diag_field('ocean_model', 'eta_hifreq', diag%axesT1, Time, & - 'High Frequency Barotropic SSH', thickness_units) + 'High Frequency Barotropic SSH', thickness_units, conversion=GV%H_to_m) CS%id_ubt_hifreq = register_diag_field('ocean_model', 'ubt_hifreq', diag%axesCu1, Time, & 'High Frequency Barotropic zonal velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_vbt_hifreq = register_diag_field('ocean_model', 'vbt_hifreq', diag%axesCv1, Time, & 'High Frequency Barotropic meridional velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_eta_pred_hifreq = register_diag_field('ocean_model', 'eta_pred_hifreq', diag%axesT1, Time, & - 'High Frequency Predictor Barotropic SSH', thickness_units) + 'High Frequency Predictor Barotropic SSH', thickness_units, & + conversion=GV%H_to_m) CS%id_uhbt_hifreq = register_diag_field('ocean_model', 'uhbt_hifreq', diag%axesCu1, Time, & - 'High Frequency Barotropic zonal transport', 'm3 s-1') + 'High Frequency Barotropic zonal transport', 'm3 s-1', & + conversion=GV%H_to_m*US%L_T_to_m_s) CS%id_vhbt_hifreq = register_diag_field('ocean_model', 'vhbt_hifreq', diag%axesCv1, Time, & - 'High Frequency Barotropic meridional transport', 'm3 s-1') + 'High Frequency Barotropic meridional transport', 'm3 s-1', & + conversion=GV%H_to_m*US%L_T_to_m_s) CS%id_frhatu = register_diag_field('ocean_model', 'frhatu', diag%axesCuL, Time, & 'Fractional thickness of layers in u-columns', 'nondim') CS%id_frhatv = register_diag_field('ocean_model', 'frhatv', diag%axesCvL, Time, & @@ -4255,9 +4258,11 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, CS%id_frhatv1 = register_diag_field('ocean_model', 'frhatv1', diag%axesCvL, Time, & 'Predictor Fractional thickness of layers in v-columns', 'nondim') CS%id_uhbt = register_diag_field('ocean_model', 'uhbt', diag%axesCu1, Time, & - 'Barotropic zonal transport averaged over a baroclinic step', 'm3 s-1') + 'Barotropic zonal transport averaged over a baroclinic step', 'm3 s-1', & + conversion=GV%H_to_m*US%L_T_to_m_s) CS%id_vhbt = register_diag_field('ocean_model', 'vhbt', diag%axesCv1, Time, & - 'Barotropic meridional transport averaged over a baroclinic step', 'm3 s-1') + 'Barotropic meridional transport averaged over a baroclinic step', 'm3 s-1', & + conversion=GV%H_to_m*US%L_T_to_m_s) if (use_BT_cont_type) then CS%id_BTC_FA_u_EE = register_diag_field('ocean_model', 'BTC_FA_u_EE', diag%axesCu1, Time, & @@ -4269,9 +4274,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, CS%id_BTC_FA_u_W0 = register_diag_field('ocean_model', 'BTC_FA_u_W0', diag%axesCu1, Time, & 'BTCont type near west face area', 'm2', conversion=US%L_to_m*GV%H_to_m) CS%id_BTC_ubt_EE = register_diag_field('ocean_model', 'BTC_ubt_EE', diag%axesCu1, Time, & - 'BTCont type far east velocity', 'm s-1') + 'BTCont type far east velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_BTC_ubt_WW = register_diag_field('ocean_model', 'BTC_ubt_WW', diag%axesCu1, Time, & - 'BTCont type far west velocity', 'm s-1') + 'BTCont type far west velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_BTC_FA_v_NN = register_diag_field('ocean_model', 'BTC_FA_v_NN', diag%axesCv1, Time, & 'BTCont type far north face area', 'm2', conversion=US%L_to_m*GV%H_to_m) CS%id_BTC_FA_v_N0 = register_diag_field('ocean_model', 'BTC_FA_v_N0', diag%axesCv1, Time, & diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 54025a0ac0..f8f3ba2a56 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1495,7 +1495,8 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag CS%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, Time, & long_name = 'Cell Thickness', standard_name='cell_thickness', units='m', v_extensive=.true.) CS%id_h_pre_sync = register_diag_field('ocean_model', 'h_pre_sync', diag%axesTL, Time, & - long_name = 'Cell thickness from the previous timestep', units='m', v_extensive=.true.) + long_name = 'Cell thickness from the previous timestep', units='m', & + v_extensive=.true., conversion=GV%H_to_m) ! Note that CS%id_volcello would normally be registered here but because it is a "cell measure" and ! must be registered first. We earlier stored the handle of volcello but need it here for posting @@ -1832,10 +1833,10 @@ subroutine register_transport_diags(Time, G, GV, US, IDs, diag) standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum') IDs%id_dynamics_h = register_diag_field('ocean_model','dynamics_h', & 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) 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.) end subroutine register_transport_diags @@ -1882,34 +1883,34 @@ subroutine write_static_fields(G, GV, US, tv, diag) 'Longitude of zonal velocity (Cu) points', 'degrees_east', interp_method='none') if (id > 0) call post_data(id, G%geoLonCu, diag, .true.) - id = register_static_field('ocean_model', 'area_t', diag%axesT1, & - 'Surface area of tracer (T) cells', 'm2', conversion=US%m_to_L**2, & + id = register_static_field('ocean_model', 'area_t', diag%axesT1, & + 'Surface area of tracer (T) cells', 'm2', conversion=US%L_to_m**2, & cmor_field_name='areacello', cmor_standard_name='cell_area', & - cmor_long_name='Ocean Grid-Cell Area', & + cmor_long_name='Ocean Grid-Cell Area', & x_cell_method='sum', y_cell_method='sum', area_cell_method='sum') if (id > 0) then call post_data(id, G%areaT, diag, .true.) call diag_register_area_ids(diag, id_area_t=id) endif - id = register_static_field('ocean_model', 'area_u', diag%axesCu1, & - 'Surface area of x-direction flow (U) cells', 'm2', conversion=US%m_to_L**2, & + id = register_static_field('ocean_model', 'area_u', diag%axesCu1, & + 'Surface area of x-direction flow (U) cells', 'm2', conversion=US%L_to_m**2, & cmor_field_name='areacello_cu', cmor_standard_name='cell_area', & - cmor_long_name='Ocean Grid-Cell Area', & + cmor_long_name='Ocean Grid-Cell Area', & x_cell_method='sum', y_cell_method='sum', area_cell_method='sum') if (id > 0) call post_data(id, G%areaCu, diag, .true.) - id = register_static_field('ocean_model', 'area_v', diag%axesCv1, & - 'Surface area of y-direction flow (V) cells', 'm2', conversion=US%m_to_L**2, & + id = register_static_field('ocean_model', 'area_v', diag%axesCv1, & + 'Surface area of y-direction flow (V) cells', 'm2', conversion=US%L_to_m**2, & cmor_field_name='areacello_cv', cmor_standard_name='cell_area', & - cmor_long_name='Ocean Grid-Cell Area', & + cmor_long_name='Ocean Grid-Cell Area', & x_cell_method='sum', y_cell_method='sum', area_cell_method='sum') if (id > 0) call post_data(id, G%areaCv, diag, .true.) - id = register_static_field('ocean_model', 'area_q', diag%axesB1, & - 'Surface area of B-grid flow (Q) cells', 'm2', conversion=US%m_to_L**2, & + id = register_static_field('ocean_model', 'area_q', diag%axesB1, & + 'Surface area of B-grid flow (Q) cells', 'm2', conversion=US%L_to_m**2, & cmor_field_name='areacello_bu', cmor_standard_name='cell_area', & - cmor_long_name='Ocean Grid-Cell Area', & + cmor_long_name='Ocean Grid-Cell Area', & x_cell_method='sum', y_cell_method='sum', area_cell_method='sum') if (id > 0) call post_data(id, G%areaBu, diag, .true.) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index cbf42d2b8b..2a17bfbd6f 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -3582,33 +3582,34 @@ subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS) Time, 'Surface mixed layer depth', 'm') CS%id_TKE_wind = register_diag_field('ocean_model', 'TKE_wind', diag%axesT1, & Time, 'Wind-stirring source of mixed layer TKE', & - 'm3 s-3', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'm3 s-3', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_TKE_RiBulk = register_diag_field('ocean_model', 'TKE_RiBulk', diag%axesT1, & Time, 'Mean kinetic energy source of mixed layer TKE', & - 'm3 s-3', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'm3 s-3', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_TKE_conv = register_diag_field('ocean_model', 'TKE_conv', diag%axesT1, & - Time, 'Convective source of mixed layer TKE', 'm3 s-3', conversion=US%Z_to_m*US%T_to_s**3) + Time, 'Convective source of mixed layer TKE', 'm3 s-3', & + conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_TKE_pen_SW = register_diag_field('ocean_model', 'TKE_pen_SW', diag%axesT1, & Time, 'TKE consumed by mixing penetrative shortwave radation through the mixed layer', & - 'm3 s-3', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'm3 s-3', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_TKE_mixing = register_diag_field('ocean_model', 'TKE_mixing', diag%axesT1, & Time, 'TKE consumed by mixing that deepens the mixed layer', & - 'm3 s-3', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'm3 s-3', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_TKE_mech_decay = register_diag_field('ocean_model', 'TKE_mech_decay', diag%axesT1, & Time, 'Mechanical energy decay sink of mixed layer TKE', & - 'm3 s-3', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'm3 s-3', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_TKE_conv_decay = register_diag_field('ocean_model', 'TKE_conv_decay', diag%axesT1, & Time, 'Convective energy decay sink of mixed layer TKE', & - 'm3 s-3', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'm3 s-3', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_TKE_conv_s2 = register_diag_field('ocean_model', 'TKE_conv_s2', diag%axesT1, & Time, 'Spurious source of mixed layer TKE from sigma2', & - 'm3 s-3', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'm3 s-3', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_PE_detrain = register_diag_field('ocean_model', 'PE_detrain', diag%axesT1, & Time, 'Spurious source of potential energy from mixed layer detrainment', & - 'W m-2', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'W m-2', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_PE_detrain2 = register_diag_field('ocean_model', 'PE_detrain2', diag%axesT1, & Time, 'Spurious source of potential energy from mixed layer only detrainment', & - 'W m-2', conversion=US%Z_to_m*US%L_to_m**2*US%T_to_s**3) + 'W m-2', conversion=US%Z_to_m*(US%L_to_m**2)*(US%s_to_T**3)) CS%id_h_mismatch = register_diag_field('ocean_model', 'h_miss_ML', diag%axesT1, & Time, 'Summed absolute mismatch in entrainment terms', 'm', conversion=US%Z_to_m) CS%id_Hsfc_used = register_diag_field('ocean_model', 'Hs_used', diag%axesT1, & diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 1bed36e75e..b282995d3f 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1712,16 +1712,20 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & 'Meridional Viscous Vertical Coupling Coefficient', 'm s-1', conversion=US%Z_to_m*US%s_to_T) CS%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, Time, & - 'Thickness at Zonal Velocity Points for Viscosity', thickness_units) + 'Thickness at Zonal Velocity Points for Viscosity', thickness_units, & + conversion=GV%H_to_m) CS%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, Time, & - 'Thickness at Meridional Velocity Points for Viscosity', thickness_units) + 'Thickness at Meridional Velocity Points for Viscosity', thickness_units, & + conversion=GV%H_to_m) CS%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, Time, & - 'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units) + 'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units, & + conversion=GV%H_to_m) CS%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, Time, & - 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units) + 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units, & + conversion=GV%H_to_m) CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, & Time, 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) From 033bca26dbbb3b79dcabf5a72b755a5d2c252409 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 3 Sep 2019 17:35:20 -0400 Subject: [PATCH 2/9] Correct tendency scaling; update KE coversions This patch fixes an issue in the time scalings of the tendencies in MOM_diagnostics, where the timestep (dt) was not being scaled, and was only corrected in some of the variables. This patch also passes through the scaling of KE and dKE_dt onto the diagnostic conversion factors. --- src/diagnostics/MOM_diagnostics.F90 | 32 +++++++++++++++-------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index f8f3ba2a56..9fa7804b26 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -92,8 +92,8 @@ module MOM_diagnostics ! The following arrays hold diagnostics in the layer-integrated energy budget. real, pointer, dimension(:,:,:) :: & - KE => NULL(), & !< KE per unit mass [m2 s-2] - dKE_dt => NULL(), & !< time derivative of the layer KE [m3 s-3] + KE => NULL(), & !< KE per unit mass [L2 T-2 ~> m2 s-2] + dKE_dt => NULL(), & !< time derivative of the layer KE [H L2 T-3 ~> m3 s-3] PE_to_KE => NULL(), & !< potential energy to KE term [m3 s-3] KE_CorAdv => NULL(), & !< KE source from the combined Coriolis and advection terms [m3 s-3]. !! The Coriolis source should be zero, but is not due to truncation @@ -254,7 +254,7 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (loc(CS)==0) call MOM_error(FATAL, & "calculate_diagnostic_fields: Module must be initialized before used.") - call calculate_derivs(dt, G, CS) + call calculate_derivs(US%s_to_T*dt, G, CS) if (dt > 0.0) then call diag_save_grids(CS%diag) @@ -916,8 +916,8 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (associated(CS%KE)) then do k=1,nz ; do j=js,je ; do i=is,ie - CS%KE(i,j,k) = US%L_T_to_m_s**2*((u(I,j,k)*u(I,j,k) + u(I-1,j,k)*u(I-1,j,k)) + & - (v(i,J,k)*v(i,J,k) + v(i,J-1,k)*v(i,J-1,k)))*0.25 + CS%KE(i,j,k) = ((u(I,j,k) * u(I,j,k) + u(I-1,j,k) * u(I-1,j,k)) & + + (v(i,J,k) * v(i,J,k) + v(i,J-1,k) * v(i,J-1,k))) * 0.25 ! DELETE THE FOLLOWING... Make this 0 to test the momentum balance, ! or a huge number to test the continuity balance. ! CS%KE(i,j,k) *= 1e20 @@ -936,19 +936,19 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (associated(CS%dKE_dt)) then do k=1,nz do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = US%L_T_to_m_s**2*US%s_to_T*uh(I,j,k)*G%dxCu(I,j)*CS%du_dt(I,j,k) + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * CS%du_dt(I,j,k) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*CS%dv_dt(i,J,k) + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * CS%dv_dt(i,J,k) enddo ; enddo do j=js,je ; do i=is,ie - KE_h(i,j) = CS%KE(i,j,k)*US%s_to_T*CS%dh_dt(i,j,k) + KE_h(i,j) = CS%KE(i,j,k) * CS%dh_dt(i,j,k) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie - CS%dKE_dt(i,j,k) = GV%H_to_m * (KE_h(i,j) + 0.5 * G%IareaT(i,j) * & - (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))) + CS%dKE_dt(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) enddo ; enddo enddo if (CS%id_dKEdt > 0) call post_data(CS%id_dKEdt, CS%dKE_dt, CS%diag) @@ -981,7 +981,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*ADp%CAv(i,J,k) enddo ; enddo do j=js,je ; do i=is,ie - KE_h(i,j) = -CS%KE(i,j,k) * G%IareaT(i,j) * & + KE_h(i,j) = -(US%L_T_to_m_s**2 * CS%KE(i,j,k)) * G%IareaT(i,j) * & US%s_to_T*(uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) enddo ; enddo if (.not.G%symmetric) & @@ -1009,7 +1009,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*ADp%gradKEv(i,J,k) enddo ; enddo do j=js,je ; do i=is,ie - KE_h(i,j) = -CS%KE(i,j,k) * G%IareaT(i,j) * & + KE_h(i,j) = -(US%L_T_to_m_s**2 * CS%KE(i,j,k)) * G%IareaT(i,j) * & US%s_to_T*(uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) enddo ; enddo if (.not.G%symmetric) & @@ -1067,7 +1067,7 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*ADp%dv_dt_dia(i,J,k) enddo ; enddo do j=js,je ; do i=is,ie - KE_h(i,j) = CS%KE(i,j,k) * & + KE_h(i,j) = (US%L_T_to_m_s**2 * CS%KE(i,j,k)) * & (CDp%diapyc_vel(i,j,k) - CDp%diapyc_vel(i,j,k+1)) enddo ; enddo if (.not.G%symmetric) & @@ -1626,11 +1626,13 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag ! terms in the kinetic energy budget CS%id_KE = register_diag_field('ocean_model', 'KE', diag%axesTL, Time, & - 'Layer kinetic energy per unit mass', 'm2 s-2') + 'Layer kinetic energy per unit mass', 'm2 s-2', & + conversion=US%L_T_to_m_s**2) if (CS%id_KE>0) call safe_alloc_ptr(CS%KE,isd,ied,jsd,jed,nz) CS%id_dKEdt = register_diag_field('ocean_model', 'dKE_dt', diag%axesTL, Time, & - 'Kinetic Energy Tendency of Layer', 'm3 s-3') + 'Kinetic Energy Tendency of Layer', 'm3 s-3', & + conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_dKEdt>0) call safe_alloc_ptr(CS%dKE_dt,isd,ied,jsd,jed,nz) CS%id_PE_to_KE = register_diag_field('ocean_model', 'PE_to_KE', diag%axesTL, Time, & From 8965e7111eea675d28b3afc9d381ea84e3f7631e Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 4 Sep 2019 10:59:52 -0400 Subject: [PATCH 3/9] More KE tendency dimensional scaling Finished the rescaling of the layer KE tendency diagnostics, which are all now passing the diagnostic checksum scaling tests. --- src/diagnostics/MOM_diagnostics.F90 | 89 ++++++++++++++++------------- 1 file changed, 48 insertions(+), 41 deletions(-) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 9fa7804b26..4aa6fad710 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -95,14 +95,15 @@ module MOM_diagnostics KE => NULL(), & !< KE per unit mass [L2 T-2 ~> m2 s-2] dKE_dt => NULL(), & !< time derivative of the layer KE [H L2 T-3 ~> m3 s-3] PE_to_KE => NULL(), & !< potential energy to KE term [m3 s-3] - KE_CorAdv => NULL(), & !< KE source from the combined Coriolis and advection terms [m3 s-3]. + KE_CorAdv => NULL(), & !< KE source from the combined Coriolis and + !! advection terms [H L2 T-3 ~> m3 s-3]. !! The Coriolis source should be zero, but is not due to truncation !! errors. There should be near-cancellation of the global integral !! of this spurious Coriolis source. - KE_adv => NULL(), & !< KE source from along-layer advection [m3 s-3] - KE_visc => NULL(), & !< KE source from vertical viscosity [m3 s-3] - KE_horvisc => NULL(), & !< KE source from horizontal viscosity [m3 s-3] - KE_dia => NULL() !< KE source from diapycnal diffusion [m3 s-3] + KE_adv => NULL(), & !< KE source from along-layer advection [H L2 T-3 ~> m3 s-3] + KE_visc => NULL(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3] + KE_horvisc => NULL(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3] + KE_dia => NULL() !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3] !>@{ Diagnostic IDs integer :: id_u = -1, id_v = -1, id_h = -1 @@ -957,16 +958,16 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (associated(CS%PE_to_KE)) then do k=1,nz do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = US%L_T_to_m_s**2*US%s_to_T*uh(I,j,k)*G%dxCu(I,j)*ADp%PFu(I,j,k) + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%PFu(I,j,k) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*ADp%PFv(i,J,k) + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%PFv(i,J,k) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie - CS%PE_to_KE(i,j,k) = GV%H_to_m * 0.5 * G%IareaT(i,j) * & - (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + CS%PE_to_KE(i,j,k) = 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) enddo ; enddo enddo if (CS%id_PE_to_KE > 0) call post_data(CS%id_PE_to_KE, CS%PE_to_KE, CS%diag) @@ -975,20 +976,20 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (associated(CS%KE_CorAdv)) then do k=1,nz do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = US%L_T_to_m_s**2*US%s_to_T*uh(I,j,k)*G%dxCu(I,j)*ADp%CAu(I,j,k) + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%CAu(I,j,k) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*ADp%CAv(i,J,k) + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%CAv(i,J,k) enddo ; enddo do j=js,je ; do i=is,ie - KE_h(i,j) = -(US%L_T_to_m_s**2 * CS%KE(i,j,k)) * G%IareaT(i,j) * & - US%s_to_T*(uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) + KE_h(i,j) = -CS%KE(i,j,k) * G%IareaT(i,j) & + * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie - CS%KE_CorAdv(i,j,k) = GV%H_to_m * (KE_h(i,j) + 0.5 * G%IareaT(i,j) * & - (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))) + CS%KE_CorAdv(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) enddo ; enddo enddo if (CS%id_KE_Coradv > 0) call post_data(CS%id_KE_Coradv, CS%KE_Coradv, CS%diag) @@ -1002,21 +1003,21 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS do k=1,nz do j=js,je ; do I=Isq,Ieq if (G%mask2dCu(i,j) /= 0.) & - KE_u(I,j) = US%L_T_to_m_s**2*US%s_to_T*uh(I,j,k)*G%dxCu(I,j)*ADp%gradKEu(I,j,k) + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%gradKEu(I,j,k) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie if (G%mask2dCv(i,j) /= 0.) & - KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*ADp%gradKEv(i,J,k) + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%gradKEv(i,J,k) enddo ; enddo do j=js,je ; do i=is,ie - KE_h(i,j) = -(US%L_T_to_m_s**2 * CS%KE(i,j,k)) * G%IareaT(i,j) * & - US%s_to_T*(uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) + KE_h(i,j) = -CS%KE(i,j,k) * G%IareaT(i,j) & + * (uh(I,j,k) - uh(I-1,j,k) + vh(i,J,k) - vh(i,J-1,k)) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie - CS%KE_adv(i,j,k) = GV%H_to_m * (KE_h(i,j) + 0.5 * G%IareaT(i,j) * & - (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))) + CS%KE_adv(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) enddo ; enddo enddo if (CS%id_KE_adv > 0) call post_data(CS%id_KE_adv, CS%KE_adv, CS%diag) @@ -1025,16 +1026,16 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (associated(CS%KE_visc)) then do k=1,nz do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = US%L_T_to_m_s**2*US%s_to_T*uh(I,j,k)*G%dxCu(I,j)*ADp%du_dt_visc(I,j,k) + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_visc(I,j,k) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*ADp%dv_dt_visc(i,J,k) + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_visc(i,J,k) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie - CS%KE_visc(i,j,k) = GV%H_to_m * (0.5 * G%IareaT(i,j) * & - (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1))) + CS%KE_visc(i,j,k) = 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) enddo ; enddo enddo if (CS%id_KE_visc > 0) call post_data(CS%id_KE_visc, CS%KE_visc, CS%diag) @@ -1043,16 +1044,16 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (associated(CS%KE_horvisc)) then do k=1,nz do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = US%s_to_T*uh(I,j,k)*US%L_to_m*G%dxCu(I,j)*US%L_T2_to_m_s2*ADp%diffu(I,j,k) + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%diffu(I,j,k) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = US%s_to_T*vh(i,J,k)*US%L_to_m*G%dyCv(i,J)*US%L_T2_to_m_s2*ADp%diffv(i,J,k) + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%diffv(i,J,k) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie - CS%KE_horvisc(i,j,k) = GV%H_to_m * 0.5 * G%IareaT(i,j) * & - (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + CS%KE_horvisc(i,j,k) = 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) enddo ; enddo enddo if (CS%id_KE_horvisc > 0) call post_data(CS%id_KE_horvisc, CS%KE_horvisc, CS%diag) @@ -1061,20 +1062,20 @@ subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS if (associated(CS%KE_dia)) then do k=1,nz do j=js,je ; do I=Isq,Ieq - KE_u(I,j) = US%L_T_to_m_s**2*US%s_to_T*uh(I,j,k)*G%dxCu(I,j)*ADp%du_dt_dia(I,j,k) + KE_u(I,j) = uh(I,j,k) * G%dxCu(I,j) * ADp%du_dt_dia(I,j,k) enddo ; enddo do J=Jsq,Jeq ; do i=is,ie - KE_v(i,J) = US%L_T_to_m_s**2*US%s_to_T*vh(i,J,k)*G%dyCv(i,J)*ADp%dv_dt_dia(i,J,k) + KE_v(i,J) = vh(i,J,k) * G%dyCv(i,J) * ADp%dv_dt_dia(i,J,k) enddo ; enddo do j=js,je ; do i=is,ie - KE_h(i,j) = (US%L_T_to_m_s**2 * CS%KE(i,j,k)) * & - (CDp%diapyc_vel(i,j,k) - CDp%diapyc_vel(i,j,k+1)) + KE_h(i,j) = CS%KE(i,j,k) & + * (US%T_to_s * (CDp%diapyc_vel(i,j,k) - CDp%diapyc_vel(i,j,k+1))) enddo ; enddo if (.not.G%symmetric) & call do_group_pass(CS%pass_KE_uv, G%domain) do j=js,je ; do i=is,ie - CS%KE_dia(i,j,k) = KE_h(i,j) + GV%H_to_m * 0.5 * G%IareaT(i,j) * & - (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) + CS%KE_dia(i,j,k) = KE_h(i,j) + 0.5 * G%IareaT(i,j) & + * (KE_u(I,j) + KE_u(I-1,j) + KE_v(i,J) + KE_v(i,J-1)) enddo ; enddo enddo if (CS%id_KE_dia > 0) call post_data(CS%id_KE_dia, CS%KE_dia, CS%diag) @@ -1636,28 +1637,34 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag if (CS%id_dKEdt>0) call safe_alloc_ptr(CS%dKE_dt,isd,ied,jsd,jed,nz) CS%id_PE_to_KE = register_diag_field('ocean_model', 'PE_to_KE', diag%axesTL, Time, & - 'Potential to Kinetic Energy Conversion of Layer', 'm3 s-3') + 'Potential to Kinetic Energy Conversion of Layer', 'm3 s-3', & + conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_PE_to_KE>0) call safe_alloc_ptr(CS%PE_to_KE,isd,ied,jsd,jed,nz) CS%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, Time, & - 'Kinetic Energy Source from Coriolis and Advection', 'm3 s-3') + 'Kinetic Energy Source from Coriolis and Advection', 'm3 s-3', & + conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_Coradv>0) call safe_alloc_ptr(CS%KE_Coradv,isd,ied,jsd,jed,nz) CS%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, Time, & - 'Kinetic Energy Source from Advection', 'm3 s-3') + 'Kinetic Energy Source from Advection', 'm3 s-3', & + conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_adv>0) call safe_alloc_ptr(CS%KE_adv,isd,ied,jsd,jed,nz) CS%id_KE_visc = register_diag_field('ocean_model', 'KE_visc', diag%axesTL, Time, & - 'Kinetic Energy Source from Vertical Viscosity and Stresses', 'm3 s-3') + 'Kinetic Energy Source from Vertical Viscosity and Stresses', 'm3 s-3', & + conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_visc>0) call safe_alloc_ptr(CS%KE_visc,isd,ied,jsd,jed,nz) CS%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, Time, & - 'Kinetic Energy Source from Horizontal Viscosity', 'm3 s-3') + 'Kinetic Energy Source from Horizontal Viscosity', 'm3 s-3', & + conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_horvisc>0) call safe_alloc_ptr(CS%KE_horvisc,isd,ied,jsd,jed,nz) if (.not. adiabatic) then CS%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, Time, & - 'Kinetic Energy Source from Diapycnal Diffusion', 'm3 s-3') + 'Kinetic Energy Source from Diapycnal Diffusion', 'm3 s-3', & + conversion=GV%H_to_m*(US%L_T_to_m_s**2)*US%s_to_T) if (CS%id_KE_dia>0) call safe_alloc_ptr(CS%KE_dia,isd,ied,jsd,jed,nz) endif From 95a677951797aee6db81657bf9ac7b9176e93a7e Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 4 Sep 2019 11:28:09 -0400 Subject: [PATCH 4/9] Barotropic transport diagnostic scaling This patch fixes the scaling of the diagnostic parameters in the [uv]hbt and [uv]hbt_hifreq transport diagnostics. --- src/core/MOM_barotropic.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 347b5a1936..7b2f367487 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -4245,10 +4245,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, conversion=GV%H_to_m) CS%id_uhbt_hifreq = register_diag_field('ocean_model', 'uhbt_hifreq', diag%axesCu1, Time, & 'High Frequency Barotropic zonal transport', 'm3 s-1', & - conversion=GV%H_to_m*US%L_T_to_m_s) + conversion=GV%H_to_m*US%L_to_m*US%L_T_to_m_s) CS%id_vhbt_hifreq = register_diag_field('ocean_model', 'vhbt_hifreq', diag%axesCv1, Time, & 'High Frequency Barotropic meridional transport', 'm3 s-1', & - conversion=GV%H_to_m*US%L_T_to_m_s) + conversion=GV%H_to_m*US%L_to_m*US%L_T_to_m_s) CS%id_frhatu = register_diag_field('ocean_model', 'frhatu', diag%axesCuL, Time, & 'Fractional thickness of layers in u-columns', 'nondim') CS%id_frhatv = register_diag_field('ocean_model', 'frhatv', diag%axesCvL, Time, & @@ -4259,10 +4259,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, 'Predictor Fractional thickness of layers in v-columns', 'nondim') CS%id_uhbt = register_diag_field('ocean_model', 'uhbt', diag%axesCu1, Time, & 'Barotropic zonal transport averaged over a baroclinic step', 'm3 s-1', & - conversion=GV%H_to_m*US%L_T_to_m_s) + conversion=GV%H_to_m*US%L_to_m*US%L_T_to_m_s) CS%id_vhbt = register_diag_field('ocean_model', 'vhbt', diag%axesCv1, Time, & 'Barotropic meridional transport averaged over a baroclinic step', 'm3 s-1', & - conversion=GV%H_to_m*US%L_T_to_m_s) + conversion=GV%H_to_m*US%L_to_m*US%L_T_to_m_s) if (use_BT_cont_type) then CS%id_BTC_FA_u_EE = register_diag_field('ocean_model', 'BTC_FA_u_EE', diag%axesCu1, Time, & From 833ad4637c0efbe3f1e0fa17bb63d3b90c5fe4b8 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 4 Sep 2019 14:02:40 -0400 Subject: [PATCH 5/9] Diagnostic global average dimensional fix Length scaling for global averaging is no longer required and has been removed. This has fixed the dimensional scaling of any diagnostics using this function. --- src/framework/MOM_spatial_means.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_spatial_means.F90 b/src/framework/MOM_spatial_means.F90 index f7084ee7ea..f6282faa52 100644 --- a/src/framework/MOM_spatial_means.F90 +++ b/src/framework/MOM_spatial_means.F90 @@ -35,8 +35,8 @@ function global_area_mean(var,G) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec 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)) ) + do j=js,je ; do i=is,ie + tmpForSumming(i,j) = var(i,j) * (G%areaT(i,j) * G%mask2dT(i,j)) enddo ; enddo global_area_mean = reproducing_sum( tmpForSumming ) * G%IareaT_global From 89b736fd28ea9ed81337204b7595213d5ef28e62 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 4 Sep 2019 16:35:18 -0400 Subject: [PATCH 6/9] (+) Tracer diffusive flux diagnostic scaling Diagnostics of diffusive tracer fluxes (*_diff[xy], *_diff[xy]_2d) were not scaling due to issues with the Coef_x term, which was in units of L2 while the tracer itself was in units of m2. Although the Coef_x term was being de-scaled to meters, the cumulative sum was failing to reproduce for many tracers, either being reverted (in temp) or failing to reproduce all bits (in salt). We resolve this by re-defining the diffusive flux terms as H L2 s-1, rather than H m2 s-1, and moving the scaling to the final conversion factor, which restored reproducibility of all tracers. This patch contains an API change. We have added the unit scaling struct (US) to the register_tracer_diagnostic function. --- src/core/MOM.F90 | 2 +- src/tracer/MOM_tracer_hor_diff.F90 | 16 ++++++++-------- src/tracer/MOM_tracer_registry.F90 | 30 +++++++++++++++++++----------- 3 files changed, 28 insertions(+), 20 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 23c11cc05b..8a5f9bfe02 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2371,7 +2371,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call register_surface_diags(Time, G, CS%sfc_IDs, CS%diag, CS%tv) call register_diags(Time, G, GV, US, CS%IDs, CS%diag) call register_transport_diags(Time, G, GV, US, CS%transport_IDs, CS%diag) - call register_tracer_diagnostics(CS%tracer_Reg, CS%h, Time, diag, G, GV, & + call register_tracer_diagnostics(CS%tracer_Reg, CS%h, Time, diag, G, GV, US, & CS%use_ALE_algorithm) if (CS%use_ALE_algorithm) then call ALE_register_diags(Time, G, GV, US, diag, CS%ALE_CSp) diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index 098a647ec8..5577115a48 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -449,20 +449,20 @@ subroutine tracer_hordiff(h, dt, MEKE, VarMix, G, GV, US, CS, Reg, tv, do_online Coef_y(i,J) * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)))) enddo ; enddo if (associated(Reg%Tr(m)%df_x)) then ; do j=js,je ; do I=G%IscB,G%IecB - Reg%Tr(m)%df_x(I,j,k) = Reg%Tr(m)%df_x(I,j,k) + US%L_to_m**2*Coef_x(I,j) * & - (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k))*Idt + Reg%Tr(m)%df_x(I,j,k) = Reg%Tr(m)%df_x(I,j,k) + Coef_x(I,j) & + * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k)) * Idt enddo ; enddo ; endif if (associated(Reg%Tr(m)%df_y)) then ; do J=G%JscB,G%JecB ; do i=is,ie - Reg%Tr(m)%df_y(i,J,k) = Reg%Tr(m)%df_y(i,J,k) + US%L_to_m**2*Coef_y(i,J) * & - (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k))*Idt + Reg%Tr(m)%df_y(i,J,k) = Reg%Tr(m)%df_y(i,J,k) + Coef_y(i,J) & + * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)) * Idt enddo ; enddo ; endif if (associated(Reg%Tr(m)%df2d_x)) then ; do j=js,je ; do I=G%IscB,G%IecB - Reg%Tr(m)%df2d_x(I,j) = Reg%Tr(m)%df2d_x(I,j) + US%L_to_m**2*Coef_x(I,j) * & - (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k))*Idt + Reg%Tr(m)%df2d_x(I,j) = Reg%Tr(m)%df2d_x(I,j) + Coef_x(I,j) & + * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i+1,j,k)) * Idt enddo ; enddo ; endif if (associated(Reg%Tr(m)%df2d_y)) then ; do J=G%JscB,G%JecB ; do i=is,ie - Reg%Tr(m)%df2d_y(i,J) = Reg%Tr(m)%df2d_y(i,J) + US%L_to_m**2*Coef_y(i,J) * & - (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k))*Idt + Reg%Tr(m)%df2d_y(i,J) = Reg%Tr(m)%df2d_y(i,J) + Coef_y(i,J) & + * (Reg%Tr(m)%t(i,j,k) - Reg%Tr(m)%t(i,j+1,k)) * Idt enddo ; enddo ; endif do j=js,je ; do i=is,ie Reg%Tr(m)%t(i,j,k) = Reg%Tr(m)%t(i,j,k) + dTr(i,j) diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index 4680c058b4..41299be3e8 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -20,6 +20,7 @@ module MOM_tracer_registry 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_verticalGrid, only : verticalGrid_type implicit none ; private @@ -53,13 +54,13 @@ module MOM_tracer_registry !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: df_x => NULL() !< diagnostic array for x-diffusive tracer flux - !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + !! [conc H L2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:,:), pointer :: df_y => NULL() !< diagnostic array for y-diffusive tracer flux - !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + !! [conc H L2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_x => NULL() !< diagnostic vertical sum x-diffusive flux - !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + !! [conc H L2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_y => NULL() !< diagnostic vertical sum y-diffusive flux - !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] + !! [conc H L2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_conc_x => NULL() !< diagnostic vertical sum x-diffusive content flux !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1] real, dimension(:,:), pointer :: df2d_conc_y => NULL() !< diagnostic vertical sum y-diffusive content flux @@ -320,9 +321,10 @@ end subroutine lock_tracer_registry !> register_tracer_diagnostics does a set of register_diag_field calls for any previously !! registered in a tracer registry with a value of registry_diags set to .true. -subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE) +subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE) 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 !< Layer thicknesses @@ -396,10 +398,12 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE) trim(flux_units), v_extensive = .true., x_cell_method = 'sum') 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., y_cell_method = 'sum') + trim(flux_units), v_extensive = .true., conversion=US%L_to_m**2, & + y_cell_method = 'sum') 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., x_cell_method = 'sum') + trim(flux_units), v_extensive = .true., conversion=US%L_to_m**2, & + x_cell_method = 'sum') 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), & @@ -409,10 +413,12 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE) flux_units, v_extensive=.true., conversion=Tr%flux_scale, x_cell_method = 'sum') Tr%id_dfx = register_diag_field("ocean_model", trim(shortnm)//"_diffx", & diag%axesCuL, Time, "Diffusive Zonal Flux of "//trim(flux_longname), & - flux_units, v_extensive=.true., conversion=Tr%flux_scale, y_cell_method = 'sum') + flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale, & + y_cell_method='sum') Tr%id_dfy = register_diag_field("ocean_model", trim(shortnm)//"_diffy", & diag%axesCvL, Time, "Diffusive Meridional Flux of "//trim(flux_longname), & - flux_units, v_extensive=.true., conversion=Tr%flux_scale, x_cell_method = 'sum') + flux_units, v_extensive=.true., conversion=(US%L_to_m**2)*Tr%flux_scale, & + x_cell_method='sum') endif 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) @@ -430,11 +436,13 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, use_ALE) Tr%id_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_diffx_2d", & diag%axesCu1, Time, & "Vertically Integrated Diffusive Zonal Flux of "//trim(flux_longname), & - flux_units, conversion=Tr%flux_scale, y_cell_method = 'sum') + flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale, & + y_cell_method='sum') Tr%id_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_diffy_2d", & diag%axesCv1, Time, & "Vertically Integrated Diffusive Meridional Flux of "//trim(flux_longname), & - flux_units, conversion=Tr%flux_scale, x_cell_method = 'sum') + flux_units, conversion=(US%L_to_m**2)*Tr%flux_scale, & + x_cell_method='sum') if (Tr%id_adx_2d > 0) call safe_alloc_ptr(Tr%ad2d_x,IsdB,IedB,jsd,jed) if (Tr%id_ady_2d > 0) call safe_alloc_ptr(Tr%ad2d_y,isd,ied,JsdB,JedB) From 70362368702df81050dd4796a2008e47edd43fcb Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 5 Sep 2019 12:01:20 -0400 Subject: [PATCH 7/9] Dimensional scaling of MEKE_src Dimensional scaling factor was added to the MEKE_src diagnostic, fixing the rescaling reproducibility. --- src/parameterizations/lateral/MOM_MEKE.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index d204db1305..6c038a4eb5 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -1148,7 +1148,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) 'MEKE derived barotropic eddy-velocity scale', 'm s-1', conversion=US%L_T_to_m_s) if (.not. associated(MEKE%MEKE)) CS%id_Ut = -1 CS%id_src = register_diag_field('ocean_model', 'MEKE_src', diag%axesT1, Time, & - 'MEKE energy source', 'm2 s-3') + 'MEKE energy source', 'm2 s-3', conversion=(US%L_T_to_m_s**2)*US%s_to_T) CS%id_decay = register_diag_field('ocean_model', 'MEKE_decay', diag%axesT1, Time, & 'MEKE decay rate', 's-1', conversion=US%s_to_T) CS%id_GM_src = register_diag_field('ocean_model', 'MEKE_GM_src', diag%axesT1, Time, & From d64b472f519611ae5f4e03b5f6d5ccc4e1af7d38 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 5 Sep 2019 15:31:03 -0400 Subject: [PATCH 8/9] Unsplit timestep diag fixes (CAu, etc) This patch fixes a bug in the registration of CAu, where the conversion argument had been included in the string of the dimensions of the diagnostic registration. The units of various accelerations have also been changed from 'meter second-2' to 'm s-2' for consistent with the rest of the model. --- src/core/MOM_dynamics_unsplit.F90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 108f4c8943..9e8be65d7a 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -678,13 +678,17 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS 'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true., & conversion=H_convert*US%L_to_m**2*US%s_to_T) CS%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, Time, & - 'Zonal Coriolis and Advective Acceleration', 'meter second-2, conversion=US%L_T2_to_m_s2') + 'Zonal Coriolis and Advective Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) CS%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, Time, & - 'Meridional Coriolis and Advective Acceleration', 'meter second-2', conversion=US%L_T2_to_m_s2) + 'Meridional Coriolis and Advective Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) CS%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, Time, & - 'Zonal Pressure Force Acceleration', 'meter second-2', conversion=US%L_T2_to_m_s2) + 'Zonal Pressure Force Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & - 'Meridional Pressure Force Acceleration', 'meter second-2', conversion=US%L_T2_to_m_s2) + 'Meridional Pressure Force Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) From c9f9787beae1a5b88656720da48c7127c4a31802 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 5 Sep 2019 17:37:56 -0400 Subject: [PATCH 9/9] ALE and mixed layer scaling diag fixes This patch fixes the dimensional scaling of the pre-ALE u and v diagnostics, and the uhml and vhml mixed layer diagnostics. --- src/ALE/MOM_ALE.F90 | 4 ++-- src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 33b498a60a..def4148f25 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -246,9 +246,9 @@ subroutine ALE_register_diags(Time, G, GV, US, diag, CS) ! These diagnostics of the state variables before ALE are useful for ! debugging the ALE code. CS%id_u_preale = register_diag_field('ocean_model', 'u_preale', diag%axesCuL, Time, & - 'Zonal velocity before remapping', 'm s-1') + 'Zonal velocity before remapping', 'm s-1', conversion=US%L_T_to_m_s) CS%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, Time, & - 'Meridional velocity before remapping', 'm s-1') + '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.) CS%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, Time, & diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index c4a2d0c38f..6543ab7af6 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -884,8 +884,8 @@ 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%Rho0*US%L_to_m**2*US%s_to_T - else ; flux_to_kg_per_s = US%L_to_m**2*US%s_to_T ; endif + 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 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, &