18 implicit none ;
private
20 #include <MOM_memory.h>
26 character(len=40) ::
mdl =
"dense_water_initialization"
37 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
40 real,
intent(in) :: max_depth
43 real,
dimension(5) :: domain_params
44 real :: sill_frac, shelf_frac
48 call get_param(param_file,
mdl,
"DENSE_WATER_DOMAIN_PARAMS", domain_params, &
49 "Fractional widths of all the domain sections for the dense water experiment.\n"//&
50 "As a 5-element vector:\n"//&
51 " - open ocean, the section at maximum depth\n"//&
52 " - downslope, the downward overflow slope\n"//&
53 " - sill separating downslope from upslope\n"//&
54 " - upslope, the upward slope accumulating dense water\n"//&
55 " - the shelf in the dense formation region.", &
56 units=
"nondim", fail_if_missing=.true.)
57 call get_param(param_file,
mdl,
"DENSE_WATER_SILL_DEPTH", sill_frac, &
58 "Depth of the sill separating downslope from upslope, as fraction of basin depth.", &
60 call get_param(param_file,
mdl,
"DENSE_WATER_SHELF_DEPTH", shelf_frac, &
61 "Depth of the shelf region accumulating dense water for overflow, as fraction of basin depth.", &
66 domain_params(i) = domain_params(i-1) + domain_params(i)
72 x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon
74 if (x <= domain_params(1))
then
77 elseif (x <= domain_params(2))
then
79 d(i,j) = max_depth - (1.0 - sill_frac) * max_depth * &
80 (x - domain_params(1)) / (domain_params(2) - domain_params(1))
81 elseif (x <= domain_params(3))
then
83 d(i,j) = sill_frac * max_depth
84 elseif (x <= domain_params(4))
then
86 d(i,j) = sill_frac * max_depth + (shelf_frac - sill_frac) * max_depth * &
87 (x - domain_params(3)) / (domain_params(4) - domain_params(3))
90 d(i,j) = shelf_frac * max_depth
102 type(
eos_type),
pointer :: eqn_of_state
103 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(out) :: t
104 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(out) :: s
105 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)),
intent(in) :: h
106 logical,
optional,
intent(in) :: just_read_params
109 real :: mld, s_ref, s_range, t_ref
112 integer :: i, j, k, nz
116 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
118 call get_param(param_file,
mdl,
"DENSE_WATER_MLD", mld, &
119 "Depth of unstratified mixed layer as a fraction of the water column.", &
120 units=
"nondim", default=
default_mld, do_not_log=just_read)
122 call get_param(param_file,
mdl,
"S_REF", s_ref, do_not_log=.true.)
123 call get_param(param_file,
mdl,
"S_RANGE", s_range, do_not_log=.true.)
124 call get_param(param_file,
mdl,
"T_REF", t_ref, do_not_log=.true.)
126 if (just_read)
return
136 zmid = zi + 0.5 * h(i,j,k) / (gv%Z_to_H * g%max_depth)
143 s(i,j,k) = s_ref + s_range * (zmid - mld) / (1.0 - mld)
146 zi = zi + h(i,j,k) / (gv%Z_to_H * g%max_depth)
158 logical,
intent(in) :: use_ale
162 real :: west_sponge_time_scale, west_sponge_width
163 real :: east_sponge_time_scale, east_sponge_width
165 real,
dimension(SZI_(G),SZJ_(G)) :: idamp
166 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h, t, s
167 real,
dimension(SZK_(GV)+1) :: e0, eta1d
169 integer :: i, j, k, nz
170 real :: x, zi, zmid, dist
171 real :: mld, s_ref, s_range, s_dense, t_ref, sill_height
175 call get_param(param_file,
mdl,
"DENSE_WATER_WEST_SPONGE_TIME_SCALE", west_sponge_time_scale, &
176 "The time scale on the west (outflow) of the domain for restoring. If zero, the sponge is disabled.", &
177 units=
"s", default=0.)
178 call get_param(param_file,
mdl,
"DENSE_WATER_WEST_SPONGE_WIDTH", west_sponge_width, &
179 "The fraction of the domain in which the western (outflow) sponge is active.", &
180 units=
"nondim", default=0.1)
181 call get_param(param_file,
mdl,
"DENSE_WATER_EAST_SPONGE_TIME_SCALE", east_sponge_time_scale, &
182 "The time scale on the east (outflow) of the domain for restoring. If zero, the sponge is disabled.", &
183 units=
"s", default=0.)
184 call get_param(param_file,
mdl,
"DENSE_WATER_EAST_SPONGE_WIDTH", east_sponge_width, &
185 "The fraction of the domain in which the eastern (outflow) sponge is active.", &
186 units=
"nondim", default=0.1)
188 call get_param(param_file,
mdl,
"DENSE_WATER_EAST_SPONGE_SALT", s_dense, &
189 "Salt anomaly of the dense water being formed in the overflow region.", &
190 units=
"1e-3", default=4.0)
195 call get_param(param_file,
mdl,
"S_REF", s_ref, do_not_log=.true.)
196 call get_param(param_file,
mdl,
"S_RANGE", s_range, do_not_log=.true.)
197 call get_param(param_file,
mdl,
"T_REF", t_ref, do_not_log=.true.)
200 if (west_sponge_time_scale <= 0. .and. east_sponge_time_scale <= 0.)
return
207 if (g%mask2dT(i,j) > 0.)
then
209 x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon
211 if (west_sponge_time_scale > 0. .and. x < west_sponge_width)
then
212 dist = 1. - x / west_sponge_width
214 idamp(i,j) = 1. / west_sponge_time_scale * max(0., min(1., dist))
215 elseif (east_sponge_time_scale > 0. .and. x > (1. - east_sponge_width))
then
216 dist = 1. - (1. - x) / east_sponge_width
217 idamp(i,j) = 1. / east_sponge_time_scale * max(0., min(1., dist))
226 e0(k) = -g%max_depth * (real(k - 1) / real(nz))
228 e0(nz+1) = -g%max_depth
232 eta1d(nz+1) = -g%bathyT(i,j)
236 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
238 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
239 h(i,j,k) = gv%Angstrom_H
241 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
257 x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon
260 zmid = zi + 0.5 * h(i,j,k) / (gv%Z_to_H * g%max_depth)
262 if (x > (1. - east_sponge_width))
then
264 s(i,j,k) = s_ref + s_dense
268 s(i,j,k) = s_ref + s_range * (zmid - mld) / (1.0 - mld)
271 zi = zi + h(i,j,k) / (gv%Z_to_H * g%max_depth)
279 call mom_error(fatal,
"dense_water_initialize_sponges: trying to use non ALE sponge")