dumbbell_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 dumbbell test case.
7
8use mom_domains, only : sum_across_pes
10use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
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
14use mom_interface_heights, only : thickness_to_dz
16use mom_tracer_registry, only : tracer_registry_type
20use regrid_consts, only : coordinatemode, default_coordinate_mode
21use regrid_consts, only : regridding_layer, regridding_zstar
22use regrid_consts, only : regridding_rho, regridding_sigma, regridding_hycom1
23use mom_ale_sponge, only : ale_sponge_cs, set_up_ale_sponge_field, initialize_ale_sponge
24
25implicit none ; private
26
27#include <MOM_memory.h>
28
29character(len=40) :: mdl = "dumbbell_initialization" !< This module's name.
30
35
36! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
37! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
38! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
39! vary with the Boussinesq approximation, the Boussinesq variant is given first.
40
41contains
42
43!> Initialization of topography.
44subroutine dumbbell_initialize_topography( D, G, param_file, max_depth )
45 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
46 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
47 intent(out) :: d !< Ocean bottom depth [Z ~> m]
48 type(param_file_type), intent(in) :: param_file !< Parameter file structure
49 real, intent(in) :: max_depth !< Maximum ocean depth [Z ~> m]
50
51 ! Local variables
52 real :: x, y ! Fractional x- and y- positions [nondim]
53 real :: dblen ! Lateral length scale for dumbbell [km] or [m]
54 real :: dbfrac ! Meridional fraction for narrow part of dumbbell [nondim]
55 logical :: dbrotate ! If true, rotate this configuration
56 integer :: i, j
57
58 call get_param(param_file, mdl, "DUMBBELL_LEN", dblen, &
59 'Lateral Length scale for dumbbell.', &
60 units=g%x_ax_unit_short, default=600., do_not_log=.false.)
61 call get_param(param_file, mdl, "DUMBBELL_FRACTION", dbfrac, &
62 'Meridional fraction for narrow part of dumbbell.', &
63 units='nondim', default=0.5, do_not_log=.false.)
64 call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
65 'Logical for rotation of dumbbell domain.', &
66 default=.false., do_not_log=.false.)
67
68 if (g%x_axis_units(1:1) == 'm') then
69 dblen = dblen*1.e3
70 endif
71
72 if (dbrotate) then
73 do j=g%jsc,g%jec ; do i=g%isc,g%iec
74 ! Compute normalized zonal coordinates (x,y=0 at center of domain)
75 x = ( g%geoLonT(i,j) ) / g%len_lon
76 y = ( g%geoLatT(i,j) ) / dblen
77 d(i,j) = g%max_depth
78 if ((y>=-0.25 .and. y<=0.25) .and. (x <= -0.5*dbfrac .or. x >= 0.5*dbfrac)) then
79 d(i,j) = 0.0
80 endif
81 enddo ; enddo
82 else
83 do j=g%jsc,g%jec ; do i=g%isc,g%iec
84 ! Compute normalized zonal coordinates (x,y=0 at center of domain)
85 x = ( g%geoLonT(i,j) ) / dblen
86 y = ( g%geoLatT(i,j) ) / g%len_lat
87 d(i,j) = g%max_depth
88 if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac)) then
89 d(i,j) = 0.0
90 endif
91 enddo ; enddo
92 endif
93
95
96!> Initializes the layer thicknesses to be uniform in the dumbbell test case
97subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read)
98 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
99 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
100 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
101 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
102 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
103 real, dimension(SZI_(G),SZJ_(G)), &
104 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
105 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
106 !! to parse for model parameter values.
107 logical, intent(in) :: just_read !< If true, this call will only read
108 !! parameters without changing h.
109
110 real :: e0(szk_(gv)+1) ! The resting interface heights [Z ~> m], usually
111 ! negative because it is positive upward.
112 real :: eta1d(szk_(gv)+1) ! Interface height relative to the sea surface
113 ! positive upward [Z ~> m].
114 real :: min_thickness ! The minimum layer thicknesses [Z ~> m].
115 real :: s_ref ! A default value for salinities [S ~> ppt].
116 real :: s_surf ! The surface salinity [S ~> ppt]
117 real :: s_range ! The range of salinities in this test case [S ~> ppt]
118 real :: s_light, s_dense ! The lightest and densest salinities in the sponges [S ~> ppt].
119 real :: eta_ic_quanta ! The granularity of quantization of initial interface heights [Z-1 ~> m-1].
120 logical :: dbrotate ! If true, rotate the domain.
121 logical :: use_ale ! True if ALE is being used, False if in layered mode
122
123 ! This include declares and sets the variable "version".
124# include "version_variable.h"
125 character(len=20) :: verticalcoordinate
126 integer :: i, j, k, is, ie, js, je, nz
127
128 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
129
130 if (.not.just_read) &
131 call mom_mesg("dumbbell_initialization.F90, dumbbell_initialize_thickness: setting thickness")
132
133 if (.not.just_read) call log_version(param_file, mdl, version, "")
134 call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
135 'Minimum thickness for layer', &
136 units='m', default=1.0e-3, scale=us%m_to_Z, do_not_log=just_read)
137 call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
138 default=default_coordinate_mode, do_not_log=just_read)
139 call get_param(param_file, mdl, "USE_REGRIDDING", use_ale, default=.false., do_not_log=.true.)
140 if (.not. use_ale) verticalcoordinate = "LAYER"
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+1
150 ! e0(k) = -G%max_depth * real(k-1) / real(nz)
151 !enddo
152
153 select case ( coordinatemode(verticalcoordinate) )
154 case ( regridding_layer) ! Initial thicknesses for isopycnal coordinates
155 call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
156 'Logical for rotation of dumbbell domain.', &
157 default=.false., do_not_log=just_read)
158 do j=js,je
159 do i=is,ie
160 ! Work relative to the center of the domain, where geoLonT and geoLatT are both 0.
161 eta1d(1) = 0.0
162 eta1d(nz+1) = -depth_tot(i,j)
163 if (((.not.dbrotate) .and. (g%geoLonT(i,j)<0.0)) .or. (dbrotate .and. (g%geoLatT(i,j)<0.0))) then
164 do k=nz,2, -1
165 eta1d(k) = eta1d(k+1) + min_thickness
166 enddo
167 else
168 do k=2,nz
169 eta1d(k) = eta1d(k-1) - min_thickness
170 enddo
171 endif
172 do k=1,nz
173 h(i,j,k) = eta1d(k) - eta1d(k+1)
174 enddo
175 enddo
176 enddo
177
178 case ( regridding_rho, regridding_hycom1) ! Initial thicknesses for isopycnal coordinates
179 call get_param(param_file, mdl, "INITIAL_SSS", s_surf, &
180 units="ppt", default=34., scale=us%ppt_to_S, do_not_log=.true.)
181 call get_param(param_file, mdl, "INITIAL_S_RANGE", s_range, &
182 units="ppt", default=2., scale=us%ppt_to_S, do_not_log=.true.)
183 call get_param(param_file, mdl, "S_REF", s_ref, &
184 units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=.true.)
185 call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", s_light, &
186 units="ppt", default=us%S_to_ppt*s_ref, scale=us%ppt_to_S, do_not_log=.true.)
187 call get_param(param_file, mdl, "TS_RANGE_S_DENSE", s_dense, &
188 units="ppt", default=us%S_to_ppt*s_ref, scale=us%ppt_to_S, do_not_log=.true.)
189 call get_param(param_file, mdl, "INTERFACE_IC_QUANTA", eta_ic_quanta, &
190 "The granularity of initial interface height values "//&
191 "per meter, to avoid sensivity to order-of-arithmetic changes.", &
192 default=2048.0, units="m-1", scale=us%Z_to_m, do_not_log=just_read)
193 if (just_read) return ! All run-time parameters have been read, so return.
194
195 do k=1,nz+1
196 ! Salinity of layer k is S_light + (k-1)/(nz-1) * (S_dense - S_light)
197 ! Salinity of interface K is S_light + (K-3/2)/(nz-1) * (S_dense - S_light)
198 ! Salinity at depth z should be S(z) = S_surf - S_range * z/max_depth
199 ! Equating: S_surf - S_range * z/max_depth = S_light + (K-3/2)/(nz-1) * (S_dense - S_light)
200 ! Equating: - S_range * z/max_depth = S_light - S_surf + (K-3/2)/(nz-1) * (S_dense - S_light)
201 ! Equating: z/max_depth = - ( S_light - S_surf + (K-3/2)/(nz-1) * (S_dense - S_light) ) / S_range
202 e0(k) = - g%max_depth * ( ( s_light - s_surf ) + ( s_dense - s_light ) * &
203 ( (real(k)-1.5) / real(nz-1) ) ) / s_range
204 ! Force round numbers ... the above expression has irrational factors ...
205 if (eta_ic_quanta > 0.0) &
206 e0(k) = nint(eta_ic_quanta*e0(k)) / eta_ic_quanta
207 e0(k) = min(real(1-k)*gv%Angstrom_Z, e0(k)) ! Bound by surface
208 e0(k) = max(-g%max_depth, e0(k)) ! Bound by bottom
209 enddo
210 do j=js,je ; do i=is,ie
211 eta1d(nz+1) = -depth_tot(i,j)
212 do k=nz,1,-1
213 eta1d(k) = e0(k)
214 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
215 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
216 h(i,j,k) = gv%Angstrom_Z
217 else
218 h(i,j,k) = eta1d(k) - eta1d(k+1)
219 endif
220 enddo
221 enddo ; enddo
222
223 case ( regridding_zstar ) ! Initial thicknesses for z coordinates
224 if (just_read) return ! All run-time parameters have been read, so return.
225 do j=js,je ; do i=is,ie
226 eta1d(nz+1) = -depth_tot(i,j)
227 do k=nz,1,-1
228 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
229 if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
230 eta1d(k) = eta1d(k+1) + min_thickness
231 h(i,j,k) = min_thickness
232 else
233 h(i,j,k) = eta1d(k) - eta1d(k+1)
234 endif
235 enddo
236 enddo ; enddo
237
238 case ( regridding_sigma ) ! Initial thicknesses for sigma coordinates
239 if (just_read) return ! All run-time parameters have been read, so return.
240 do j=js,je ; do i=is,ie
241 h(i,j,:) = depth_tot(i,j) / real(nz)
242 enddo ; enddo
243
244end select
245
247
248!> Initial values for temperature and salinity for the dumbbell test case
249subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, just_read)
250 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
251 type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
252 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t !< Potential temperature [C ~> degC]
253 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s !< Salinity [S ~> ppt]
254 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m]
255 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
256 type(param_file_type), intent(in) :: param_file !< Parameter file structure
257 logical, intent(in) :: just_read !< If true, this call will
258 !! only read parameters without changing h.
259
260 ! Local variables
261 integer :: i, j, k, is, ie, js, je, nz
262 real :: s_surf ! The surface salinity [S ~> ppt]
263 real :: s_range ! The range of salinities in this test case [S ~> ppt]
264 real :: t_surf ! The surface temperature [C ~> degC]
265 real :: x ! The fractional position in the domain [nondim]
266 real :: dblen ! The size of the dumbbell test case [km] or [m]
267 logical :: dbrotate ! If true, rotate the domain.
268 logical :: use_ale ! If false, use layer mode.
269 character(len=20) :: verticalcoordinate, density_profile
270
271 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
272
273 ! layer mode
274 call get_param(param_file, mdl, "USE_REGRIDDING", use_ale, default=.false., do_not_log=.true.)
275 if (.not. use_ale) call mom_error(fatal, "dumbbell_initialize_temperature_salinity: "//&
276 "Please use 'fit' for 'TS_CONFIG' in the LAYER mode.")
277
278 call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
279 default=default_coordinate_mode, do_not_log=just_read)
280 call get_param(param_file, mdl, "INITIAL_DENSITY_PROFILE", density_profile, &
281 'Initial profile shape. Valid values are "linear", "parabolic" '// &
282 'and "exponential".', default='linear', do_not_log=just_read)
283 call get_param(param_file, mdl, "DUMBBELL_T_SURF", t_surf, &
284 'Initial surface temperature in the DUMBBELL configuration', &
285 units='degC', default=20., scale=us%degC_to_C, do_not_log=just_read)
286 call get_param(param_file, mdl, "DUMBBELL_SREF", s_surf, &
287 'DUMBBELL REFERENCE SALINITY', &
288 units="ppt", default=34., scale=us%ppt_to_S, do_not_log=just_read)
289 call get_param(param_file, mdl, "DUMBBELL_S_RANGE", s_range, &
290 'DUMBBELL salinity range (right-left)', &
291 units="ppt", default=2., scale=us%ppt_to_S, do_not_log=just_read)
292 call get_param(param_file, mdl, "DUMBBELL_LEN", dblen, &
293 'Lateral Length scale for dumbbell ', &
294 units=g%x_ax_unit_short, default=600., do_not_log=just_read)
295 call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
296 'Logical for rotation of dumbbell domain.', &
297 default=.false., do_not_log=just_read)
298
299 if (g%x_axis_units(1:1) == 'm') then
300 dblen = dblen*1.e3
301 endif
302
303 do j=g%jsc,g%jec
304 do i=g%isc,g%iec
305 ! Compute normalized zonal coordinates (x,y=0 at center of domain)
306 if (dbrotate) then
307 ! This is really y in the rotated case
308 x = ( g%geoLatT(i,j) ) / dblen
309 else
310 x = ( g%geoLonT(i,j) ) / dblen
311 endif
312 do k=1,nz
313 t(i,j,k) = t_surf
314 enddo
315 if (x>=0. ) then
316 do k=1,nz
317 s(i,j,k) = s_surf + 0.5*s_range
318 enddo
319 endif
320 if (x<0. ) then
321 do k=1,nz
322 s(i,j,k) = s_surf - 0.5*s_range
323 enddo
324 endif
325
326 enddo
327 enddo
328
330
331!> Initialize the restoring sponges for the dumbbell test case
332subroutine dumbbell_initialize_sponges(G, GV, US, tv, h_in, depth_tot, param_file, use_ALE, CSp, ACSp)
333 type(ocean_grid_type), intent(in) :: g !< Horizontal grid control structure
334 type(verticalgrid_type), intent(in) :: gv !< Vertical grid control structure
335 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
336 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
337 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2]
338 real, dimension(SZI_(G),SZJ_(G)), &
339 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
340 type(param_file_type), intent(in) :: param_file !< Parameter file structure
341 logical, intent(in) :: use_ale !< ALE flag
342 type(sponge_cs), pointer :: csp !< Layered sponge control structure pointer
343 type(ale_sponge_cs), pointer :: acsp !< ALE sponge control structure pointer
344
345 real :: sponge_time_scale ! The damping time scale [T ~> s]
346
347 real, dimension(SZI_(G),SZJ_(G)) :: idamp ! inverse damping timescale [T-1 ~> s-1]
348 real :: dz(szi_(g),szj_(g),szk_(gv)) ! Sponge thicknesses in height units [Z ~> m]
349 real :: s(szi_(g),szj_(g),szk_(gv)) ! Sponge salinities [S ~> ppt]
350 real :: t(szi_(g),szj_(g),szk_(gv)) ! Sponge tempertures [C ~> degC], used only to convert thicknesses
351 ! in non-Boussinesq mode
352 real, dimension(SZK_(GV)+1) :: eta1d ! Interface positions for ALE sponge [Z ~> m]
353 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta ! A temporary array for interface heights [Z ~> m].
354
355 integer :: i, j, k, nz
356 real :: x ! The fractional position in the domain [nondim]
357 real :: dblen ! The size of the dumbbell test case [km] or [m]
358 real :: min_thickness ! The minimum layer thickness [Z ~> m]
359 real :: s_ref, s_range ! A reference salinity and the range of salinities in this test case [S ~> ppt]
360 real :: t_surf ! The surface temperature [C ~> degC]
361 logical :: dbrotate ! If true, rotate the domain.
362
363 call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, &
364 'Lateral Length scale for dumbbell ', &
365 units='km', default=600., do_not_log=.true.)
366 call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
367 'Logical for rotation of dumbbell domain.', &
368 default=.false., do_not_log=.true.)
369
370 if (g%x_axis_units(1:1) == 'm') then
371 dblen = dblen*1.e3
372 endif
373
374 nz = gv%ke
375
376 call get_param(param_file, mdl, "DUMBBELL_SPONGE_TIME_SCALE", sponge_time_scale, &
377 "The time scale in the reservoir for restoring. If zero, the sponge is disabled.", &
378 units="s", default=0., scale=us%s_to_T)
379 call get_param(param_file, mdl, "DUMBBELL_T_SURF", t_surf, &
380 'Initial surface temperature in the DUMBBELL configuration', &
381 units='degC', default=20., scale=us%degC_to_C, do_not_log=.true.)
382 call get_param(param_file, mdl, "DUMBBELL_SREF", s_ref, &
383 'DUMBBELL REFERENCE SALINITY', &
384 units="ppt", default=34., scale=us%ppt_to_S, do_not_log=.true.)
385 call get_param(param_file, mdl, "DUMBBELL_S_RANGE", s_range, &
386 'DUMBBELL salinity range (right-left)', &
387 units="ppt", default=2., scale=us%ppt_to_S, do_not_log=.true.)
388 call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
389 'Minimum thickness for layer', &
390 units='m', default=1.0e-3, scale=us%m_to_Z, do_not_log=.true.)
391
392 ! no active sponges
393 if (sponge_time_scale <= 0.) return
394
395 ! everywhere is initially unsponged
396 idamp(:,:) = 0.0
397
398 do j = g%jsc, g%jec
399 do i = g%isc,g%iec
400 if (g%mask2dT(i,j) > 0.) then
401 ! nondimensional x position
402 if (dbrotate) then
403 ! This is really y in the rotated case
404 x = ( g%geoLatT(i,j) ) / dblen
405 else
406 x = ( g%geoLonT(i,j) ) / dblen
407 endif
408 if (x > 0.25 .or. x < -0.25) then
409 ! scale restoring by depth into sponge
410 idamp(i,j) = 1. / sponge_time_scale
411 endif
412 endif
413 enddo
414 enddo
415
416 if (use_ale) then
417 ! construct a uniform grid for the sponge
418 do j=g%jsc,g%jec ; do i=g%isc,g%iec
419 eta1d(nz+1) = depth_tot(i,j)
420 do k=nz,1,-1
421 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
422 if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
423 eta1d(k) = eta1d(k+1) + min_thickness
424 dz(i,j,k) = min_thickness
425 else
426 dz(i,j,k) = eta1d(k) - eta1d(k+1)
427 endif
428 enddo
429 enddo ; enddo
430
431 ! construct temperature and salinity for the sponge
432 ! start with initial condition
433 s(:,:,:) = 0.0
434 t(:,:,:) = t_surf
435
436 do j=g%jsc,g%jec ; do i=g%isc,g%iec
437 ! Compute normalized zonal coordinates (x,y=0 at center of domain)
438 if (dbrotate) then
439 ! This is really y in the rotated case
440 x = ( g%geoLatT(i,j) ) / dblen
441 else
442 x = ( g%geoLonT(i,j) ) / dblen
443 endif
444 if (x>=0.25 ) then
445 do k=1,nz
446 s(i,j,k) = s_ref + 0.5*s_range
447 enddo
448 endif
449 if (x<=-0.25 ) then
450 do k=1,nz
451 s(i,j,k) = s_ref - 0.5*s_range
452 enddo
453 endif
454 enddo ; enddo
455
456 ! Store damping rates and the grid on which the T/S sponge data will reside
457 call initialize_ale_sponge(idamp, g, gv, param_file, acsp, dz, nz, data_h_is_z=.true.)
458
459 if (associated(tv%S)) call set_up_ale_sponge_field(s, g, gv, tv%S, acsp, 'salt', &
460 sp_long_name='salinity', sp_unit='g kg-1 s-1')
461 else
462 ! Convert thicknesses from thickness units to height units
463 call thickness_to_dz(h_in, tv, dz, g, gv, us)
464
465 do j=g%jsc,g%jec ; do i=g%isc,g%iec
466 eta(i,j,1) = 0.0
467 do k=2,nz
468 eta(i,j,k) = eta(i,j,k-1) - dz(i,j,k-1)
469 enddo
470 eta(i,j,nz+1) = -depth_tot(i,j)
471 do k=1,nz
472 s(i,j,k)= tv%S(i,j,k)
473 enddo
474 enddo ; enddo
475
476 ! This call sets up the damping rates and interface heights.
477 ! This sets the inverse damping timescale fields in the sponges. !
478 call initialize_sponge(idamp, eta, g, param_file, csp, gv)
479
480 ! The remaining calls to set_up_sponge_field can be in any order. !
481 if ( associated(tv%S) ) call set_up_sponge_field(s, tv%S, g, gv, nz, csp)
482 endif
483
484end subroutine dumbbell_initialize_sponges
485