Skip to content

Commit

Permalink
*Revise BFB_set_coord and BFB_buoyancy_forcing
Browse files Browse the repository at this point in the history 
  Revised BFB_set_coord and BFB_buoyancy_forcing to be consistent with other
instances of a linear equation of state, and to set g_prime in non-Boussinesq
mode similarly to how it is set in other places.   Also use RESTORE_FLUX_RHO in
place of RHO_0 when to specify the densities that are used to convert the piston
velocities into restoring heat or salt fluxes with BFB_buoyancy_forcing, like in
other places in the MOM6 code.  This change will change answers by default when
BFB_set_coord is used, but the old Boussinesq answers can be recovered by
setting the reference salinity S_REF to 38.75 ppt or by setting RHO_T0_S0 to
1003.0 kg m-3.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Jul 7, 2023
1 parent 2a4b1f6 commit 835ca9d
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 18 deletions.
31 changes: 22 additions & 9 deletions src/user/BFB_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,32 +38,45 @@ subroutine BFB_set_coord(Rlay, g_prime, GV, US, param_file)
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters

real :: Rho_T0_S0 ! The density at T=0, S=0 [R ~> kg m-3]
real :: dRho_dT ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
real :: dRho_dS ! The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
real :: SST_s, T_bot ! Temperatures at the surface and seafloor [C ~> degC]
real :: S_ref ! Reference salinity [S ~> ppt]
real :: rho_top, rho_bot ! Densities at the surface and seafloor [R ~> kg m-3]
integer :: k, nz
! This include declares and sets the variable "version".
# include "version_variable.h"
character(len=40) :: mdl = "BFB_initialization" ! This module's name.

call log_version(param_file, mdl, version, "")
call get_param(param_file, mdl, "DRHO_DT", drho_dt, &
"Rate of change of density with temperature.", &
units="kg m-3 K-1", default=-0.2, scale=US%kg_m3_to_R*US%C_to_degC)
call get_param(param_file, mdl, "DRHO_DT", dRho_dT, &
"The partial derivative of density with temperature.", &
units="kg m-3 K-1", default=-0.2, scale=US%kg_m3_to_R*US%C_to_degC)
call get_param(param_file, mdl, "DRHO_DS", dRho_dS, &
"The partial derivative of density with salinity.", &
units="kg m-3 PSU-1", default=0.8, scale=US%kg_m3_to_R*US%S_to_ppt)
call get_param(param_file, mdl, "RHO_T0_S0", Rho_T0_S0, &
"The density at T=0, S=0.", units="kg m-3", default=1000.0, scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "SST_S", SST_s, &
"SST at the southern edge of the domain.", units="degC", default=20.0, scale=US%degC_to_C)
"SST at the southern edge of the domain.", &
units="degC", default=20.0, scale=US%degC_to_C)
call get_param(param_file, mdl, "T_BOT", T_bot, &
"Bottom temperature", units="degC", default=5.0, scale=US%degC_to_C)
rho_top = GV%Rho0 + drho_dt*SST_s
rho_bot = GV%Rho0 + drho_dt*T_bot
call get_param(param_file, mdl, "S_REF", S_ref, &
"The initial salinities.", units="PSU", default=35.0, scale=US%ppt_to_S)
rho_top = (Rho_T0_S0 + dRho_dS*S_ref) + dRho_dT*SST_s
rho_bot = (Rho_T0_S0 + dRho_dS*S_ref) + dRho_dT*T_bot
nz = GV%ke

do k = 1,nz
Rlay(k) = (rho_bot - rho_top)/(nz-1)*real(k-1) + rho_top
if (k >1) then
g_prime(k) = (Rlay(k) - Rlay(k-1)) * GV%g_Earth / (GV%Rho0)
else
if (k==1) then
g_prime(k) = GV%g_Earth
elseif (GV%Boussinesq) then
g_prime(k) = (Rlay(k) - Rlay(k-1)) * GV%g_Earth / GV%Rho0
else
g_prime(k) = (Rlay(k) - Rlay(k-1)) * GV%g_Earth / (0.5*(Rlay(k) + Rlay(k-1)))
endif
enddo

Expand Down
35 changes: 26 additions & 9 deletions src/user/BFB_surface_forcing.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,12 +29,17 @@ module BFB_surface_forcing
real :: Rho0 !< The density used in the Boussinesq approximation [R ~> kg m-3].
real :: G_Earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
real :: Flux_const !< The restoring rate at the surface [Z T-1 ~> m s-1].
real :: rho_restore !< The density that is used to convert piston velocities into salt
!! or heat fluxes with salinity or temperature restoring [R ~> kg m-3]
real :: SST_s !< SST at the southern edge of the linear forcing ramp [C ~> degC]
real :: SST_n !< SST at the northern edge of the linear forcing ramp [C ~> degC]
real :: S_ref !< Reference salinity used throughout the domain [S ~> ppt]
real :: lfrslat !< Southern latitude where the linear forcing ramp begins [degrees_N] or [km]
real :: lfrnlat !< Northern latitude where the linear forcing ramp ends [degrees_N] or [km]
real :: drho_dt !< Rate of change of density with temperature [R C-1 ~> kg m-3 degC-1].
!! Note that temperature is being used as a dummy variable here.
real :: Rho_T0_S0 !< The density at T=0, S=0 [R ~> kg m-3]
real :: dRho_dT !< The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
real :: dRho_dS !< The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1]
!! Note that temperature and salinity are being used as dummy variables here.
!! All temperatures are converted into density.

