18 implicit none ;
private
20 #include <MOM_memory.h>
22 public set_up_sponge_field, set_up_sponge_ml_density
23 public initialize_sponge, apply_sponge, sponge_end, init_sponge_diags
32 real,
dimension(:,:,:),
pointer :: p => null()
36 real,
dimension(:,:),
pointer :: p => null()
41 logical :: bulkmixedlayer
55 integer,
pointer :: col_i(:) => null()
56 integer,
pointer :: col_j(:) => null()
57 real,
pointer :: iresttime_col(:) => null()
58 real,
pointer :: rcv_ml_ref(:) => null()
60 real,
pointer :: ref_eta(:,:) => null()
62 type(
p3d) :: var(max_fields_)
63 type(
p2d) :: ref_val(max_fields_)
65 logical :: do_i_mean_sponge
66 real,
pointer :: iresttime_im(:) => null()
68 real,
pointer :: rcv_ml_ref_im(:) => null()
70 real,
pointer :: ref_eta_im(:,:) => null()
72 type(
p2d) :: ref_val_im(max_fields_)
77 integer :: id_w_sponge = -1
87 subroutine initialize_sponge(Iresttime, int_height, G, param_file, CS, GV, &
88 Iresttime_i_mean, int_height_i_mean)
90 real,
dimension(SZI_(G),SZJ_(G)), &
91 intent(in) :: iresttime
92 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
93 intent(in) :: int_height
98 real,
dimension(SZJ_(G)), &
99 optional,
intent(in) :: iresttime_i_mean
101 real,
dimension(SZJ_(G),SZK_(G)+1), &
102 optional,
intent(in) :: int_height_i_mean
107 #include "version_variable.h"
108 character(len=40) :: mdl =
"MOM_sponge"
109 logical :: use_sponge
110 integer :: i, j, k, col, total_sponge_cols
112 if (
associated(cs))
then
113 call mom_error(warning,
"initialize_sponge called with an associated "// &
114 "control structure.")
120 call get_param(param_file, mdl,
"SPONGE", use_sponge, &
121 "If true, sponges may be applied anywhere in the domain. "//&
122 "The exact location and properties of those sponges are "//&
123 "specified from MOM_initialization.F90.", default=.false.)
125 if (.not.use_sponge)
return
128 if (
present(iresttime_i_mean) .neqv.
present(int_height_i_mean)) &
129 call mom_error(fatal,
"initialize_sponge: The optional arguments \n"//&
130 "Iresttime_i_mean and int_height_i_mean must both be present \n"//&
133 cs%do_i_mean_sponge =
present(iresttime_i_mean)
139 cs%bulkmixedlayer = .false.
141 cs%num_col = 0 ; cs%fldno = 0
142 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
143 if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) &
144 cs%num_col = cs%num_col + 1
147 if (cs%num_col > 0)
then
149 allocate(cs%Iresttime_col(cs%num_col)) ; cs%Iresttime_col = 0.0
150 allocate(cs%col_i(cs%num_col)) ; cs%col_i = 0
151 allocate(cs%col_j(cs%num_col)) ; cs%col_j = 0
154 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
155 if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0))
then
156 cs%col_i(col) = i ; cs%col_j(col) = j
157 cs%Iresttime_col(col) = iresttime(i,j)
162 allocate(cs%Ref_eta(cs%nz+1,cs%num_col))
163 do col=1,cs%num_col ;
do k=1,cs%nz+1
164 cs%Ref_eta(k,col) = int_height(cs%col_i(col),cs%col_j(col),k)
169 if (cs%do_i_mean_sponge)
then
170 allocate(cs%Iresttime_im(g%jsd:g%jed)) ; cs%Iresttime_im(:) = 0.0
171 allocate(cs%Ref_eta_im(g%jsd:g%jed,g%ke+1)) ; cs%Ref_eta_im(:,:) = 0.0
174 cs%Iresttime_im(j) = iresttime_i_mean(j)
176 do k=1,cs%nz+1 ;
do j=g%jsc,g%jec
177 cs%Ref_eta_im(j,k) = int_height_i_mean(j,k)
181 total_sponge_cols = cs%num_col
182 call sum_across_pes(total_sponge_cols)
184 call log_param(param_file, mdl,
"!Total sponge columns", total_sponge_cols, &
185 "The total number of columns where sponges are applied.")
187 end subroutine initialize_sponge
192 subroutine init_sponge_diags(Time, G, diag, CS)
193 type(time_type),
target,
intent(in) :: time
195 type(
diag_ctrl),
target,
intent(inout) :: diag
199 if (.not.
associated(cs))
return
202 cs%id_w_sponge = register_diag_field(
'ocean_model',
'w_sponge', diag%axesTi, &
203 time,
'The diapycnal motion due to the sponges',
'm s-1')
205 end subroutine init_sponge_diags
210 subroutine set_up_sponge_field(sp_val, f_ptr, G, nlay, CS, sp_val_i_mean)
212 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
214 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
215 target,
intent(in) :: f_ptr
216 integer,
intent(in) :: nlay
219 real,
dimension(SZJ_(G),SZK_(G)),&
220 optional,
intent(in) :: sp_val_i_mean
224 character(len=256) :: mesg
226 if (.not.
associated(cs))
return
228 cs%fldno = cs%fldno + 1
230 if (cs%fldno > max_fields_)
then
231 write(mesg,
'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease &
232 &the number of fields to be damped in the call to &
233 &initialize_sponge." )') cs%fldno
234 call mom_error(fatal,
"set_up_sponge_field: "//mesg)
237 allocate(cs%Ref_val(cs%fldno)%p(cs%nz,cs%num_col))
238 cs%Ref_val(cs%fldno)%p(:,:) = 0.0
241 cs%Ref_val(cs%fldno)%p(k,col) = sp_val(cs%col_i(col),cs%col_j(col),k)
244 cs%Ref_val(cs%fldno)%p(k,col) = 0.0
248 cs%var(cs%fldno)%p => f_ptr
250 if (nlay/=cs%nz)
then
251 write(mesg,
'("Danger: Sponge reference fields require nz (",I3,") layers.&
252 & A field with ",I3," layers was passed to set_up_sponge_field.")') &
254 if (is_root_pe())
call mom_error(warning,
"set_up_sponge_field: "//mesg)
257 if (cs%do_i_mean_sponge)
then
258 if (.not.
present(sp_val_i_mean))
call mom_error(fatal, &
259 "set_up_sponge_field: sp_val_i_mean must be present with i-mean sponges.")
261 allocate(cs%Ref_val_im(cs%fldno)%p(cs%jsd:cs%jed,cs%nz))
262 cs%Ref_val(cs%fldno)%p(:,:) = 0.0
263 do k=1,cs%nz ;
do j=cs%jsc,cs%jec
264 cs%Ref_val_im(cs%fldno)%p(j,k) = sp_val_i_mean(j,k)
268 end subroutine set_up_sponge_field
273 subroutine set_up_sponge_ml_density(sp_val, G, CS, sp_val_i_mean)
275 real,
dimension(SZI_(G),SZJ_(G)), &
279 real,
dimension(SZJ_(G)), &
280 optional,
intent(in) :: sp_val_i_mean
291 character(len=256) :: mesg
293 if (.not.
associated(cs))
return
295 if (
associated(cs%Rcv_ml_ref))
then
296 call mom_error(fatal,
"set_up_sponge_ML_density appears to have been "//&
300 cs%bulkmixedlayer = .true.
301 allocate(cs%Rcv_ml_ref(cs%num_col)) ; cs%Rcv_ml_ref(:) = 0.0
303 cs%Rcv_ml_ref(col) = sp_val(cs%col_i(col),cs%col_j(col))
306 if (cs%do_i_mean_sponge)
then
307 if (.not.
present(sp_val_i_mean))
call mom_error(fatal, &
308 "set_up_sponge_field: sp_val_i_mean must be present with i-mean sponges.")
310 allocate(cs%Rcv_ml_ref_im(cs%jsd:cs%jed)) ; cs%Rcv_ml_ref_im(:) = 0.0
312 cs%Rcv_ml_ref_im(j) = sp_val_i_mean(j)
316 end subroutine set_up_sponge_ml_density
320 subroutine apply_sponge(h, dt, G, GV, ea, eb, CS, Rcv_ml)
323 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
325 real,
intent(in) :: dt
326 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
330 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
336 real,
dimension(SZI_(G),SZJ_(G)), &
337 optional,
intent(inout) :: rcv_ml
344 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
349 real,
dimension(SZI_(G), SZJ_(G)) :: &
354 real,
dimension(SZJ_(G), SZK_(G)+1) :: &
356 real,
allocatable,
dimension(:,:,:) :: &
358 real,
dimension(SZI_(G), SZK_(G)+1) :: &
361 real,
dimension(SZI_(G)) :: &
380 integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz
381 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
383 if (.not.
associated(cs))
return
384 if (cs%bulkmixedlayer) nkmb = gv%nk_rho_varies
385 if (cs%bulkmixedlayer .and. (.not.
present(rcv_ml))) &
386 call mom_error(fatal,
"Rml must be provided to apply_sponge when using "//&
387 "a bulk mixed layer.")
389 if ((cs%id_w_sponge > 0) .or. cs%do_i_mean_sponge)
then
390 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
392 enddo ;
enddo ;
enddo
395 if (cs%do_i_mean_sponge)
then
398 if (cs%bulkmixedlayer)
call mom_error(fatal,
"apply_sponge is not yet set up to "//&
399 "work properly with i-mean sponges and a bulk mixed layer.")
401 do j=js,je ;
do i=is,ie ; e_d(i,j,nz+1) = -g%bathyT(i,j) ;
enddo ;
enddo
402 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
403 e_d(i,j,k) = e_d(i,j,k+1) + h(i,j,k)*gv%H_to_Z
404 enddo ;
enddo ;
enddo
407 dilate(i) = g%bathyT(i,j) / (e_d(i,j,1) + g%bathyT(i,j))
409 do k=1,nz+1 ;
do i=is,ie
410 e_d(i,j,k) = dilate(i) * (e_d(i,j,k) + g%bathyT(i,j)) - g%bathyT(i,j)
415 do j=js,je ;
do i=is,ie
416 eta_anom(i,j) = e_d(i,j,k) - cs%Ref_eta_im(j,k)
417 if (cs%Ref_eta_im(j,k) < -g%bathyT(i,j)) eta_anom(i,j) = 0.0
419 call global_i_mean(eta_anom(:,:), eta_mean_anom(:,k), g)
422 if (cs%fldno > 0)
allocate(fld_mean_anom(g%isd:g%ied,nz,cs%fldno))
424 do j=js,je ;
do i=is,ie
425 fld_anom(i,j) = cs%var(m)%p(i,j,k) - cs%Ref_val_im(m)%p(j,k)
427 call global_i_mean(fld_anom(:,:), fld_mean_anom(:,k,m), g, h(:,:,k))
430 do j=js,je ;
if (cs%Iresttime_im(j) > 0.0)
then
431 damp = dt*cs%Iresttime_im(j) ; damp_1pdamp = damp / (1.0 + damp)
434 h_above(i,1) = 0.0 ; h_below(i,nz+1) = 0.0
436 do k=nz,1,-1 ;
do i=is,ie
437 h_below(i,k) = h_below(i,k+1) + max(h(i,j,k)-gv%Angstrom_H, 0.0)
439 do k=2,nz+1 ;
do i=is,ie
440 h_above(i,k) = h_above(i,k-1) + max(h(i,j,k-1)-gv%Angstrom_H, 0.0)
445 w = damp_1pdamp * eta_mean_anom(j,k) * gv%Z_to_H
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)
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)
456 do k=1,nz ;
do i=is,ie
457 ea_k = max(0.0, -w_int(i,j,k))
458 eb_k = max(0.0, w_int(i,j,k+1))
460 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
461 cs%Ref_val_im(m)%p(j,k) * (ea_k + eb_k)) / &
462 (h(i,j,k) + (ea_k + eb_k)) - &
463 damp_1pdamp * fld_mean_anom(j,k,m)
466 h(i,j,k) = max(h(i,j,k) + (w_int(i,j,k+1) - w_int(i,j,k)), &
467 min(h(i,j,k), gv%Angstrom_H))
471 if (cs%fldno > 0)
deallocate(fld_mean_anom)
478 i = cs%col_i(c) ; j = cs%col_j(c)
479 damp = dt*cs%Iresttime_col(c)
481 e(1) = 0.0 ; e0 = 0.0
483 e(k+1) = e(k) - h(i,j,k)*gv%H_to_Z
485 e_str = e(nz+1) / cs%Ref_eta(nz+1,c)
487 if ( cs%bulkmixedlayer )
then
488 i1pdamp = 1.0 / (1.0 + damp)
489 if (
associated(cs%Rcv_ml_ref)) &
490 rcv_ml(i,j) = i1pdamp * (rcv_ml(i,j) + cs%Rcv_ml_ref(c)*damp)
493 cs%var(m)%p(i,j,k) = i1pdamp * &
494 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(k,c)*damp)
500 if (gv%Rlay(k) > rcv_ml(i,j))
then
501 w = min((((e(k)-e0) - e_str*cs%Ref_eta(k,c)) * damp)*gv%Z_to_H, &
502 ((wb + h(i,j,k)) - gv%Angstrom_H))
505 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
506 cs%Ref_val(m)%p(k,c)*(damp*h(i,j,k) + (wpb - wm))) / &
507 (h(i,j,k)*(1.0 + damp) + (wpb - wm))
511 cs%var(m)%p(i,j,k) = i1pdamp * &
512 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(k,c)*damp)
514 w = wb + (h(i,j,k) - gv%Angstrom_H)
517 eb(i,j,k) = eb(i,j,k) + wpb
518 ea(i,j,k) = ea(i,j,k) - wm
519 h(i,j,k) = h(i,j,k) + (wb - w)
526 w = min((wb + (h(i,j,k) - gv%Angstrom_H)),0.0)
527 h(i,j,k) = h(i,j,k) + (wb - w)
528 ea(i,j,k) = ea(i,j,k) - w
534 eb(i,j,k) = eb(i,j,k) + w
538 h(i,j,k) = h(i,j,k) + w
540 cs%var(m)%p(i,j,k) = (cs%var(m)%p(i,j,k)*h(i,j,k) + &
541 cs%Ref_val(m)%p(k,c)*w) / (h(i,j,k) + w)
547 cs%var(m)%p(i,j,k) = i1pdamp * &
548 (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(gv%nkml,c)*damp)
557 w = min((((e(k)-e0) - e_str*cs%Ref_eta(k,c)) * damp)*gv%Z_to_H, &
558 ((wb + h(i,j,k)) - gv%Angstrom_H))
559 wm = 0.5*(w - abs(w))
561 cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
562 cs%Ref_val(m)%p(k,c) * (damp*h(i,j,k) + (wpb - wm))) / &
563 (h(i,j,k)*(1.0 + damp) + (wpb - wm))
565 eb(i,j,k) = eb(i,j,k) + wpb
566 ea(i,j,k) = ea(i,j,k) - wm
567 h(i,j,k) = h(i,j,k) + (wb - w)
575 if (
associated(cs%diag))
then ;
if (query_averaging_enabled(cs%diag))
then
576 if (cs%id_w_sponge > 0)
then
578 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
579 w_int(i,j,k) = w_int(i,j,k) * idt
580 enddo ;
enddo ;
enddo
581 call post_data(cs%id_w_sponge, w_int(:,:,:), cs%diag)
585 end subroutine apply_sponge
588 subroutine sponge_end(CS)
593 if (.not.
associated(cs))
return
595 if (
associated(cs%col_i))
deallocate(cs%col_i)
596 if (
associated(cs%col_j))
deallocate(cs%col_j)
598 if (
associated(cs%Iresttime_col))
deallocate(cs%Iresttime_col)
599 if (
associated(cs%Rcv_ml_ref))
deallocate(cs%Rcv_ml_ref)
600 if (
associated(cs%Ref_eta))
deallocate(cs%Ref_eta)
602 if (
associated(cs%Iresttime_im))
deallocate(cs%Iresttime_im)
603 if (
associated(cs%Rcv_ml_ref_im))
deallocate(cs%Rcv_ml_ref_im)
604 if (
associated(cs%Ref_eta_im))
deallocate(cs%Ref_eta_im)
607 if (
associated(cs%Ref_val(cs%fldno)%p))
deallocate(cs%Ref_val(cs%fldno)%p)
608 if (
associated(cs%Ref_val_im(cs%fldno)%p)) &
609 deallocate(cs%Ref_val_im(cs%fldno)%p)
614 end subroutine sponge_end