Set up sponges in 2d DOME configuration.
357 type(ocean_grid_type),
intent(in) :: G
358 type(verticalGrid_type),
intent(in) :: GV
359 type(thermo_var_ptrs),
intent(in) :: tv
360 type(param_file_type),
intent(in) :: param_file
361 logical,
intent(in) :: use_ALE
362 type(sponge_CS),
pointer :: CSp
363 type(ALE_sponge_CS),
pointer :: ACSp
365 real :: T(SZI_(G),SZJ_(G),SZK_(G))
366 real :: S(SZI_(G),SZJ_(G),SZK_(G))
367 real :: RHO(SZI_(G),SZJ_(G),SZK_(G))
368 real :: h(SZI_(G),SZJ_(G),SZK_(G))
369 real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1)
370 real :: Idamp(SZI_(G),SZJ_(G))
372 real :: S_range, T_range
373 real :: e0(SZK_(G)+1)
375 real :: eta1D(SZK_(G)+1)
377 real :: d_eta(SZK_(G))
378 real :: dome2d_width_bay, dome2d_width_bottom, dome2d_depth_bay
379 real :: dome2d_west_sponge_time_scale, dome2d_east_sponge_time_scale
380 real :: dome2d_west_sponge_width, dome2d_east_sponge_width
382 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
384 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
385 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
387 call get_param(param_file, mdl,
"DOME2D_WEST_SPONGE_TIME_SCALE", dome2d_west_sponge_time_scale, &
388 'The time-scale on the west edge of the domain for restoring T/S '//&
389 'in the sponge. If zero, the western sponge is disabled', &
390 units=
's', default=0.)
391 call get_param(param_file, mdl,
"DOME2D_EAST_SPONGE_TIME_SCALE", dome2d_east_sponge_time_scale, &
392 'The time-scale on the east edge of the domain for restoring T/S '//&
393 'in the sponge. If zero, the eastern sponge is disabled', &
394 units=
's', default=0.)
395 call get_param(param_file, mdl,
"DOME2D_WEST_SPONGE_WIDTH", dome2d_west_sponge_width, &
396 'The fraction of the domain in which the western sponge for restoring T/S '//&
398 units=
'nondim', default=0.1)
399 call get_param(param_file, mdl,
"DOME2D_EAST_SPONGE_WIDTH", dome2d_east_sponge_width, &
400 'The fraction of the domain in which the eastern sponge for restoring T/S '//&
402 units=
'nondim', default=0.1)
405 if (dome2d_west_sponge_time_scale <= 0. .and. dome2d_east_sponge_time_scale <= 0.)
return
407 if (
associated(csp))
call mom_error(fatal, &
408 "DOME2d_initialize_sponges called with an associated control structure.")
409 if (
associated(acsp))
call mom_error(fatal, &
410 "DOME2d_initialize_sponges called with an associated ALE-sponge control structure.")
412 call get_param(param_file, mdl,
"DOME2D_SHELF_WIDTH", dome2d_width_bay, &
413 default=0.1, do_not_log=.true.)
414 call get_param(param_file, mdl,
"DOME2D_BASIN_WIDTH", dome2d_width_bottom, &
415 default=0.3, do_not_log=.true.)
416 call get_param(param_file, mdl,
"DOME2D_SHELF_DEPTH", dome2d_depth_bay, &
417 default=0.2, do_not_log=.true.)
418 call get_param(param_file, mdl,
"S_REF",s_ref)
419 call get_param(param_file, mdl,
"T_REF",t_ref)
420 call get_param(param_file, mdl,
"S_RANGE",s_range,default=2.0)
421 call get_param(param_file, mdl,
"T_RANGE",t_range,default=0.0)
426 do j=js,je ;
do i=is,ie
427 if (g%mask2dT(i,j) > 0.)
then
428 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
429 if ( dome2d_west_sponge_time_scale > 0. .and. x < dome2d_west_sponge_width )
then
431 dummy1 = 1. - x / dome2d_west_sponge_width
432 idamp(i,j) = 1./dome2d_west_sponge_time_scale * max(0., min(1., dummy1))
433 elseif ( dome2d_east_sponge_time_scale > 0. .and. x > ( 1. - dome2d_east_sponge_width ) )
then
435 dummy1 = 1. - ( 1. - x ) / dome2d_east_sponge_width
436 idamp(i,j) = 1./dome2d_east_sponge_time_scale * max(0., min(1., dummy1))
450 e0(k) = -g%max_depth * ( real(k-1) / real(nz) )
452 e0(nz+1) = -g%max_depth
453 do j=js,je ;
do i=is,ie
454 eta1d(nz+1) = -g%bathyT(i,j)
457 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
458 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
459 h(i,j,k) = gv%Angstrom_H
461 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
466 call initialize_ale_sponge(idamp, g, param_file, acsp, h, nz)
469 t(:,:,:) = 0.0 ; s(:,:,:) = 0.0
470 do j=js,je ;
do i=is,ie
473 z = z + 0.5 * gv%H_to_Z * h(i,j,k)
474 s(i,j,k) = 34.0 - 1.0 * (z/g%max_depth)
475 if ( ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon < dome2d_west_sponge_width ) &
476 s(i,j,k) = s_ref + s_range
477 z = z + 0.5 * gv%H_to_Z * h(i,j,k)
481 if (
associated(tv%T) )
then
482 call set_up_ale_sponge_field(t, g, tv%T, acsp)
484 if (
associated(tv%S) )
then
485 call set_up_ale_sponge_field(s, g, tv%S, acsp)
491 do j=js,je ;
do i=is,ie
492 eta1d(nz+1) = -g%bathyT(i,j)
494 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
495 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
496 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
497 d_eta(k) = gv%Angstrom_Z
499 d_eta(k) = (eta1d(k) - eta1d(k+1))
503 x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon
504 if ( x <= dome2d_width_bay )
then
505 do k=1,nz-1 ; d_eta(k) = gv%Angstrom_Z ;
enddo
506 d_eta(nz) = dome2d_depth_bay * g%max_depth - (nz-1) * gv%Angstrom_Z
509 eta(i,j,nz+1) = -g%bathyT(i,j)
511 eta(i,j,k) = eta(i,j,k+1) + d_eta(k)
514 call initialize_sponge(idamp, eta, g, param_file, csp, gv)