Skip to content

Commit

Permalink
+Dimensional rescaling of user OBC test cases
Browse files Browse the repository at this point in the history 
  Add dimensional rescaling of user OBC test cases, including documentation of
the units of variables in the Kelvin, shelfwave, tidal_bay and dyed_channel
initialization and rescaling parameters parameters via optional scale arguments
to get_param calls.  These changes also incorporate the answer-changing
correction to the Kelvin wave OBC test case in PR mom-ocean#1406, with a comment noting
what seems like an additional bug in this test case.  This commit includes
adding a unit_scale_type argument to call_OBC_register, register_file_OBC,
register_tidal_bay_OBC, register_Kelvin_OBC, register_shelfwave_OBC and
register_dyed_channel_OBC.  These Kelvin wave OBC configuration from
ESMG-configs now passes the dimensional rescaling tests.  All answers are
bitwise identical in the MOM6_examples test cases, but there are interface
changes.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed May 23, 2021
1 parent 1562e8f commit 29061de
Show file tree
Hide file tree
Showing 7 changed files with 145 additions and 116 deletions.
2 changes: 1 addition & 1 deletion src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2168,7 +2168,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &

call MOM_timing_init(CS)

if (associated(CS%OBC)) call call_OBC_register(param_file, CS%update_OBC_CSp, CS%OBC)
if (associated(CS%OBC)) call call_OBC_register(param_file, CS%update_OBC_CSp, US, CS%OBC)

call tracer_registry_init(param_file, CS%tracer_Reg)

Expand Down
15 changes: 8 additions & 7 deletions src/core/MOM_boundary_update.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,10 @@ module MOM_boundary_update
!> The following subroutines and associated definitions provide the
!! machinery to register and call the subroutines that initialize
!! open boundary conditions.
subroutine call_OBC_register(param_file, CS, OBC)
subroutine call_OBC_register(param_file, CS, US, OBC)
type(param_file_type), intent(in) :: param_file !< Parameter file to parse
type(update_OBC_CS), pointer :: CS !< Control structure for OBCs
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(ocean_OBC_type), pointer :: OBC !< Open boundary structure

! Local variables
Expand Down Expand Up @@ -92,19 +93,19 @@ subroutine call_OBC_register(param_file, CS, OBC)
default=.false.)