type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
Expand Down Expand Up @@ -125,7 +130,7 @@ subroutine BFB_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
call MOM_error(FATAL, "User_buoyancy_surface_forcing: " // &
"Temperature and salinity restoring used without modification." )

rhoXcp = CS%Rho0 * fluxes%C_p
rhoXcp = CS%rho_restore * fluxes%C_p
do j=js,je ; do i=is,ie
! Set Temp_restore and Salin_restore to the temperature (in [C ~> degC]) and
! salinity (in [S ~> ppt]) that are being restored toward.
Expand All @@ -134,7 +139,7 @@ subroutine BFB_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)

fluxes%heat_added(i,j) = (G%mask2dT(i,j) * (rhoXcp * CS%Flux_const)) * &
(Temp_restore - sfc_state%SST(i,j))
fluxes%vprec(i,j) = - (G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)) * &
fluxes%vprec(i,j) = - (G%mask2dT(i,j) * (CS%rho_restore*CS%Flux_const)) * &
((Salin_restore - sfc_state%SSS(i,j)) / (0.5 * (Salin_restore + sfc_state%SSS(i,j))))
enddo ; enddo
else
Expand All @@ -144,7 +149,7 @@ subroutine BFB_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
! "Buoyancy restoring used without modification." )

! The -1 is because density has the opposite sign to buoyancy.
buoy_rest_const = -1.0 * (CS%G_Earth * CS%Flux_const) / CS%Rho0
buoy_rest_const = -1.0 * (CS%G_Earth * CS%Flux_const) / CS%rho_restore
Temp_restore = 0.0
do j=js,je ; do i=is,ie
! Set density_restore to an expression for the surface potential
Expand All @@ -158,8 +163,7 @@ subroutine BFB_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
(G%geoLatT(i,j) - CS%lfrslat) + CS%SST_s
endif

density_restore = Temp_restore*CS%drho_dt + CS%Rho0

density_restore = (CS%Rho_T0_S0 + CS%dRho_dS*CS%S_ref) + CS%dRho_dT*Temp_restore
fluxes%buoy(i,j) = G%mask2dT(i,j) * buoy_rest_const * &
(density_restore - sfc_state%sfc_density(i,j))
enddo ; enddo
Expand Down Expand Up @@ -216,9 +220,17 @@ subroutine BFB_surface_forcing_init(Time, G, US, param_file, diag, CS)
call get_param(param_file, mdl, "SST_N", CS%SST_n, &
"SST at the northern edge of the linear forcing ramp.", &
units="degC", default=10.0, scale=US%degC_to_C)
call get_param(param_file, mdl, "DRHO_DT", CS%drho_dt, &
"The rate of change of density with temperature.", &
call get_param(param_file, mdl, "DRHO_DT", CS%dRho_dT, &
"The partial derivative of density with temperature.", &
units="kg m-3 K-1", default=-0.2, scale=US%kg_m3_to_R*US%C_to_degC)
call get_param(param_file, mdl, "DRHO_DS", CS%dRho_dS, &
"The partial derivative of density with salinity.", &
units="kg m-3 PSU-1", default=0.8, scale=US%kg_m3_to_R*US%S_to_ppt)
call get_param(param_file, mdl, "RHO_T0_S0", CS%Rho_T0_S0, &
"The density at T=0, S=0.", units="kg m-3", default=1000.0, scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "S_REF", CS%S_ref, &
"The reference salinity used here throughout the domain.", &
units="PSU", default=35.0, scale=US%ppt_to_S)

call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, &
"If true, the buoyancy fluxes drive the model back "//&
Expand All @@ -231,6 +243,11 @@ subroutine BFB_surface_forcing_init(Time, G, US, param_file, diag, CS)
default=0.0, units="m day-1", scale=US%m_to_Z*US%T_to_s)
! Convert CS%Flux_const from m day-1 to m s-1.
CS%Flux_const = CS%Flux_const / 86400.0
call get_param(param_file, mdl, "RESTORE_FLUX_RHO", CS%rho_restore, &
"The density that is used to convert piston velocities into salt or heat "//&
"fluxes with RESTORE_SALINITY or RESTORE_TEMPERATURE.", &
units="kg m-3", default=CS%Rho0*US%R_to_kg_m3, scale=US%kg_m3_to_R, &
do_not_log=(CS%Flux_const==0.0))
endif

end subroutine BFB_surface_forcing_init
Expand Down

0 comments on commit 835ca9d

Please sign in to comment.