From ee6baaa23b5bd0179ee41c59932150eacdae96d2 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 27 Sep 2019 08:23:16 -0400 Subject: [PATCH] +Changed units of GV%Rho0 to [R] Changed the units of GV%Rho0 from [kg m-3] to [R] for dimensional consistency testing. This required the addition of unit_scale_type arguments to several interfaces. All answers are bitwise identical, but new arguments have been added to several public interfaces and the units of an element in a public type have changed. --- config_src/mct_driver/mom_ocean_model_mct.F90 | 2 +- .../nuopc_driver/mom_ocean_model_nuopc.F90 | 2 +- src/core/MOM.F90 | 4 ++-- src/core/MOM_PressureForce_Montgomery.F90 | 4 ++-- src/core/MOM_PressureForce_analytic_FV.F90 | 4 ++-- src/core/MOM_PressureForce_blocked_AFV.F90 | 4 ++-- src/core/MOM_barotropic.F90 | 8 +++---- src/core/MOM_forcing_type.F90 | 14 ++++++----- src/core/MOM_isopycnal_slopes.F90 | 2 +- src/core/MOM_verticalGrid.F90 | 8 +++---- src/diagnostics/MOM_diagnostics.F90 | 4 ++-- src/diagnostics/MOM_sum_output.F90 | 8 +++---- src/diagnostics/MOM_wave_speed.F90 | 4 ++-- src/diagnostics/MOM_wave_structure.F90 | 6 ++--- .../MOM_coord_initialization.F90 | 24 +++++++++---------- .../MOM_state_initialization.F90 | 23 +++++++++--------- src/parameterizations/lateral/MOM_MEKE.F90 | 2 +- .../lateral/MOM_internal_tides.F90 | 4 ++-- .../lateral/MOM_mixed_layer_restrat.F90 | 4 ++-- .../lateral/MOM_thickness_diffuse.F90 | 2 +- .../vertical/MOM_CVMix_KPP.F90 | 2 +- .../vertical/MOM_CVMix_conv.F90 | 2 +- .../vertical/MOM_CVMix_shear.F90 | 2 +- .../vertical/MOM_bulk_mixed_layer.F90 | 18 +++++++------- .../vertical/MOM_diabatic_aux.F90 | 8 +++---- .../vertical/MOM_energetic_PBL.F90 | 8 +++---- .../vertical/MOM_internal_tide_input.F90 | 4 ++-- .../vertical/MOM_kappa_shear.F90 | 2 +- .../vertical/MOM_set_diffusivity.F90 | 24 +++++++++---------- .../vertical/MOM_set_viscosity.F90 | 6 ++--- .../vertical/MOM_tidal_mixing.F90 | 16 ++++++------- .../vertical/MOM_vert_friction.F90 | 4 ++-- src/tracer/MOM_OCMIP2_CFC.F90 | 4 ++-- src/user/BFB_initialization.F90 | 6 ++--- src/user/DOME_initialization.F90 | 2 +- src/user/MOM_wave_interface.F90 | 2 +- src/user/Rossby_front_2d_initialization.F90 | 13 +++++----- 37 files changed, 130 insertions(+), 126 deletions(-) diff --git a/config_src/mct_driver/mom_ocean_model_mct.F90 b/config_src/mct_driver/mom_ocean_model_mct.F90 index 4f1c7d963a..8873f283ff 100644 --- a/config_src/mct_driver/mom_ocean_model_mct.F90 +++ b/config_src/mct_driver/mom_ocean_model_mct.F90 @@ -582,7 +582,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & #endif endif - call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%GV%Rho0) + call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%US%R_to_kg_m3*OS%GV%Rho0) call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid) if (OS%use_waves) then diff --git a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 index e04064f672..db475754c9 100644 --- a/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 +++ b/config_src/nuopc_driver/mom_ocean_model_nuopc.F90 @@ -570,7 +570,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) !weight of the current flux in the running average #endif endif - call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%GV%Rho0) + call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%US, OS%US%R_to_kg_m3*OS%GV%Rho0) call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid) if (OS%use_waves) then diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 775a15b427..7837f72b3b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2678,7 +2678,7 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS) call calculate_density(tv%T(i,j,1), tv%S(i,j,1), p_atm(i,j)/2.0, & Rho_conv, tv%eqn_of_state) else - Rho_conv=GV%Rho0 + Rho_conv = US%R_to_kg_m3*GV%Rho0 endif IgR0 = 1.0 / (Rho_conv * GV%mks_g_Earth) ssh(i,j) = ssh(i,j) + p_atm(i,j) * IgR0 @@ -2914,7 +2914,7 @@ subroutine extract_surface_state(CS, sfc_state) if (G%mask2dT(i,j)>0.) then ! instantaneous melt_potential [J m-2] - sfc_state%melt_potential(i,j) = CS%tv%C_p * CS%GV%Rho0 * delT(i) + sfc_state%melt_potential(i,j) = CS%tv%C_p * CS%US%R_to_kg_m3*GV%Rho0 * delT(i) endif enddo enddo ! end of j loop diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index ebcc3e4afc..e627cba724 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -435,7 +435,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, h_neglect = GV%H_subroundoff * GV%H_to_Z I_Rho0 = 1.0/CS%Rho0 - G_Rho0 = GV%g_Earth/GV%Rho0 + G_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) if (CS%tides) then ! Determine the surface height anomaly for calculating self attraction @@ -640,7 +640,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke Rho0xG = Rho0*US%L_T_to_m_s**2 * GV%g_Earth - G_Rho0 = GV%g_Earth / GV%Rho0 + G_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) use_EOS = associated(tv%eqn_of_state) z_neglect = GV%H_subroundoff*GV%H_to_Z diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index 0d56722825..3e1e2f72e1 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -531,9 +531,9 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at h_neglect = GV%H_subroundoff dz_neglect = GV%H_subroundoff * GV%H_to_Z - I_Rho0 = US%m_s_to_L_T**2 / GV%Rho0 + I_Rho0 = US%m_s_to_L_T**2 / (US%R_to_kg_m3*GV%Rho0) g_Earth_z = US%L_T_to_m_s**2 * GV%g_Earth - G_Rho0 = GV%g_Earth/GV%Rho0 + G_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) rho_ref = CS%Rho0 if (CS%tides) then diff --git a/src/core/MOM_PressureForce_blocked_AFV.F90 b/src/core/MOM_PressureForce_blocked_AFV.F90 index 87b325ef15..87d8d0fc8f 100644 --- a/src/core/MOM_PressureForce_blocked_AFV.F90 +++ b/src/core/MOM_PressureForce_blocked_AFV.F90 @@ -515,9 +515,9 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, h_neglect = GV%H_subroundoff dz_neglect = GV%H_subroundoff * GV%H_to_Z - I_Rho0 = US%m_s_to_L_T**2 / GV%Rho0 + I_Rho0 = US%m_s_to_L_T**2 / (US%R_to_kg_m3*GV%Rho0) g_Earth_z = US%L_T_to_m_s**2 * GV%g_Earth - G_Rho0 = GV%g_Earth / GV%Rho0 + G_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) rho_ref = CS%Rho0 if (CS%tides) then diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 7b2f367487..fd2d6560be 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -724,8 +724,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, dtbt = dt_in_T * Instep bebt = CS%bebt be_proj = CS%bebt - mass_accel_to_Z = US%m_to_L*US%T_to_s**2 * US%m_to_Z / GV%Rho0 - mass_to_Z = US%m_to_Z / GV%Rho0 + mass_accel_to_Z = US%m_to_L*US%T_to_s**2 * US%m_to_Z / (US%R_to_kg_m3*GV%Rho0) + mass_to_Z = US%m_to_Z / (US%R_to_kg_m3*GV%Rho0) !--- setup the weight when computing vbt_trans and ubt_trans if (project_velocity) then @@ -4345,10 +4345,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, enddo ; enddo ! else ! do j=js,je ; do I=is-1,ie -! CS%IDatu(I,j) = G%mask2dCu(I,j) * 2.0 / (GV%Rho0*(G%bathyT(i+1,j) + G%bathyT(i,j))) +! CS%IDatu(I,j) = G%mask2dCu(I,j) * 2.0 / (US%R_to_kg_m3*GV%Rho0*(G%bathyT(i+1,j) + G%bathyT(i,j))) ! enddo ; enddo ! do J=js-1,je ; do i=is,ie -! CS%IDatv(i,J) = G%mask2dCv(i,J) * 2.0 / (GV%Rho0*(G%bathyT(i,j+1) + G%bathyT(i,j))) +! CS%IDatv(i,J) = G%mask2dCv(i,J) * 2.0 / (US%R_to_kg_m3*GV%Rho0*(G%bathyT(i,j+1) + G%bathyT(i,j))) ! enddo ; enddo ! endif diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index a8c6f7bf1a..ececc6d1e7 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -336,7 +336,7 @@ module MOM_forcing_type !! for optimization purposes. The 2d (i,j) wrapper is the next subroutine below. !! This routine multiplies fluxes by dt, so that the result is an accumulation of fluxes !! over a time step. -subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, & +subroutine extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & FluxRescaleDepth, useRiverHeatContent, useCalvingHeatContent, & h, T, netMassInOut, netMassOut, net_heat, net_salt, pen_SW_bnd, tv, & aggregate_FW, nonpenSW, netmassInOut_rate,net_Heat_Rate, & @@ -344,6 +344,7 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(forcing), intent(inout) :: fluxes !< structure containing pointers to possible !! forcing fields. NULL unused fields. type(optics_type), pointer :: optics !< pointer to optics @@ -433,7 +434,7 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, !}BGR Ih_limit = 1.0 / FluxRescaleDepth - Irho0 = 1.0 / GV%Rho0 + Irho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) I_Cp = 1.0 / fluxes%C_p J_m2_to_H = 1.0 / (GV%H_to_kg_m2 * fluxes%C_p) @@ -804,13 +805,14 @@ end subroutine extractFluxes1d !> 2d wrapper for 1d extract fluxes from surface fluxes type. !! This subroutine extracts fluxes from the surface fluxes type. It multiplies the !! fluxes by dt, so that the result is an accumulation of the fluxes over a time step. -subroutine extractFluxes2d(G, GV, fluxes, optics, nsw, dt, FluxRescaleDepth, & +subroutine extractFluxes2d(G, GV, US, fluxes, optics, nsw, dt, FluxRescaleDepth, & useRiverHeatContent, useCalvingHeatContent, h, T, & netMassInOut, netMassOut, net_heat, Net_salt, Pen_SW_bnd, tv, & aggregate_FW) type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(forcing), intent(inout) :: fluxes !< structure containing pointers to forcing. type(optics_type), pointer :: optics !< pointer to optics integer, intent(in) :: nsw !< number of bands of penetrating SW @@ -854,7 +856,7 @@ subroutine extractFluxes2d(G, GV, fluxes, optics, nsw, dt, FluxRescaleDepth, & !$OMP h,T,netMassInOut,netMassOut,Net_heat,Net_salt,Pen_SW_bnd,tv, & !$OMP aggregate_FW) do j=G%jsc, G%jec - call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, & + call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & FluxRescaleDepth, useRiverHeatContent, useCalvingHeatContent,& h(:,j,:), T(:,j,:), netMassInOut(:,j), netMassOut(:,j), & net_heat(:,j), net_salt(:,j), pen_SW_bnd(:,:,j), tv, aggregate_FW) @@ -916,7 +918,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt depthBeforeScalingFluxes = max( GV%Angstrom_H, 1.e-30*GV%m_to_H ) pressure(:) = 0. ! Ignore atmospheric pressure - GoRho = (GV%g_Earth*US%m_to_Z * GV%H_to_m*US%T_to_s) / GV%Rho0 + GoRho = (GV%g_Earth*US%m_to_Z * GV%H_to_m*US%T_to_s) / (US%R_to_kg_m3*GV%Rho0) start = 1 + G%isc - G%isd npts = 1 + G%iec - G%isc @@ -929,7 +931,7 @@ subroutine calculateBuoyancyFlux1d(G, GV, US, fluxes, optics, nsw, h, Temp, Salt ! netSalt = salt via surface fluxes [ppt H s-1 ~> ppt m s-1 or gSalt m-2 s-1] ! Note that unlike other calls to extractFLuxes1d() that return the time-integrated flux ! this call returns the rate because dt=1 - call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, & + call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & depthBeforeScalingFluxes, useRiverHeatContent, useCalvingHeatContent, & h(:,j,:), Temp(:,j,:), netH, netEvap, netHeatMinusSW, & netSalt, penSWbnd, tv, .false., skip_diags=skip_diags) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 61118074fd..ae06413e90 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -121,7 +121,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & present_N2_u = PRESENT(N2_u) present_N2_v = PRESENT(N2_v) - G_Rho0 = (US%L_to_Z*L_to_Z*GV%g_Earth) / GV%Rho0 + G_Rho0 = (US%L_to_Z*L_to_Z*GV%g_Earth) / (US%R_to_kg_m3*GV%Rho0) if (present_N2_u) then do j=js,je ; do I=is-1,ie N2_u(I,j,1) = 0. diff --git a/src/core/MOM_verticalGrid.F90 b/src/core/MOM_verticalGrid.F90 index 2d313c5148..093db28c07 100644 --- a/src/core/MOM_verticalGrid.F90 +++ b/src/core/MOM_verticalGrid.F90 @@ -101,7 +101,7 @@ subroutine verticalGridInit( param_file, GV, US ) "calculate accelerations and the mass for conservation "//& "properties, or with BOUSSINSEQ false to convert some "//& "parameters from vertical units of m to kg m-2.", & - units="kg m-3", default=1035.0) + units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "BOUSSINESQ", GV%Boussinesq, & "If true, make the Boussinesq approximation.", default=.true.) call get_param(param_file, mdl, "ANGSTROM", GV%Angstrom_m, & @@ -143,15 +143,15 @@ subroutine verticalGridInit( param_file, GV, US ) GV%ke = nk if (GV%Boussinesq) then - GV%H_to_kg_m2 = GV%Rho0 * GV%H_to_m + GV%H_to_kg_m2 = US%R_to_kg_m3*GV%Rho0 * GV%H_to_m 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%m_to_H = US%R_to_kg_m3*GV%Rho0 * GV%kg_m2_to_H + GV%H_to_m = GV%H_to_kg_m2 / (US%R_to_kg_m3*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 diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 1d9e7f39b7..47aeaf547e 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -842,7 +842,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) z_bot(i,j) = z_top(i,j) - GV%H_to_Z*h(i,j,k) enddo ; enddo call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), & - z_top, z_bot, 0.0, GV%Rho0, GV%mks_g_Earth*US%Z_to_m, & + z_top, z_bot, 0.0, US%R_to_kg_m3*GV%Rho0, GV%mks_g_Earth*US%Z_to_m, & G%HI, G%HI, tv%eqn_of_state, dpress) do j=js,je ; do i=is,ie mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth @@ -2006,7 +2006,7 @@ subroutine write_static_fields(G, GV, US, tv, diag) id = register_static_field('ocean_model','Rho_0', diag%axesNull, & 'mean ocean density used with the Boussinesq approximation', & - 'kg m-3', cmor_field_name='rhozero', & + 'kg m-3', cmor_field_name='rhozero', conversion=US%R_to_kg_m3, & cmor_standard_name='reference_sea_water_density_for_boussinesq_approximation', & cmor_long_name='reference sea water density for boussinesq approximation') if (id > 0) call post_data(id, GV%Rho0, diag, .true.) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index d6f495faa5..9d8cff542f 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -542,7 +542,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ tmp1(i,j,k) = H_to_kg_m2 * h(i,j,k) * areaTm(i,j) enddo ; enddo ; enddo mass_tot = reproducing_sum(tmp1, sums=mass_lay, EFP_sum=mass_EFP) - do k=1,nz ; vol_lay(k) = US%m_to_Z * (mass_lay(k) / GV%Rho0) ; enddo + do k=1,nz ; vol_lay(k) = US%m_to_Z * (mass_lay(k) / (US%R_to_kg_m3*GV%Rho0)) ; enddo endif endif ! Boussinesq @@ -666,7 +666,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ hint = Z_0APE(K) + (hbelow - G%bathyT(i,j)) hbot = Z_0APE(K) - G%bathyT(i,j) hbot = (hbot + ABS(hbot)) * 0.5 - PE_pt(i,j,K) = 0.5 * areaTm(i,j) * US%Z_to_m*US%L_T_to_m_s**2*(GV%Rho0*GV%g_prime(K)) * & + PE_pt(i,j,K) = 0.5 * areaTm(i,j) * US%Z_to_m*US%L_T_to_m_s**2*(US%R_to_kg_m3*GV%Rho0*GV%g_prime(K)) * & (hint * hint - hbot * hbot) enddo enddo ; enddo @@ -675,7 +675,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ do k=nz,1,-1 hint = Z_0APE(K) + eta(i,j,K) ! eta and H_0 have opposite signs. hbot = max(Z_0APE(K) - G%bathyT(i,j), 0.0) - PE_pt(i,j,K) = 0.5 * (areaTm(i,j) * US%Z_to_m*US%L_T_to_m_s**2*(GV%Rho0*GV%g_prime(K))) * & + PE_pt(i,j,K) = 0.5 * (areaTm(i,j) * US%Z_to_m*US%L_T_to_m_s**2*(US%R_to_kg_m3*GV%Rho0*GV%g_prime(K))) * & (hint * hint - hbot * hbot) enddo enddo ; enddo @@ -750,7 +750,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_ CS%salt_prev_EFP = salt_EFP ; CS%net_salt_in_EFP = real_to_EFP(0.0) CS%heat_prev_EFP = heat_EFP ; CS%net_heat_in_EFP = real_to_EFP(0.0) endif - Irho0 = 1.0/GV%Rho0 + Irho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) if (CS%use_temperature) then Salt_chg_EFP = Salt_EFP - CS%salt_prev_EFP diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index cd9dd9dbb8..f8ac508a28 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -132,7 +132,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, & endif S => tv%S ; T => tv%T - g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / GV%Rho0 + g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) Z_to_Pa = GV%Z_to_H * GV%H_to_Pa use_EOS = associated(tv%eqn_of_state) @@ -600,7 +600,7 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos) endif ; endif S => tv%S ; T => tv%T - g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / GV%Rho0 + g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) use_EOS = associated(tv%eqn_of_state) Z_to_Pa = GV%Z_to_H * GV%H_to_Pa diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90 index d8c7cc5a02..e282b0e43a 100644 --- a/src/diagnostics/MOM_wave_structure.F90 +++ b/src/diagnostics/MOM_wave_structure.F90 @@ -178,7 +178,7 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo Pi = (4.0*atan(1.0)) S => tv%S ; T => tv%T - g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth /GV%Rho0 + g_Rho0 = US%L_T_to_m_s**2 * GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) cg_subRO = 1e-100*US%m_s_to_L_T ! The hard-coded value here might need to increase. use_EOS = associated(tv%eqn_of_state) @@ -479,8 +479,8 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo ! Back-calculate amplitude from energy equation if (Kmag2 > 0.0) then !### This should be simpified to use a single division. - KE_term = 0.25*GV%Rho0*( ((1.0 + f2/freq**2) / Kmag2)*int_dwdz2 + int_w2 ) - PE_term = 0.25*GV%Rho0*( int_N2w2/(US%s_to_T*freq)**2 ) + KE_term = 0.25*US%R_to_kg_m3*GV%Rho0*( ((1.0 + f2/freq**2) / Kmag2)*int_dwdz2 + int_w2 ) + PE_term = 0.25*US%R_to_kg_m3*GV%Rho0*( int_N2w2/(US%s_to_T*freq)**2 ) if (En(i,j) >= 0.0) then W0 = sqrt( En(i,j)/(KE_term + PE_term) ) else diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index 19cb9774f0..b2519d47ad 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -148,8 +148,8 @@ subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file) g_prime(1) = g_fs do k=2,nz ; g_prime(k) = g_int ; enddo - Rlay(1) = US%kg_m3_to_R*GV%Rho0 - do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(US%kg_m3_to_R*GV%Rho0/GV%g_Earth) ; enddo + Rlay(1) = GV%Rho0 + do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo call callTree_leave(trim(mdl)//'()') @@ -179,7 +179,7 @@ subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file) default=GV%mks_g_Earth, scale=US%m_s_to_L_T**2*US%Z_to_m) call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, & "The reference potential density used for layer 1.", & - units="kg m-3", default=GV%Rho0, scale=US%kg_m3_to_R) + units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "DENSITY_RANGE", Rlay_range, & "The range of reference potential densities in the layers.", & units="kg m-3", default=2.0, scale=US%kg_m3_to_R) @@ -191,7 +191,7 @@ subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file) enddo ! These statements set the interface reduced gravities. ! do k=2,nz - g_prime(k) = (GV%g_Earth/(US%kg_m3_to_R*GV%Rho0)) * (Rlay(k) - Rlay(k-1)) + g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) enddo call callTree_leave(trim(mdl)//'()') @@ -243,7 +243,7 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state, scale=US%kg_m3_to_R) ! These statements set the layer densities. ! - do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(US%kg_m3_to_R*GV%Rho0/GV%g_Earth) ; enddo + do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo call callTree_leave(trim(mdl)//'()') end subroutine set_coord_from_TS_ref @@ -291,7 +291,7 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, & g_prime(1) = g_fs do k=1,nz ; Pref(k) = P_ref ; enddo call calculate_density(T0, S0, Pref, Rlay, 1, nz, eqn_of_state, scale=US%kg_m3_to_R) - do k=2,nz; g_prime(k) = (GV%g_Earth/(US%kg_m3_to_R*GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo + do k=2,nz; g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo call callTree_leave(trim(mdl)//'()') end subroutine set_coord_from_TS_profile @@ -375,7 +375,7 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, & do k=k_light-1,1,-1 Rlay(k) = 2.0*Rlay(k+1) - Rlay(k+2) enddo - do k=2,nz ; g_prime(k) = (GV%g_Earth/(US%kg_m3_to_R*GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo + do k=2,nz ; g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo call callTree_leave(trim(mdl)//'()') end subroutine set_coord_from_TS_range @@ -418,7 +418,7 @@ subroutine set_coord_from_file(Rlay, g_prime, GV, US, param_file) call read_axis_data(filename, coord_var, Rlay) do k=1,nz ; Rlay(k) = US%kg_m3_to_R*Rlay(k) ; enddo g_prime(1) = g_fs - do k=2,nz ; g_prime(k) = (GV%g_Earth/(US%kg_m3_to_R*GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo + do k=2,nz ; g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) ; enddo do k=1,nz ; if (g_prime(k) <= 0.0) then call MOM_error(FATAL, "MOM_initialization set_coord_from_file: "//& "Zero or negative g_primes read from variable "//"Layer"//" in file "//& @@ -451,7 +451,7 @@ subroutine set_coord_linear(Rlay, g_prime, GV, US, param_file) call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, & "The reference potential density used for the surface interface.", & - units="kg m-3", default=GV%Rho0, scale=US%kg_m3_to_R) + units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "DENSITY_RANGE", Rlay_range, & "The range of reference potential densities across all interfaces.", & units="kg m-3", default=2.0, scale=US%kg_m3_to_R) @@ -468,7 +468,7 @@ subroutine set_coord_linear(Rlay, g_prime, GV, US, param_file) ! These statements set the interface reduced gravities. g_prime(1) = g_fs do k=2,nz - g_prime(k) = (GV%g_Earth/(US%kg_m3_to_R*GV%Rho0)) * (Rlay(k) - Rlay(k-1)) + g_prime(k) = (GV%g_Earth/(GV%Rho0)) * (Rlay(k) - Rlay(k-1)) enddo call callTree_leave(trim(mdl)//'()') @@ -499,8 +499,8 @@ subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file) g_prime(1) = g_fs do k=2,nz ; g_prime(k) = 0. ; enddo - Rlay(1) = US%kg_m3_to_R*GV%Rho0 - do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(US%kg_m3_to_R*GV%Rho0/GV%g_Earth) ; enddo + Rlay(1) = GV%Rho0 + do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo call callTree_leave(trim(mdl)//'()') diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 9210da72da..c061169854 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -318,7 +318,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & case ("soliton"); call soliton_initialize_thickness(h, G, GV, US) case ("phillips"); call Phillips_initialize_thickness(h, G, GV, US, PF, & just_read_params=just_read) - case ("rossby_front"); call Rossby_front_initialize_thickness(h, G, GV, & + case ("rossby_front"); call Rossby_front_initialize_thickness(h, G, GV, US, & PF, just_read_params=just_read) case ("USER"); call user_initialize_thickness(h, G, GV, PF, & just_read_params=just_read) @@ -952,7 +952,7 @@ subroutine convert_thickness(h, G, GV, US, tv) max_itt = 10 Boussinesq = GV%Boussinesq I_gEarth = 1.0 / (GV%mks_g_Earth) - Hm_rho_to_Pa = GV%mks_g_Earth * GV%H_to_m ! = GV%H_to_Pa / GV%Rho0 + Hm_rho_to_Pa = GV%mks_g_Earth * GV%H_to_m ! = GV%H_to_Pa / (US%R_to_kg_m3*GV%Rho0) if (Boussinesq) then call MOM_error(FATAL,"Not yet converting thickness with Boussinesq approx.") @@ -995,7 +995,7 @@ subroutine convert_thickness(h, G, GV, US, tv) do k=1,nz ; do j=js,je ; do i=is,ie h(i,j,k) = (h(i,j,k) * US%R_to_kg_m3*GV%Rlay(k)) * Hm_rho_to_Pa * GV%kg_m2_to_H**2 ! This is mathematically equivalent to - ! h(i,j,k) = h(i,j,k) * (US%R_to_kg_m3*GV%Rlay(k) / GV%Rho0) + ! h(i,j,k) = h(i,j,k) * (GV%Rlay(k) / GV%Rho0) enddo ; enddo ; enddo endif endif @@ -1154,7 +1154,7 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read_params) endif do j=G%jsc,G%jec ; do i=G%isc,G%iec - call cut_off_column_top(GV%ke, tv, GV, GV%mks_g_Earth*US%Z_to_m, G%bathyT(i,j), & + call cut_off_column_top(GV%ke, tv, GV, US, GV%mks_g_Earth*US%Z_to_m, G%bathyT(i,j), & min_thickness, tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), & tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), p_surf(i,j), h(i,j,:), remap_CS, & z_tol=1.0e-5*US%m_to_Z) @@ -1165,11 +1165,12 @@ end subroutine trim_for_ice !> Adjust the layer thicknesses by removing the top of the water column above the !! depth where the hydrostatic pressure matches p_surf -subroutine cut_off_column_top(nk, tv, GV, G_earth, depth, min_thickness, & +subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, & T, T_t, T_b, S, S_t, S_b, p_surf, h, remap_CS, z_tol) integer, intent(in) :: nk !< Number of layers type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics 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 real, intent(in) :: G_earth !< Gravitational acceleration [m2 Z-1 s-2 ~> m s-2] real, intent(in) :: depth !< Depth of ocean column [Z ~> m]. real, intent(in) :: min_thickness !< Smallest thickness allowed [Z ~> m]. @@ -1203,7 +1204,7 @@ subroutine cut_off_column_top(nk, tv, GV, G_earth, depth, min_thickness, & e_top = e(1) do k=1,nk call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), & - P_t, p_surf, GV%Rho0, G_earth, tv%eqn_of_state, & + P_t, p_surf, US%R_to_kg_m3*GV%Rho0, G_earth, tv%eqn_of_state, & P_b, z_out, z_tol=z_tol) if (z_out>=e(K)) then ! Imposed pressure was less that pressure at top of cell @@ -2406,15 +2407,15 @@ subroutine MOM_state_init_tests(G, GV, US, tv) S_t(k) = 35.-(0./500.)*e(k) S(k) = 35.+(0./500.)*z(k) S_b(k) = 35.-(0./500.)*e(k+1) - call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -GV%Rho0*GV%mks_g_Earth*z(k), & + call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -US%R_to_kg_m3*GV%Rho0*GV%mks_g_Earth*z(k), & rho(k), tv%eqn_of_state) P_tot = P_tot + GV%mks_g_Earth * rho(k) * h(k) enddo P_t = 0. do k = 1, nk - call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), & - P_t, 0.5*P_tot, GV%Rho0, GV%mks_g_Earth, tv%eqn_of_state, P_b, z_out) + call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), P_t, 0.5*P_tot, & + US%R_to_kg_m3*GV%Rho0, GV%mks_g_Earth, tv%eqn_of_state, P_b, z_out) write(0,*) k,P_t,P_b,0.5*P_tot,e(K),e(K+1),z_out P_t = P_b enddo @@ -2424,8 +2425,8 @@ subroutine MOM_state_init_tests(G, GV, US, tv) write(0,*) ' ==================================================================== ' write(0,*) '' write(0,*) h - call cut_off_column_top(nk, tv, GV, GV%mks_g_Earth, -e(nk+1), GV%Angstrom_H, & - T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS) + call cut_off_column_top(nk, tv, GV, US, GV%mks_g_Earth, -e(nk+1), GV%Angstrom_H, & + T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS) write(0,*) h end subroutine MOM_state_init_tests diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index dc44601f71..cdaa8151c9 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -675,7 +675,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & (G%dF_dy(i,j) + beta_topo_y)**2 ) - I_H = US%L_to_m*GV%Rho0 * I_mass(i,j) + I_H = US%L_to_m*US%R_to_kg_m3*GV%Rho0 * I_mass(i,j) if (KhCoeff*SN*I_H>0.) then ! Solve resid(E) = 0, where resid = Kh(E) * (SN)^2 - damp_rate(E) E diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 9014cb1dbb..4f91cd7ea5 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -202,7 +202,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & if (.not.associated(CS)) return is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nAngle = CS%NAngle - I_rho0 = 1.0 / GV%Rho0 + I_rho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) dt_in_T = US%s_to_T*dt cn_subRO = 1e-100*US%m_s_to_L_T ! The hard-coded value here might need to increase. @@ -2307,7 +2307,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) h2(i,j) = min(0.01*(G%bathyT(i,j))**2, h2(i,j)) ! Compute the fixed part; units are [kg m-2] here ! will be multiplied by N and En to get into [W m-2] - CS%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*GV%Rho0 * US%L_to_Z*kappa_itides * h2(i,j) + CS%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*US%R_to_kg_m3*GV%Rho0 * US%L_to_Z*kappa_itides * h2(i,j) enddo ; enddo deallocate(h2) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index ba241ea4b1..ca62160bc1 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -282,7 +282,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt_in_T, MLD_in uDml(:) = 0.0 ; vDml(:) = 0.0 uDml_slow(:) = 0.0 ; vDml_slow(:) = 0.0 I4dt = 0.25 / (dt_in_T) - g_Rho0 = GV%g_Earth / GV%Rho0 + g_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) h_neglect = GV%H_subroundoff dz_neglect = GV%H_subroundoff*GV%H_to_Z proper_averaging = .not. CS%MLE_use_MLD_ave_bug @@ -616,7 +616,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt_in_T, G, GV, US, uDml(:) = 0.0 ; vDml(:) = 0.0 I4dt = 0.25 / (dt_in_T) - g_Rho0 = GV%g_Earth / GV%Rho0 + g_Rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) use_EOS = associated(tv%eqn_of_state) h_neglect = GV%H_subroundoff dz_neglect = GV%H_subroundoff*GV%H_to_Z diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index dc235a369e..63385733ec 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -649,7 +649,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt_in_T, G_scale = GV%g_Earth*US%L_to_m**2*US%s_to_T**3 * GV%H_to_m h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 dz_neglect = GV%H_subroundoff*GV%H_to_Z - G_rho0 = GV%g_Earth / GV%Rho0 + G_rho0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) N2_floor = CS%N2_floor*US%Z_to_L**2 use_EOS = associated(tv%eqn_of_state) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 2ff0b3efe1..f5ee25c743 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -953,7 +953,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyF #endif ! some constants - GoRho = GV%mks_g_Earth / GV%Rho0 + GoRho = GV%mks_g_Earth / (US%R_to_kg_m3*GV%Rho0) buoy_scale = US%L_to_m**2*US%s_to_T**3 ! loop over horizontal points on processor diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 1fbbc15120..19a71116f3 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -172,7 +172,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl) real :: pref, rhok, rhokm1, dz, dh, hcorr integer :: i, j, k - g_o_rho0 = GV%mks_g_Earth / GV%Rho0 + g_o_rho0 = GV%mks_g_Earth / (US%R_to_kg_m3*GV%Rho0) ! initialize dummy variables rho_lwr(:) = 0.0; rho_1d(:) = 0.0 diff --git a/src/parameterizations/vertical/MOM_CVMix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 index 3ab0567db1..68081a97d9 100644 --- a/src/parameterizations/vertical/MOM_CVMix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -82,7 +82,7 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS ) real, parameter :: epsln = 1.e-10 !< Threshold to identify vanished layers ! some constants - GoRho = GV%mks_g_Earth / GV%Rho0 + GoRho = GV%mks_g_Earth / (US%R_to_kg_m3*GV%Rho0) do j = G%jsc, G%jec do i = G%isc, G%iec diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 1dfb1c36e4..aa101fb9f1 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -372,7 +372,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt_in_T, ea, eb, G, GV, ! dt_in_T = dt * US%s_to_T - Irho0 = 1.0 / (US%kg_m3_to_R*GV%Rho0) + Irho0 = 1.0 / (GV%Rho0) dt__diag = dt_in_T ; if (present(dt_diag)) dt__diag = dt_diag Idt_diag = 1.0 / (dt__diag) write_diags = .true. ; if (present(last_call)) write_diags = last_call @@ -533,7 +533,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt_in_T, ea, eb, G, GV, ! net_heat = heat via surface fluxes [degC H ~> degC m or degC kg m-2] ! net_salt = salt via surface fluxes [ppt H ~> dppt m or gSalt m-2] ! Pen_SW_bnd = components to penetrative shortwave radiation - call extractFluxes1d(G, GV, fluxes, optics, nsw, j, US%T_to_s*dt_in_T, & + call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, US%T_to_s*dt_in_T, & CS%H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, & h(:,1:), T(:,1:), netMassInOut, netMassOut, Net_heat, Net_salt, Pen_SW_bnd,& tv, aggregate_FW_forcing) @@ -865,7 +865,7 @@ subroutine convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, & integer :: is, ie, nz, i, k, k1, nzc, nkmb is = G%isc ; ie = G%iec ; nz = GV%ke - g_H2_2Rho0 = (GV%g_Earth * GV%H_to_Z**2) / (2.0 * US%kg_m3_to_R*GV%Rho0) + g_H2_2Rho0 = (GV%g_Earth * GV%H_to_Z**2) / (2.0 * GV%Rho0) nzc = nz ; if (present(nz_conv)) nzc = nz_conv nkmb = CS%nkml+CS%nkbl @@ -1068,7 +1068,7 @@ subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & Angstrom = GV%Angstrom_H C1_3 = 1.0/3.0 ; C1_6 = 1.0/6.0 - g_H2_2Rho0 = (GV%g_Earth * GV%H_to_Z**2) / (2.0 * US%kg_m3_to_R*GV%Rho0) + g_H2_2Rho0 = (GV%g_Earth * GV%H_to_Z**2) / (2.0 * GV%Rho0) Idt = 1.0 / dt_in_T is = G%isc ; ie = G%iec ; nz = GV%ke @@ -1611,7 +1611,7 @@ subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, & integer :: is, ie, nz, i, k, ks, itt, n C1_3 = 1.0/3.0 ; C1_6 = 1.0/6.0 ; C1_24 = 1.0/24.0 - g_H_2Rho0 = (GV%g_Earth * GV%H_to_Z) / (2.0 * US%kg_m3_to_R*GV%Rho0) + g_H_2Rho0 = (GV%g_Earth * GV%H_to_Z) / (2.0 * GV%Rho0) Hmix_min = CS%Hmix_min h_neglect = GV%H_subroundoff is = G%isc ; ie = G%iec ; nz = GV%ke @@ -2362,9 +2362,9 @@ subroutine mixedlayer_detrain_2(h, T, S, R0, Rcv, RcvTgt, dt_in_T, dt_diag, d_ea nkmb = CS%nkml+CS%nkbl h_neglect = GV%H_subroundoff g_2 = 0.5 * GV%g_Earth - Rho0xG = US%kg_m3_to_R*GV%Rho0 * GV%g_Earth + Rho0xG = GV%Rho0 * GV%g_Earth Idt_H2 = GV%H_to_Z**2 / dt_diag - I2Rho0 = 0.5 / (US%kg_m3_to_R*GV%Rho0) + I2Rho0 = 0.5 / (GV%Rho0) Angstrom = GV%Angstrom_H ! This is hard coding of arbitrary and dimensional numbers. @@ -2802,7 +2802,7 @@ subroutine mixedlayer_detrain_2(h, T, S, R0, Rcv, RcvTgt, dt_in_T, dt_diag, d_ea h_det_to_h2*( (R0(i,kb1)-R0_det)*h1 + (R0(i,kb2)-R0_det)*h2 ) + & h_ml_to_h2*( (R0(i,kb2)-R0(i,0))*h2 + (R0(i,kb1)-R0(i,0))*h1 + & (R0_det-R0(i,0))*h_det_to_h2 ) + & - h_det_to_h1*h_ml_to_h1*(R0_det-R0(i,0))) - 2.0*US%kg_m3_to_R*GV%Rho0*dPE_extrap ) + h_det_to_h1*h_ml_to_h1*(R0_det-R0(i,0))) - 2.0*GV%Rho0*dPE_extrap ) if (allocated(CS%diag_PE_detrain)) & CS%diag_PE_detrain(i,j) = CS%diag_PE_detrain(i,j) + s1en @@ -3163,7 +3163,7 @@ subroutine mixedlayer_detrain_1(h, T, S, R0, Rcv, RcvTgt, dt_in_T, dt_diag, d_ea "CS%nkbl must be 1 in mixedlayer_detrain_1.") dt_Time = dt_in_T / CS%BL_detrain_time - g_H2_2Rho0dt = (GV%g_Earth * GV%H_to_Z**2) / (2.0 * US%kg_m3_to_R*GV%Rho0 * dt_diag) + g_H2_2Rho0dt = (GV%g_Earth * GV%H_to_Z**2) / (2.0 * GV%Rho0 * dt_diag) g_H2_2dt = (GV%g_Earth * GV%H_to_Z**2) / (2.0 * dt_diag) ! Move detrained water into the buffer layer. diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 96652a9f45..b50011efed 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -750,7 +750,7 @@ subroutine diagnoseMLDbyDensityDifference(id_MLD, h, tv, densityDiff, G, GV, US, id_SQ = -1 ; if (PRESENT(id_MLDsq)) id_SQ = id_MLDsq - gE_rho0 = US%L_to_Z**2*GV%g_Earth / GV%Rho0 + gE_rho0 = US%L_to_Z**2*GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) dH_subML = 50.*GV%m_to_H ; if (present(dz_subML)) dH_subML = GV%Z_to_H*dz_subML is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -946,7 +946,7 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t if (present(cTKE)) cTKE(:,:,:) = 0.0 if (calculate_buoyancy) then SurfPressure(:) = 0.0 - GoRho = US%L_to_Z**2*GV%g_Earth / GV%Rho0 + GoRho = US%L_to_Z**2*GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) start = 1 + G%isc - G%isd npts = 1 + G%iec - G%isc endif @@ -1053,14 +1053,14 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t ! but do change answers. !----------------------------------------------------------------------------------------- if (calculate_buoyancy) then - call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, & + call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, & h2d, T2d, netMassInOut, netMassOut, netHeat, netSalt, & Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW=nonpenSW, & net_Heat_rate=netheat_rate, net_salt_rate=netsalt_rate, & netmassinout_rate=netmassinout_rate, pen_sw_bnd_rate=pen_sw_bnd_rate) else - call extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt, & + call extractFluxes1d(G, GV, US, fluxes, optics, nsw, j, dt, & H_limit_fluxes, CS%use_river_heat_content, CS%use_calving_heat_content, & h2d, T2d, netMassInOut, netMassOut, netHeat, netSalt, & Pen_SW_bnd, tv, aggregate_FW_forcing, nonpenSW=nonpenSW) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 5527866793..a99aa7c1e2 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -789,7 +789,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs max_itt = 20 h_tt_min = 0.0 - I_dtrho = 0.0 ; if (dt*GV%Rho0 > 0.0) I_dtrho = (US%Z_to_m**3*US%s_to_T**3) / (dt*GV%Rho0) + I_dtrho = 0.0 ; if (dt*GV%Rho0 > 0.0) I_dtrho = (US%Z_to_m**3*US%s_to_T**3) / (dt*US%R_to_kg_m3*GV%Rho0) vstar_unit_scale = US%m_to_Z * US%T_to_s MLD_guess = MLD_io @@ -863,9 +863,9 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs !/ Apply MStar to get mech_TKE if ((CS%answers_2018) .and. (CS%mstar_scheme==Use_Fixed_MStar)) then - mech_TKE = (dt*MSTAR_total*GV%Rho0) * u_star**3 + mech_TKE = (dt*MSTAR_total*US%R_to_kg_m3*GV%Rho0) * u_star**3 else - mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3) + mech_TKE = MSTAR_total * (dt*US%R_to_kg_m3*GV%Rho0* u_star**3) endif if (CS%TKE_diagnostics) then @@ -970,7 +970,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs ! on a curve fit from the data of Wang (GRL, 2003). ! Note: Ro = 1.0 / sqrt(0.5 * dt * Rho0 * (absf*htot)**3 / conv_PErel) nstar_FC = CS%nstar * conv_PErel / (conv_PErel + 0.2 * & - sqrt(0.5 * dt * GV%Rho0 * (absf*(htot*GV%H_to_Z))**3 * conv_PErel)) + sqrt(0.5 * dt * US%R_to_kg_m3*GV%Rho0 * (absf*(htot*GV%H_to_Z))**3 * conv_PErel)) endif if (debug) nstar_k(K) = nstar_FC diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index 5a9b9b5bbd..36066a20fb 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -185,7 +185,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) logical :: do_i(SZI_(G)), do_any integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke - G_Rho0 = (US%L_to_Z**2*GV%g_Earth) / GV%Rho0 + G_Rho0 = (US%L_to_Z**2*GV%g_Earth) / (US%R_to_kg_m3*GV%Rho0) ! Find the (limited) density jump across each interface. do i=is,ie @@ -403,7 +403,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) itide%h2(i,j) = min((max_frac_rough*G%bathyT(i,j))**2, itide%h2(i,j)) ! Compute the fixed part of internal tidal forcing; units are [J m-2] here. - CS%TKE_itidal_coef(i,j) = 0.5*kappa_h2_factor*GV%Rho0*& + CS%TKE_itidal_coef(i,j) = 0.5*kappa_h2_factor*US%R_to_kg_m3*GV%Rho0*& kappa_itides * US%Z_to_m**2*itide%h2(i,j) * itide%tideamp(i,j)**2 enddo ; enddo diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 437c52bd6d..3cc1e3b34d 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -793,7 +793,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & Ri_crit = CS%Rino_crit gR0 = GV%z_to_H*GV%H_to_Pa - g_R0 = (US%L_to_Z**2 * GV%g_Earth) / GV%Rho0 + g_R0 = (US%L_to_Z**2 * GV%g_Earth) / (US%R_to_kg_m3*GV%Rho0) k0dt = dt*CS%kappa_0 ! These are hard-coded for now. Perhaps these could be made dynamic later? ! tol_dksrc = 0.5*tol_ksrc_chg ; tol_dksrc_low = 1.0 - 1.0/tol_ksrc_chg ? diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 9a73801b1b..de312ce1c0 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -279,7 +279,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt_in_T, if (.not.associated(CS)) call MOM_error(FATAL,"set_diffusivity: "//& "Module must be initialized before it is used.") - I_Rho0 = 1.0 / GV%Rho0 + I_Rho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) ! ### Dimensional parameters if (CS%answers_2018) then kappa_dt_fill = US%m_to_Z**2 * 1.e-3 * 7200. @@ -509,7 +509,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt_in_T, CS%dissip_N0 + CS%dissip_N1 * sqrt(N2_lay(i,k)), & ! Floor aka Gargett CS%dissip_N2 * N2_lay(i,k)) ! Floor of Kd_min*rho0/F_Ri Kd_lay(i,j,k) = max(Kd_lay(i,j,k) , & ! Apply floor to Kd - dissip * (CS%FluxRi_max / (GV%Rho0 * (N2_lay(i,k) + Omega2)))) + dissip * (CS%FluxRi_max / (US%R_to_kg_m3*GV%Rho0 * (N2_lay(i,k) + Omega2)))) enddo ; enddo if (present(Kd_int)) then ; do K=2,nz ; do i=is,ie @@ -517,13 +517,13 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt_in_T, CS%dissip_N0 + CS%dissip_N1 * sqrt(N2_int(i,K)), & ! Floor aka Gargett CS%dissip_N2 * N2_int(i,K)) ! Floor of Kd_min*rho0/F_Ri Kd_int(i,j,K) = max(Kd_int(i,j,K) , & ! Apply floor to Kd - dissip * (CS%FluxRi_max / (GV%Rho0 * (N2_int(i,K) + Omega2)))) + dissip * (CS%FluxRi_max / (US%R_to_kg_m3*GV%Rho0 * (N2_int(i,K) + Omega2)))) enddo ; enddo ; endif endif if (associated(dd%Kd_work)) then do k=1,nz ; do i=is,ie - dd%Kd_Work(i,j,k) = GV%Rho0 * Kd_lay(i,j,k) * N2_lay(i,k) * & + dd%Kd_Work(i,j,k) = US%R_to_kg_m3*GV%Rho0 * Kd_lay(i,j,k) * N2_lay(i,k) * & GV%H_to_Z*h(i,j,k) ! Watt m-2 s = kg s-3 enddo ; enddo endif @@ -690,9 +690,9 @@ subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, & I_dt = 1.0 / dt Omega2 = CS%omega**2 H_neglect = GV%H_subroundoff - G_Rho0 = (US%L_to_Z**2 * GV%g_Earth) / GV%Rho0 + G_Rho0 = (US%L_to_Z**2 * GV%g_Earth) / (US%R_to_kg_m3*GV%Rho0) if (CS%answers_2018) then - I_Rho0 = 1.0 / GV%Rho0 + I_Rho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) G_IRho0 = (US%L_to_Z**2 * GV%g_Earth) * I_Rho0 else G_IRho0 = G_Rho0 @@ -890,7 +890,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, & integer :: i, k, is, ie, nz is = G%isc ; ie = G%iec ; nz = G%ke - G_Rho0 = (US%L_to_Z**2 * GV%g_Earth) / GV%Rho0 + G_Rho0 = (US%L_to_Z**2 * GV%g_Earth) / (US%R_to_kg_m3*GV%Rho0) H_neglect = GV%H_subroundoff ! Find the (limited) density jump across each interface. @@ -1177,8 +1177,8 @@ subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & TKE_Ray = 0.0 ; Rayleigh_drag = .false. if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) Rayleigh_drag = .true. - I_Rho0 = 1.0/GV%Rho0 - R0_g = GV%Rho0 / (US%L_to_Z**2 * GV%g_Earth) + I_Rho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) + R0_g = US%R_to_kg_m3*GV%Rho0 / (US%L_to_Z**2 * GV%g_Earth) do K=2,nz ; Rint(K) = 0.5*(US%R_to_kg_m3*GV%Rlay(k-1)+US%R_to_kg_m3*GV%Rlay(k)) ; enddo @@ -1394,7 +1394,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, & ! Determine whether to add Rayleigh drag contribution to TKE Rayleigh_drag = .false. if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) Rayleigh_drag = .true. - I_Rho0 = 1.0/GV%Rho0 + I_Rho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) cdrag_sqrt = sqrt(CS%cdrag) do i=G%isc,G%iec ! Developed in single-column mode @@ -1818,7 +1818,7 @@ subroutine set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0) enddo if (CS%bulkmixedlayer) then - g_R0 = GV%g_Earth / GV%Rho0 + g_R0 = GV%g_Earth / (US%R_to_kg_m3*GV%Rho0) kmb = GV%nk_rho_varies eps = 0.1 do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo @@ -2122,7 +2122,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ (CS%dissip_N0>0.) .or. (CS%dissip_Kd_min>0.) CS%dissip_N2 = 0.0 if (CS%FluxRi_max > 0.0) & - CS%dissip_N2 = CS%dissip_Kd_min * GV%Rho0 / CS%FluxRi_max + CS%dissip_N2 = CS%dissip_Kd_min * US%R_to_kg_m3*GV%Rho0 / CS%FluxRi_max CS%id_Kd_layer = register_diag_field('ocean_model', 'Kd_layer', diag%axesTL, Time, & 'Diapycnal diffusivity of layers (as set)', 'm2 s-1', & diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 02b5c9691d..00d964106d 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -269,7 +269,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize) Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB nkmb = GV%nk_rho_varies ; nkml = GV%nkml h_neglect = GV%H_subroundoff - Rho0x400_G = 400.0*(GV%Rho0 / (US%L_to_Z**2 * GV%g_Earth)) * GV%Z_to_H + Rho0x400_G = 400.0*(US%R_to_kg_m3*GV%Rho0 / (US%L_to_Z**2 * GV%g_Earth)) * GV%Z_to_H Vol_quit = 0.9*GV%Angstrom_H + h_neglect C2pi_3 = 8.0*atan(1.0)/3.0 @@ -1134,7 +1134,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt_in_T, G, GV, US, CS, sym Jsq = js-1 ; Isq = is-1 endif ; endif - Rho0x400_G = 400.0*(GV%Rho0/(US%L_to_Z**2 * GV%g_Earth)) * GV%Z_to_H + Rho0x400_G = 400.0*(US%R_to_kg_m3*GV%Rho0/(US%L_to_Z**2 * GV%g_Earth)) * GV%Z_to_H U_bg_sq = CS%drag_bg_vel * CS%drag_bg_vel cdrag_sqrt = sqrt(CS%cdrag) cdrag_sqrt_Z = US%L_to_Z * sqrt(CS%cdrag) @@ -1144,7 +1144,7 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt_in_T, G, GV, US, CS, sym dt_Rho0 = US%T_to_s*dt_in_T / GV%H_to_kg_m2 h_neglect = GV%H_subroundoff h_tiny = 2.0*GV%Angstrom_H + h_neglect - g_H_Rho0 = (GV%g_Earth*GV%H_to_Z) / GV%Rho0 + g_H_Rho0 = (GV%g_Earth*GV%H_to_Z) / (US%R_to_kg_m3*GV%Rho0) if (associated(forces%frac_shelf_u) .neqv. associated(forces%frac_shelf_v)) & call MOM_error(FATAL, "set_viscous_ML: one of forces%frac_shelf_u and "//& diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index fd910697af..aa158581fc 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -480,7 +480,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, CS) utide = CS%tideamp(i,j) ! Compute the fixed part of internal tidal forcing. ! The units here are [kg Z3 m-3 T-2 ~> J m-2 = kg s-2] here. - CS%TKE_itidal(i,j) = 0.5 * CS%kappa_h2_factor * GV%Rho0 * & + CS%TKE_itidal(i,j) = 0.5 * CS%kappa_h2_factor * US%R_to_kg_m3*GV%Rho0 * & CS%kappa_itides * CS%h2(i,j) * utide*utide enddo ; enddo @@ -1021,7 +1021,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, htot(i) = htot(i) + GV%H_to_Z*h(i,j,k) enddo ; enddo - I_Rho0 = 1.0/GV%Rho0 + I_Rho0 = 1.0 / (US%R_to_kg_m3*GV%Rho0) use_Polzin = ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09)) .or. & (CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09)) .or. & @@ -1255,7 +1255,7 @@ subroutine add_int_tide_diffusivity(h, N2_bot, j, TKE_to_Kd, max_TKE, G, GV, US, if (k 0) then !$OMP parallel do default(shared) private(trunc_any,CFL) diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90 index 0268c04f17..a5fc04fc06 100644 --- a/src/tracer/MOM_OCMIP2_CFC.F90 +++ b/src/tracer/MOM_OCMIP2_CFC.F90 @@ -460,9 +460,9 @@ subroutine OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS ! The -GV%Rho0 changes the sign convention of the flux and changes the units ! of the flux from [Conc. m s-1] to [Conc. kg m-2 s-1]. call coupler_type_extract_data(fluxes%tr_fluxes, CS%ind_cfc_11_flux, ind_flux, & - CFC11_flux, -GV%Rho0, idim=idim, jdim=jdim) + CFC11_flux, -G%US%R_to_kg_m3*GV%Rho0, idim=idim, jdim=jdim) call coupler_type_extract_data(fluxes%tr_fluxes, CS%ind_cfc_12_flux, ind_flux, & - CFC12_flux, -GV%Rho0, idim=idim, jdim=jdim) + CFC12_flux, -G%US%R_to_kg_m3*GV%Rho0, idim=idim, jdim=jdim) ! Use a tridiagonal solver to determine the concentrations after the ! surface source is applied and diapycnal advection and diffusion occurs. diff --git a/src/user/BFB_initialization.F90 b/src/user/BFB_initialization.F90 index fcfca47d50..546efcf0b9 100644 --- a/src/user/BFB_initialization.F90 +++ b/src/user/BFB_initialization.F90 @@ -56,14 +56,14 @@ subroutine BFB_set_coord(Rlay, g_prime, GV, US, param_file, eqn_of_state) "SST at the suothern edge of the domain.", units="C", default=20.0) call get_param(param_file, mdl, "T_BOT", T_bot, & "Bottom Temp", units="C", default=5.0) - rho_top = US%kg_m3_to_R*GV%rho0 + drho_dt*SST_s - rho_bot = US%kg_m3_to_R*GV%rho0 + drho_dt*T_bot + rho_top = GV%Rho0 + drho_dt*SST_s + rho_bot = GV%Rho0 + drho_dt*T_bot nz = GV%ke do k = 1,nz Rlay(k) = (rho_bot - rho_top)/(nz-1)*real(k-1) + rho_top if (k >1) then - g_prime(k) = (Rlay(k) - Rlay(k-1)) * GV%g_Earth / (US%kg_m3_to_R*GV%rho0) + g_prime(k) = (Rlay(k) - Rlay(k-1)) * GV%g_Earth / (GV%Rho0) else g_prime(k) = GV%g_Earth endif diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 77e0cb44c8..fa3a18b411 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -290,7 +290,7 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, US, param_file, tr_Reg) if (.not.associated(OBC)) return - g_prime_tot = (GV%g_Earth / GV%Rho0)*2.0 + g_prime_tot = (GV%g_Earth / (US%R_to_kg_m3*GV%Rho0))*2.0 Def_Rad = US%L_to_m*sqrt(D_edge*g_prime_tot) / (1.0e-4*US%T_to_s * 1000.0) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5e3*US%m_to_L*Def_Rad) * GV%Z_to_H diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 0da6285f37..a048d85491 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -1024,7 +1024,7 @@ subroutine get_StokesSL_LiFoxKemper(ustar, hbl, GV, US, UStokes_SL, LA) if (ustar > 0.0) then ! Computing u10 based on u_star and COARE 3.5 relationships - call ust_2_u10_coare3p5(US%Z_to_m*US%s_to_T*ustar*sqrt(GV%Rho0/1.225), u10, GV, US) + call ust_2_u10_coare3p5(US%Z_to_m*US%s_to_T*ustar*sqrt(US%R_to_kg_m3*GV%Rho0/1.225), u10, GV, US) ! surface Stokes drift UStokes = us_to_u10*u10 ! diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index b991fa95bc..2ef4dbd644 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -36,9 +36,10 @@ module Rossby_front_2d_initialization contains !> Initialization of thicknesses in 2D Rossby front test -subroutine Rossby_front_initialize_thickness(h, G, GV, param_file, just_read_params) - type(ocean_grid_type), intent(in) :: G !< Grid structure - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure +subroutine Rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read_params) + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -78,7 +79,7 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, param_file, just_read_par case (REGRIDDING_LAYER, REGRIDDING_RHO) do j = G%jsc,G%jec ; do i = G%isc,G%iec Dml = Hml( G, G%geoLatT(i,j) ) - eta = -( -dRho_DT / GV%Rho0 ) * Tz * 0.5 * ( Dml * Dml ) + eta = -( -dRho_DT / (US%R_to_kg_m3*GV%Rho0) ) * Tz * 0.5 * ( Dml * Dml ) stretch = ( ( G%max_depth + eta ) / G%max_depth ) h0 = ( G%max_depth / real(nz) ) * stretch do k = 1, nz @@ -89,7 +90,7 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, param_file, just_read_par case (REGRIDDING_ZSTAR, REGRIDDING_SIGMA) do j = G%jsc,G%jec ; do i = G%isc,G%iec Dml = Hml( G, G%geoLatT(i,j) ) - eta = -( -dRho_DT / GV%Rho0 ) * Tz * 0.5 * ( Dml * Dml ) + eta = -( -dRho_DT / (US%R_to_kg_m3*GV%Rho0) ) * Tz * 0.5 * ( Dml * Dml ) stretch = ( ( G%max_depth + eta ) / G%max_depth ) h0 = ( G%max_depth / real(nz) ) * stretch do k = 1, nz @@ -205,7 +206,7 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just do j = G%jsc,G%jec ; do I = G%isc-1,G%iec+1 f = 0.5* (G%CoriolisBu(I,j) + G%CoriolisBu(I,j-1) ) dUdT = 0.0 ; if (abs(f) > 0.0) & - dUdT = ( GV%g_Earth*dRho_dT ) / ( f * GV%Rho0 ) + dUdT = ( GV%g_Earth*dRho_dT ) / ( f * US%R_to_kg_m3*GV%Rho0 ) Dml = Hml( G, G%geoLatT(i,j) ) Ty = US%L_to_m*dTdy( G, T_range, G%geoLatT(i,j) ) zi = 0.