if (CS%use_files) CS%use_files = &
register_file_OBC(param_file, CS%file_OBC_CSp, &
register_file_OBC(param_file, CS%file_OBC_CSp, US, &
OBC%OBC_Reg)
if (CS%use_tidal_bay) CS%use_tidal_bay = &
register_tidal_bay_OBC(param_file, CS%tidal_bay_OBC_CSp, &
register_tidal_bay_OBC(param_file, CS%tidal_bay_OBC_CSp, US, &
OBC%OBC_Reg)
if (CS%use_Kelvin) CS%use_Kelvin = &
register_Kelvin_OBC(param_file, CS%Kelvin_OBC_CSp, &
register_Kelvin_OBC(param_file, CS%Kelvin_OBC_CSp, US, &
OBC%OBC_Reg)
if (CS%use_shelfwave) CS%use_shelfwave = &
register_shelfwave_OBC(param_file, CS%shelfwave_OBC_CSp, &
register_shelfwave_OBC(param_file, CS%shelfwave_OBC_CSp, US, &
OBC%OBC_Reg)
if (CS%use_dyed_channel) CS%use_dyed_channel = &
register_dyed_channel_OBC(param_file, CS%dyed_channel_OBC_CSp, &
register_dyed_channel_OBC(param_file, CS%dyed_channel_OBC_CSp, US, &
OBC%OBC_Reg)

end subroutine call_OBC_register
Expand All @@ -128,7 +129,7 @@ subroutine update_OBC_data(OBC, G, GV, US, tv, h, CS, Time)
if (CS%use_Kelvin) &
call Kelvin_set_OBC_data(OBC, CS%Kelvin_OBC_CSp, G, GV, US, h, Time)
if (CS%use_shelfwave) &
call shelfwave_set_OBC_data(OBC, CS%shelfwave_OBC_CSp, G, GV, h, Time)
call shelfwave_set_OBC_data(OBC, CS%shelfwave_OBC_CSp, G, GV, US, h, Time)
if (CS%use_dyed_channel) &
call dyed_channel_update_flow(OBC, CS%dyed_channel_OBC_CSp, G, GV, Time)
if (OBC%needs_IO_for_data .or. OBC%add_tide_constituents) &
Expand Down
6 changes: 4 additions & 2 deletions src/core/MOM_open_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,8 @@ module MOM_open_boundary
!! segment [H L2 T-1 ~> m3 s-1].
real, pointer, dimension(:,:) :: normal_vel_bt=>NULL() !< The barotropic velocity normal to
!! the OB segment [L T-1 ~> m s-1].
real, pointer, dimension(:,:) :: eta=>NULL() !< The sea-surface elevation along the segment [m].
real, pointer, dimension(:,:) :: eta=>NULL() !< The sea-surface elevation along the
!! segment [H ~> m or kg m-2].
real, pointer, dimension(:,:,:) :: grad_normal=>NULL() !< The gradient of the normal flow along the
!! segment times the grid spacing [L T-1 ~> m s-1]
real, pointer, dimension(:,:,:) :: grad_tan=>NULL() !< The gradient of the tangential flow along the
Expand Down Expand Up @@ -4464,9 +4465,10 @@ subroutine OBC_registry_init(param_file, Reg)
end subroutine OBC_registry_init

!> Add file to OBC registry.
function register_file_OBC(param_file, CS, OBC_Reg)
function register_file_OBC(param_file, CS, US, OBC_Reg)
type(param_file_type), intent(in) :: param_file !< parameter file.
type(file_OBC_CS), pointer :: CS !< file control structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(OBC_registry_type), pointer :: OBC_Reg !< OBC registry.
logical :: register_file_OBC
character(len=32) :: casename = "OBC file" !< This case's name.
Expand Down
130 changes: 66 additions & 64 deletions src/user/Kelvin_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,13 +35,13 @@ module Kelvin_initialization
!> Control structure for Kelvin wave open boundaries.
type, public :: Kelvin_OBC_CS ; private
integer :: mode = 0 !< Vertical mode
real :: coast_angle = 0 !< Angle of coastline
real :: coast_offset1 = 0 !< Longshore distance to coastal angle
real :: coast_offset2 = 0 !< Longshore distance to coastal angle
real :: H0 = 0 !< Bottom depth
real :: F_0 !< Coriolis parameter
real :: rho_range !< Density range
real :: rho_0 !< Mean density
real :: coast_angle = 0 !< Angle of coastline [rad]
real :: coast_offset1 = 0 !< Longshore distance to coastal angle [L ~> m]
real :: coast_offset2 = 0 !< Longshore distance to coastal angle [L ~> m]
real :: H0 = 0 !< Bottom depth [Z ~> m]f
real :: F_0 !< Coriolis parameter [T-1 ~> s-1]
real :: rho_range !< Density range [R ~> kg m-3]
real :: rho_0 !< Mean density [R ~> kg m-3]
end type Kelvin_OBC_CS

! This include declares and sets the variable "version".
Expand All @@ -50,9 +50,10 @@ module Kelvin_initialization
contains

!> Add Kelvin wave to OBC registry.
function register_Kelvin_OBC(param_file, CS, OBC_Reg)
function register_Kelvin_OBC(param_file, CS, US, OBC_Reg)
type(param_file_type), intent(in) :: param_file !< parameter file.
type(Kelvin_OBC_CS), pointer :: CS !< Kelvin wave control structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(OBC_registry_type), pointer :: OBC_Reg !< OBC registry.

! Local variables
Expand All @@ -73,31 +74,29 @@ function register_Kelvin_OBC(param_file, CS, OBC_Reg)
"Vertical Kelvin wave mode imposed at upstream open boundary.", &
default=0)
call get_param(param_file, mdl, "F_0", CS%F_0, &
default=0.0, do_not_log=.true.)
default=0.0, units="s-1", scale=US%T_to_s, do_not_log=.true.)
call get_param(param_file, mdl, "TOPO_CONFIG", config, do_not_log=.true.)
if (trim(config) == "Kelvin") then
call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", CS%coast_offset1, &
"The distance along the southern and northern boundaries "//&
"at which the coasts angle in.", &
units="km", default=100.0)
units="km", default=100.0, scale=1.0e3*US%m_to_L)
call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", CS%coast_offset2, &
"The distance from the southern and northern boundaries "//&
"at which the coasts angle in.", &
units="km", default=10.0)
units="km", default=10.0, scale=1.0e3*US%m_to_L)
call get_param(param_file, mdl, "ROTATED_COAST_ANGLE", CS%coast_angle, &
"The angle of the southern bondary beyond X=ROTATED_COAST_OFFSET.", &
units="degrees", default=11.3)
CS%coast_angle = CS%coast_angle * (atan(1.0)/45.) ! Convert to radians
CS%coast_offset1 = CS%coast_offset1 * 1.e3 ! Convert to m
CS%coast_offset2 = CS%coast_offset2 * 1.e3 ! Convert to m
endif
if (CS%mode /= 0) then
call get_param(param_file, mdl, "DENSITY_RANGE", CS%rho_range, &
default=2.0, do_not_log=.true.)
default=2.0, do_not_log=.true., scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "RHO_0", CS%rho_0, &
default=1035.0, do_not_log=.true.)
default=1035.0, do_not_log=.true., scale=US%kg_m3_to_R)
call get_param(param_file, mdl, "MAXIMUM_DEPTH", CS%H0, &
default=1000.0, do_not_log=.true.)
default=1000.0, do_not_log=.true., scale=US%m_to_Z)
endif

