dense_water_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!> Initialization routines for the dense water formation
6!! and overflow experiment.
8
9use mom_ale_sponge, only : ale_sponge_cs, set_up_ale_sponge_field, initialize_ale_sponge
11use mom_eos, only : eos_type
12use mom_error_handler, only : mom_error, fatal
13use mom_file_parser, only : get_param, param_file_type
14use mom_grid, only : ocean_grid_type
15use mom_sponge, only : sponge_cs
19
20implicit none ; private
21
22#include <MOM_memory.h>
23
27
28character(len=40) :: mdl = "dense_water_initialization" !< Module name
29
30real, parameter :: default_sill = 0.2 !< Default depth of the sill [nondim]
31real, parameter :: default_shelf = 0.4 !< Default depth of the shelf [nondim]
32real, parameter :: default_mld = 0.25 !< Default depth of the mixed layer [nondim]
33
34contains
35
36!> Initialize the topography field for the dense water experiment
37subroutine dense_water_initialize_topography(D, G, param_file, max_depth)
38 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
39 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
40 intent(out) :: d !< Ocean bottom depth [Z ~> m]
41 type(param_file_type), intent(in) :: param_file !< Parameter file structure
42 real, intent(in) :: max_depth !< Maximum ocean depth [Z ~> m]
43
44 ! Local variables
45 real, dimension(5) :: domain_params ! nondimensional widths of all domain sections [nondim]
46 real :: sill_frac ! Depth of the sill separating downslope from upslope, as a fraction of
47 ! the basin depth [nondim]
48 real :: shelf_frac ! Depth of the shelf region accumulating dense water for overflow,
49 ! as a fraction the basin depth [nondim]
50 real :: x ! Horizontal position normalized by the domain width [nondim]
51 integer :: i, j
52
53 call get_param(param_file, mdl, "DENSE_WATER_DOMAIN_PARAMS", domain_params, &
54 "Fractional widths of all the domain sections for the dense water experiment.\n"//&
55 "As a 5-element vector:\n"//&
56 " - open ocean, the section at maximum depth\n"//&
57 " - downslope, the downward overflow slope\n"//&
58 " - sill separating downslope from upslope\n"//&
59 " - upslope, the upward slope accumulating dense water\n"//&
60 " - the shelf in the dense formation region.", &
61 units="nondim", fail_if_missing=.true.)
62 call get_param(param_file, mdl, "DENSE_WATER_SILL_DEPTH", sill_frac, &
63 "Depth of the sill separating downslope from upslope, as fraction of basin depth.", &
64 units="nondim", default=default_sill)
65 call get_param(param_file, mdl, "DENSE_WATER_SHELF_DEPTH", shelf_frac, &
66 "Depth of the shelf region accumulating dense water for overflow, as fraction of basin depth.", &
67 units="nondim", default=default_shelf)
68
69 do i = 2, 5
70 ! turn widths into positions
71 domain_params(i) = domain_params(i-1) + domain_params(i)
72 enddo
73
74 do j = g%jsc,g%jec
75 do i = g%isc,g%iec
76 ! compute normalised zonal coordinate
77 x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon
78
79 if (x <= domain_params(1)) then
80 ! open ocean region
81 d(i,j) = max_depth
82 elseif (x <= domain_params(2)) then
83 ! downslope region, linear
84 d(i,j) = max_depth - (1.0 - sill_frac) * max_depth * &
85 (x - domain_params(1)) / (domain_params(2) - domain_params(1))
86 elseif (x <= domain_params(3)) then
87 ! sill region
88 d(i,j) = sill_frac * max_depth
89 elseif (x <= domain_params(4)) then
90 ! upslope region
91 d(i,j) = sill_frac * max_depth + (shelf_frac - sill_frac) * max_depth * &
92 (x - domain_params(3)) / (domain_params(4) - domain_params(3))
93 else
94 ! shelf region
95 d(i,j) = shelf_frac * max_depth
96 endif
97 enddo
98 enddo
99
101
102!> Initialize the temperature and salinity for the dense water experiment
103subroutine dense_water_initialize_ts(G, GV, US, param_file, T, S, h, just_read)
104 type(ocean_grid_type), intent(in) :: g !< Horizontal grid control structure
105 type(verticalgrid_type), intent(in) :: gv !< Vertical grid control structure
106 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
107 type(param_file_type), intent(in) :: param_file !< Parameter file structure
108 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t !< Output temperature [C ~> degC]
109 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s !< Output salinity [S ~> ppt]
110 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [Z ~> m]
111 logical, intent(in) :: just_read !< If true, this call will
112 !! only read parameters without changing T & S.
113 ! Local variables
114 real :: mld ! The initial mixed layer depth as a fraction of the maximum depth [nondim]
115 real :: s_ref, s_range ! The reference salinity and its range in the initial conditions [S ~> ppt]
116 real :: t_ref ! The reference temperature [C ~> degC]
117 real :: zi, zmid ! Depths from the surface nondimensionalized by the maximum depth [nondim]
118 integer :: i, j, k, nz
119
120 nz = gv%ke
121
122 call get_param(param_file, mdl, "DENSE_WATER_MLD", mld, &
123 "Depth of unstratified mixed layer as a fraction of the water column.", &
124 units="nondim", default=default_mld, do_not_log=just_read)
125 call get_param(param_file, mdl, "S_REF", s_ref, 'Reference salinity', &
126 default=35.0, units="ppt", scale=us%ppt_to_S, do_not_log=just_read)
127 call get_param(param_file, mdl,"T_REF", t_ref, 'Reference temperature', &
128 units='degC', scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
129 call get_param(param_file, mdl,"S_RANGE", s_range, 'Initial salinity range', &
130 units="ppt", default=2.0, scale=us%ppt_to_S, do_not_log=just_read)
131
132 if (just_read) return ! All run-time parameters have been read, so return.
133
134 ! uniform temperature everywhere
135 t(:,:,:) = t_ref
136
137 do j = g%jsc,g%jec
138 do i = g%isc,g%iec
139 zi = 0.
140 do k = 1,nz
141 ! nondimensional middle of layer
142 zmid = zi + 0.5 * h(i,j,k) / g%max_depth
143
144 if (zmid < mld) then
145 ! use reference salinity in the mixed layer
146 s(i,j,k) = s_ref
147 else
148 ! linear between bottom of mixed layer and bottom
149 s(i,j,k) = s_ref + s_range * (zmid - mld) / (1.0 - mld)
150 endif
151
152 zi = zi + h(i,j,k) / g%max_depth
153 enddo
154 enddo
155 enddo
156end subroutine dense_water_initialize_ts
157
158!> Initialize the restoring sponges for the dense water experiment
159subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp)
160 type(ocean_grid_type), intent(in) :: g !< Horizontal grid control structure
161 type(verticalgrid_type), intent(in) :: gv !< Vertical grid control structure
162 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
163 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
164 real, dimension(SZI_(G),SZJ_(G)), &
165 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
166 type(param_file_type), intent(in) :: param_file !< Parameter file structure
167 logical, intent(in) :: use_ale !< ALE flag
168 type(sponge_cs), pointer :: csp !< Layered sponge control structure pointer
169 type(ale_sponge_cs), pointer :: acsp !< ALE sponge control structure pointer
170
171 ! Local variables
172 real :: west_sponge_time_scale, east_sponge_time_scale ! Sponge timescales [T ~> s]
173 real :: west_sponge_width ! The fraction of the domain in which the western (outflow) sponge is active [nondim]
174 real :: east_sponge_width ! The fraction of the domain in which the eastern (outflow) sponge is active [nondim]
175
176 real, dimension(SZI_(G),SZJ_(G)) :: idamp ! inverse damping timescale [T-1 ~> s-1]
177 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! sponge layer thicknesses in height units [Z ~> m]
178 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: t ! sponge temperature [C ~> degC]
179 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: s ! sponge salinity [S ~> ppt]
180 real, dimension(SZK_(GV)+1) :: e0, eta1d ! interface positions for ALE sponge [Z ~> m]
181 real :: x ! Horizontal position normalized by the domain width [nondim]
182 real :: zi, zmid ! Depths from the surface nondimensionalized by the maximum depth [nondim]
183 real :: dist ! Distance from the edge of a sponge normalized by the width of that sponge [nondim]
184 real :: mld ! The initial mixed layer depth as a fraction of the maximum depth [nondim]
185 real :: s_ref, s_range ! The reference salinity and its range in the initial conditions [S ~> ppt]
186 real :: s_dense ! The salinity of the dense water being formed on the shelf [S ~> ppt]
187 real :: t_ref ! The reference temperature [C ~> degC]
188 real :: sill_frac ! Fractional depths of the sill, relative to the maximum depth [nondim]
189 integer :: i, j, k, nz
190
191 nz = gv%ke
192
193 call get_param(param_file, mdl, "DENSE_WATER_WEST_SPONGE_TIME_SCALE", west_sponge_time_scale, &
194 "The time scale on the west (outflow) of the domain for restoring. "//&
195 "If zero, the sponge is disabled.", units="s", default=0., scale=us%s_to_T)
196 call get_param(param_file, mdl, "DENSE_WATER_WEST_SPONGE_WIDTH", west_sponge_width, &
197 "The fraction of the domain in which the western (outflow) sponge is active.", &
198 units="nondim", default=0.1)
199 call get_param(param_file, mdl, "DENSE_WATER_EAST_SPONGE_TIME_SCALE", east_sponge_time_scale, &
200 "The time scale on the east (outflow) of the domain for restoring. "//&
201 "If zero, the sponge is disabled.", units="s", default=0., scale=us%s_to_T)
202 call get_param(param_file, mdl, "DENSE_WATER_EAST_SPONGE_WIDTH", east_sponge_width, &
203 "The fraction of the domain in which the eastern (outflow) sponge is active.", &
204 units="nondim", default=0.1)
205 call get_param(param_file, mdl, "DENSE_WATER_EAST_SPONGE_SALT", s_dense, &
206 "Salt anomaly of the dense water being formed in the overflow region.", &
207 units="ppt", default=4.0, scale=us%ppt_to_S)
208
209 call get_param(param_file, mdl, "DENSE_WATER_MLD", mld, &
210 units="nondim", default=default_mld, do_not_log=.true.)
211 call get_param(param_file, mdl, "DENSE_WATER_SILL_DEPTH", sill_frac, &
212 units="nondim", default=default_sill, do_not_log=.true.)
213
214 call get_param(param_file, mdl, "S_REF", s_ref, &
215 units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=.true.)
216 call get_param(param_file, mdl, "S_RANGE", s_range, &
217 units="ppt", default=2.0, scale=us%ppt_to_S, do_not_log=.true.)
218 call get_param(param_file, mdl, "T_REF", t_ref, &
219 units='degC', scale=us%degC_to_C, fail_if_missing=.true., do_not_log=.true.)
220
221 ! no active sponges
222 if (west_sponge_time_scale <= 0. .and. east_sponge_time_scale <= 0.) return
223
224 ! everywhere is initially unsponged
225 idamp(:,:) = 0.0
226
227 do j = g%jsc, g%jec
228 do i = g%isc,g%iec
229 if (g%mask2dT(i,j) > 0.) then
230 ! nondimensional x position
231 x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon
232
233 if (west_sponge_time_scale > 0. .and. x < west_sponge_width) then
234 dist = 1. - x / west_sponge_width
235 ! scale restoring by depth into sponge
236 idamp(i,j) = 1. / west_sponge_time_scale * max(0., min(1., dist))
237 elseif (east_sponge_time_scale > 0. .and. x > (1. - east_sponge_width)) then
238 dist = 1. - (1. - x) / east_sponge_width
239 idamp(i,j) = 1. / east_sponge_time_scale * max(0., min(1., dist))
240 endif
241 endif
242 enddo
243 enddo
244
245 if (use_ale) then
246 ! construct a uniform grid for the sponge
247 do k = 1,nz
248 e0(k) = -g%max_depth * (real(k - 1) / real(nz))
249 enddo
250 e0(nz+1) = -g%max_depth
251
252 do j = g%jsc,g%jec
253 do i = g%isc,g%iec
254 eta1d(nz+1) = -depth_tot(i,j)
255 do k = nz,1,-1
256 eta1d(k) = e0(k)
257
258 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
259 ! is this layer vanished?
260 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
261 dz(i,j,k) = gv%Angstrom_Z
262 else
263 dz(i,j,k) = eta1d(k) - eta1d(k+1)
264 endif
265 enddo
266 enddo
267 enddo
268
269 ! construct temperature and salinity for the sponge
270 ! start with initial condition
271 t(:,:,:) = t_ref
272 s(:,:,:) = s_ref
273
274 do j = g%jsc,g%jec
275 do i = g%isc,g%iec
276 zi = 0.
277 x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon
278 do k = 1,nz
279 ! nondimensional middle of layer
280 zmid = zi + 0.5 * dz(i,j,k) / g%max_depth
281
282 if (x > (1. - east_sponge_width)) then
283 !if (zmid >= 0.9 * sill_frac) &
284 s(i,j,k) = s_ref + s_dense
285 else
286 ! linear between bottom of mixed layer and bottom
287 if (zmid >= mld) &
288 s(i,j,k) = s_ref + s_range * (zmid - mld) / (1.0 - mld)
289 endif
290
291 zi = zi + dz(i,j,k) / g%max_depth
292 enddo
293 enddo
294 enddo
295
296 ! This call sets up the damping rates and interface heights in the sponges.
297 call initialize_ale_sponge(idamp, g, gv, param_file, acsp, dz, nz, data_h_is_z=.true.)
298
299 if ( associated(tv%T) ) call set_up_ale_sponge_field(t, g, gv, tv%T, acsp, 'temp', &
300 sp_long_name='temperature', sp_unit='degC s-1')
301 if ( associated(tv%S) ) call set_up_ale_sponge_field(s, g, gv, tv%S, acsp, 'salt', &
302 sp_long_name='salinity', sp_unit='g kg-1 s-1')
303 else
304 call mom_error(fatal, "dense_water_initialize_sponges: trying to use non ALE sponge")
305 endif
307
309
310!> \namespace dense_water_initialization
311!!
312!! This experiment consists of a shelf accumulating dense water, which spills
313!! over an upward slope and a sill, before flowing down a slope into an open
314!! ocean region. It's intended as a test of one of the motivating situations
315!! for the adaptive coordinate.
316!!
317!! The nondimensional widths of the 5 regions are controlled by the
318!! <code>DENSE_WATER_DOMAIN_PARAMS</code>, and the heights of the sill and shelf
319!! as a fraction of the total domain depth are controlled by
320!! <code>DENSE_WATER_SILL_DEPTH</code> and <code>DENSE_WATER_SHELF_DEPTH</code>.
321!!
322!! The density in the domain is governed by a linear equation of state, and
323!! is set up with a mixed layer of non-dimensional depth <code>DENSE_WATER_MLD</code>
324!! below which there is a linear stratification from <code>S_REF</code>, increasing
325!! by <code>S_RANGE</code> to the bottom.
326!!
327!! To force the experiment, there are sponges on the inflow and outflow of the
328!! domain. The inflow sponge has a salinity anomaly of
329!! <code>DENSE_WATER_EAST_SPONGE_SALT</code> through the entire depth. The outflow
330!! sponge simply restores to the initial condition. Both sponges have controllable
331!! widths and restoring timescales.