shelfwave_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Configures the model for the idealized shelfwave test case.
7
8use mom_domains, only : sum_across_pes
10use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
11use mom_file_parser, only : get_param, log_version, param_file_type
12use mom_grid, only : ocean_grid_type
13use mom_open_boundary, only : ocean_obc_type, obc_none, obc_direction_w
16use mom_time_manager, only : time_type, time_to_real
19
20implicit none ; private
21
22#include <MOM_memory.h>
23
24character(len=40) :: mdl = "shelfwave_initialization" !< This module's name.
25
26! The following routines are visible to the outside world
30
31!> Control structure for shelfwave open boundaries.
32type, public :: shelfwave_obc_cs ; private
33 real :: my_amp !< Amplitude of the open boundary current inflows [L T-1 ~> m s-1]
34 real :: kk !< Cross-shore wavenumber [km-1] or [m-1]
35 real :: ll !< Longshore wavenumber [km-1] or [m-1]
36 real :: alpha !< Exponential decay rate in the y-direction [km-1] or [m-1]
37 real :: omega !< Frequency of the shelf wave [T-1 ~> s-1]
38 logical :: shelfwave_correct_amplitude !< If true, SHELFWAVE_AMPLITUDE gives the actual inflow
39 !! velocity, rather than giving an overall scaling factor for the flow.
40end type shelfwave_obc_cs
41
42contains
43
44!> Add shelfwave to OBC registry.
45function register_shelfwave_obc(param_file, CS, G, US, OBC_Reg)
46 type(param_file_type), intent(in) :: param_file !< parameter file.
47 type(shelfwave_obc_cs), pointer :: cs !< shelfwave control structure.
48 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
49 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
50 type(obc_registry_type), pointer :: obc_reg !< Open boundary condition registry.
51 logical :: register_shelfwave_obc
52
53 ! Local variables
54 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
55 character(len=32) :: casename = "shelfwave" !< This case's name.
56 real :: jj ! Cross-shore wave mode [nondim]
57 real :: f0 ! Coriolis parameter [T-1 ~> s-1]
58 real :: lx ! Long-shore length scale of bathymetry [km] or [m]
59 real :: ly ! Cross-shore length scale [km] or [m]
60 real :: default_amp ! The default velocity amplitude [m s-1] or amplitude scaling factor [nondim]
61
62 pi = 4.0*atan(1.0)
63
64 if (associated(cs)) then
65 call mom_error(warning, "register_shelfwave_OBC called with an "// &
66 "associated control structure.")
67 return
68 endif
69 allocate(cs)
70
71 ! Register the tracer for horizontal advection & diffusion.
72 call register_obc(casename, param_file, obc_reg)
73 call get_param(param_file, mdl, "F_0", f0, &
74 default=0.0, units="s-1", scale=us%T_to_s, do_not_log=.true.)
75 call get_param(param_file, mdl,"SHELFWAVE_X_WAVELENGTH", lx, &
76 "Length scale of shelfwave in x-direction.",&
77 units=g%x_ax_unit_short, default=100.)
78 call get_param(param_file, mdl, "SHELFWAVE_Y_LENGTH_SCALE", ly, &
79 "Length scale of exponential dropoff of topography in the y-direction.", &
80 units=g%y_ax_unit_short, default=50.)
81 call get_param(param_file, mdl, "SHELFWAVE_Y_MODE", jj, &
82 "Cross-shore wave mode.", &
83 units="nondim", default=1.)
84 call get_param(param_file, mdl, "SHELFWAVE_CORRECT_AMPLITUDE", cs%shelfwave_correct_amplitude, &
85 "If true, SHELFWAVE_AMPLITUDE gives the actual inflow velocity, rather than giving "//&
86 "an overall scaling factor for the flow.", default=.true.)
87 default_amp = 1.0 ; if (cs%shelfwave_correct_amplitude) default_amp = 0.1
88 call get_param(param_file, mdl, "SHELFWAVE_AMPLITUDE", cs%my_amp, &
89 "Amplitude of the open boundary current inflows in the shelfwave configuration.", &
90 units="m s-1", default=default_amp, scale=us%m_s_to_L_T)
91
92 cs%alpha = 1. / ly
93 cs%ll = 2. * pi / lx
94 cs%kk = jj * pi / g%len_lat
95 cs%omega = 2 * cs%alpha * f0 * cs%ll / &
96 (cs%kk*cs%kk + cs%alpha*cs%alpha + cs%ll*cs%ll)
98
99end function register_shelfwave_obc
100
101!> Clean up the shelfwave OBC from registry.
102subroutine shelfwave_obc_end(CS)
103 type(shelfwave_obc_cs), pointer :: cs !< shelfwave control structure.
104
105 if (associated(cs)) then
106 deallocate(cs)
107 endif
108end subroutine shelfwave_obc_end
109
110!> Initialization of topography.
111subroutine shelfwave_initialize_topography( D, G, param_file, max_depth, US )
112 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
113 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
114 intent(out) :: d !< Ocean bottom depth [Z ~> m]
115 type(param_file_type), intent(in) :: param_file !< Parameter file structure
116 real, intent(in) :: max_depth !< Maximum model depth [Z ~> m]
117 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
118
119 ! Local variables
120 real :: y ! Position relative to the southern boundary [km] or [m] or [degrees_N]
121 real :: rly ! Exponential decay rate of the topography [km-1] or [m-1] or [degrees_N-1]
122 real :: ly ! Exponential decay lengthscale of the topography [km] or [m] or [degrees_N]
123 real :: h0 ! The minimum depth of the ocean [Z ~> m]
124 integer :: i, j
125
126 call get_param(param_file, mdl,"SHELFWAVE_Y_LENGTH_SCALE", ly, &
127 units=g%y_ax_unit_short, default=50., do_not_log=.true.)
128 call get_param(param_file, mdl,"MINIMUM_DEPTH", h0, &
129 units="m", default=10., scale=us%m_to_Z, do_not_log=.true.)
130
131 rly = 0. ; if (ly>0.) rly = 1. / ly
132
133 do j=g%jsc,g%jec ; do i=g%isc,g%iec
134 ! Compute normalized zonal coordinates (x,y=0 at center of domain)
135 y = ( g%geoLatT(i,j) - g%south_lat )
136 d(i,j) = h0 * exp(2 * rly * y)
137 enddo ; enddo
138
140
141!> This subroutine sets the properties of flow at open boundary conditions.
142subroutine shelfwave_set_obc_data(OBC, CS, G, GV, US, h, Time)
143 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
144 !! whether, where, and what open boundary
145 !! conditions are used.
146 type(shelfwave_obc_cs), pointer :: cs !< tidal bay control structure.
147 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
148 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
149 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
150 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2]
151 type(time_type), intent(in) :: time !< model time.
152
153 ! The following variables are used to set up the transport in the shelfwave example.
154 real :: time_sec ! The time in the run [T ~> s]
155 real :: cos_wt, sin_wt ! Cosine and sine associated with the propagating x-direction structure [nondim]
156 real :: cos_ky, sin_ky ! Cosine and sine associated with the y-direction structure [nondim]
157 real :: x ! Position relative to the western boundary [km] or [m] or [degrees_E]
158 real :: y ! Position relative to the southern boundary [km] or [m] or [degrees_N]
159 real :: i_yscale ! A factor to give the correct inflow velocity [km-1] or [m-1] or [degrees_N-1] or
160 ! to compensate for the variable units of the y-coordinate [km axis_unit-1], usually 1 [nondim]
161 real :: my_amp ! Amplitude of the open boundary current inflows, including sign changes
162 ! to account for grid rotation [L T-1 ~> m s-1]
163 integer :: i, j, is, ie, js, je, n
164 integer :: turns ! Number of index quarter turns
165 type(obc_segment_type), pointer :: segment => null()
166
167 if (.not.associated(obc)) return
168
169 turns = modulo(g%HI%turns, 4)
170 my_amp = cs%my_amp ; if ((turns==2) .or. (turns==3)) my_amp = -cs%my_amp
171
172 time_sec = time_to_real(time, scale=us%s_to_T)
173 if (cs%shelfwave_correct_amplitude) then
174 ! This makes the units and edge value of normal_vel_bt the same as my_amp.
175 i_yscale = 1.0 / cs%kk
176 else ! This preserves the previous answers.
177 if (g%grid_unit_to_L == 0.0) call mom_error(fatal, &
178 "shelfwave_set_OBC_data requires the use of Cartesian coordinates.")
179 i_yscale = (1.0e3 * us%m_to_L) / g%grid_unit_to_L
180 endif
181 do n = 1, obc%number_of_segments
182 segment => obc%segment(n)
183 if (.not. segment%on_pe) cycle
184 if (rotate_obc_segment_direction(segment%direction, -turns) /= obc_direction_w) cycle
185
186 if (segment%is_E_or_W) then
187 ! segment thicknesses are defined at cell face centers.
188 is = segment%HI%isdB ; ie = segment%HI%iedB
189 js = segment%HI%jsd ; je = segment%HI%jed
190 else
191 is = segment%HI%isd ; ie = segment%HI%ied
192 js = segment%HI%jsdB ; je = segment%HI%jedB
193 endif
194
195 do j=js,je ; do i=is,ie
196 if (segment%is_E_or_W) then
197 x = g%geoLonCu(i,j) - g%west_lon
198 y = g%geoLatCu(i,j) - g%south_lat
199 else
200 x = g%geoLonCv(i,j) - g%west_lon
201 y = g%geoLatCv(i,j) - g%south_lat
202 endif
203 sin_wt = sin(cs%ll*x - cs%omega*time_sec)
204 cos_wt = cos(cs%ll*x - cs%omega*time_sec)
205 sin_ky = sin(cs%kk * y)
206 cos_ky = cos(cs%kk * y)
207 segment%normal_vel_bt(i,j) = my_amp * exp(- cs%alpha * y) * cos_wt * &
208 ((cs%alpha * sin_ky + cs%kk * cos_ky) * i_yscale)
209! segment%tangential_vel_bt(I,j) = my_amp * (CS%ll * I_yscale) * exp(- CS%alpha * y) * sin_wt * sin_ky
210! segment%vorticity_bt(I,j) = my_amp * exp(- CS%alpha * y) * cos_wt * sin_ky * &
211! ((CS%ll**2 + CS%kk**2 + CS%alpha**2) * (I_yscale / G%grid_unit_to_L))
212 enddo ; enddo
213 enddo
214
215end subroutine shelfwave_set_obc_data
216