! Register the Kelvin open boundary.
Expand All @@ -122,7 +121,7 @@ subroutine Kelvin_initialize_topography(D, G, param_file, max_depth, US)
real, dimension(G%isd:G%ied,G%jsd:G%jed), &
intent(out) :: D !< Ocean bottom depth in m or Z if US is present
type(param_file_type), intent(in) :: param_file !< Parameter file structure
real, intent(in) :: max_depth !< Maximum model depth in the units of D
real, intent(in) :: max_depth !< Maximum model depth in the units of D [Z ~> m or m]
type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type

! Local variables
Expand Down Expand Up @@ -176,22 +175,27 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
type(Kelvin_OBC_CS), pointer :: CS !< Kelvin wave control structure.
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
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2].
type(time_type), intent(in) :: Time !< model time.

! The following variables are used to set up the transport in the Kelvin example.
real :: time_sec, cff
real :: N0 ! Brunt-Vaisala frequency [s-1]
real :: plx !< Longshore wave parameter
real :: pmz !< Vertical wave parameter
real :: lambda !< Offshore decay scale
real :: omega !< Wave frequency [s-1]
real :: time_sec ! The time in the run [T ~> s]
real :: cff ! The wave speed [L T-1 ~> m s-1]
real :: N0 ! Brunt-Vaisala frequency times a rescaling of slopes [L Z-1 T-1 ~> s-1]
real :: lambda ! Offshore decay scale [L-1 ~> m-1]
real :: omega ! Wave frequency [T-1 ~> s-1]
real :: PI
integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz
integer :: IsdB, IedB, JsdB, JedB
real :: fac, x, y, x1, y1
real :: val1, val2, sina, cosa
real :: mag_SSH ! An overall magnitude of the external wave sea surface height at the coastline [Z ~> m]
real :: mag_int ! An overall magnitude of the internal wave at the coastline [L2 T-2 ~> m2 s-2]
real :: x1, y1 ! Various positions [L ~> m]
real :: x, y ! Various positions [L ~> m]
real :: val1 ! The periodicity factor [nondim]
real :: val2 ! The local wave amplitude [Z ~> m]
real :: km_to_L_scale ! A scaling factor from longitudes in km to L [L km-1 ~> 1e3]
real :: sina, cosa ! The sine and cosine of the coast angle [nondim]
type(OBC_segment_type), pointer :: segment => NULL()

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Expand All @@ -201,23 +205,20 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
if (.not.associated(OBC)) call MOM_error(FATAL, 'Kelvin_initialization.F90: '// &
'Kelvin_set_OBC_data() was called but OBC type was not initialized!')

time_sec = time_type_to_real(Time)
time_sec = US%s_to_T*time_type_to_real(Time)
PI = 4.0*atan(1.0)
fac = 1.0
km_to_L_scale = 1000.0*US%m_to_L

