Skip to content

Commit

Permalink
+Rescaled pressure arguments to multiple routines
Browse files Browse the repository at this point in the history 
  Rescaled the reference pressure arguments to 3 set_coord routines, 5
initialization routines, and kappa_shear_column.  Also removed the unused pres
argument to convert_temp_salt_for_TEOS10 and replaced its ocean_grid_type
argument with a hor_index_type.  All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 12, 2020
1 parent 3223eb6 commit fb820c1
Show file tree
Hide file tree
Showing 7 changed files with 101 additions and 83 deletions.
27 changes: 13 additions & 14 deletions src/equation_of_state/MOM_EOS.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3129,19 +3129,16 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref,
end subroutine int_spec_vol_dp_generic_plm

!> Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10
subroutine convert_temp_salt_for_TEOS10(T, S, press, G, kd, mask_z, EOS)
use MOM_grid, only : ocean_grid_type

type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
subroutine convert_temp_salt_for_TEOS10(T, S, HI, kd, mask_z, EOS)
integer, intent(in) :: kd !< The number of layers to work on
type(hor_index_type), intent(in) :: HI !< The horizontal index structure
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
intent(inout) :: T !< Potential temperature referenced to the surface [degC]
real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
intent(inout) :: S !< Salinity [ppt]
real, dimension(:), intent(in) :: press !< Pressure at the top of the layer [Pa].
type(EOS_type), pointer :: EOS !< Equation of state structure
real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed,kd), &
intent(in) :: mask_z !< 3d mask regulating which points to convert.
integer, intent(in) :: kd !< The number of layers to work on
type(EOS_type), pointer :: EOS !< Equation of state structure

integer :: i,j,k
real :: gsw_sr_from_sp, gsw_ct_from_pt, gsw_sa_from_sp
Expand All @@ -3152,12 +3149,14 @@ subroutine convert_temp_salt_for_TEOS10(T, S, press, G, kd, mask_z, EOS)

if ((EOS%form_of_EOS /= EOS_TEOS10) .and. (EOS%form_of_EOS /= EOS_NEMO)) return

do k=1,kd ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
do k=1,kd ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec
if (mask_z(i,j,k) >= 1.0) then
S(i,j,k) = gsw_sr_from_sp(S(i,j,k))
! p=press(k)/10000. !convert pascal to dbar
! S(i,j,k) = gsw_sa_from_sp(S(i,j,k),p,G%geoLonT(i,j),G%geoLatT(i,j))
T(i,j,k) = gsw_ct_from_pt(S(i,j,k),T(i,j,k))
! Get absolute salnity from practical salinity, converting pressures from Pascal to dbar.
! If this option is activated, pressure will need to be added as an argument, and it should be
! moved out into module that is not shared between components, where the ocean_grid can be used.
! S(i,j,k) = gsw_sa_from_sp(S(i,j,k),pres(i,j,k)*1.0e-4,G%geoLonT(i,j),G%geoLatT(i,j))
T(i,j,k) = gsw_ct_from_pt(S(i,j,k), T(i,j,k))
endif
enddo ; enddo ; enddo
end subroutine convert_temp_salt_for_TEOS10
Expand Down
37 changes: 20 additions & 17 deletions src/initialization/MOM_coord_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,11 +89,11 @@ subroutine MOM_initialize_coord(GV, US, PF, write_geom, output_dir, tv, max_dept
case ("linear")
call set_coord_linear(GV%Rlay, GV%g_prime, GV, US, PF)
case ("ts_ref")
call set_coord_from_ts_ref(GV%Rlay, GV%g_prime, GV, US, PF, eos, US%RL2_T2_to_Pa*tv%P_Ref)
call set_coord_from_TS_ref(GV%Rlay, GV%g_prime, GV, US, PF, eos, tv%P_Ref)
case ("ts_profile")
call set_coord_from_TS_profile(GV%Rlay, GV%g_prime, GV, US, PF, eos, US%RL2_T2_to_Pa*tv%P_Ref)
call set_coord_from_TS_profile(GV%Rlay, GV%g_prime, GV, US, PF, eos, tv%P_Ref)
case ("ts_range")
call set_coord_from_TS_range(GV%Rlay, GV%g_prime, GV, US, PF, eos, US%RL2_T2_to_Pa*tv%P_Ref)
call set_coord_from_TS_range(GV%Rlay, GV%g_prime, GV, US, PF, eos, tv%P_Ref)
case ("file")
call set_coord_from_file(GV%Rlay, GV%g_prime, GV, US, PF)
case ("USER")
Expand Down Expand Up @@ -198,8 +198,7 @@ subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file)
end subroutine set_coord_from_layer_density

!> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a profile of g'.
subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state, &
P_Ref)
subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref)
real, dimension(:), intent(out) :: Rlay !< The layers' target coordinate values
!! (potential density) [R ~> kg m-3].
real, dimension(:), intent(out) :: g_prime !< The reduced gravity across the interfaces
Expand All @@ -209,7 +208,8 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
!! parameters
type(EOS_type), pointer :: eqn_of_state !< integer selecting the equation of state.
real, intent(in) :: P_Ref !< The coordinate-density reference pressure [Pa].
real, intent(in) :: P_Ref !< The coordinate-density reference pressure
!! [R L2 T-2 ~> Pa].
! Local variables
real :: T_ref ! Reference temperature
real :: S_ref ! Reference salinity
Expand Down Expand Up @@ -240,7 +240,8 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state
! The uppermost layer's density is set here. Subsequent layers' !
! densities are determined from this value and the g values. !
! T0 = 28.228 ; S0 = 34.5848 ; Pref = P_Ref
call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)

! These statements set the layer densities. !
do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo
Expand All @@ -249,8 +250,7 @@ subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state
end subroutine set_coord_from_TS_ref

!> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a T-S profile.
subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, &
eqn_of_state, P_Ref)
subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref)
real, dimension(:), intent(out) :: Rlay !< The layers' target coordinate values
!! (potential density) [R ~> kg m-3].
real, dimension(:), intent(out) :: g_prime !< The reduced gravity across the interfaces
Expand All @@ -260,7 +260,9 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, &
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
!! parameters
type(EOS_type), pointer :: eqn_of_state !< integer that selects equation of state.
real, intent(in) :: P_Ref !< The coordinate-density reference pressure [Pa].
real, intent(in) :: P_Ref !< The coordinate-density reference pressure
!! [R L2 T-2 ~> Pa].

! Local variables
real, dimension(GV%ke) :: T0, S0, Pref
real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2].
Expand Down Expand Up @@ -289,16 +291,15 @@ subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, &
" set_coord_from_TS_profile: Unable to open " //trim(filename))
! These statements set the interface reduced gravities. !
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=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, pres_scale=US%RL2_T2_to_Pa)
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

!> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a linear T-S profile.
subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, &
eqn_of_state, P_Ref)
subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref)
real, dimension(:), intent(out) :: Rlay !< The layers' target coordinate values
!! (potential density) [R ~> kg m-3].
real, dimension(:), intent(out) :: g_prime !< The reduced gravity across the interfaces
Expand All @@ -308,7 +309,8 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, &
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
!! parameters
type(EOS_type), pointer :: eqn_of_state !< integer that selects equation of state
real, intent(in) :: P_Ref !< The coordinate-density reference pressure [Pa]
real, intent(in) :: P_Ref !< The coordinate-density reference pressure
!! [R L2 T-2 ~> Pa].

! Local variables
real, dimension(GV%ke) :: T0, S0, Pref
Expand Down Expand Up @@ -369,8 +371,9 @@ subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, &
enddo

