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

Jump to bottom

+Add 6 runtime parameters for OBC testcase initialization modules #290

Merged
merged 4 commits into from
Dec 29, 2022
Merged

+Add 6 runtime parameters for OBC testcase initialization modules #290

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9b3f1ba
+Add runtime parameter for dumbbell_initialization
Hallberg-NOAA Dec 28, 2022
780ff97
+Add runtime parameters for Kelvin_initialization
Hallberg-NOAA Dec 28, 2022
7c12ead
+Add runtime parameters for tidal_bay_initialization
Hallberg-NOAA Dec 28, 2022
b4d72ae
Merge branch 'dev/gfdl' into explicit_OBC_user_params
marshallward Dec 29, 2022
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
64 changes: 39 additions & 25 deletions src/user/Kelvin_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ module Kelvin_initialization
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_E
use MOM_open_boundary, only : OBC_DIRECTION_S, OBC_DIRECTION_W
use MOM_open_boundary, only : OBC_registry_type
use MOM_unit_scaling, only : unit_scale_type
use MOM_unit_scaling, only : unit_scale_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_time_manager, only : time_type, time_type_to_real

Expand All @@ -35,13 +35,16 @@ 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 [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]
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]
real :: coast_angle = 0 !< Angle of coastline [rad]
real :: coast_offset1 = 0 !< Longshore distance to coastal angle [L ~> m]
real :: coast_offset2 = 0 !< Offshore distance to coastal angle [L ~> m]
real :: H0 = 0 !< Bottom depth [Z ~> m]
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]
real :: wave_period !< Period of the mode-0 waves [T ~> s]
real :: ssh_amp !< Amplitude of the sea surface height forcing for mode-0 waves [Z ~> m]
real :: inflow_amp !< Amplitude of the boundary velocity forcing for internal waves [L T-1 ~> m s-1]
end type Kelvin_OBC_CS

! This include declares and sets the variable "version".
Expand Down Expand Up @@ -87,16 +90,28 @@ function register_Kelvin_OBC(param_file, CS, US, OBC_Reg)
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
units="degrees", default=11.3, scale=atan(1.0)/45.) ! Convert to radians
else
CS%coast_offset1 = 0.0 ; CS%coast_offset2 = 0.0 ; CS%coast_angle = 0.0
endif
if (CS%mode /= 0) then
if (CS%mode == 0) then
call get_param(param_file, mdl, "KELVIN_WAVE_PERIOD", CS%wave_period, &
"The period of the Kelvin wave forcing at the open boundaries. "//&
"The default value is the M2 tide period.", &
units="s", default=12.42*3600.0, scale=US%s_to_T)
call get_param(param_file, mdl, "KELVIN_WAVE_SSH_AMP", CS%ssh_amp, &
"The amplitude of the Kelvin wave sea surface height anomaly forcing "//&
"at the open boundaries.", units="m", default=1.0, scale=US%m_to_Z)
else
call get_param(param_file, mdl, "DENSITY_RANGE", CS%rho_range, &
units="kg m-3", default=2.0, scale=US%kg_m3_to_R, do_not_log=.true.)
call get_param(param_file, mdl, "RHO_0", CS%rho_0, &
units="kg m-3", default=1035.0, scale=US%kg_m3_to_R, do_not_log=.true.)
call get_param(param_file, mdl, "MAXIMUM_DEPTH", CS%H0, &
units="m", default=1000.0, scale=US%m_to_Z, do_not_log=.true.)
call get_param(param_file, mdl, "KELVIN_WAVE_INFLOW_AMP", CS%inflow_amp, &
"The amplitude of the Kelvin wave sea surface inflow velocity forcing "//&
"at the open boundaries.", units="m s-1", default=1.0, scale=US%m_s_to_L_T)
endif

! Register the Kelvin open boundary.
Expand Down Expand Up @@ -126,8 +141,10 @@ subroutine Kelvin_initialize_topography(D, G, param_file, max_depth, US)

! Local variables
character(len=40) :: mdl = "Kelvin_initialize_topography" ! This subroutine's name.
real :: min_depth ! The minimum and maximum depths [Z ~> m].
real :: coast_offset1, coast_offset2, coast_angle, right_angle
real :: min_depth ! The minimum and maximum depths [Z ~> m].
real :: coast_angle ! Angle of coastline [rad]
real :: coast_offset1 ! Longshore distance to coastal angle [L ~> m]
real :: coast_offset2 ! Offshore distance to coastal angle [L ~> m]
integer :: i, j

call MOM_mesg(" Kelvin_initialization.F90, Kelvin_initialize_topography: setting topography", 5)
Expand All @@ -139,22 +156,19 @@ subroutine Kelvin_initialize_topography(D, G, param_file, max_depth, US)
call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", coast_offset2, &
units="km", default=10.0, do_not_log=.true.)
call get_param(param_file, mdl, "ROTATED_COAST_ANGLE", coast_angle, &
units="degrees", default=11.3, do_not_log=.true.)