if (CS%mode == 0) then
omega = 2.0 * PI / (12.42 * 3600.0) ! M2 Tide period
val1 = US%m_to_Z * sin(omega * time_sec)
mag_SSH = 1.0*US%m_to_Z
omega = 2.0 * PI / (12.42 * 3600.0*US%s_to_T) ! M2 Tide period
val1 = sin(omega * time_sec)
else
N0 = US%L_to_m*US%s_to_T * sqrt((CS%rho_range / CS%rho_0) * GV%g_Earth * (US%m_to_Z * CS%H0))
mag_int = 1.0*US%m_s_to_L_T**2
N0 = sqrt((CS%rho_range / CS%rho_0) * GV%g_Earth / CS%H0)
lambda = PI * CS%mode * CS%F_0 / (CS%H0 * N0)
! Two wavelengths in domain
plx = 4.0 * PI / G%len_lon
pmz = PI * CS%mode / CS%H0
lambda = pmz * CS%F_0 / N0
omega = CS%F_0 * plx / lambda

! lambda = PI * CS%mode * CS%F_0 / (CS%H0 * N0)
! omega = (4.0 * CS%H0 * N0) / (CS%mode * G%len_lon)
omega = (4.0 * CS%H0 * N0) / (CS%mode * US%m_to_L*G%len_lon)
endif

sina = sin(CS%coast_angle)
Expand All @@ -230,22 +231,23 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
if (segment%direction == OBC_DIRECTION_N) cycle

! This should be somewhere else...
segment%Velocity_nudging_timescale_in = 1.0/(0.3*86400)
!### This is supposed to be a timescale [T ~> s] but appears to be a rate in [s-1].
segment%Velocity_nudging_timescale_in = US%s_to_T * 1.0/(0.3*86400)

