Skip to content

Commit

Permalink
+Add optional US arg to calculate_compress
Browse files Browse the repository at this point in the history 
  Added a new optional unit_scale_type argument to the calculate_compress to
trigger dimensional rescaling of their input and output variables.  Also use
this new argument in calls to calculate_compress from build_slight_column and
rescaled internal variables in the same routine.  All answers are bitwise
identical, but there is a new optional argument to a public interface.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 15, 2020
1 parent bb73bb8 commit 4182ad2
Show file tree
Hide file tree
Showing 2 changed files with 54 additions and 28 deletions.
30 changes: 15 additions & 15 deletions src/ALE/coord_slight.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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
Expand Down Expand Up @@ -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)) + &
Expand Down
52 changes: 39 additions & 13 deletions src/equation_of_state/MOM_EOS.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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, &
Expand All @@ -1026,26 +1041,37 @@ 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

if (.not.associated(EOS)) call MOM_error(FATAL, &
"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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 4182ad2

Please sign in to comment.