Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

*Revise BFB_set_coord and BFB_buoyancy_forcing #454

Merged
merged 2 commits into from
Aug 25, 2023
Merged

*Revise BFB_set_coord and BFB_buoyancy_forcing #454

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
fa1088c
*Revise BFB_set_coord and BFB_buoyancy_forcing
Hallberg-NOAA Jul 7, 2023
5a5d6fa
Merge branch 'dev/gfdl' into nonBous_BFB_test_case
marshallward Aug 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 22 additions & 9 deletions src/user/BFB_initialization.F90
View file
Open in desktop
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
Open in desktop
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