diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index cfd286450e..fb2a1b1ca6 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -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 @@ -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 diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index d4afa115af..691ca4a60c 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -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") @@ -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 @@ -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 @@ -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 @@ -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 @@ -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]. @@ -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 @@ -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 @@ -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) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 6339211d3e..d58a6b4704 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -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) @@ -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, & @@ -1561,7 +1561,7 @@ 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 @@ -1569,7 +1569,7 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P 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]. @@ -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. @@ -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 @@ -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 @@ -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 @@ -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) diff --git a/src/parameterizations/vertical/MOM_kappa_shear.F90 b/src/parameterizations/vertical/MOM_kappa_shear.F90 index 6d773c67a0..781085e794 100644 --- a/src/parameterizations/vertical/MOM_kappa_shear.F90 +++ b/src/parameterizations/vertical/MOM_kappa_shear.F90 @@ -160,7 +160,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & kappa_avg, & ! The time-weighted average of kappa [Z2 T-1 ~> m2 s-1]. tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2]. real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2]. - real :: surface_pres ! The top surface pressure [Pa]. + real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa]. real :: dz_in_lay ! The running sum of the thickness in a layer [Z ~> m]. real :: k0dt ! The background diffusivity times the timestep [Z2 ~> m2]. @@ -283,7 +283,7 @@ subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, & endif f2 = 0.25 * ((G%CoriolisBu(I,j)**2 + G%CoriolisBu(I-1,J-1)**2) + & (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2)) - surface_pres = 0.0 ; if (associated(p_surf)) surface_pres = US%RL2_T2_to_Pa*p_surf(i,j) + surface_pres = 0.0 ; if (associated(p_surf)) surface_pres = p_surf(i,j) ! ---------------------------------------------------- I_Ld2_1d, dz_Int_1d @@ -430,7 +430,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ kappa_avg, & ! The time-weighted average of kappa [Z2 T-1 ~> m2 s-1]. tke_avg ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2]. real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2]. - real :: surface_pres ! The top surface pressure [Pa]. + real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa]. real :: dz_in_lay ! The running sum of the thickness in a layer [Z ~> m]. real :: k0dt ! The background diffusivity times the timestep [Z2 ~> m2]. @@ -585,8 +585,8 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ endif f2 = G%CoriolisBu(I,J)**2 surface_pres = 0.0 ; if (associated(p_surf)) & - surface_pres = 0.25 * US%RL2_T2_to_Pa*((p_surf(i,j) + p_surf(i+1,j+1)) + & - (p_surf(i+1,j) + p_surf(i,j+1))) + surface_pres = 0.25 * ((p_surf(i,j) + p_surf(i+1,j+1)) + & + (p_surf(i+1,j) + p_surf(i,j+1))) ! ---------------------------------------------------- ! Set the initial guess for kappa, here defined at interfaces. @@ -661,7 +661,7 @@ subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_ if (CS%debug) then call hchksum(kappa_io, "kappa", G%HI, scale=US%Z2_T_to_m2_s) - call Bchksum(tke_io, "tke", G%HI) + call Bchksum(tke_io, "tke", G%HI, scale=US%Z_to_m**2*US%s_to_T**2) endif if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag) @@ -679,7 +679,6 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & dz, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, & tke_avg, tv, CS, GV, US, I_Ld2_1d, dz_Int_1d) 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, dimension(SZK_(GV)+1), & intent(inout) :: kappa !< The time-weighted average of kappa [Z2 T-1 ~> m2 s-1]. real, dimension(SZK_(GV)+1), & @@ -687,7 +686,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & !! an interface [Z2 T-2 ~> m2 s-2]. integer, intent(in) :: nzc !< The number of active layers in the column. real, intent(in) :: f2 !< The square of the Coriolis parameter [T-2 ~> s-2]. - real, intent(in) :: surface_pres !< The surface pressure [Pa]. + real, intent(in) :: surface_pres !< The surface pressure [R L2 T-2 ~> Pa]. real, dimension(SZK_(GV)), & intent(in) :: dz !< The layer thickness [Z ~> m]. real, dimension(SZK_(GV)), & @@ -708,6 +707,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & !! have NULL ptrs. type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous !! call to kappa_shear_init. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZK_(GV)+1), & optional, intent(out) :: I_Ld2_1d !< The inverse of the squared mixing length [Z-2 ~> m-2]. real, dimension(SZK_(GV)+1), & @@ -741,15 +741,15 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & kappa_mid, & ! The average of the initial and predictor estimates of kappa [Z2 T-1 ~> m2 s-1]. tke_pred, & ! The value of TKE from a predictor step [Z2 T-2 ~> m2 s-2]. kappa_pred, & ! The value of kappa from a predictor step [Z2 T-1 ~> m2 s-1]. - pressure, & ! The pressure at an interface [Pa]. + pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa]. T_int, & ! The temperature interpolated to an interface [degC]. Sal_int, & ! The salinity interpolated to an interface [ppt]. dbuoy_dT, & ! The partial derivatives of buoyancy with changes in temperature dbuoy_dS, & ! and salinity, [Z T-2 degC-1 ~> m s-2 degC-1] and [Z T-2 ppt-1 ~> m s-2 ppt-1]. I_L2_bdry, & ! The inverse of the square of twice the harmonic mean ! distance to the top and bottom boundaries [Z-2 ~> m-2]. - K_Q, & ! Diffusivity divided by TKE [Z2 m-2 s2 T-1 ~> s]. - K_Q_tmp, & ! A temporary copy of diffusivity divided by TKE [Z2 m-2 s2 T-1 ~> s]. + K_Q, & ! Diffusivity divided by TKE [T ~> s]. + K_Q_tmp, & ! A temporary copy of diffusivity divided by TKE [T ~> s]. local_src_avg, & ! The time-integral of the local source [nondim]. tol_min, & ! Minimum tolerated ksrc for the corrector step [T-1 ~> s-1]. tol_max, & ! Maximum tolerated ksrc for the corrector step [T-1 ~> s-1]. @@ -762,8 +762,8 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & real :: b1 ! The inverse of the pivot in the tridiagonal equations. real :: bd1 ! A term in the denominator of b1. real :: d1 ! 1 - c1 in the tridiagonal equations. - real :: gR0 ! A conversion factor from Z to Pa equal to Rho_0 times g - ! [Pa Z-1 = kg m-1 s-2 Z-1 ~> kg m-2 s-2]. + real :: gR0 ! A conversion factor from Z to pressure, given by Rho_0 times g + ! [R L2 T-2 Z-1 ~> kg m-2 s-2]. real :: g_R0 ! g_R0 is a rescaled version of g/Rho [Z R-1 T-2 ~> m4 kg-1 s-2]. real :: Norm ! A factor that normalizes two weights to 1 [Z-2 ~> m-2]. real :: tol_dksrc ! Tolerance for the change in the kappa source within an iteration @@ -813,7 +813,7 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & #endif Ri_crit = CS%Rino_crit - gR0 = GV%z_to_H*GV%H_to_Pa + gR0 = GV%Rho0 * GV%g_Earth g_R0 = (US%L_to_Z**2 * GV%g_Earth) / (GV%Rho0) k0dt = dt*CS%kappa_0 @@ -910,8 +910,8 @@ subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, & T_int(K) = 0.5*(T(k-1) + T(k)) Sal_int(K) = 0.5*(Sal(k-1) + Sal(k)) enddo - call calculate_density_derivs(T_int, Sal_int, pressure, dbuoy_dT, & - dbuoy_dS, 2, nzc-1, tv%eqn_of_state, scale=-g_R0*US%kg_m3_to_R) + call calculate_density_derivs(T_int, Sal_int, pressure, dbuoy_dT, dbuoy_dS, 2, nzc-1, & + tv%eqn_of_state, scale=-g_R0*US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) else do K=1,nzc+1 ; dbuoy_dT(K) = -g_R0 ; dbuoy_dS(K) = 0.0 ; enddo endif @@ -1388,7 +1388,7 @@ subroutine find_kappa_tke(N2, S2, kappa_in, Idz, dz_Int, I_L2_bdry, f2, & type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(nz+1), intent(inout) :: K_Q !< The shear-driven diapycnal diffusivity divided by !! the turbulent kinetic energy per unit mass at - !! interfaces [Z2 m-2 s2 T-1 ~> s]. + !! interfaces [T ~> s]. real, dimension(nz+1), intent(out) :: tke !< The turbulent kinetic energy per unit mass at !! interfaces [Z2 T-2 ~> m2 s-2]. real, dimension(nz+1), intent(out) :: kappa !< The diapycnal diffusivity at interfaces diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index aaa670070b..948705f4b3 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -744,7 +744,7 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h, real, dimension(:,:,:), intent(inout) :: temp !< potential temperature [degC] real, dimension(:,:,:), intent(inout) :: salt !< salinity [PSU] real, dimension(size(temp,3)), intent(in) :: R_tgt !< desired potential density [R ~> kg m-3]. - real, intent(in) :: p_ref !< reference pressure [Pa]. + real, intent(in) :: p_ref !< reference pressure [R L2 T-2 ~> Pa]. integer, intent(in) :: niter !< maximum number of iterations integer, intent(in) :: k_start !< starting index (i.e. below the buffer layer) real, intent(in) :: land_fill !< land fill value @@ -763,7 +763,7 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h, hin, & ! Input layer thicknesses [H ~> m or kg m-2] drho_dT, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1] drho_dS ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1] - real, dimension(size(temp,1)) :: press + real, dimension(size(temp,1)) :: press ! Reference pressures [R L2 T-2 ~> Pa] integer :: nx, ny, nz, nt, i, j, k, n, itt real :: dT_dS_gauge ! The relative penalizing of temperature to salinity changes when ! minimizing property changes while correcting density [degC ppt-1]. @@ -801,9 +801,10 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h, adjust_salt = .true. iter_loop: do itt = 1,niter do k=1, nz - call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, scale=US%kg_m3_to_R) + call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, & + scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, & - eos, scale=US%kg_m3_to_R) + eos, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) enddo do k=k_start,nz ; do i=1,nx @@ -831,9 +832,10 @@ subroutine determine_temperature(temp, salt, R_tgt, p_ref, niter, land_fill, h, if (adjust_salt .and. old_fit) then ; do itt = 1,niter do k=1, nz - call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, scale=US%kg_m3_to_R) + call calculate_density(T(:,k), S(:,k), press, rho(:,k), 1, nx, eos, & + scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) call calculate_density_derivs(T(:,k), S(:,k), press, drho_dT(:,k), drho_dS(:,k), 1, nx, & - eos, scale=US%kg_m3_to_R) + eos, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) enddo do k=k_start,nz ; do i=1,nx ! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln ) then diff --git a/src/user/Neverland_initialization.F90 b/src/user/Neverland_initialization.F90 index 949530e773..64afe85ab5 100644 --- a/src/user/Neverland_initialization.F90 +++ b/src/user/Neverland_initialization.F90 @@ -122,7 +122,7 @@ subroutine Neverland_initialize_thickness(h, G, GV, US, param_file, eqn_of_state type(EOS_type), pointer :: eqn_of_state !< integer that selects the !! equation of state. real, intent(in) :: P_Ref !< The coordinate-density - !! reference pressure [Pa]. + !! reference pressure [R L2 T-2 ~> Pa]. ! Local variables real :: e0(SZK_(G)+1) ! The resting interface heights, in depth units [Z ~> m], ! usually negative because it is positive upward. diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index 5641035ded..a173150b9d 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -83,7 +83,7 @@ end subroutine benchmark_initialize_topography !! temperature profile with an exponentially decaying thermocline on top of a !! linear stratification. subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, & - P_ref, just_read_params) + P_Ref, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid 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 @@ -94,7 +94,7 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state type(EOS_type), pointer :: eqn_of_state !< integer that selects the !! equation of state. real, intent(in) :: P_Ref !< The coordinate-density - !! reference pressure [Pa]. + !! reference pressure [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 @@ -109,10 +109,11 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state real :: ML_depth ! The specified initial mixed layer depth, in depth units [Z ~> m]. real :: thermocline_scale ! The e-folding scale of the thermocline, in depth units [Z ~> m]. real, dimension(SZK_(GV)) :: & - T0, pres, S0, & ! drho + T0, S0, & ! Profiles of temperature [degC] and salinity [ppt] rho_guess, & ! Potential density at T0 & S0 [R ~> kg m-3]. drho_dT, & ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]. drho_dS ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]. + real :: pres(SZK_(GV)) ! Reference pressure [R L2 T-2 ~> Pa]. real :: a_exp ! The fraction of the overall stratification that is exponential. real :: I_ts, I_md ! Inverse lengthscales [Z-1 ~> m-1]. real :: T_frac ! A ratio of the interface temperature to the range @@ -151,8 +152,10 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state pres(k) = P_Ref ; S0(k) = 35.0 enddo T0(k1) = 29.0 - call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), eqn_of_state, scale=US%kg_m3_to_R) - call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, k1, 1, eqn_of_state, scale=US%kg_m3_to_R) + call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), 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, k1, 1, eqn_of_state, & + scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) ! A first guess of the layers' temperatures. do k=1,nz @@ -161,8 +164,10 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state ! 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 T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) enddo @@ -227,12 +232,12 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, & type(EOS_type), pointer :: eqn_of_state !< integer that selects the !! equation of state. real, intent(in) :: P_Ref !< The coordinate-density - !! reference pressure [Pa]. + !! reference pressure [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)), S0(SZK_(G)) - real :: pres(SZK_(G)) ! Reference pressure [Pa]. + real :: pres(SZK_(G)) ! 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]. @@ -256,8 +261,10 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, & enddo T0(k1) = 29.0 - call calculate_density(T0(k1),S0(k1),pres(k1),rho_guess(k1),eqn_of_state, scale=US%kg_m3_to_R) - call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,k1,1,eqn_of_state, scale=US%kg_m3_to_R) + call calculate_density(T0(k1), S0(k1), pres(k1), rho_guess(k1), 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, k1, 1, eqn_of_state, & + scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) ! A first guess of the layers' temperatures. ! do k=1,nz @@ -266,8 +273,10 @@ subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, & ! 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 T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) enddo