17 implicit none ;
private
19 #include <MOM_memory.h>
21 character(len=40) ::
mdl =
"shelfwave_initialization"
49 real :: kk, ll, pi, len_lat
51 character(len=32) :: casename =
"shelfwave"
55 if (
associated(cs))
then
56 call mom_error(warning,
"register_shelfwave_OBC called with an "// &
57 "associated control structure.")
63 call register_obc(casename, param_file, obc_reg)
68 call get_param(param_file,
mdl,
"SHELFWAVE_X_WAVELENGTH",cs%Lx, &
69 "Length scale of shelfwave in x-direction.",&
70 units=
"Same as x,y", default=100.)
71 call get_param(param_file,
mdl,
"SHELFWAVE_Y_LENGTH_SCALE",cs%Ly, &
72 "Length scale of exponential dropoff of topography "//&
73 "in the y-direction.", &
74 units=
"Same as x,y", default=50.)
75 call get_param(param_file,
mdl,
"SHELFWAVE_Y_MODE",cs%jj, &
76 "Cross-shore wave mode.", &
77 units=
"nondim", default=1.)
79 cs%ll = 2. * pi / cs%Lx
80 cs%kk = cs%jj * pi / len_lat
81 cs%omega = 2 * cs%alpha * cs%f0 * cs%ll / &
82 (cs%kk*cs%kk + cs%alpha*cs%alpha + cs%ll*cs%ll)
91 if (
associated(cs))
then
99 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
102 real,
intent(in) :: max_depth
108 real :: y, rly, ly, h0
110 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
112 call get_param(param_file,
mdl,
"SHELFWAVE_Y_LENGTH_SCALE",ly, &
113 default=50., do_not_log=.true.)
115 default=10., units=
"m", scale=m_to_z, do_not_log=.true.)
117 rly = 0. ;
if (ly>0.) rly = 1. / ly
119 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
121 y = ( g%geoLatT(i,j) - g%south_lat )
122 d(i,j) = h0 * exp(2 * rly * y)
134 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
135 type(time_type),
intent(in) :: time
138 real :: my_amp, time_sec
139 real :: cos_wt, cos_ky, sin_wt, sin_ky, omega, alpha
140 real :: x, y, jj, kk, ll
141 character(len=40) ::
mdl =
"shelfwave_set_OBC_data"
142 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, n
143 integer :: isdb, iedb, jsdb, jedb
146 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
147 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
148 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
150 if (.not.
associated(obc))
return
152 time_sec = time_type_to_real(time)
159 do n = 1, obc%number_of_segments
160 segment => obc%segment(n)
161 if (.not. segment%on_pe) cycle
162 if (segment%direction /= obc_direction_w) cycle
164 isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
165 jsd = segment%HI%jsd ; jed = segment%HI%jed
166 do j=jsd,jed ;
do i=isdb,iedb
167 x = g%geoLonCu(i,j) - g%west_lon
168 y = g%geoLatCu(i,j) - g%south_lat
169 sin_wt = sin(ll*x - omega*time_sec)
170 cos_wt = cos(ll*x - omega*time_sec)
173 segment%normal_vel_bt(i,j) = g%US%m_s_to_L_T*my_amp * exp(- alpha * y) * cos_wt * &
174 (alpha * sin_ky + kk * cos_ky)