DOME2d_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 of the 2D DOME experiment with density water initialized on a coastal shelf.
7
8use mom_ale_sponge, only : ale_sponge_cs, set_up_ale_sponge_field, initialize_ale_sponge
10use mom_error_handler, only : mom_mesg, mom_error, fatal
11use mom_file_parser, only : get_param, log_version, param_file_type
12use mom_get_input, only : directories
13use mom_grid, only : ocean_grid_type
18use regrid_consts, only : coordinatemode, default_coordinate_mode
19use regrid_consts, only : regridding_layer, regridding_zstar
20use regrid_consts, only : regridding_rho, regridding_sigma
21
22implicit none ; private
23
24#include <MOM_memory.h>
25
26! Public functions
31
32! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
33! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
34! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
35! vary with the Boussinesq approximation, the Boussinesq variant is given first.
36
37character(len=40) :: mdl = "DOME2D_initialization" !< This module's name.
38
39contains
40
41!> Initialize topography with a shelf and slope in a 2D domain
42subroutine dome2d_initialize_topography( D, G, param_file, max_depth )
43 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
44 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
45 intent(out) :: d !< Ocean bottom depth [Z ~> m]
46 type(param_file_type), intent(in) :: param_file !< Parameter file structure
47 real, intent(in) :: max_depth !< Maximum ocean depth [Z ~> m]
48
49 ! Local variables
50 real :: bay_depth ! Depth of shelf, as fraction of basin depth [nondim]
51 real :: l1, l2 ! Fractional horizontal positions where the slope changes [nondim]
52 real :: x ! Fractional horizontal positions [nondim]
53 real :: dome2d_width_bay ! Width of shelf, as fraction of domain [nondim]
54 real :: dome2d_width_bottom ! Width of deep ocean basin, as fraction of domain [nondim]
55 real :: dome2d_depth_bay ! Depth of shelf, as fraction of basin depth [nondim]
56 integer :: i, j
57 ! This include declares and sets the variable "version".
58# include "version_variable.h"
59
60 call log_version(param_file, mdl, version, "")
61 call get_param(param_file, mdl, "DOME2D_SHELF_WIDTH", dome2d_width_bay, &
62 'Width of shelf, as fraction of domain, in 2d DOME configuration.', &
63 units='nondim', default=0.1)
64 call get_param(param_file, mdl, "DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
65 'Width of deep ocean basin, as fraction of domain, in 2d DOME configuration.', &
66 units='nondim', default=0.3)
67 call get_param(param_file, mdl, "DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
68 'Depth of shelf, as fraction of basin depth, in 2d DOME configuration.', &
69 units='nondim', default=0.2)
70
71 ! location where downslope starts
72 l1 = dome2d_width_bay
73
74 ! location where downslope reaches maximum depth
75 l2 = 1.0 - dome2d_width_bottom
76
77 bay_depth = dome2d_depth_bay
78
79 do j=g%jsc,g%jec ; do i=g%isc,g%iec
80
81 ! Compute normalized zonal coordinate
82 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
83
84 if ( x <= l1 ) then
85 d(i,j) = bay_depth * max_depth
86 elseif (( x > l1 ) .and. ( x < l2 )) then
87 d(i,j) = bay_depth * max_depth + (1.0-bay_depth) * max_depth * &
88 ( x - l1 ) / (l2 - l1)
89 else
90 d(i,j) = max_depth
91 endif
92
93 enddo ; enddo
94
96
97!> Initialize thicknesses according to coordinate mode
98subroutine dome2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read )
99 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
100 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
101 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
102 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
103 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
104 real, dimension(SZI_(G),SZJ_(G)), &
105 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
106 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
107 !! to parse for model parameter values.
108 logical, intent(in) :: just_read !< If true, this call will only read
109 !! parameters without changing h.
110
111 ! Local variables
112 real :: e0(szk_(gv)) ! The resting interface heights, in depth units [Z ~> m], usually
113 ! negative because it is positive upward.
114 real :: eta1d(szk_(gv)+1)! Interface height relative to the sea surface
115 ! positive upward, in depth units [Z ~> m]
116 real :: x ! Fractional horizontal positions [nondim]
117 real :: min_thickness ! Minimum layer thicknesses [Z ~> m]
118 real :: dome2d_width_bay ! Width of shelf, as fraction of domain [nondim]
119 real :: dome2d_width_bottom ! Width of deep ocean basin, as fraction of domain [nondim]
120 real :: dome2d_depth_bay ! Depth of shelf, as fraction of basin depth [nondim]
121 character(len=40) :: verticalcoordinate
122 integer :: i, j, k, is, ie, js, je, nz
123
124 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
125
126 if (.not.just_read) &
127 call mom_mesg("MOM_initialization.F90, DOME2d_initialize_thickness: setting thickness")
128
129 call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
130 default=1.e-3, units="m", do_not_log=.true., scale=us%m_to_Z)
131 call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
132 default=default_coordinate_mode, do_not_log=.true.)
133 call get_param(param_file, mdl, "DOME2D_SHELF_WIDTH", dome2d_width_bay, &
134 units="nondim", default=0.1, do_not_log=.true.)
135 call get_param(param_file, mdl, "DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
136 units="nondim", default=0.3, do_not_log=.true.)
137 call get_param(param_file, mdl, "DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
138 units="nondim", default=0.2, do_not_log=.true.)
139
140 if (just_read) return ! All run-time parameters have been read, so return.
141
142 ! WARNING: this routine specifies the interface heights so that the last layer
143 ! is vanished, even at maximum depth. In order to have a uniform
144 ! layer distribution, use this line of code within the loop:
145 ! e0(k) = -G%max_depth * real(k-1) / real(nz)
146 ! To obtain a thickness distribution where the last layer is
147 ! vanished and the other thicknesses uniformly distributed, use:
148 ! e0(k) = -G%max_depth * real(k-1) / real(nz-1)
149 do k=1,nz
150 e0(k) = -g%max_depth * real(k-1) / real(nz)
151 enddo
152
153 select case ( coordinatemode(verticalcoordinate) )
154
155 case ( regridding_layer, regridding_rho )
156
157 do j=js,je ; do i=is,ie
158 eta1d(nz+1) = -depth_tot(i,j)
159 do k=nz,1,-1
160 eta1d(k) = e0(k)
161 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
162 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
163 h(i,j,k) = gv%Angstrom_Z
164 else
165 h(i,j,k) = eta1d(k) - eta1d(k+1)
166 endif
167 enddo
168
169 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
170 if ( x <= dome2d_width_bay ) then
171 h(i,j,1:nz-1) = gv%Angstrom_Z
172 h(i,j,nz) = dome2d_depth_bay * g%max_depth - (nz-1) * gv%Angstrom_Z
173 endif
174
175 enddo ; enddo
176
177 ! case ( IC_RHO_C )
178 !
179 ! do j=js,je ; do i=is,ie
180 ! eta1D(nz+1) = -depth_tot(i,j)
181 ! do k=nz,1,-1
182 ! eta1D(k) = e0(k)
183 ! if (eta1D(k) < (eta1D(k+1) + min_thickness)) then
184 ! eta1D(k) = eta1D(k+1) + min_thickness
185 ! h(i,j,k) = min_thickness
186 ! else
187 ! h(i,j,k) = eta1D(k) - eta1D(k+1)
188 ! endif
189 ! enddo
190 !
191 ! x = G%geoLonT(i,j) / G%len_lon
192 ! if ( x <= dome2d_width_bay ) then
193 ! h(i,j,1:nz-1) = min_thickness
194 ! h(i,j,nz) = dome2d_depth_bay * G%max_depth - (nz-1) * min_thickness
195 ! endif
196 !
197 ! enddo ; enddo
198
199 case ( regridding_zstar )
200
201 do j=js,je ; do i=is,ie
202 eta1d(nz+1) = -depth_tot(i,j)
203 do k=nz,1,-1
204 eta1d(k) = e0(k)
205 if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
206 eta1d(k) = eta1d(k+1) + min_thickness
207 h(i,j,k) = min_thickness
208 else
209 h(i,j,k) = eta1d(k) - eta1d(k+1)
210 endif
211 enddo
212 enddo ; enddo
213
214 case ( regridding_sigma )
215 do j=js,je ; do i=is,ie
216 h(i,j,:) = depth_tot(i,j) / nz
217 enddo ; enddo
218
219 case default
220 call mom_error(fatal,"dome2d_initialize: "// &
221 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
222
223 end select
224
225end subroutine dome2d_initialize_thickness
226
227
228!> Initialize temperature and salinity in the 2d DOME configuration
229subroutine dome2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, just_read)
230 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
231 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
232 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t !< Potential temperature [C ~> degC]
233 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s !< Salinity [S ~> ppt]
234 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m]
235 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
236 type(param_file_type), intent(in) :: param_file !< Parameter file structure
237 logical, intent(in) :: just_read !< If true, this call will
238 !! only read parameters without changing T & S.
239
240 real :: x ! Fractional horizontal positions [nondim]
241 real :: delta_s ! Change in salinity between layers [S ~> ppt]
242 real :: s_ref, t_ref ! Reference salinity [S ~> ppt] and temperature [C ~> degC] within surface layer
243 real :: s_range, t_range ! Range of salinities [S ~> ppt] and temperatures [C ~> degC] over the vertical
244 real :: s_surf ! Initial surface salinity [S ~> ppt]
245 real :: t_bay ! Temperature in the inflow embayment [C ~> degC]
246 real :: xi0, xi1 ! Fractional vertical positions [nondim]
247 real :: dome2d_width_bay ! Width of shelf, as fraction of domain [nondim]
248 real :: dome2d_width_bottom ! Width of deep ocean basin, as fraction of domain [nondim]
249 real :: dome2d_depth_bay ! Depth of shelf, as fraction of basin depth [nondim]
250 character(len=40) :: verticalcoordinate
251 integer :: index_bay_z
252 integer :: i, j, k, is, ie, js, je, nz
253
254 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
255
256 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
257 default=default_coordinate_mode, do_not_log=.true.)
258 call get_param(param_file, mdl, "DOME2D_SHELF_WIDTH", dome2d_width_bay, &
259 units="nondim", default=0.1, do_not_log=.true.)
260 call get_param(param_file, mdl, "DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
261 units="nondim", default=0.3, do_not_log=.true.)
262 call get_param(param_file, mdl, "DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
263 units="nondim", default=0.2, do_not_log=.true.)
264 call get_param(param_file, mdl, "S_REF", s_ref, 'Reference salinity', &
265 units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=just_read)
266 call get_param(param_file, mdl, "T_REF", t_ref, 'Reference temperature', &
267 units='degC', scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
268 call get_param(param_file, mdl, "S_RANGE", s_range,' Initial salinity range', &
269 units="ppt", default=2.0, scale=us%ppt_to_S, do_not_log=just_read)
270 call get_param(param_file, mdl, "T_RANGE", t_range, 'Initial temperature range', &
271 units='degC', default=0.0, scale=us%degC_to_C, do_not_log=just_read)
272 call get_param(param_file, mdl, "INITIAL_SSS", s_surf, "Initial surface salinity", &
273 units="ppt", default=34.0, scale=us%ppt_to_S, do_not_log=just_read)
274 call get_param(param_file, mdl, "DOME2D_T_BAY", t_bay, &
275 "Temperature in the inflow embayment in the DOME2d test case", &
276 units="degC", default=1.0, scale=us%degC_to_C, do_not_log=just_read)
277
278 if (just_read) return ! All run-time parameters have been read, so return.
279
280 t(:,:,:) = 0.0
281 s(:,:,:) = 0.0
282
283 ! Linear salinity profile
284
285 select case ( coordinatemode(verticalcoordinate) )
286
287 case ( regridding_zstar, regridding_sigma )
288
289 do j=js,je ; do i=is,ie
290 xi0 = 0.0
291 do k = 1,nz
292 xi1 = xi0 + h(i,j,k) / g%max_depth
293 s(i,j,k) = s_surf + 0.5 * s_range * (xi0 + xi1)
294 xi0 = xi1
295 enddo
296 enddo ; enddo
297
298 case ( regridding_rho )
299
300 do j=js,je ; do i=is,ie
301 xi0 = 0.0
302 do k = 1,nz
303 xi1 = xi0 + h(i,j,k) / g%max_depth
304 s(i,j,k) = s_surf + 0.5 * s_range * (xi0 + xi1)
305 xi0 = xi1
306 enddo
307 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
308 if ( x <= dome2d_width_bay ) then
309 s(i,j,nz) = s_surf + s_range
310 endif
311 enddo ; enddo
312
313 case ( regridding_layer )
314
315 delta_s = s_range / ( gv%ke - 1.0 )
316 s(:,:,1) = s_ref
317 do k = 2,gv%ke
318 s(:,:,k) = s(:,:,k-1) + delta_s
319 enddo
320
321 case default
322 call mom_error(fatal,"dome2d_initialize: "// &
323 "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
324
325 end select
326
327 ! Modify salinity and temperature when z coordinates are used
328 if ( coordinatemode(verticalcoordinate) == regridding_zstar ) then
329 index_bay_z = nint( dome2d_depth_bay * gv%ke )
330 do j = g%jsc,g%jec ; do i = g%isc,g%iec
331 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
332 if ( x <= dome2d_width_bay ) then
333 s(i,j,1:index_bay_z) = s_ref + s_range ! Use for z coordinates
334 t(i,j,1:index_bay_z) = t_bay ! Use for z coordinates
335 endif
336 enddo ; enddo ! i and j loops
337 endif ! Z initial conditions
338
339 ! Modify salinity and temperature when sigma coordinates are used
340 if ( coordinatemode(verticalcoordinate) == regridding_sigma ) then
341 do i = g%isc,g%iec ; do j = g%jsc,g%jec
342 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
343 if ( x <= dome2d_width_bay ) then
344 s(i,j,1:gv%ke) = s_ref + s_range ! Use for sigma coordinates
345 t(i,j,1:gv%ke) = t_bay ! Use for sigma coordinates
346 endif
347 enddo ; enddo
348 endif
349
350 ! Modify temperature when rho coordinates are used
351 if (( coordinatemode(verticalcoordinate) == regridding_rho ) .or. &
352 ( coordinatemode(verticalcoordinate) == regridding_layer )) then
353 do i = g%isc,g%iec ; do j = g%jsc,g%jec
354 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
355 if ( x <= dome2d_width_bay ) then
356 t(i,j,gv%ke) = t_bay
357 endif
358 enddo ; enddo
359 endif
360
362
363!> Set up sponges in 2d DOME configuration
364subroutine dome2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp)
365 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
366 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
367 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
368 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
369 real, dimension(SZI_(G),SZJ_(G)), &
370 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
371 type(param_file_type), intent(in) :: param_file !< Parameter file structure
372 logical, intent(in) :: use_ale !< If true, indicates model is in ALE mode
373 type(sponge_cs), pointer :: csp !< Layer-mode sponge structure
374 type(ale_sponge_cs), pointer :: acsp !< ALE-mode sponge structure
375 ! Local variables
376 real :: t(szi_(g),szj_(g),szk_(gv)) ! A temporary array for temp [C ~> degC]
377 real :: s(szi_(g),szj_(g),szk_(gv)) ! A temporary array for salt [S ~> ppt]
378 real :: dz(szi_(g),szj_(g),szk_(gv)) ! A temporary array for thickness in height units [Z ~> m]
379 real :: eta(szi_(g),szj_(g),szk_(gv)+1) ! A temporary array for interface heights [Z ~> m]
380 real :: idamp(szi_(g),szj_(g)) ! The sponge damping rate [T-1 ~> s-1]
381 real :: s_ref ! Reference salinity within the surface layer [S ~> ppt]
382 real :: t_ref ! Reference temperature within the surface layer [C ~> degC]
383 real :: s_range ! Range of salinities in the vertical [S ~> ppt]
384 real :: t_range ! Range of temperatures in the vertical [C ~> degC]
385 real :: s_range_sponge ! Range of salinities in the vertical in the east sponge [S ~> ppt]
386 real :: s_surf ! Initial surface salinity [S ~> ppt]
387 real :: e0(szk_(gv)+1) ! The resting interface heights [Z ~> m],
388 ! usually negative because it is positive upward.
389 real :: eta1d(szk_(gv)+1) ! Interface height relative to the sea surface
390 ! positive upward [Z ~> m].
391 real :: d_eta(szk_(gv)) ! The layer thickness in a column [Z ~> m].
392 real :: dome2d_width_bay ! Width of shelf, as fraction of domain [nondim]
393 real :: dome2d_width_bottom ! Width of deep ocean basin, as fraction of domain [nondim]
394 real :: dome2d_depth_bay ! Depth of shelf, as fraction of basin depth [nondim]
395 real :: dome2d_west_sponge_time_scale, dome2d_east_sponge_time_scale ! Sponge timescales [T ~> s]
396 real :: dome2d_west_sponge_width ! The fraction of the domain in which the western sponge for
397 ! restoring T/S is active [nondim]
398 real :: dome2d_east_sponge_width ! The fraction of the domain in which the eastern sponge for
399 ! restoring T/S is active [nondim]
400 real :: dummy1, x ! Nondimensional local variables indicating horizontal positions [nondim]
401 real :: z ! Vertical positions [Z ~> m]
402 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
403
404 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
405 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
406
407 call get_param(param_file, mdl, "DOME2D_WEST_SPONGE_TIME_SCALE", dome2d_west_sponge_time_scale, &
408 'The time-scale on the west edge of the domain for restoring T/S '//&
409 'in the sponge. If zero, the western sponge is disabled', &
410 units='s', default=0., scale=us%s_to_T)
411 call get_param(param_file, mdl, "DOME2D_EAST_SPONGE_TIME_SCALE", dome2d_east_sponge_time_scale, &
412 'The time-scale on the east edge of the domain for restoring T/S '//&
413 'in the sponge. If zero, the eastern sponge is disabled', &
414 units='s', default=0., scale=us%s_to_T)
415 call get_param(param_file, mdl, "DOME2D_WEST_SPONGE_WIDTH", dome2d_west_sponge_width, &
416 'The fraction of the domain in which the western sponge for restoring T/S '//&
417 'is active.', &
418 units='nondim', default=0.1)
419 call get_param(param_file, mdl, "DOME2D_EAST_SPONGE_WIDTH", dome2d_east_sponge_width, &
420 'The fraction of the domain in which the eastern sponge for restoring T/S '//&
421 'is active.', &
422 units='nondim', default=0.1)
423
424 ! Return if sponges are not in use
425 if (dome2d_west_sponge_time_scale <= 0. .and. dome2d_east_sponge_time_scale <= 0.) return
426
427 if (associated(csp)) call mom_error(fatal, &
428 "DOME2d_initialize_sponges called with an associated control structure.")
429 if (associated(acsp)) call mom_error(fatal, &
430 "DOME2d_initialize_sponges called with an associated ALE-sponge control structure.")
431
432 call get_param(param_file, mdl, "DOME2D_SHELF_WIDTH", dome2d_width_bay, &
433 units="nondim", default=0.1, do_not_log=.true.)
434 call get_param(param_file, mdl, "DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
435 units="nondim", default=0.3, do_not_log=.true.)
436 call get_param(param_file, mdl, "DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
437 units="nondim", default=0.2, do_not_log=.true.)
438 call get_param(param_file, mdl, "S_REF", s_ref, units="ppt", default=35.0, scale=us%ppt_to_S)
439 call get_param(param_file, mdl, "T_REF", t_ref, units="degC", scale=us%degC_to_C, fail_if_missing=.false.)
440 call get_param(param_file, mdl, "S_RANGE", s_range, units="ppt", default=2.0, scale=us%ppt_to_S)
441 call get_param(param_file, mdl, "T_RANGE", t_range, units="degC", default=0.0, scale=us%degC_to_C)
442 call get_param(param_file, mdl, "INITIAL_SSS", s_surf, "Initial surface salinity", &
443 units="ppt", default=34.0, scale=us%ppt_to_S, do_not_log=.true.)
444 call get_param(param_file, mdl, "DOME2D_EAST_SPONGE_S_RANGE", s_range_sponge, &
445 "Range of salinities in the eastern sponge region in the DOME2D configuration", &
446 units="ppt", default=1.0, scale=us%ppt_to_S)
447
448 ! Set the sponge damping rate as a function of position
449 idamp(:,:) = 0.0
450 do j=js,je ; do i=is,ie
451 if (g%mask2dT(i,j) > 0.) then ! Only set damping rate for wet points
452 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon ! Non-dimensional position within domain (0,1)
453 if ( dome2d_west_sponge_time_scale > 0. .and. x < dome2d_west_sponge_width ) then
454 ! Within half the shelf width from the left edge
455 dummy1 = 1. - x / dome2d_west_sponge_width
456 idamp(i,j) = 1./dome2d_west_sponge_time_scale * max(0., min(1., dummy1))
457 elseif ( dome2d_east_sponge_time_scale > 0. .and. x > ( 1. - dome2d_east_sponge_width ) ) then
458 ! Within a quarter of the basin width from the right
459 dummy1 = 1. - ( 1. - x ) / dome2d_east_sponge_width
460 idamp(i,j) = 1./dome2d_east_sponge_time_scale * max(0., min(1., dummy1))
461 else
462 idamp(i,j) = 0.
463 endif
464 else
465 idamp(i,j) = 0.
466 endif
467 enddo ; enddo
468
469 if (use_ale) then
470
471 ! Construct a grid (somewhat arbitrarily) to describe the sponge T/S on
472 do k=1,nz
473 e0(k) = -g%max_depth * ( real(k-1) / real(nz) )
474 enddo
475 e0(nz+1) = -g%max_depth
476 do j=js,je ; do i=is,ie
477 eta1d(nz+1) = -depth_tot(i,j)
478 do k=nz,1,-1
479 eta1d(k) = e0(k)
480 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
481 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
482 dz(i,j,k) = gv%Angstrom_Z
483 else
484 dz(i,j,k) = eta1d(k) - eta1d(k+1)
485 endif
486 enddo
487 enddo ; enddo
488
489 ! Construct temperature and salinity on the arbitrary grid
490 t(:,:,:) = 0.0 ; s(:,:,:) = 0.0
491 do j=js,je ; do i=is,ie
492 z = -depth_tot(i,j)
493 do k = nz,1,-1
494 z = z + 0.5 * dz(i,j,k) ! Position of the center of layer k
495 ! Use salinity stratification in the eastern sponge.
496 s(i,j,k) = s_surf - s_range_sponge * (z / g%max_depth)
497 ! Use a constant salinity in the western sponge.
498 if ( ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon < dome2d_west_sponge_width ) &
499 s(i,j,k) = s_ref + s_range
500 z = z + 0.5 * dz(i,j,k) ! Position of the interface k
501 enddo
502 enddo ; enddo
503
504 ! Store damping rates and the grid on which the T/S sponge data will reside
505 call initialize_ale_sponge(idamp, g, gv, param_file, acsp, dz, nz, data_h_is_z=.true.)
506
507 if ( associated(tv%T) ) call set_up_ale_sponge_field(t, g, gv, tv%T, acsp, 'temp', &
508 sp_long_name='temperature', sp_unit='degC s-1')
509 if ( associated(tv%S) ) call set_up_ale_sponge_field(s, g, gv, tv%S, acsp, 'salt', &
510 sp_long_name='salinity', sp_unit='g kg-1 s-1')
511
512 else
513
514 ! Construct interface heights to restore toward
515 do j=js,je ; do i=is,ie
516 eta1d(nz+1) = -depth_tot(i,j)
517 do k=nz,1,-1
518 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
519 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
520 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
521 d_eta(k) = gv%Angstrom_Z
522 else
523 d_eta(k) = (eta1d(k) - eta1d(k+1))
524 endif
525 enddo
526
527 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
528 if ( x <= dome2d_width_bay ) then
529 do k=1,nz-1 ; d_eta(k) = gv%Angstrom_Z ; enddo
530 d_eta(nz) = dome2d_depth_bay * g%max_depth - (nz-1) * gv%Angstrom_Z
531 endif
532
533 eta(i,j,nz+1) = -depth_tot(i,j)
534 do k=nz,1,-1
535 eta(i,j,k) = eta(i,j,k+1) + d_eta(k)
536 enddo
537 enddo ; enddo
538 call initialize_sponge(idamp, eta, g, param_file, csp, gv)
539
540 endif
541
542end subroutine dome2d_initialize_sponges
543
544end module dome2d_initialization