diff --git a/src/ALE/coord_slight.F90 b/src/ALE/coord_slight.F90 index 89c78d6c64..000315bae8 100644 --- a/src/ALE/coord_slight.F90 +++ b/src/ALE/coord_slight.F90 @@ -205,7 +205,7 @@ subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, & logical, dimension(nz+1) :: reliable ! If true, this interface is in a reliable position. real, dimension(nz+1) :: T_int, S_int ! Temperature [degC] and salinity [ppt] interpolated to interfaces. real, dimension(nz+1) :: rho_tmp ! A temporary density [R ~> kg m-3] - real, dimension(nz+1) :: drho_dp ! The partial derivative of density with pressure [kg m-3 Pa-1] + real, dimension(nz+1) :: drho_dp ! The partial derivative of density with pressure [T2 L-2 ~> kg m-3 Pa-1] real, dimension(nz+1) :: p_IS, p_R ! Pressures [R L2 T-2 ~> Pa] real, dimension(nz+1) :: drhoIS_dT ! The partial derivative of in situ density with temperature ! in [R degC-1 ~> kg m-3 degC-1] @@ -216,19 +216,20 @@ subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, & real, dimension(nz+1) :: drhoR_dS ! The partial derivative of reference density with salinity ! in [R ppt-1 ~> kg m-3 ppt-1] real, dimension(nz+1) :: strat_rat - real :: H_to_cPa + real :: H_to_cPa ! A conversion factor from thicknesses to the compressibility fraction times + ! the units of pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] real :: drIS, drR ! In situ and reference density differences [R ~> kg m-3] - real :: Fn_now, I_HStol, Fn_zero_val - real :: z_int_unst - real :: dz ! A uniform layer thickness in very shallow water [H ~> m or kg m-2]. - real :: dz_ur ! The total thickness of an unstable region [H ~> m or kg m-2]. + real :: Fn_now, I_HStol, Fn_zero_val ! Nondimensional variables [nondim] + real :: z_int_unst ! The depth where the stratification allows the interior grid to start [H ~> m or kg m-2] + real :: dz ! A uniform layer thickness in very shallow water [H ~> m or kg m-2]. + real :: dz_ur ! The total thickness of an unstable region [H ~> m or kg m-2]. real :: wgt, cowgt ! A weight and its complement [nondim]. - real :: rho_ml_av ! The average potential density in a near-surface region [R ~> kg m-3]. - real :: H_ml_av ! A thickness to try to use in taking the near-surface average [H ~> m or kg m-2]. - real :: rho_x_z ! A cumulative integral of a density [R H ~> kg m-2 or kg2 m-5]. - real :: z_wt ! The thickness actually used in taking the near-surface average [H ~> m or kg m-2]. - real :: k_interior ! The (real) value of k where the interior grid starts. - real :: k_int2 ! The (real) value of k where the interior grid starts. + real :: rho_ml_av ! The average potential density in a near-surface region [R ~> kg m-3]. + real :: H_ml_av ! A thickness to try to use in taking the near-surface average [H ~> m or kg m-2]. + real :: rho_x_z ! A cumulative integral of a density [R H ~> kg m-2 or kg2 m-5]. + real :: z_wt ! The thickness actually used in taking the near-surface average [H ~> m or kg m-2]. + real :: k_interior ! The (real) value of k where the interior grid starts [nondim]. + real :: k_int2 ! The (real) value of k where the interior grid starts [nondim]. real :: z_interior ! The depth where the interior grid starts [H ~> m or kg m-2]. real :: z_ml_fix ! The depth at which the fixed-thickness near-surface layers end [H ~> m or kg m-2]. real :: dz_dk ! The thickness of layers between the fixed-thickness @@ -378,13 +379,12 @@ subroutine build_slight_column(CS, US, eqn_of_state, H_to_pres, H_subroundoff, & call calculate_density_derivs(T_int, S_int, p_R, drhoR_dT, drhoR_dS, 2, nz-1, & eqn_of_state, US) if (CS%compressibility_fraction > 0.0) then - call calculate_compress(T_int, S_int, US%RL2_T2_to_Pa*p_R(:), rho_tmp, drho_dp, 2, nz-1, & - eqn_of_state) + call calculate_compress(T_int, S_int, p_R(:), rho_tmp, drho_dp, 2, nz-1, eqn_of_state, US) else do K=2,nz ; drho_dp(K) = 0.0 ; enddo endif - H_to_cPa = CS%compressibility_fraction * H_to_pres * US%L_T_to_m_s**2 + H_to_cPa = CS%compressibility_fraction * H_to_pres strat_rat(1) = 1.0 do K=2,nz drIS = drhoIS_dT(K) * (T_f(k) - T_f(k-1)) + & diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index c1fd5fd42f..4730d92807 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -394,6 +394,7 @@ subroutine calc_spec_vol_US(T, S, pressure, specvol, start, npts, EOS, US, spv_r real, optional, intent(in) :: scale !< A multiplicative factor by which to scale specific !! volume in combination with scaling given by US [various] + ! Local variables real, dimension(size(pressure)) :: pres ! Pressure converted to [Pa] real :: spv_reference ! spv_ref converted to [m3 kg-1] real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg R-1 m-3 ~> 1] @@ -995,21 +996,35 @@ subroutine calc_spec_vol_derivs_HI_1d(T, S, pressure, dSV_dT, dSV_dS, HI, EOS, U end subroutine calc_spec_vol_derivs_HI_1d -!> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs. -subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, EOS) +!> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array +!! inputs. If US is present, the units of the inputs and outputs are rescaled. +subroutine calculate_compress_array(T, S, press, rho, drho_dp, start, npts, EOS, US) real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] real, dimension(:), intent(in) :: S !< Salinity [PSU] - real, dimension(:), intent(in) :: pressure !< Pressure [Pa] - real, dimension(:), intent(inout) :: rho !< In situ density [kg m-3] + real, dimension(:), intent(in) :: press !< Pressure [Pa] or [R L2 T-2 ~> Pa] + real, dimension(:), intent(inout) :: rho !< In situ density [kg m-3] or [R ~> kg m-3] real, dimension(:), intent(inout) :: drho_dp !< The partial derivative of density with pressure - !! (also the inverse of the square of sound speed) [s2 m-2] + !! (also the inverse of the square of sound speed) + !! [s2 m-2] or [T2 L-2] integer, intent(in) :: start !< Starting index within the array integer, intent(in) :: npts !< The number of values to calculate type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + + ! Local variables + real, dimension(size(press)) :: pressure ! Pressure converted to [Pa] + integer :: i, is, ie if (.not.associated(EOS)) call MOM_error(FATAL, & "calculate_compress called with an unassociated EOS_type EOS.") + is = start ; ie = is + npts - 1 + if (present(US)) then + do i=is,ie ; pressure(i) = US%RL2_T2_to_Pa * press(i) ; enddo + else + do i=is,ie ; pressure(i) = press(i) ; enddo + endif + select case (EOS%form_of_EOS) case (EOS_LINEAR) call calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts, & @@ -1026,18 +1041,29 @@ subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, E call MOM_error(FATAL, "calculate_compress: EOS%form_of_EOS is not valid.") end select + if (present(US)) then + if (US%kg_m3_to_R /= 1.0) then ; do i=is,ie + rho(i) = US%kg_m3_to_R * rho(i) + enddo ; endif + if (US%L_T_to_m_s /= 1.0) then ; do i=is,ie + drho_dp(i) = US%L_T_to_m_s**2 * drho_dp(i) + enddo ; endif + endif + end subroutine calculate_compress_array !> Calculate density and compressibility for a scalar. This just promotes the scalar to an array -!! with a singleton dimension and calls calculate_compress_array -subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS) +!! with a singleton dimension and calls calculate_compress_array. If US is present, the units of +!! the inputs and outputs are rescaled. +subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS, US) real, intent(in) :: T !< Potential temperature referenced to the surface [degC] real, intent(in) :: S !< Salinity [ppt] - real, intent(in) :: pressure !< Pressure [Pa] - real, intent(out) :: rho !< In situ density [kg m-3] - real, intent(out) :: drho_dp !< The partial derivative of density with pressure - !! (also the inverse of the square of sound speed) [s2 m-2] + real, intent(in) :: pressure !< Pressure [Pa] or [R L2 T-2 ~> Pa] + real, intent(out) :: rho !< In situ density [kg m-3] or [R ~> kg m-3] + real, intent(out) :: drho_dp !< The partial derivative of density with pressure (also the + !! inverse of the square of sound speed) [s2 m-2] or [T2 L-2] type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type real, dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa @@ -1045,7 +1071,7 @@ subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS) "calculate_compress called with an unassociated EOS_type EOS.") Ta(1) = T ; Sa(1) = S; pa(1) = pressure - call calculate_compress_array(Ta, Sa, pa, rhoa, drho_dpa, 1, 1, EOS) + call calculate_compress_array(Ta, Sa, pa, rhoa, drho_dpa, 1, 1, EOS, US) rho = rhoa(1) ; drho_dp = drho_dpa(1) end subroutine calculate_compress_scalar @@ -1076,7 +1102,7 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(inout) :: dza !< The change in the geopotential anomaly across - !! the layer [T-2 ~> m2 s-2] or [m2 s-2] + !! the layer [L2 T-2 ~> m2 s-2] or [m2 s-2] real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & optional, intent(inout) :: intp_dza !< The integral in pressure through the layer of the !! geopotential anomaly relative to the anomaly at the bottom of the