coast_angle = coast_angle * (atan(1.0)/45.) ! Convert to radians
right_angle = 2 * atan(1.0)
units="degrees", default=11.3, scale=(atan(1.0)/45.), do_not_log=.true.) ! Convert to radians

do j=G%jsc,G%jec ; do i=G%isc,G%iec
D(i,j) = max_depth
! Southern side
if ((G%geoLonT(i,j) - G%west_lon > coast_offset1) .AND. &
(atan2(G%geoLatT(i,j) - G%south_lat + coast_offset2, &
G%geoLonT(i,j) - G%west_lon - coast_offset1) < coast_angle)) &
G%geoLonT(i,j) - G%west_lon - coast_offset1) < coast_angle)) &
D(i,j) = 0.5*min_depth
! Northern side
if ((G%geoLonT(i,j) - G%west_lon < G%len_lon - coast_offset1) .AND. &
(atan2(G%len_lat + G%south_lat + coast_offset2 - G%geoLatT(i,j), &
G%len_lon + G%west_lon - coast_offset1 - G%geoLonT(i,j)) < coast_angle)) &
G%len_lon + G%west_lon - coast_offset1 - G%geoLonT(i,j)) < coast_angle)) &
D(i,j) = 0.5*min_depth

if (D(i,j) > max_depth) D(i,j) = max_depth
Expand All @@ -181,10 +195,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
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
real :: PI ! The ratio of the circumference of a circle to its diameter [nondim]
real :: depth_tot(SZI_(G),SZJ_(G)) ! The total depth of the ocean [Z ~> m]
integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz
integer :: IsdB, IedB, JsdB, JedB
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]
Expand All @@ -194,6 +206,8 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
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()
integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz
integer :: IsdB, IedB, JsdB, JedB

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
Expand All @@ -214,11 +228,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time)
enddo ; enddo ; enddo

if (CS%mode == 0) then
mag_SSH = 1.0*US%m_to_Z
omega = 2.0 * PI / (12.42 * 3600.0*US%s_to_T) ! M2 Tide period
mag_SSH = CS%ssh_amp
omega = 2.0 * PI / CS%wave_period
val1 = sin(omega * time_sec)
else
mag_int = 1.0*US%m_s_to_L_T**2
mag_int = CS%inflow_amp**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
Expand Down
21 changes: 10 additions & 11 deletions src/user/dumbbell_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,11 +52,10 @@ subroutine dumbbell_initialize_topography( D, G, param_file, max_depth )
logical :: dbrotate ! If true, rotate this configuration
integer :: i, j

call get_param(param_file, mdl, "DUMBBELL_LEN",dblen, &
call get_param(param_file, mdl, "DUMBBELL_LEN", dblen, &
'Lateral Length scale for dumbbell.', &
units='km', default=600., do_not_log=.false.)
! units=G%x_ax_unit_short, default=600., do_not_log=.false.)
call get_param(param_file, mdl, "DUMBBELL_FRACTION",dbfrac, &
units=G%x_ax_unit_short, default=600., do_not_log=.false.)
call get_param(param_file, mdl, "DUMBBELL_FRACTION", dbfrac, &
'Meridional fraction for narrow part of dumbbell.', &
units='nondim', default=0.5, do_not_log=.false.)
call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
Expand Down Expand Up @@ -275,8 +274,6 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, US, param_

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

T_surf = 20.0*US%degC_to_C

! layer mode
call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.)
if (.not. use_ALE) call MOM_error(FATAL, "dumbbell_initialize_temperature_salinity: "//&
Expand All @@ -287,16 +284,18 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, US, param_
call get_param(param_file, mdl, "INITIAL_DENSITY_PROFILE", density_profile, &
'Initial profile shape. Valid values are "linear", "parabolic" '// &
'and "exponential".', default='linear', do_not_log=just_read)
call get_param(param_file, mdl, "DUMBBELL_T_SURF", T_surf, &
'Initial surface temperature in the DUMBBELL configuration', &
units='degC', default=20., scale=US%degC_to_C, do_not_log=just_read)
call get_param(param_file, mdl, "DUMBBELL_SREF", S_surf, &
'DUMBBELL REFERENCE SALINITY', &
units='1e-3', default=34., scale=US%ppt_to_S, do_not_log=just_read)
call get_param(param_file, mdl, "DUMBBELL_S_RANGE", S_range, &
'DUMBBELL salinity range (right-left)', &
units='1e-3', default=2., scale=US%ppt_to_S, do_not_log=just_read)
call get_param(param_file, mdl, "DUMBBELL_LEN", dblen, &
'Lateral Length scale for dumbbell ', &
units='km', default=600., do_not_log=just_read)
! units=G%x_ax_unit_short, default=600., do_not_log=.false.)
'Lateral Length scale for dumbbell ', &
units=G%x_ax_unit_short, default=600., do_not_log=just_read)
call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
'Logical for rotation of dumbbell domain.', &
default=.false., do_not_log=just_read)
Expand Down Expand Up @@ -376,8 +375,8 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, h_in, depth_tot, param_fil
nz = GV%ke

