Merge pull request #1002 from marshallward/diag_scaling_p2

Diag scaling p2

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, &
@@ -259,11 +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')
-  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.)
+      '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
 

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')

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

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
 

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, &

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,15 +742,16 @@ 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
       enddo
     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

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
 

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)