From 1252e3b82594d3241c4b34c4830f3cc951840e1f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 9 Apr 2020 09:01:46 -0400 Subject: [PATCH] Pass rescaled pressures to calculate_TFreeze Use the new pres_scale argument to TFreeze and pass rescaled pressures to calculate_TFreeze and several instances of calculate_density. All answers are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 26 ++++++++++--------- .../vertical/MOM_diabatic_aux.F90 | 17 +++++++----- .../vertical/MOM_entrain_diffusive.F90 | 10 +++---- .../vertical/MOM_full_convection.F90 | 20 +++++++------- .../vertical/MOM_internal_tide_input.F90 | 10 +++---- .../vertical/MOM_set_diffusivity.F90 | 10 +++---- 6 files changed, 49 insertions(+), 44 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 36f97a65fa..298fbc4507 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -158,7 +158,8 @@ module MOM_ice_shelf logical :: find_salt_root !< If true, if true find Sbdry using a quadratic eq. real :: TFr_0_0 !< The freezing point at 0 pressure and 0 salinity [degC] real :: dTFr_dS !< Partial derivative of freezing temperature with salinity [degC ppt-1] - real :: dTFr_dp !< Partial derivative of freezing temperature with pressure [degC Pa-1] + real :: dTFr_dp !< Partial derivative of freezing temperature with + !! pressure [degC T2 R-1 L-2 ~> degC Pa-1] !>@{ Diagnostic handles integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, & id_tfreeze = -1, id_tfl_shelf = -1, & @@ -217,7 +218,7 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) !< with temperature [kg m-3 degC-1]. dR0_dS, & !< Partial derivative of the mixed layer density !< with salinity [kg m-3 ppt-1]. - p_int !< The pressure at the ice-ocean interface [Pa]. + p_int !< The pressure at the ice-ocean interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: & exch_vel_t, & !< Sub-shelf thermal exchange velocity [Z T-1 ~> m s-1] @@ -371,13 +372,13 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) do j=js,je ! Find the pressure at the ice-ocean interface, averaged only over the ! part of the cell covered by ice shelf. - do i=is,ie ; p_int(i) = US%RL2_T2_to_Pa*CS%g_Earth * ISS%mass_shelf(i,j) ; enddo + do i=is,ie ; p_int(i) = CS%g_Earth * ISS%mass_shelf(i,j) ; enddo ! Calculate insitu densities and expansion coefficients call calculate_density(state%sst(:,j), state%sss(:,j), p_int, & - Rhoml(:), is, ie-is+1, CS%eqn_of_state) + Rhoml(:), is, ie-is+1, CS%eqn_of_state, pres_scale=US%RL2_T2_to_Pa) call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, & - dR0_dT, dR0_dS, is, ie-is+1, CS%eqn_of_state) + dR0_dT, dR0_dS, is, ie-is+1, CS%eqn_of_state, pres_scale=US%RL2_T2_to_Pa) do i=is,ie if ((state%ocean_mass(i,j) > US%RZ_to_kg_m2*CS%col_mass_melt_threshold) .and. & @@ -445,7 +446,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) do it1 = 1,20 ! Determine the potential temperature at the ice-ocean interface. - call calculate_TFreeze(Sbdry(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state) + call calculate_TFreeze(Sbdry(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, & + pres_scale=US%RL2_T2_to_Pa) dT_ustar = (ISS%tfreeze(i,j) - state%sst(i,j)) * ustar_h dS_ustar = (Sbdry(i,j) - state%sss(i,j)) * ustar_h @@ -588,7 +590,8 @@ subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces) ! is specified and large enough that the ocean salinity at the interface ! is about the same as the boundary layer salinity. - call calculate_TFreeze(state%sss(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state) + call calculate_TFreeze(state%sss(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, & + pres_scale=US%RL2_T2_to_Pa) exch_vel_t(i,j) = CS%gamma_t ISS%tflux_ocn(i,j) = RhoCp * exch_vel_t(i,j) * (ISS%tfreeze(i,j) - state%sst(i,j)) @@ -1272,12 +1275,11 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl "this is the freezing potential temperature at "//& "S=0, P=0.", units="degC", default=0.0, do_not_log=.true.) call get_param(param_file, mdl, "DTFREEZE_DS", CS%dTFr_dS, & - "this is the derivative of the freezing potential "//& - "temperature with salinity.", units="degC psu-1", default=-0.054, do_not_log=.true.) + "this is the derivative of the freezing potential temperature with salinity.", & + units="degC psu-1", default=-0.054, do_not_log=.true.) call get_param(param_file, mdl, "DTFREEZE_DP", CS%dTFr_dp, & - "this is the derivative of the freezing potential "//& - "temperature with pressure.", & - units="degC Pa-1", default=0.0, do_not_log=.true.) + "this is the derivative of the freezing potential temperature with pressure.", & + units="degC Pa-1", default=0.0, scale=US%RL2_T2_to_Pa, do_not_log=.true.) endif call get_param(param_file, mdl, "G_EARTH", CS%g_Earth, & diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 8fec9a4ca2..1d079d451b 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -118,9 +118,10 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) real, dimension(SZI_(G)) :: & fraz_col, & ! The accumulated heat requirement due to frazil [Q R Z ~> J m-2]. T_freeze, & ! The freezing potential temperature at the current salinity [degC]. - ps ! pressure + ps ! Surface pressure [R L2 T-2 ~> Pa] real, dimension(SZI_(G),SZK_(G)) :: & - pressure ! The pressure at the middle of each layer [Pa]. + pressure ! The pressure at the middle of each layer [R L2 T-2 ~> Pa]. + real :: H_to_RL2_T2 ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1] real :: hc ! A layer's heat capacity [Q R Z degC-1 ~> J m-2 degC-1]. logical :: T_fr_set ! True if the freezing point has been calculated for a ! row of points. @@ -135,6 +136,8 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) if (.not.CS%pressure_dependent_frazil) then do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo + else + H_to_RL2_T2 = GV%H_to_RZ * GV%g_Earth endif !$OMP parallel do default(none) shared(is,ie,js,je,CS,G,GV,US,h,nz,tv,p_surf) & !$OMP private(fraz_col,T_fr_set,T_freeze,hc,ps) & @@ -142,18 +145,18 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) do j=js,je ps(:) = 0.0 if (PRESENT(p_surf)) then ; do i=is,ie - ps(i) = US%RL2_T2_to_Pa*p_surf(i,j) + ps(i) = p_surf(i,j) enddo ; endif do i=is,ie ; fraz_col(i) = 0.0 ; enddo if (CS%pressure_dependent_frazil) then do i=is,ie - pressure(i,1) = ps(i) + (0.5*GV%H_to_Pa)*h(i,j,1) + pressure(i,1) = ps(i) + (0.5*H_to_RL2_T2)*h(i,j,1) enddo do k=2,nz ; do i=is,ie pressure(i,k) = pressure(i,k-1) + & - (0.5*GV%H_to_Pa) * (h(i,j,k) + h(i,j,k-1)) + (0.5*H_to_RL2_T2) * (h(i,j,k) + h(i,j,k-1)) enddo ; enddo endif @@ -162,7 +165,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) do i=is,ie ; if (tv%frazil(i,j) > 0.0) then if (.not.T_fr_set) then call calculate_TFreeze(tv%S(i:,j,1), pressure(i:,1), T_freeze(i:), & - 1, ie-i+1, tv%eqn_of_state) + 1, ie-i+1, tv%eqn_of_state, pres_scale=US%RL2_T2_to_Pa) T_fr_set = .true. endif @@ -188,7 +191,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then if (.not.T_fr_set) then call calculate_TFreeze(tv%S(i:,j,k), pressure(i:,k), T_freeze(i:), & - 1, ie-i+1, tv%eqn_of_state) + 1, ie-i+1, tv%eqn_of_state, pres_scale=US%RL2_T2_to_Pa) T_fr_set = .true. endif diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index 3413f41389..e5366897e4 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -174,7 +174,7 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & real :: g_2dt ! 0.5 * G_Earth / dt, times unit conversion factors ! [m3 H-2 s-2 T-1 ~> m s-3 or m7 kg-2 s-3]. real, dimension(SZI_(G)) :: & - pressure, & ! The pressure at an interface [Pa]. + pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa]. T_eos, S_eos, & ! The potential temperature and salinity at which to ! evaluate dRho_dT and dRho_dS [degC] and [ppt]. dRho_dT, dRho_dS ! The partial derivatives of potential density with temperature and @@ -836,12 +836,12 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ; enddo if (associated(tv%eqn_of_state)) then if (associated(fluxes%p_surf)) then - do i=is,ie ; pressure(i) = US%RL2_T2_to_Pa*fluxes%p_surf(i,j) ; enddo + do i=is,ie ; pressure(i) = fluxes%p_surf(i,j) ; enddo else do i=is,ie ; pressure(i) = 0.0 ; enddo endif do K=2,nz - do i=is,ie ; pressure(i) = pressure(i) + GV%H_to_Pa*h(i,j,k-1) ; enddo + do i=is,ie ; pressure(i) = pressure(i) + (GV%g_Earth*GV%H_to_RZ)*h(i,j,k-1) ; enddo do i=is,ie if (k==kb(i)) then T_eos(i) = 0.5*(tv%T(i,j,kmb) + tv%T(i,j,k)) @@ -851,8 +851,8 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, & S_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k)) endif enddo - call calculate_density_derivs(T_eos, S_eos, pressure, & - dRho_dT, dRho_dS, is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) + call calculate_density_derivs(T_eos, S_eos, pressure, dRho_dT, dRho_dS, is, ie-is+1, & + tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) do i=is,ie if ((k>kmb) .and. (k Pa]. real :: T_EOS(SZI_(G)) ! Filtered and vertically averaged temperatures [degC] real :: S_EOS(SZI_(G)) ! Filtered and vertically averaged salinities [ppt] real :: kap_dt_x2 ! The product of 2*kappa*dt [H2 ~> m2 or kg2 m-4]. @@ -403,24 +403,24 @@ subroutine smoothed_dRdT_dRdS(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, h endif if (associated(p_surf)) then - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa*p_surf(i,j) ; enddo + do i=is,ie ; pres(i) = p_surf(i,j) ; enddo else do i=is,ie ; pres(i) = 0.0 ; enddo endif - call calculate_density_derivs(T_f(:,1), S_f(:,1), pres, dR_dT(:,1), dR_dS(:,1), & - is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) - do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*GV%H_to_Pa ; enddo + call calculate_density_derivs(T_f(:,1), S_f(:,1), pres, dR_dT(:,1), dR_dS(:,1), is-G%isd+1, ie-is+1, & + tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*(GV%H_to_RZ*GV%g_Earth) ; enddo do K=2,nz do i=is,ie T_EOS(i) = 0.5*(T_f(i,k-1) + T_f(i,k)) S_EOS(i) = 0.5*(S_f(i,k-1) + S_f(i,k)) enddo - call calculate_density_derivs(T_EOS, S_EOS, pres, dR_dT(:,K), dR_dS(:,K), & - is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) - do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*GV%H_to_Pa ; enddo + call calculate_density_derivs(T_EOS, S_EOS, pres, dR_dT(:,K), dR_dS(:,K), is-G%isd+1, ie-is+1, & + tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*(GV%H_to_RZ*GV%g_Earth) ; enddo enddo - call calculate_density_derivs(T_f(:,nz), S_f(:,nz), pres, dR_dT(:,nz+1), dR_dS(:,nz+1), & - is-G%isd+1, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) + call calculate_density_derivs(T_f(:,nz), S_f(:,nz), pres, dR_dT(:,nz+1), dR_dS(:,nz+1), is-G%isd+1, & + ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) end subroutine smoothed_dRdT_dRdS diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index a9275c1ccc..6c366e5ff9 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -167,7 +167,7 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) real, dimension(SZI_(G),SZK_(G)+1) :: & dRho_int ! The unfiltered density differences across interfaces [R ~> kg m-3]. real, dimension(SZI_(G)) :: & - pres, & ! The pressure at each interface [Pa]. + pres, & ! The pressure at each interface [R L2 T-2 ~> Pa]. Temp_int, & ! The temperature at each interface [degC]. Salin_int, & ! The salinity at each interface [ppt]. drho_bot, & ! The density difference at the bottom of a layer [R ~> kg m-3] @@ -199,18 +199,18 @@ subroutine find_N2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot) do j=js,je if (associated(tv%eqn_of_state)) then if (associated(fluxes%p_surf)) then - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa*fluxes%p_surf(i,j) ; enddo + do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo else do i=is,ie ; pres(i) = 0.0 ; enddo endif do K=2,nz do i=is,ie - pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1) + pres(i) = pres(i) + (GV%g_Earth*GV%H_to_RZ)*h(i,j,k-1) Temp_Int(i) = 0.5 * (T_f(i,j,k) + T_f(i,j,k-1)) Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) enddo - call calculate_density_derivs(Temp_int, Salin_int, pres, & - dRho_dT(:), dRho_dS(:), is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) + call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT(:), dRho_dS(:), & + is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) do i=is,ie dRho_int(i,K) = max(dRho_dT(i)*(T_f(i,j,k) - T_f(i,j,k-1)) + & dRho_dS(i)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 4f9a6bf478..9c47841748 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -868,7 +868,7 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, & dRho_dS ! partial derivative of density wrt saln [R ppt-1 ~> kg m-3 ppt-1] real, dimension(SZI_(G)) :: & - pres, & ! pressure at each interface [Pa] + pres, & ! pressure at each interface [R L2 T-2 ~> Pa] Temp_int, & ! temperature at each interface [degC] Salin_int, & ! salinity at each interface [ppt] drho_bot, & ! A density difference [R ~> kg m-3] @@ -896,18 +896,18 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, & enddo if (associated(tv%eqn_of_state)) then if (associated(fluxes%p_surf)) then - do i=is,ie ; pres(i) = US%RL2_T2_to_Pa*fluxes%p_surf(i,j) ; enddo + do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ; enddo else do i=is,ie ; pres(i) = 0.0 ; enddo endif do K=2,nz do i=is,ie - pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1) + pres(i) = pres(i) + (GV%g_Earth*GV%H_to_RZ)*h(i,j,k-1) Temp_Int(i) = 0.5 * (T_f(i,j,k) + T_f(i,j,k-1)) Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1)) enddo - call calculate_density_derivs(Temp_int, Salin_int, pres, & - dRho_dT(:,K), dRho_dS(:,K), is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R) + call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT(:,K), dRho_dS(:,K), & + is, ie-is+1, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) do i=is,ie dRho_int(i,K) = max(dRho_dT(i,K)*(T_f(i,j,k) - T_f(i,j,k-1)) + & dRho_dS(i,K)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0)