call get_param(param_file, mdl, "DUMBBELL_SPONGE_TIME_SCALE", sponge_time_scale, &
"The time scale in the reservoir for restoring. If zero, the sponge is disabled.", &
units="s", default=0., scale=US%s_to_T)
"The time scale in the reservoir for restoring. If zero, the sponge is disabled.", &
units="s", default=0., scale=US%s_to_T)
call get_param(param_file, mdl, "DUMBBELL_SREF", S_ref, &
'DUMBBELL REFERENCE SALINITY', &
units='1e-3', default=34., scale=US%ppt_to_S, do_not_log=.true.)
Expand Down
26 changes: 18 additions & 8 deletions src/user/tidal_bay_initialization.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,10 @@ module tidal_bay_initialization

!> Control structure for tidal bay open boundaries.
type, public :: tidal_bay_OBC_CS ; private
real :: tide_flow = 3.0e6 !< Maximum tidal flux [L2 Z T-1 ~> m3 s-1]
real :: tide_flow = 3.0e6 !< Maximum tidal flux with the tidal bay configuration [L2 Z T-1 ~> m3 s-1]
real :: tide_period !< The period associated with the tidal bay configuration [T ~> s-1]
real :: tide_ssh_amp !< The magnitude of the sea surface height anomalies at the inflow
!! with the tidal bay configuration [Z ~> m]
end type tidal_bay_OBC_CS

contains
Expand All @@ -43,6 +46,13 @@ function register_tidal_bay_OBC(param_file, CS, US, OBC_Reg)
call get_param(param_file, mdl, "TIDAL_BAY_FLOW", CS%tide_flow, &
"Maximum total tidal volume flux.", &
units="m3 s-1", default=3.0e6, scale=US%m_s_to_L_T*US%m_to_L*US%m_to_Z)
call get_param(param_file, mdl, "TIDAL_BAY_PERIOD", CS%tide_period, &
"Period of the inflow in the tidal bay configuration.", &
units="s", default=12.0*3600.0, scale=US%s_to_T)
call get_param(param_file, mdl, "TIDAL_BAY_SSH_ANOM", CS%tide_ssh_amp, &
"Magnitude of the sea surface height anomalies at the inflow with the "//&
"tidal bay configuration.", &
units="m", default=0.1, scale=US%m_to_Z)

! Register the open boundaries.
call register_OBC(casename, param_file, OBC_Reg)
Expand All @@ -63,11 +73,11 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time)
type(time_type), intent(in) :: Time !< model time.

! The following variables are used to set up the transport in the tidal_bay example.
real :: time_sec
real :: time_sec ! Elapsed model time [T ~> s]
real :: cff_eta ! The total column thickness anomalies associated with the inflow [H ~> m or kg m-2]
real :: my_flux ! The vlume flux through the face [L2 Z T-1 ~> m3 s-1]
real :: total_area ! The total face area of the OBCs [L Z ~> m2]
real :: PI
real :: PI ! The ratio of the circumference of a circle to its diameter [nondim]
real :: flux_scale ! A scaling factor for the areas [m2 H-1 L-1 ~> nondim or m3 kg-1]
real, allocatable :: my_area(:,:) ! The total OBC inflow area [m2]
integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, n
Expand All @@ -86,10 +96,10 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time)

flux_scale = GV%H_to_m*US%L_to_m

time_sec = time_type_to_real(Time)
cff_eta = 0.1*GV%m_to_H * sin(2.0*PI*time_sec/(12.0*3600.0))
my_area=0.0
my_flux=0.0
time_sec = US%s_to_T*time_type_to_real(Time)
cff_eta = CS%tide_ssh_amp*GV%Z_to_H * sin(2.0*PI*time_sec / CS%tide_period)
my_area = 0.0
my_flux = 0.0
segment => OBC%segment(1)

do j=segment%HI%jsc,segment%HI%jec ; do I=segment%HI%IscB,segment%HI%IecB
Expand All @@ -101,7 +111,7 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, US, h, Time)
endif
enddo ; enddo
total_area = reproducing_sum(my_area)
my_flux = - CS%tide_flow*SIN(2.0*PI*time_sec/(12.0*3600.0))
my_flux = - CS%tide_flow * SIN(2.0*PI*time_sec / CS%tide_period)

do n = 1, OBC%number_of_segments
segment => OBC%segment(n)
Expand Down