MOM_sponge.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!> Implements sponge regions in isopycnal mode
6module mom_sponge
7
8use mom_coms, only : sum_across_pes
11use mom_error_handler, only : mom_error, fatal, note, warning, is_root_pe
12use mom_file_parser, only : get_param, log_param, log_version, param_file_type
13use mom_grid, only : ocean_grid_type
15use mom_time_manager, only : time_type
19
20! Planned extension: Support for time varying sponge targets.
21
22implicit none ; private
23
24#include <MOM_memory.h>
25
28
29! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
30! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
31! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
32! vary with the Boussinesq approximation, the Boussinesq variant is given first.
33
34!> A structure for creating arrays of pointers to 3D arrays
35type, public :: p3d
36 real, dimension(:,:,:), pointer :: p => null() !< A pointer to a 3D array [various]
37end type p3d
38!> A structure for creating arrays of pointers to 2D arrays
39type, public :: p2d
40 real, dimension(:,:), pointer :: p => null() !< A pointer to a 2D array [various]
41end type p2d
42
43!> This control structure holds memory and parameters for the MOM_sponge module
44type, public :: sponge_cs ; private
45 logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with
46 !! nkml sublayers and nkbl buffer layer.
47 integer :: nz !< The total number of layers.
48 integer :: num_col !< The number of sponge points within the computational domain.
49 integer :: fldno = 0 !< The number of fields which have already been
50 !! registered by calls to set_up_sponge_field
51 integer, pointer :: col_i(:) => null() !< Array of the i-indicies of each of the columns being damped.
52 integer, pointer :: col_j(:) => null() !< Array of the j-indicies of each of the columns being damped.
53 real, pointer :: iresttime_col(:) => null() !< The inverse restoring time of each column [T-1 ~> s-1].
54 real, pointer :: rcv_ml_ref(:) => null() !< The value toward which the mixed layer
55 !! coordinate-density is being damped [R ~> kg m-3].
56 real, pointer :: ref_eta(:,:) => null() !< The value toward which the interface
57 !! heights are being damped [Z ~> m].
58 type(p3d) :: var(max_fields_) !< Pointers to the fields that are being damped.
59 type(p2d) :: ref_val(max_fields_) !< The values to which the fields are damped.
60
61 logical :: do_i_mean_sponge !< If true, apply sponges to the i-mean fields.
62 real, pointer :: iresttime_im(:) => null() !< The inverse restoring time of
63 !! each row for i-mean sponges [T-1 ~> s-1].
64 real, pointer :: rcv_ml_ref_im(:) => null() !! The value toward which the i-mean
65 !< mixed layer coordinate-density is being damped [R ~> kg m-3].
66 real, pointer :: ref_eta_im(:,:) => null() !< The value toward which the i-mean
67 !! interface heights are being damped [Z ~> m].
68 type(p2d) :: ref_val_im(max_fields_) !< The values toward which the i-means of
69 !! fields are damped.
70
71 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
72 !! regulate the timing of diagnostic output.
73 integer :: id_w_sponge = -1 !< A diagnostic ID
74end type sponge_cs
75
76contains
77
78!> This subroutine determines the number of points which are within sponges in
79!! this computational domain. Only points that have positive values of
80!! Iresttime and which mask2dT indicates are ocean points are included in the
81!! sponges. It also stores the target interface heights.
82subroutine initialize_sponge(Iresttime, int_height, G, param_file, CS, GV, &
83 Iresttime_i_mean, int_height_i_mean)
84 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
85 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
86 real, dimension(SZI_(G),SZJ_(G)), &
87 intent(in) :: iresttime !< The inverse of the restoring time [T-1 ~> s-1].
88 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
89 intent(in) :: int_height !< The interface heights to damp back toward [Z ~> m].
90 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
91 type(sponge_cs), pointer :: cs !< A pointer that is set to point to the control
92 !! structure for this module
93 real, dimension(SZJ_(G)), &
94 optional, intent(in) :: iresttime_i_mean !< The inverse of the restoring time for
95 !! the zonal mean properties [T-1 ~> s-1].
96 real, dimension(SZJ_(G),SZK_(GV)+1), &
97 optional, intent(in) :: int_height_i_mean !< The interface heights toward which to
98 !! damp the zonal mean heights [Z ~> m].
99
100
101 ! This include declares and sets the variable "version".
102# include "version_variable.h"
103 character(len=40) :: mdl = "MOM_sponge" ! This module's name.
104 logical :: use_sponge
105 integer :: i, j, k, col, total_sponge_cols
106
107 if (associated(cs)) then
108 call mom_error(warning, "initialize_sponge called with an associated "// &
109 "control structure.")
110 return
111 endif
112
113! Set default, read and log parameters
114 call log_version(param_file, mdl, version)
115 call get_param(param_file, mdl, "SPONGE", use_sponge, &
116 "If true, sponges may be applied anywhere in the domain. "//&
117 "The exact location and properties of those sponges are "//&
118 "specified from MOM_initialization.F90.", default=.false.)
119
120 if (.not.use_sponge) return
121 allocate(cs)
122
123 if (present(iresttime_i_mean) .neqv. present(int_height_i_mean)) &
124 call mom_error(fatal, "initialize_sponge: The optional arguments \n"//&
125 "Iresttime_i_mean and int_height_i_mean must both be present \n"//&
126 "if either one is.")
127
128 cs%do_i_mean_sponge = present(iresttime_i_mean)
129
130 cs%nz = gv%ke
131
132 ! CS%bulkmixedlayer may be set later via a call to set_up_sponge_ML_density.
133 cs%bulkmixedlayer = .false.
134
135 cs%num_col = 0 ; cs%fldno = 0
136 do j=g%jsc,g%jec ; do i=g%isc,g%iec
137 if ((iresttime(i,j) > 0.0) .and. (g%mask2dT(i,j) > 0.0)) &
138 cs%num_col = cs%num_col + 1
139 enddo ; enddo
140
141 if (cs%num_col > 0) then
142
143 allocate(cs%Iresttime_col(cs%num_col), source=0.0)
144 allocate(cs%col_i(cs%num_col), source=0)
145 allocate(cs%col_j(cs%num_col), source=0)
146
147 col = 1
148 do j=g%jsc,g%jec ; do i=g%isc,g%iec
149 if ((iresttime(i,j) > 0.0) .and. (g%mask2dT(i,j) > 0.0)) then
150 cs%col_i(col) = i ; cs%col_j(col) = j
151 cs%Iresttime_col(col) = iresttime(i,j)
152 col = col +1
153 endif
154 enddo ; enddo
155
156 allocate(cs%Ref_eta(cs%nz+1,cs%num_col))
157 do col=1,cs%num_col ; do k=1,cs%nz+1
158 cs%Ref_eta(k,col) = int_height(cs%col_i(col),cs%col_j(col),k)
159 enddo ; enddo
160
161 endif
162
163 if (cs%do_i_mean_sponge) then
164 allocate(cs%Iresttime_im(g%jsd:g%jed), source=0.0)
165 allocate(cs%Ref_eta_im(g%jsd:g%jed,gv%ke+1), source=0.0)
166
167 do j=g%jsc,g%jec
168 cs%Iresttime_im(j) = iresttime_i_mean(j)
169 enddo
170 do k=1,cs%nz+1 ; do j=g%jsc,g%jec
171 cs%Ref_eta_im(j,k) = int_height_i_mean(j,k)
172 enddo ; enddo
173 endif
174
175 total_sponge_cols = cs%num_col
176 call sum_across_pes(total_sponge_cols)
177
178 call log_param(param_file, mdl, "!Total sponge columns", total_sponge_cols, &
179 "The total number of columns where sponges are applied.")
180
181end subroutine initialize_sponge
182
183!> This subroutine sets up diagnostics for the sponges. It is separate
184!! from initialize_sponge because it requires fields that are not readily
185!! available where initialize_sponge is called.
186subroutine init_sponge_diags(Time, G, GV, US, diag, CS)
187 type(time_type), target, intent(in) :: time !< The current model time
188 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
189 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
190 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
191 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic output
192 type(sponge_cs), pointer :: cs !< A pointer to the control structure for this module that
193 !! is set by a previous call to initialize_sponge.
194
195 if (.not.associated(cs)) return
196
197 cs%diag => diag
198 cs%id_w_sponge = register_diag_field('ocean_model', 'w_sponge', diag%axesTi, &
199 time, 'The diapycnal motion due to the sponges', 'm s-1', conversion=gv%H_to_m*us%s_to_T)
200
201end subroutine init_sponge_diags
202
203!> This subroutine stores the reference profile for the variable whose
204!! address is given by f_ptr. nlay is the number of layers in this variable.
205subroutine set_up_sponge_field(sp_val, f_ptr, G, GV, nlay, CS, sp_val_i_mean)
206 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
207 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
208 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
209 intent(in) :: sp_val !< The reference profiles of the quantity being registered [various]
210 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
211 target, intent(in) :: f_ptr !< a pointer to the field which will be damped [various]
212 integer, intent(in) :: nlay !< the number of layers in this quantity
213 type(sponge_cs), pointer :: cs !< A pointer to the control structure for this module that
214 !! is set by a previous call to initialize_sponge.
215 real, dimension(SZJ_(G),SZK_(GV)),&
216 optional, intent(in) :: sp_val_i_mean !< The i-mean reference value for
217 !! this field with i-mean sponges [various]
218
219 integer :: j, k, col
220 character(len=256) :: mesg ! String for error messages
221
222 if (.not.associated(cs)) return
223
224 cs%fldno = cs%fldno + 1
225
226 if (cs%fldno > max_fields_) then
227 write(mesg,'("Increase MAX_FIELDS_ to at least ",I0," in MOM_memory.h or decrease &
228 &the number of fields to be damped in the call to &
229 &initialize_sponge." )') cs%fldno
230 call mom_error(fatal,"set_up_sponge_field: "//mesg)
231 endif
232
233 allocate(cs%Ref_val(cs%fldno)%p(cs%nz,cs%num_col), source=0.0)
234 do col=1,cs%num_col
235 do k=1,nlay
236 cs%Ref_val(cs%fldno)%p(k,col) = sp_val(cs%col_i(col),cs%col_j(col),k)
237 enddo
238 do k=nlay+1,cs%nz
239 cs%Ref_val(cs%fldno)%p(k,col) = 0.0
240 enddo
241 enddo
242
243 cs%var(cs%fldno)%p => f_ptr
244
245 if (nlay/=cs%nz) then
246 write(mesg,'("Danger: Sponge reference fields require nz (",I0,") layers.&
247 & A field with ",I0," layers was passed to set_up_sponge_field.")') &
248 cs%nz, nlay
249 if (is_root_pe()) call mom_error(warning, "set_up_sponge_field: "//mesg)
250 endif
251
252 if (cs%do_i_mean_sponge) then
253 if (.not.present(sp_val_i_mean)) call mom_error(fatal, &
254 "set_up_sponge_field: sp_val_i_mean must be present with i-mean sponges.")
255
256 allocate(cs%Ref_val_im(cs%fldno)%p(g%jsd:g%jed,cs%nz), source=0.0)
257 do k=1,cs%nz ; do j=g%jsc,g%jec
258 cs%Ref_val_im(cs%fldno)%p(j,k) = sp_val_i_mean(j,k)
259 enddo ; enddo
260 endif
261
262end subroutine set_up_sponge_field
263
264
265!> This subroutine stores the reference value for mixed layer density. It is handled differently
266!! from other values because it is only used in determining which layers can be inflated.
267subroutine set_up_sponge_ml_density(sp_val, G, CS, sp_val_i_mean)
268 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
269 real, dimension(SZI_(G),SZJ_(G)), &
270 intent(in) :: sp_val !< The reference values of the mixed layer density [R ~> kg m-3]
271 type(sponge_cs), pointer :: cs !< A pointer to the control structure for this module that is
272 !! set by a previous call to initialize_sponge.
273 ! The contents of this structure are intent(inout) here.
274 real, dimension(SZJ_(G)), &
275 optional, intent(in) :: sp_val_i_mean !< the reference values of the zonal mean mixed
276 !! layer density [R ~> kg m-3], for use if Iresttime_i_mean > 0.
277
278 integer :: j, col
279
280 if (.not.associated(cs)) return
281
282 if (associated(cs%Rcv_ml_ref)) then
283 call mom_error(fatal, "set_up_sponge_ML_density appears to have been "//&
284 "called twice.")
285 endif
286
287 cs%bulkmixedlayer = .true.
288 allocate(cs%Rcv_ml_ref(cs%num_col), source=0.0)
289 do col=1,cs%num_col
290 cs%Rcv_ml_ref(col) = sp_val(cs%col_i(col),cs%col_j(col))
291 enddo
292
293 if (cs%do_i_mean_sponge) then
294 if (.not.present(sp_val_i_mean)) call mom_error(fatal, &
295 "set_up_sponge_field: sp_val_i_mean must be present with i-mean sponges.")
296
297 allocate(cs%Rcv_ml_ref_im(g%jsd:g%jed), source=0.0)
298 do j=g%jsc,g%jec
299 cs%Rcv_ml_ref_im(j) = sp_val_i_mean(j)
300 enddo
301 endif
302
303end subroutine set_up_sponge_ml_density
304
305!> This subroutine applies damping to the layers thicknesses, mixed layer buoyancy, and a variety of
306!! tracers for every column where there is damping.
307subroutine apply_sponge(h, tv, dt, G, GV, US, ea, eb, CS, Rcv_ml)
308 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
309 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
310 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
311 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
312 intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
313 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
314 !! thermodynamic variables
315 real, intent(in) :: dt !< The amount of time covered by this call [T ~> s].
316 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
317 intent(inout) :: ea !< An array to which the amount of fluid entrained
318 !! from the layer above during this call will be
319 !! added [H ~> m or kg m-2].
320 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
321 intent(inout) :: eb !< An array to which the amount of fluid entrained
322 !! from the layer below during this call will be
323 !! added [H ~> m or kg m-2].
324 type(sponge_cs), pointer :: cs !< A pointer to the control structure for this module
325 !! that is set by a previous call to initialize_sponge.
326 real, dimension(SZI_(G),SZJ_(G)), &
327 optional, intent(inout) :: rcv_ml !< The coordinate density of the mixed layer [R ~> kg m-3].
328
329 ! Local variables
330 real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: &
331 w_int, & ! Water moved upward across an interface within a timestep,
332 ! [H ~> m or kg m-2].
333 e_d ! Interface heights that are dilated to have a value of 0
334 ! at the surface [Z ~> m].
335 real, dimension(SZI_(G), SZJ_(G)) :: &
336 eta_anom, & ! Anomalies in the interface height, relative to the i-mean
337 ! target value [Z ~> m].
338 fld_anom ! Anomalies in a tracer concentration, relative to the
339 ! i-mean target value [various]
340 real, dimension(SZJ_(G), SZK_(GV)+1) :: &
341 eta_mean_anom ! The i-mean interface height anomalies [Z ~> m].
342 real, allocatable, dimension(:,:,:) :: &
343 fld_mean_anom ! The i-mean tracer concentration anomalies [various]
344 real, dimension(SZI_(G), SZK_(GV)+1) :: &
345 h_above, & ! The total thickness above an interface [H ~> m or kg m-2].
346 h_below ! The total thickness below an interface [H ~> m or kg m-2].
347 real, dimension(SZI_(G)) :: &
348 dilate ! A nondimensional factor by which to dilate layers to
349 ! give 0 at the surface [nondim].
350
351 real :: e(szk_(gv)+1) ! The interface heights [Z ~> m], usually negative.
352 real :: dz_to_h(szk_(gv)+1) ! Factors used to convert interface height movement
353 ! to thickness fluxes [H Z-1 ~> nondim or kg m-3]
354 real :: e0 ! The height of the free surface [Z ~> m].
355 real :: e_str ! A nondimensional amount by which the reference
356 ! profile must be stretched for the free surfaces
357 ! heights in the two profiles to agree [nondim].
358 real :: w_mean ! The vertical displacement of water moving upward through an
359 ! interface within 1 timestep [Z ~> m].
360 real :: w ! The thickness of water moving upward through an
361 ! interface within 1 timestep [H ~> m or kg m-2].
362 real :: wm ! wm is w if w is negative and 0 otherwise [H ~> m or kg m-2].
363 real :: wb ! w at the interface below a layer [H ~> m or kg m-2].
364 real :: wpb ! wpb is wb if wb is positive and 0 otherwise [H ~> m or kg m-2].
365 real :: ea_k ! Water entrained from above within a timestep [H ~> m or kg m-2]
366 real :: eb_k ! Water entrained from below within a timestep [H ~> m or kg m-2]
367 real :: damp ! The timestep times the local damping coefficient [nondim].
368 real :: i1pdamp ! I1pdamp is 1/(1 + damp). [nondim]
369 real :: damp_1pdamp ! damp_1pdamp is damp/(1 + damp). [nondim]
370 real :: idt ! The inverse of the timestep [T-1 ~> s-1]
371 integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz
372 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
373
374 if (.not.associated(cs)) return
375 if (cs%bulkmixedlayer) nkmb = gv%nk_rho_varies
376 if (cs%bulkmixedlayer .and. (.not.present(rcv_ml))) &
377 call mom_error(fatal, "Rml must be provided to apply_sponge when using "//&
378 "a bulk mixed layer.")
379
380 if ((cs%id_w_sponge > 0) .or. cs%do_i_mean_sponge) then
381 do k=1,nz+1 ; do j=js,je ; do i=is,ie
382 w_int(i,j,k) = 0.0
383 enddo ; enddo ; enddo
384 endif
385
386 if (cs%do_i_mean_sponge) then
387 ! Apply forcing to restore the zonal-mean properties to prescribed values.
388
389 if (cs%bulkmixedlayer) call mom_error(fatal, "apply_sponge is not yet set up to "//&
390 "work properly with i-mean sponges and a bulk mixed layer.")
391
392 do j=js,je ; do i=is,ie ; e_d(i,j,nz+1) = -g%bathyT(i,j) ; enddo ; enddo
393 if ((.not.gv%Boussinesq) .and. allocated(tv%SpV_avg)) then
394 do k=nz,1,-1 ; do j=js,je ; do i=is,ie
395 e_d(i,j,k) = e_d(i,j,k+1) + gv%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
396 enddo ; enddo ; enddo
397 else
398 do k=nz,1,-1 ; do j=js,je ; do i=is,ie
399 e_d(i,j,k) = e_d(i,j,k+1) + h(i,j,k)*gv%H_to_Z
400 enddo ; enddo ; enddo
401 endif
402 do j=js,je
403 do i=is,ie
404 dilate(i) = (g%bathyT(i,j) + g%Z_ref) / (e_d(i,j,1) + g%bathyT(i,j))
405 enddo
406 do k=1,nz+1 ; do i=is,ie
407 e_d(i,j,k) = dilate(i) * (e_d(i,j,k) + g%bathyT(i,j)) - (g%bathyT(i,j) + g%Z_ref)
408 enddo ; enddo
409 enddo
410
411 do k=2,nz
412 do j=js,je ; do i=is,ie
413 eta_anom(i,j) = e_d(i,j,k) - cs%Ref_eta_im(j,k)
414 if (cs%Ref_eta_im(j,k) < -(g%bathyT(i,j) + g%Z_ref)) eta_anom(i,j) = 0.0
415 enddo ; enddo
416 call global_i_mean(eta_anom(:,:), eta_mean_anom(:,k), g, tmp_scale=us%Z_to_m)
417 enddo
418
419 if (cs%fldno > 0) allocate(fld_mean_anom(g%isd:g%ied,nz,cs%fldno))
420 do m=1,cs%fldno
421 do j=js,je ; do i=is,ie
422 fld_anom(i,j) = cs%var(m)%p(i,j,k) - cs%Ref_val_im(m)%p(j,k)
423 enddo ; enddo
424 call global_i_mean(fld_anom(:,:), fld_mean_anom(:,k,m), g, h(:,:,k))
425 enddo
426
427 do j=js,je ; if (cs%Iresttime_im(j) > 0.0) then
428 damp = dt * cs%Iresttime_im(j) ; damp_1pdamp = damp / (1.0 + damp)
429
430 do i=is,ie
431 h_above(i,1) = 0.0 ; h_below(i,nz+1) = 0.0
432 enddo
433 do k=nz,1,-1 ; do i=is,ie
434 h_below(i,k) = h_below(i,k+1) + max(h(i,j,k)-gv%Angstrom_H, 0.0)
435 enddo ; enddo
436 do k=2,nz+1 ; do i=is,ie
437 h_above(i,k) = h_above(i,k-1) + max(h(i,j,k-1)-gv%Angstrom_H, 0.0)
438 enddo ; enddo
439
440 ! In both blocks below, w is positive for an upward (lightward) flux of mass,
441 ! resulting in the downward movement of an interface.
442 if ((.not.gv%Boussinesq) .and. allocated(tv%SpV_avg)) then
443 do k=2,nz
444 w_mean = damp_1pdamp * eta_mean_anom(j,k)
445 do i=is,ie
446 w = w_mean * 2.0*gv%RZ_to_H / (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k))
447 if (w > 0.0) then
448 w_int(i,j,k) = min(w, h_below(i,k))
449 eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,k)
450 else
451 w_int(i,j,k) = max(w, -h_above(i,k))
452 ea(i,j,k) = ea(i,j,k) - w_int(i,j,k)
453 endif
454 enddo
455 enddo
456 else
457 do k=2,nz
458 w = damp_1pdamp * eta_mean_anom(j,k) * gv%Z_to_H
459 if (w > 0.0) then
460 do i=is,ie
461 w_int(i,j,k) = min(w, h_below(i,k))
462 eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,k)
463 enddo
464 else
465 do i=is,ie
466 w_int(i,j,k) = max(w, -h_above(i,k))
467 ea(i,j,k) = ea(i,j,k) - w_int(i,j,k)
468 enddo
469 endif
470 enddo
471 endif
472 do k=1,nz ; do i=is,ie
473 ea_k = max(0.0, -w_int(i,j,k))
474 eb_k = max(0.0, w_int(i,j,k+1))
475 do m=1,cs%fldno
476 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
477 cs%Ref_val_im(m)%p(j,k) * (ea_k + eb_k)) / &
478 (h(i,j,k) + (ea_k + eb_k)) - &
479 damp_1pdamp * fld_mean_anom(j,k,m)
480 enddo
481
482 h(i,j,k) = max(h(i,j,k) + (w_int(i,j,k+1) - w_int(i,j,k)), &
483 min(h(i,j,k), gv%Angstrom_H))
484 enddo ; enddo
485 endif ; enddo
486
487 if (cs%fldno > 0) deallocate(fld_mean_anom)
488
489 endif
490
491 do c=1,cs%num_col
492 i = cs%col_i(c) ; j = cs%col_j(c)
493 damp = dt * cs%Iresttime_col(c)
494
495 e(1) = 0.0 ; e0 = 0.0
496 if ((.not.gv%Boussinesq) .and. allocated(tv%SpV_avg)) then
497 do k=1,nz
498 e(k+1) = e(k) - gv%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k)
499 enddo
500 dz_to_h(1) = gv%RZ_to_H / tv%SpV_avg(i,j,1)
501 do k=2,nz
502 dz_to_h(k) = 2.0*gv%RZ_to_H / (tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k))
503 enddo
504 else
505 do k=1,nz
506 e(k+1) = e(k) - h(i,j,k)*gv%H_to_Z
507 dz_to_h(k) = gv%Z_to_H
508 enddo
509 endif
510 e_str = e(nz+1) / cs%Ref_eta(nz+1,c)
511
512 if ( cs%bulkmixedlayer ) then
513 i1pdamp = 1.0 / (1.0 + damp)
514 if (associated(cs%Rcv_ml_ref)) &
515 rcv_ml(i,j) = i1pdamp * (rcv_ml(i,j) + cs%Rcv_ml_ref(c)*damp)
516 do k=1,nkmb
517 do m=1,cs%fldno
518 cs%var(m)%p(i,j,k) = i1pdamp * &
519 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(k,c)*damp)
520 enddo
521 enddo
522
523 wpb = 0.0 ; wb = 0.0
524 do k=nz,nkmb+1,-1
525 if (gv%Rlay(k) > rcv_ml(i,j)) then
526 w = min((((e(k)-e0) - e_str*cs%Ref_eta(k,c)) * damp)*dz_to_h(k), &
527 ((wb + h(i,j,k)) - gv%Angstrom_H))
528 wm = 0.5*(w-abs(w))
529 do m=1,cs%fldno
530 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
531 cs%Ref_val(m)%p(k,c)*(damp*h(i,j,k) + (wpb - wm))) / &
532 (h(i,j,k)*(1.0 + damp) + (wpb - wm))
533 enddo
534 else
535 do m=1,cs%fldno
536 cs%var(m)%p(i,j,k) = i1pdamp * &
537 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(k,c)*damp)
538 enddo
539 w = wb + (h(i,j,k) - gv%Angstrom_H)
540 wm = 0.5*(w-abs(w))
541 endif
542 eb(i,j,k) = eb(i,j,k) + wpb
543 ea(i,j,k) = ea(i,j,k) - wm
544 h(i,j,k) = h(i,j,k) + (wb - w)
545 wb = w
546 wpb = w - wm
547 enddo
548
549 if (wb < 0) then
550 do k=nkmb,1,-1
551 w = min((wb + (h(i,j,k) - gv%Angstrom_H)),0.0)
552 h(i,j,k) = h(i,j,k) + (wb - w)
553 ea(i,j,k) = ea(i,j,k) - w
554 wb = w
555 enddo
556 else
557 w = wb
558 do k=gv%nkml,nkmb
559 eb(i,j,k) = eb(i,j,k) + w
560 enddo
561
562 k = gv%nkml
563 h(i,j,k) = h(i,j,k) + w
564 do m=1,cs%fldno
565 cs%var(m)%p(i,j,k) = (cs%var(m)%p(i,j,k)*h(i,j,k) + &
566 cs%Ref_val(m)%p(k,c)*w) / (h(i,j,k) + w)
567 enddo
568 endif
569
570 do k=1,nkmb
571 do m=1,cs%fldno
572 cs%var(m)%p(i,j,k) = i1pdamp * &
573 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(gv%nkml,c)*damp)
574 enddo
575 enddo
576
577 else ! not BULKMIXEDLAYER
578
579 wpb = 0.0
580 wb = 0.0
581 do k=nz,1,-1
582 w = min((((e(k)-e0) - e_str*cs%Ref_eta(k,c)) * damp)*dz_to_h(k), &
583 ((wb + h(i,j,k)) - gv%Angstrom_H))
584 wm = 0.5*(w - abs(w))
585 do m=1,cs%fldno
586 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
587 cs%Ref_val(m)%p(k,c) * (damp*h(i,j,k) + (wpb - wm))) / &
588 (h(i,j,k)*(1.0 + damp) + (wpb - wm))
589 enddo
590 eb(i,j,k) = eb(i,j,k) + wpb
591 ea(i,j,k) = ea(i,j,k) - wm
592 h(i,j,k) = h(i,j,k) + (wb - w)
593 wb = w
594 wpb = w - wm
595 enddo
596
597 endif ! end BULKMIXEDLAYER
598 enddo ! end of c loop
599
600 if (associated(cs%diag)) then ; if (query_averaging_enabled(cs%diag)) then
601 if (cs%id_w_sponge > 0) then
602 idt = 1.0 / dt
603 do k=1,nz+1 ; do j=js,je ; do i=is,ie
604 w_int(i,j,k) = w_int(i,j,k) * idt ! Scale values by clobbering array since it is local
605 enddo ; enddo ; enddo
606 call post_data(cs%id_w_sponge, w_int(:,:,:), cs%diag)
607 endif
608 endif ; endif
609
610end subroutine apply_sponge
611
612!> This call deallocates any memory in the sponge control structure.
613subroutine sponge_end(CS)
614 type(sponge_cs), pointer :: cs !< A pointer to the control structure for this module
615 !! that is set by a previous call to initialize_sponge.
616 integer :: m
617
618 if (.not.associated(cs)) return
619
620 if (associated(cs%col_i)) deallocate(cs%col_i)
621 if (associated(cs%col_j)) deallocate(cs%col_j)
622
623 if (associated(cs%Iresttime_col)) deallocate(cs%Iresttime_col)
624 if (associated(cs%Rcv_ml_ref)) deallocate(cs%Rcv_ml_ref)
625 if (associated(cs%Ref_eta)) deallocate(cs%Ref_eta)
626
627 if (associated(cs%Iresttime_im)) deallocate(cs%Iresttime_im)
628 if (associated(cs%Rcv_ml_ref_im)) deallocate(cs%Rcv_ml_ref_im)
629 if (associated(cs%Ref_eta_im)) deallocate(cs%Ref_eta_im)
630
631 do m=1,cs%fldno
632 if (associated(cs%Ref_val(cs%fldno)%p)) deallocate(cs%Ref_val(cs%fldno)%p)
633 if (associated(cs%Ref_val_im(cs%fldno)%p)) &
634 deallocate(cs%Ref_val_im(cs%fldno)%p)
635 enddo
636
637 deallocate(cs)
638
639end subroutine sponge_end
640
641!> \namespace mom_sponge
642!!
643!! By Robert Hallberg, March 1999-June 2000
644!!
645!! This program contains the subroutines that implement sponge
646!! regions, in which the stratification and water mass properties
647!! are damped toward some profiles. There are three externally
648!! callable subroutines in this file.
649!!
650!! initialize_sponge determines the mapping from the model
651!! variables into the arrays of damped columns. This remapping is
652!! done for efficiency and to conserve memory. Only columns which
653!! have positive inverse damping times and which are deeper than a
654!! supplied depth are placed in sponges. The inverse damping
655!! time is also stored in this subroutine, and memory is allocated
656!! for all of the reference profiles which will subsequently be
657!! provided through calls to set_up_sponge_field. The first two
658!! arguments are a two-dimensional array containing the damping
659!! rates, and the interface heights to damp towards.
660!!
661!! set_up_sponge_field is called to provide a reference profile
662!! and the location of the field that will be damped back toward
663!! that reference profile. A third argument, the number of layers
664!! in the field is also provided, but this should always be nz.
665!!
666!! Apply_sponge damps all of the fields that have been registered
667!! with set_up_sponge_field toward their reference profiles. The
668!! four arguments are the thickness to be damped, the amount of time
669!! over which the damping occurs, and arrays to which the movement
670!! of fluid into a layer from above and below will be added. The
671!! effect on momentum of the sponge may be accounted for later using
672!! the movement of water recorded in these later arrays.
673
674end module mom_sponge