diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4659b685e5..4d16260f7c 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -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) diff --git a/src/core/MOM_boundary_update.F90 b/src/core/MOM_boundary_update.F90 index 17712491c4..2e25af2460 100644 --- a/src/core/MOM_boundary_update.F90 +++ b/src/core/MOM_boundary_update.F90 @@ -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 @@ -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 @@ -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) & diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 0cb81e9978..bd76f5a9aa 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -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 @@ -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. diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index d7cee50c99..3eba7c81bc 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -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". @@ -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 @@ -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. @@ -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 @@ -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 @@ -201,18 +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 / 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 - 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) @@ -225,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 @@ -261,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 @@ -277,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) ) ) @@ -294,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 @@ -322,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 @@ -336,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))) ) diff --git a/src/user/dyed_channel_initialization.F90 b/src/user/dyed_channel_initialization.F90 index 4c633ebdc9..317ed4ac21 100644 --- a/src/user/dyed_channel_initialization.F90 +++ b/src/user/dyed_channel_initialization.F90 @@ -14,6 +14,7 @@ module dyed_channel_initialization use MOM_time_manager, only : time_type, time_type_to_real use MOM_tracer_registry, only : tracer_registry_type, tracer_name_lookup use MOM_tracer_registry, only : tracer_type +use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -26,9 +27,9 @@ module dyed_channel_initialization !> Control structure for dyed-channel open boundaries. type, public :: dyed_channel_OBC_CS ; private - real :: zonal_flow = 8.57 !< Mean inflow - real :: tidal_amp = 0.0 !< Sloshing amplitude - real :: frequency = 0.0 !< Sloshing frequency + real :: zonal_flow = 8.57 !< Mean inflow [L T-1 ~> m s-1] + real :: tidal_amp = 0.0 !< Sloshing amplitude [L T-1 ~> m s-1] + real :: frequency = 0.0 !< Sloshing frequency [T-1 ~> s-1] end type dyed_channel_OBC_CS integer :: ntr = 0 !< Number of dye tracers @@ -37,9 +38,10 @@ module dyed_channel_initialization contains !> Add dyed channel to OBC registry. -function register_dyed_channel_OBC(param_file, CS, OBC_Reg) +function register_dyed_channel_OBC(param_file, CS, US, OBC_Reg) type(param_file_type), intent(in) :: param_file !< parameter file. type(dyed_channel_OBC_CS), pointer :: CS !< Dyed channel 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 logical :: register_dyed_channel_OBC @@ -55,13 +57,13 @@ function register_dyed_channel_OBC(param_file, CS, OBC_Reg) call get_param(param_file, mdl, "CHANNEL_MEAN_FLOW", CS%zonal_flow, & "Mean zonal flow imposed at upstream open boundary.", & - units="m/s", default=8.57) + units="m/s", default=8.57, scale=US%m_s_to_L_T) call get_param(param_file, mdl, "CHANNEL_TIDAL_AMP", CS%tidal_amp, & "Sloshing amplitude imposed at upstream open boundary.", & - units="m/s", default=0.0) + units="m/s", default=0.0, scale=US%m_s_to_L_T) call get_param(param_file, mdl, "CHANNEL_FLOW_FREQUENCY", CS%frequency, & "Frequency of oscillating zonal flow.", & - units="s-1", default=0.0) + units="s-1", default=0.0, scale=US%T_to_s) ! Register the open boundaries. call register_OBC(casename, param_file, OBC_Reg) @@ -142,7 +144,9 @@ subroutine dyed_channel_update_flow(OBC, CS, G, GV, Time) ! Local variables character(len=40) :: mdl = "dyed_channel_update_flow" ! This subroutine's name. character(len=80) :: name - real :: flow, time_sec, PI + real :: flow ! The OBC velocity [L T-1 ~> m s-1] + real :: PI ! 3.1415926535... + real :: time_sec ! The elapsed time since the start of the calendar [T ~> s] integer :: i, j, k, l, itt, isd, ied, jsd, jed, m, n integer :: IsdB, IedB, JsdB, JedB type(OBC_segment_type), pointer :: segment => NULL() @@ -150,7 +154,7 @@ subroutine dyed_channel_update_flow(OBC, CS, G, GV, Time) if (.not.associated(OBC)) call MOM_error(FATAL, 'dyed_channel_initialization.F90: '// & 'dyed_channel_update_flow() was called but OBC type was not initialized!') - time_sec = time_type_to_real(Time) + time_sec = G%US%s_to_T * time_type_to_real(Time) PI = 4.0*atan(1.0) do l=1, OBC%number_of_segments @@ -163,9 +167,9 @@ subroutine dyed_channel_update_flow(OBC, CS, G, GV, Time) jsd = segment%HI%jsd ; jed = segment%HI%jed IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB if (CS%frequency == 0.0) then - flow = G%US%m_s_to_L_T*CS%zonal_flow + flow = CS%zonal_flow else - flow = G%US%m_s_to_L_T*CS%zonal_flow + CS%tidal_amp * cos(2 * PI * CS%frequency * time_sec) + flow = CS%zonal_flow + CS%tidal_amp * cos(2 * PI * CS%frequency * time_sec) endif do k=1,GV%ke do j=jsd,jed ; do I=IsdB,IedB diff --git a/src/user/shelfwave_initialization.F90 b/src/user/shelfwave_initialization.F90 index 7bf6aebf59..041d77d9f9 100644 --- a/src/user/shelfwave_initialization.F90 +++ b/src/user/shelfwave_initialization.F90 @@ -30,20 +30,21 @@ module shelfwave_initialization type, public :: shelfwave_OBC_CS ; private real :: Lx = 100.0 !< Long-shore length scale of bathymetry. real :: Ly = 50.0 !< Cross-shore length scale. - real :: f0 = 1.e-4 !< Coriolis parameter. + real :: f0 = 1.e-4 !< Coriolis parameter [T-1 ~> s-1] real :: jj = 1 !< Cross-shore wave mode. real :: kk !< Parameter. real :: ll !< Longshore wavenumber. real :: alpha !< 1/Ly. - real :: omega !< Frequency. + real :: omega !< Frequency of the shelf wave [T-1 ~> s-1] end type shelfwave_OBC_CS contains !> Add shelfwave to OBC registry. -function register_shelfwave_OBC(param_file, CS, OBC_Reg) +function register_shelfwave_OBC(param_file, CS, US, OBC_Reg) type(param_file_type), intent(in) :: param_file !< parameter file. type(shelfwave_OBC_CS), pointer :: CS !< shelfwave 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_shelfwave_OBC ! Local variables @@ -62,18 +63,20 @@ function register_shelfwave_OBC(param_file, CS, OBC_Reg) ! Register the tracer for horizontal advection & diffusion. call register_OBC(casename, param_file, OBC_Reg) - call get_param(param_file, mdl,"F_0",CS%f0, & - do_not_log=.true.) - call get_param(param_file, mdl,"LENLAT",len_lat, & + call get_param(param_file, mdl, "F_0", CS%f0, & + default=0.0, units="s-1", scale=US%T_to_s, do_not_log=.true.) + call get_param(param_file, mdl, "LENLAT", len_lat, & do_not_log=.true.) call get_param(param_file, mdl,"SHELFWAVE_X_WAVELENGTH",CS%Lx, & "Length scale of shelfwave in x-direction.",& units="Same as x,y", default=100.) - call get_param(param_file, mdl,"SHELFWAVE_Y_LENGTH_SCALE",CS%Ly, & +! units="km", default=100.0, scale=1.0e3*US%m_to_L) + call get_param(param_file, mdl, "SHELFWAVE_Y_LENGTH_SCALE", CS%Ly, & "Length scale of exponential dropoff of topography "//& "in the y-direction.", & units="Same as x,y", default=50.) - call get_param(param_file, mdl,"SHELFWAVE_Y_MODE",CS%jj, & +! units="km", default=50.0, scale=1.0e3*US%m_to_L) + call get_param(param_file, mdl, "SHELFWAVE_Y_MODE", CS%jj, & "Cross-shore wave mode.", & units="nondim", default=1.) CS%alpha = 1. / CS%Ly @@ -126,19 +129,23 @@ subroutine shelfwave_initialize_topography( D, G, param_file, max_depth, US ) end subroutine shelfwave_initialize_topography !> This subroutine sets the properties of flow at open boundary conditions. -subroutine shelfwave_set_OBC_data(OBC, CS, G, GV, h, Time) - type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies - !! whether, where, and what open boundary - !! conditions are used. - type(shelfwave_OBC_CS), pointer :: CS !< tidal bay 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 +subroutine shelfwave_set_OBC_data(OBC, CS, G, GV, US, h, Time) + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies + !! whether, where, and what open boundary + !! conditions are used. + type(shelfwave_OBC_CS), pointer :: CS !< tidal bay 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 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness. type(time_type), intent(in) :: Time !< model time. ! The following variables are used to set up the transport in the shelfwave example. - real :: my_amp, time_sec - real :: cos_wt, cos_ky, sin_wt, sin_ky, omega, alpha + real :: my_amp ! Amplitude of the open boundary current inflows [L T-1 ~> m s-1] + real :: time_sec ! The time in the run [T ~> s] + real :: cos_wt, cos_ky, sin_wt, sin_ky + real :: omega ! Frequency of the shelf wave [T-1 ~> s-1] + real :: alpha real :: x, y, jj, kk, ll character(len=40) :: mdl = "shelfwave_set_OBC_data" ! This subroutine's name. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, n @@ -151,10 +158,10 @@ subroutine shelfwave_set_OBC_data(OBC, CS, G, GV, h, Time) if (.not.associated(OBC)) return - time_sec = time_type_to_real(Time) + time_sec = US%s_to_T*time_type_to_real(Time) omega = CS%omega alpha = CS%alpha - my_amp = 1.0 + my_amp = 1.0*G%US%m_s_to_L_T jj = CS%jj kk = CS%kk ll = CS%ll @@ -172,9 +179,9 @@ subroutine shelfwave_set_OBC_data(OBC, CS, G, GV, h, Time) cos_wt = cos(ll*x - omega*time_sec) sin_ky = sin(kk * y) cos_ky = cos(kk * y) - segment%normal_vel_bt(I,j) = G%US%m_s_to_L_T*my_amp * exp(- alpha * y) * cos_wt * & + segment%normal_vel_bt(I,j) = my_amp * exp(- alpha * y) * cos_wt * & (alpha * sin_ky + kk * cos_ky) -! segment%tangential_vel_bt(I,j) = G%US%m_s_to_L_T*my_amp * ll * exp(- alpha * y) * sin_wt * sin_ky +! segment%tangential_vel_bt(I,j) = my_amp * ll * exp(- alpha * y) * sin_wt * sin_ky ! segment%vorticity_bt(I,j) = my_amp * exp(- alpha * y) * cos_wt * sin_ky& ! (ll*ll + kk*kk + alpha*alpha) enddo ; enddo diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90 index e6db433f60..b3c8f45843 100644 --- a/src/user/tidal_bay_initialization.F90 +++ b/src/user/tidal_bay_initialization.F90 @@ -12,6 +12,7 @@ module tidal_bay_initialization use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE use MOM_open_boundary, only : OBC_segment_type, register_OBC use MOM_open_boundary, only : OBC_registry_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 @@ -24,18 +25,20 @@ 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. + real :: tide_flow = 3.0e6 !< Maximum tidal flux [L2 Z T-1 ~> m3 s-1] end type tidal_bay_OBC_CS contains !> Add tidal bay to OBC registry. -function register_tidal_bay_OBC(param_file, CS, OBC_Reg) +function register_tidal_bay_OBC(param_file, CS, US, OBC_Reg) type(param_file_type), intent(in) :: param_file !< parameter file. type(tidal_bay_OBC_CS), pointer :: CS !< tidal bay 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_tidal_bay_OBC character(len=32) :: casename = "tidal bay" !< This case's name. + character(len=40) :: mdl = "tidal_bay_initialization" ! This module's name. if (associated(CS)) then call MOM_error(WARNING, "register_tidal_bay_OBC called with an "// & @@ -44,6 +47,10 @@ function register_tidal_bay_OBC(param_file, CS, OBC_Reg) endif allocate(CS) + call get_param(param_file, mdl, "TIDAL_BAY_FLOW", CS%tide_flow, & + "Maximum total tidal volume flux.", & + units="m3 s-1", default=3.0d6, scale=US%m_s_to_L_T*US%m_to_L*US%m_to_Z) + ! Register the open boundaries. call register_OBC(casename, param_file, OBC_Reg) register_tidal_bay_OBC = .true. @@ -67,14 +74,17 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, h, Time) type(tidal_bay_OBC_CS), pointer :: CS !< tidal bay 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 - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness. + 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 tidal_bay example. - real :: time_sec, cff - real :: my_flux, total_area + real :: time_sec + 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, allocatable :: my_area(:,:) + 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] character(len=40) :: mdl = "tidal_bay_set_OBC_data" ! This subroutine's name. integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, n integer :: IsdB, IedB, JsdB, JedB @@ -90,8 +100,10 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, h, Time) allocate(my_area(1:1,js:je)) + flux_scale = GV%H_to_m*G%US%L_to_m + time_sec = time_type_to_real(Time) - cff = 0.1*sin(2.0*PI*time_sec/(12.0*3600.0)) + 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 segment => OBC%segment(1) @@ -99,7 +111,8 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, h, Time) do j=segment%HI%jsc,segment%HI%jec ; do I=segment%HI%IscB,segment%HI%IecB if (OBC%segnum_u(I,j) /= OBC_NONE) then do k=1,nz - my_area(1,j) = my_area(1,j) + h(I,j,k)*G%US%L_to_m*G%dyCu(I,j) + ! This area has to be in MKS units to work with reproducing_sum. + my_area(1,j) = my_area(1,j) + h(I,j,k)*flux_scale*G%dyCu(I,j) enddo endif enddo ; enddo @@ -111,8 +124,8 @@ subroutine tidal_bay_set_OBC_data(OBC, CS, G, GV, h, Time) if (.not. segment%on_pe) cycle - segment%normal_vel_bt(:,:) = G%US%m_s_to_L_T*my_flux/total_area - segment%eta(:,:) = cff + segment%normal_vel_bt(:,:) = my_flux / (G%US%m_to_Z*G%US%m_to_L*total_area) + segment%eta(:,:) = cff_eta enddo ! end segment loop