g_prime(1) = g_fs
do k=1,nz ; Pref(k) = P_ref ; enddo
call calculate_density(T0, S0, Pref, Rlay, k_light, nz-k_light+1, eqn_of_state, scale=US%kg_m3_to_R)
do k=1,nz ; Pref(k) = P_Ref ; enddo
call calculate_density(T0, S0, Pref, Rlay, k_light, nz-k_light+1, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
! Extrapolate target densities for the variable density mixed and buffer layers.
do k=k_light-1,1,-1
Rlay(k) = 2.0*Rlay(k+1) - Rlay(k+2)
Expand Down
35 changes: 20 additions & 15 deletions src/initialization/MOM_state_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -291,9 +291,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
case ("ISOMIP"); call ISOMIP_initialize_thickness(h, G, GV, US, PF, tv, &
just_read_params=just_read)
case ("benchmark"); call benchmark_initialize_thickness(h, G, GV, US, PF, &
tv%eqn_of_state, US%RL2_T2_to_Pa*tv%P_Ref, just_read_params=just_read)
tv%eqn_of_state, tv%P_Ref, just_read_params=just_read)
case ("Neverland"); call Neverland_initialize_thickness(h, G, GV, US, PF, &
tv%eqn_of_state, US%RL2_T2_to_Pa*tv%P_Ref)
tv%eqn_of_state, tv%P_Ref)
case ("search"); call initialize_thickness_search
case ("circle_obcs"); call circle_obcs_initialize_thickness(h, G, GV, PF, &
just_read_params=just_read)
Expand Down Expand Up @@ -348,11 +348,11 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
! " \t baroclinic_zone - an analytic baroclinic zone. \n"//&
select case (trim(config))
case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, G, GV, US, PF, &
eos, US%RL2_T2_to_Pa*tv%P_Ref, just_read_params=just_read)
eos, tv%P_Ref, just_read_params=just_read)
case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, G, &
PF, just_read_params=just_read)
case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, &
G, GV, US, PF, eos, US%RL2_T2_to_Pa*tv%P_Ref, just_read_params=just_read)
G, GV, US, PF, eos, tv%P_Ref, just_read_params=just_read)
case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, &
G, PF, just_read_params=just_read)
case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, G, PF, &
Expand Down Expand Up @@ -1561,15 +1561,15 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
!! parameters.
type(EOS_type), pointer :: eqn_of_state !< Integer that selects the equatio of state.
real, intent(in) :: P_Ref !< The coordinate-density reference pressure
!! [Pa].
!! [R L2 T-2 ~> Pa].
logical, optional, intent(in) :: just_read_params !< If present and true, this call will
!! only read parameters without changing h.
! Local variables
real :: T0(SZK_(G)) ! Layer potential temperatures [degC]
real :: S0(SZK_(G)) ! Layer salinities [degC]
real :: T_Ref ! Reference Temperature [degC]
real :: S_Ref ! Reference Salinity [ppt]
real :: pres(SZK_(G)) ! An array of the reference pressure [Pa].
real :: pres(SZK_(G)) ! An array of the reference pressure [R L2 T-2 ~> Pa].
real :: drho_dT(SZK_(G)) ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
real :: drho_dS(SZK_(G)) ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 [R ~> kg m-3].
Expand Down Expand Up @@ -1601,8 +1601,10 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
T0(k) = T_Ref
enddo

call calculate_density(T0(1),S0(1),pres(1),rho_guess(1),eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,1,eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, 1, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)

if (fit_salin) then
! A first guess of the layers' temperatures.
Expand All @@ -1611,8 +1613,10 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
enddo
! Refine the guesses for each layer.
do itt=1,6
call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
do k=1,nz
S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS(k))
enddo
Expand All @@ -1623,8 +1627,10 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P
T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT(1)
enddo
do itt=1,6
call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state, scale=US%kg_m3_to_R)
call calculate_density(T0, S0, pres, rho_guess, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, 1, nz, eqn_of_state, &
scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa)
do k=1,nz
T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k)
enddo
Expand Down Expand Up @@ -2185,8 +2191,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param
allocate(frac_shelf_h(isd:ied,jsd:jed))

! Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO
press(:) = US%RL2_T2_to_Pa*tv%P_Ref
call convert_temp_salt_for_TEOS10(temp_z, salt_z, press, G, kd, mask_z, eos)
call convert_temp_salt_for_TEOS10(temp_z, salt_z, G%HI, kd, mask_z, eos)

press(:) = tv%P_Ref
do k=1,kd ; do j=js,je
Expand Down Expand Up @@ -2397,7 +2402,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param

if (adjust_temperature .and. .not. useALEremapping) then
call determine_temperature(tv%T(is:ie,js:je,:), tv%S(is:ie,js:je,:), &
GV%Rlay(1:nz), US%RL2_T2_to_Pa*tv%P_Ref, niter, missing_value, h(is:ie,js:je,:), ks, US, eos)
GV%Rlay(1:nz), tv%P_Ref, niter, missing_value, h(is:ie,js:je,:), ks, US, eos)
endif

deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z)
Expand Down
Loading

0 comments on commit fb820c1

Please sign in to comment.