if (segment%direction == OBC_DIRECTION_W) then
IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB
jsd = segment%HI%jsd ; jed = segment%HI%jed
JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB
do j=jsd,jed ; do I=IsdB,IedB
x1 = 1000. * G%geoLonCu(I,j)
y1 = 1000. * G%geoLatCu(I,j)
x1 = km_to_L_scale * G%geoLonCu(I,j)
y1 = km_to_L_scale * G%geoLatCu(I,j)
x = (x1 - CS%coast_offset1) * cosa + y1 * sina
y = - (x1 - CS%coast_offset1) * sina + y1 * cosa
y = -(x1 - CS%coast_offset1) * sina + y1 * cosa
if (CS%mode == 0) then
! Use inside bathymetry
cff = sqrt(GV%g_Earth * G%bathyT(i+1,j) )
val2 = fac * exp(- US%T_to_s*CS%F_0 * US%m_to_L*y / cff)
segment%eta(I,j) = val2 * cos(omega * time_sec)
val2 = mag_SSH * exp(- CS%F_0 * y / cff)
segment%eta(I,j) = GV%Z_to_H*val2 * cos(omega * time_sec)
segment%normal_vel_bt(I,j) = (val2 * (val1 * cff * cosa / &
(G%bathyT(i+1,j) )) )
if (segment%nudged) then
Expand All @@ -266,14 +268,14 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
segment%normal_vel_bt(I,j) = 0.0
if (segment%nudged) then
do k=1,nz
segment%nudged_normal_vel(I,j,k) = US%m_s_to_L_T * fac * lambda / CS%F_0 * &
exp(- lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * &
segment%nudged_normal_vel(I,j,k) = mag_int * lambda / CS%F_0 * &
exp(-lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * &
cos(omega * time_sec)
enddo
elseif (segment%specified) then
do k=1,nz
segment%normal_vel(I,j,k) = US%m_s_to_L_T * fac * lambda / CS%F_0 * &
exp(- lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * &
segment%normal_vel(I,j,k) = mag_int * lambda / CS%F_0 * &
exp(-lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * &
cos(omega * time_sec)
segment%normal_trans(I,j,k) = segment%normal_vel(I,j,k) * h(i+1,j,k) * G%dyCu(I,j)
enddo
Expand All @@ -282,12 +284,12 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
enddo ; enddo
if (associated(segment%tangential_vel)) then
do J=JsdB+1,JedB-1 ; do I=IsdB,IedB
x1 = 1000. * G%geoLonBu(I,J)
y1 = 1000. * G%geoLatBu(I,J)
x1 = km_to_L_scale * G%geoLonBu(I,J)
y1 = km_to_L_scale * G%geoLatBu(I,J)
x = (x1 - CS%coast_offset1) * cosa + y1 * sina
y = - (x1 - CS%coast_offset1) * sina + y1 * cosa
cff =sqrt(GV%g_Earth * G%bathyT(i+1,j) )
val2 = fac * exp(- US%T_to_s*CS%F_0 * US%m_to_L*y / cff)
cff = sqrt(GV%g_Earth * G%bathyT(i+1,j) )
val2 = mag_SSH * exp(- CS%F_0 * y / cff)
if (CS%mode == 0) then ; do k=1,nz
segment%tangential_vel(I,J,k) = (val1 * val2 * cff * sina) / &
( 0.5*(G%bathyT(i+1,j+1) + G%bathyT(i+1,j) ) )
Expand All @@ -299,24 +301,24 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
isd = segment%HI%isd ; ied = segment%HI%ied
JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB
do J=JsdB,JedB ; do i=isd,ied
x1 = 1000. * G%geoLonCv(i,J)
y1 = 1000. * G%geoLatCv(i,J)
x1 = km_to_L_scale * G%geoLonCv(i,J)
y1 = km_to_L_scale * G%geoLatCv(i,J)
x = (x1 - CS%coast_offset1) * cosa + y1 * sina
y = - (x1 - CS%coast_offset1) * sina + y1 * cosa
if (CS%mode == 0) then
cff = sqrt(GV%g_Earth * G%bathyT(i,j+1) )
val2 = fac * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * US%m_to_L*y / cff)
segment%eta(I,j) = val2 * cos(omega * time_sec)
segment%normal_vel_bt(I,j) = US%L_T_to_m_s * (val1 * cff * sina / &
val2 = mag_SSH * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * y / cff)
segment%eta(I,j) = GV%Z_to_H*val2 * cos(omega * time_sec)
segment%normal_vel_bt(I,j) = (val1 * cff * sina / &
(G%bathyT(i,j+1) )) * val2
if (segment%nudged) then
do k=1,nz
segment%nudged_normal_vel(I,j,k) = US%L_T_to_m_s * (val1 * cff * sina / &
segment%nudged_normal_vel(I,j,k) = (val1 * cff * sina / &
(G%bathyT(i,j+1) )) * val2
enddo
elseif (segment%specified) then
do k=1,nz
segment%normal_vel(I,j,k) = US%L_T_to_m_s * (val1 * cff * sina / &
segment%normal_vel(I,j,k) = (val1 * cff * sina / &
(G%bathyT(i,j+1) )) * val2
segment%normal_trans(i,J,k) = segment%normal_vel(i,J,k) * h(i,j+1,k) * G%dxCv(i,J)
enddo
Expand All @@ -327,12 +329,12 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
segment%normal_vel_bt(i,J) = 0.0
if (segment%nudged) then
do k=1,nz
segment%nudged_normal_vel(i,J,k) = US%m_s_to_L_T*fac * lambda / CS%F_0 * &
segment%nudged_normal_vel(i,J,k) = mag_int * lambda / CS%F_0 * &
exp(- lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * cosa
enddo
elseif (segment%specified) then
do k=1,nz
segment%normal_vel(i,J,k) = US%m_s_to_L_T*fac * lambda / CS%F_0 * &
segment%normal_vel(i,J,k) = mag_int * lambda / CS%F_0 * &
exp(- lambda * y) * cos(PI * CS%mode * (k - 0.5) / nz) * cosa
segment%normal_trans(i,J,k) = segment%normal_vel(i,J,k) * h(i,j+1,k) * G%dxCv(i,J)
enddo
Expand All @@ -341,12 +343,12 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
enddo ; enddo
if (associated(segment%tangential_vel)) then
do J=JsdB,JedB ; do I=IsdB+1,IedB-1
x1 = 1000. * G%geoLonBu(I,J)
y1 = 1000. * G%geoLatBu(I,J)
x1 = km_to_L_scale * G%geoLonBu(I,J)
y1 = km_to_L_scale * G%geoLatBu(I,J)
x = (x1 - CS%coast_offset1) * cosa + y1 * sina
y = - (x1 - CS%coast_offset1) * sina + y1 * cosa
cff = sqrt(GV%g_Earth * G%bathyT(i,j+1) )
val2 = fac * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * US%m_to_L*y / cff)
val2 = mag_SSH * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * y / cff)
if (CS%mode == 0) then ; do k=1,nz
segment%tangential_vel(I,J,k) = ((val1 * val2 * cff * sina) / &
( 0.5*((G%bathyT(i+1,j+1)) + G%bathyT(i,j+1))) )
Expand Down
Loading

0 comments on commit 29061de

Please sign in to comment.