19 implicit none ;
private
21 #include <MOM_memory.h>
33 real,
dimension(:,:,:),
pointer :: p => null()
37 real,
dimension(:,:),
pointer :: p => null()
42 logical :: bulkmixedlayer
56 integer,
pointer :: col_i(:) => null()
57 integer,
pointer :: col_j(:) => null()
58 real,
pointer :: iresttime_col(:) => null()
59 real,
pointer :: rcv_ml_ref(:) => null()
61 real,
pointer :: ref_eta(:,:) => null()
63 type(
p3d) :: var(max_fields_)
64 type(
p2d) :: ref_val(max_fields_)
66 logical :: do_i_mean_sponge
67 real,
pointer :: iresttime_im(:) => null()
69 real,
pointer :: rcv_ml_ref_im(:) => null()
71 real,
pointer :: ref_eta_im(:,:) => null()
73 type(
p2d) :: ref_val_im(max_fields_)
78 integer :: id_w_sponge = -1
89 Iresttime_i_mean, int_height_i_mean)
91 real,
dimension(SZI_(G),SZJ_(G)), &
92 intent(in) :: iresttime
93 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
94 intent(in) :: int_height
99 real,
dimension(SZJ_(G)), &
100 optional,
intent(in) :: iresttime_i_mean
102 real,
dimension(SZJ_(G),SZK_(G)+1), &
103 optional,
intent(in) :: int_height_i_mean
108 #include "version_variable.h"
109 character(len=40) :: mdl =
"MOM_sponge"
110 logical :: use_sponge
111 integer :: i, j, k, col, total_sponge_cols
113 if (
associated(cs))
then
114 call mom_error(warning,
"initialize_sponge called with an associated "// &
115 "control structure.")
121 call get_param(param_file, mdl,
"SPONGE", use_sponge, &
122 "If true, sponges may be applied anywhere in the domain. "//&
123 "The exact location and properties of those sponges are "//&
124 "specified from MOM_initialization.F90.", default=.false.)
126 if (.not.use_sponge)
return
129 if (
present(iresttime_i_mean) .neqv.
present(int_height_i_mean)) &
130 call mom_error(fatal,
"initialize_sponge: The optional arguments \n"//&
131 "Iresttime_i_mean and int_height_i_mean must both be present \n"//&
134 cs%do_i_mean_sponge =
present(iresttime_i_mean)
140 cs%bulkmixedlayer = .false.
142 cs%num_col = 0 ; cs%fldno = 0
143 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
144 if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) &
145 cs%num_col = cs%num_col + 1
148 if (cs%num_col > 0)
then
150 allocate(cs%Iresttime_col(cs%num_col)) ; cs%Iresttime_col = 0.0
151 allocate(cs%col_i(cs%num_col)) ; cs%col_i = 0
152 allocate(cs%col_j(cs%num_col)) ; cs%col_j = 0
155 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
156 if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0))
then
157 cs%col_i(col) = i ; cs%col_j(col) = j
158 cs%Iresttime_col(col) = g%US%T_to_s*iresttime(i,j)
163 allocate(cs%Ref_eta(cs%nz+1,cs%num_col))
164 do col=1,cs%num_col ;
do k=1,cs%nz+1
165 cs%Ref_eta(k,col) = int_height(cs%col_i(col),cs%col_j(col),k)
170 if (cs%do_i_mean_sponge)
then
171 allocate(cs%Iresttime_im(g%jsd:g%jed)) ; cs%Iresttime_im(:) = 0.0
172 allocate(cs%Ref_eta_im(g%jsd:g%jed,g%ke+1)) ; cs%Ref_eta_im(:,:) = 0.0
175 cs%Iresttime_im(j) = g%US%T_to_s*iresttime_i_mean(j)
177 do k=1,cs%nz+1 ;
do j=g%jsc,g%jec
178 cs%Ref_eta_im(j,k) = int_height_i_mean(j,k)
182 total_sponge_cols = cs%num_col
183 call sum_across_pes(total_sponge_cols)
185 call log_param(param_file, mdl,
"!Total sponge columns", total_sponge_cols, &
186 "The total number of columns where sponges are applied.")
194 type(time_type),
target,
intent(in) :: time
198 type(
diag_ctrl),
target,
intent(inout) :: diag
202 if (.not.
associated(cs))
return
205 cs%id_w_sponge = register_diag_field(
'ocean_model',
'w_sponge', diag%axesTi, &
206 time,
'The diapycnal motion due to the sponges',
'm s-1', conversion=us%s_to_T)
215 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
217 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
218 target,
intent(in) :: f_ptr
219 integer,
intent(in) :: nlay
222 real,
dimension(SZJ_(G),SZK_(G)),&
223 optional,
intent(in) :: sp_val_i_mean
227 character(len=256) :: mesg
229 if (.not.
associated(cs))
return
231 cs%fldno = cs%fldno + 1
233 if (cs%fldno > max_fields_)
then
234 write(mesg,
'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease &
235 &the number of fields to be damped in the call to &
236 &initialize_sponge." )') cs%fldno
237 call mom_error(fatal,
"set_up_sponge_field: "//mesg)
240 allocate(cs%Ref_val(cs%fldno)%p(cs%nz,cs%num_col))
241 cs%Ref_val(cs%fldno)%p(:,:) = 0.0
244 cs%Ref_val(cs%fldno)%p(k,col) = sp_val(cs%col_i(col),cs%col_j(col),k)
247 cs%Ref_val(cs%fldno)%p(k,col) = 0.0
251 cs%var(cs%fldno)%p => f_ptr
253 if (nlay/=cs%nz)
then
254 write(mesg,
'("Danger: Sponge reference fields require nz (",I3,") layers.&
255 & A field with ",I3," layers was passed to set_up_sponge_field.")') &
260 if (cs%do_i_mean_sponge)
then
261 if (.not.
present(sp_val_i_mean))
call mom_error(fatal, &
262 "set_up_sponge_field: sp_val_i_mean must be present with i-mean sponges.")
264 allocate(cs%Ref_val_im(cs%fldno)%p(cs%jsd:cs%jed,cs%nz))
265 cs%Ref_val(cs%fldno)%p(:,:) = 0.0
266 do k=1,cs%nz ;
do j=cs%jsc,cs%jec
267 cs%Ref_val_im(cs%fldno)%p(j,k) = sp_val_i_mean(j,k)
278 real,
dimension(SZI_(G),SZJ_(G)), &
282 real,
dimension(SZJ_(G)), &
283 optional,
intent(in) :: sp_val_i_mean
294 character(len=256) :: mesg
296 if (.not.
associated(cs))
return
298 if (
associated(cs%Rcv_ml_ref))
then
299 call mom_error(fatal,
"set_up_sponge_ML_density appears to have been "//&
303 cs%bulkmixedlayer = .true.
304 allocate(cs%Rcv_ml_ref(cs%num_col)) ; cs%Rcv_ml_ref(:) = 0.0
306 cs%Rcv_ml_ref(col) = sp_val(cs%col_i(col),cs%col_j(col))
309 if (cs%do_i_mean_sponge)
then
310 if (.not.
present(sp_val_i_mean))
call mom_error(fatal, &
311 "set_up_sponge_field: sp_val_i_mean must be present with i-mean sponges.")
313 allocate(cs%Rcv_ml_ref_im(cs%jsd:cs%jed)) ; cs%Rcv_ml_ref_im(:) = 0.0
315 cs%Rcv_ml_ref_im(j) = sp_val_i_mean(j)
323 subroutine apply_sponge(h, dt, G, GV, US, ea, eb, CS, Rcv_ml)
327 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
329 real,
intent(in) :: dt
330 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
334 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
340 real,
dimension(SZI_(G),SZJ_(G)), &
341 optional,
intent(inout) :: rcv_ml
348 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
353 real,
dimension(SZI_(G), SZJ_(G)) :: &
358 real,
dimension(SZJ_(G), SZK_(G)+1) :: &
360 real,
allocatable,
dimension(:,:,:) :: &
362 real,
dimension(SZI_(G), SZK_(G)+1) :: &
365 real,
dimension(SZI_(G)) :: &
384 integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz
385 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
387 if (.not.
associated(cs))
return
388 if (cs%bulkmixedlayer) nkmb = gv%nk_rho_varies
389 if (cs%bulkmixedlayer .and. (.not.
present(rcv_ml))) &
390 call mom_error(fatal,
"Rml must be provided to apply_sponge when using "//&
391 "a bulk mixed layer.")
393 if ((cs%id_w_sponge > 0) .or. cs%do_i_mean_sponge)
then
394 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
396 enddo ;
enddo ;
enddo
399 if (cs%do_i_mean_sponge)
then
402 if (cs%bulkmixedlayer)
call mom_error(fatal,
"apply_sponge is not yet set up to "//&
403 "work properly with i-mean sponges and a bulk mixed layer.")
405 do j=js,je ;
do i=is,ie ; e_d(i,j,nz+1) = -g%bathyT(i,j) ;
enddo ;
enddo
406 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
407 e_d(i,j,k) = e_d(i,j,k+1) + h(i,j,k)*gv%H_to_Z
408 enddo ;
enddo ;
enddo
411 dilate(i) = g%bathyT(i,j) / (e_d(i,j,1) + g%bathyT(i,j))
413 do k=1,nz+1 ;
do i=is,ie
414 e_d(i,j,k) = dilate(i) * (e_d(i,j,k) + g%bathyT(i,j)) - g%bathyT(i,j)
419 do j=js,je ;
do i=is,ie
420 eta_anom(i,j) = e_d(i,j,k) - cs%Ref_eta_im(j,k)
421 if (cs%Ref_eta_im(j,k) < -g%bathyT(i,j)) eta_anom(i,j) = 0.0
423 call global_i_mean(eta_anom(:,:), eta_mean_anom(:,k), g, tmp_scale=us%Z_to_m)
426 if (cs%fldno > 0)
allocate(fld_mean_anom(g%isd:g%ied,nz,cs%fldno))
428 do j=js,je ;
do i=is,ie
429 fld_anom(i,j) = cs%var(m)%p(i,j,k) - cs%Ref_val_im(m)%p(j,k)
431 call global_i_mean(fld_anom(:,:), fld_mean_anom(:,k,m), g, h(:,:,k))
434 do j=js,je ;
if (cs%Iresttime_im(j) > 0.0)
then
435 damp = dt * cs%Iresttime_im(j) ; damp_1pdamp = damp / (1.0 + damp)
438 h_above(i,1) = 0.0 ; h_below(i,nz+1) = 0.0
440 do k=nz,1,-1 ;
do i=is,ie
441 h_below(i,k) = h_below(i,k+1) + max(h(i,j,k)-gv%Angstrom_H, 0.0)
443 do k=2,nz+1 ;
do i=is,ie
444 h_above(i,k) = h_above(i,k-1) + max(h(i,j,k-1)-gv%Angstrom_H, 0.0)
449 w = damp_1pdamp * eta_mean_anom(j,k) * gv%Z_to_H
452 w_int(i,j,k) = min(w, h_below(i,k))
453 eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,k)
455 w_int(i,j,k) = max(w, -h_above(i,k))
456 ea(i,j,k) = ea(i,j,k) - w_int(i,j,k)
460 do k=1,nz ;
do i=is,ie
461 ea_k = max(0.0, -w_int(i,j,k))
462 eb_k = max(0.0, w_int(i,j,k+1))
464 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
465 cs%Ref_val_im(m)%p(j,k) * (ea_k + eb_k)) / &
466 (h(i,j,k) + (ea_k + eb_k)) - &
467 damp_1pdamp * fld_mean_anom(j,k,m)
470 h(i,j,k) = max(h(i,j,k) + (w_int(i,j,k+1) - w_int(i,j,k)), &
471 min(h(i,j,k), gv%Angstrom_H))
475 if (cs%fldno > 0)
deallocate(fld_mean_anom)
482 i = cs%col_i(c) ; j = cs%col_j(c)
483 damp = dt * cs%Iresttime_col(c)
485 e(1) = 0.0 ; e0 = 0.0
487 e(k+1) = e(k) - h(i,j,k)*gv%H_to_Z
489 e_str = e(nz+1) / cs%Ref_eta(nz+1,c)
491 if ( cs%bulkmixedlayer )
then
492 i1pdamp = 1.0 / (1.0 + damp)
493 if (
associated(cs%Rcv_ml_ref)) &
494 rcv_ml(i,j) = i1pdamp * (rcv_ml(i,j) + cs%Rcv_ml_ref(c)*damp)
497 cs%var(m)%p(i,j,k) = i1pdamp * &
498 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(k,c)*damp)
504 if (gv%Rlay(k) > rcv_ml(i,j))
then
505 w = min((((e(k)-e0) - e_str*cs%Ref_eta(k,c)) * damp)*gv%Z_to_H, &
506 ((wb + h(i,j,k)) - gv%Angstrom_H))
509 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
510 cs%Ref_val(m)%p(k,c)*(damp*h(i,j,k) + (wpb - wm))) / &
511 (h(i,j,k)*(1.0 + damp) + (wpb - wm))
515 cs%var(m)%p(i,j,k) = i1pdamp * &
516 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(k,c)*damp)
518 w = wb + (h(i,j,k) - gv%Angstrom_H)
521 eb(i,j,k) = eb(i,j,k) + wpb
522 ea(i,j,k) = ea(i,j,k) - wm
523 h(i,j,k) = h(i,j,k) + (wb - w)
530 w = min((wb + (h(i,j,k) - gv%Angstrom_H)),0.0)
531 h(i,j,k) = h(i,j,k) + (wb - w)
532 ea(i,j,k) = ea(i,j,k) - w
538 eb(i,j,k) = eb(i,j,k) + w
542 h(i,j,k) = h(i,j,k) + w
544 cs%var(m)%p(i,j,k) = (cs%var(m)%p(i,j,k)*h(i,j,k) + &
545 cs%Ref_val(m)%p(k,c)*w) / (h(i,j,k) + w)
551 cs%var(m)%p(i,j,k) = i1pdamp * &
552 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(gv%nkml,c)*damp)
561 w = min((((e(k)-e0) - e_str*cs%Ref_eta(k,c)) * damp)*gv%Z_to_H, &
562 ((wb + h(i,j,k)) - gv%Angstrom_H))
563 wm = 0.5*(w - abs(w))
565 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
566 cs%Ref_val(m)%p(k,c) * (damp*h(i,j,k) + (wpb - wm))) / &
567 (h(i,j,k)*(1.0 + damp) + (wpb - wm))
569 eb(i,j,k) = eb(i,j,k) + wpb
570 ea(i,j,k) = ea(i,j,k) - wm
571 h(i,j,k) = h(i,j,k) + (wb - w)
579 if (
associated(cs%diag))
then ;
if (query_averaging_enabled(cs%diag))
then
580 if (cs%id_w_sponge > 0)
then
582 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
583 w_int(i,j,k) = w_int(i,j,k) * idt
584 enddo ;
enddo ;
enddo
585 call post_data(cs%id_w_sponge, w_int(:,:,:), cs%diag)
597 if (.not.
associated(cs))
return
599 if (
associated(cs%col_i))
deallocate(cs%col_i)
600 if (
associated(cs%col_j))
deallocate(cs%col_j)
602 if (
associated(cs%Iresttime_col))
deallocate(cs%Iresttime_col)
603 if (
associated(cs%Rcv_ml_ref))
deallocate(cs%Rcv_ml_ref)
604 if (
associated(cs%Ref_eta))
deallocate(cs%Ref_eta)
606 if (
associated(cs%Iresttime_im))
deallocate(cs%Iresttime_im)
607 if (
associated(cs%Rcv_ml_ref_im))
deallocate(cs%Rcv_ml_ref_im)
608 if (
associated(cs%Ref_eta_im))
deallocate(cs%Ref_eta_im)
611 if (
associated(cs%Ref_val(cs%fldno)%p))
deallocate(cs%Ref_val(cs%fldno)%p)
612 if (
associated(cs%Ref_val_im(cs%fldno)%p)) &
613 deallocate(cs%Ref_val_im(cs%fldno)%p)