MOM_ALE_sponge.F90

Go to the documentation of this file.

!> This module contains the routines used to apply sponge layers when using
!! the ALE mode.
!! Applying sponges requires the following:
!! (1) initialize_ALE_sponge
!! (2) set_up_ALE_sponge_field (tracers) and set_up_ALE_sponge_vel_field (vel)
!! (3) apply_ALE_sponge
!! (4) init_ALE_sponge_diags (not being used for now)
!! (5) ALE_sponge_end (not being used for now)

module mom_ale_sponge

! This file is part of MOM6. See LICENSE.md for the license.

use mom_coms, only : sum_across_pes
use mom_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field
use mom_diag_mediator, only : diag_ctrl
use mom_error_handler, only : mom_error, fatal, note, warning, is_root_pe
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_grid, only : ocean_grid_type
use mom_spatial_means, only : global_i_mean
use mom_time_manager, only : time_type
use mom_remapping, only : remapping_cs, remapping_core_h, initialize_remapping

! GMM - Planned extension:  Support for time varying sponge targets.

implicit none ; private

#include <MOM_memory.h>

!< Publicly available functions
public set_up_ale_sponge_field, set_up_ale_sponge_vel_field
public initialize_ale_sponge, apply_ale_sponge, ale_sponge_end, init_ale_sponge_diags

type :: p3d
  real, dimension(:,:,:), pointer :: p => null()
end type p3d
type :: p2d
  real, dimension(:,:), pointer :: p => null()
end type p2d

!> SPONGE control structure
type, public :: ale_sponge_cs ; private
  integer :: nz                      !< The total number of layers.
  integer :: nz_data                 !< The total number of arbritary layers.
  integer :: isc, iec, jsc, jec      !< The index ranges of the computational domain at h.
  integer :: iscb, iecb, jscb, jecb  !< The index ranges of the computational domain at u/v.
  integer :: isd, ied, jsd, jed      !< The index ranges of the data domain.

  integer :: num_col, num_col_u, num_col_v  !< The number of sponge points within the
                             !! computational domain.
  integer :: fldno = 0       !< The number of fields which have already been
                             !! registered by calls to set_up_sponge_field
  logical :: sponge_uv      !< Control whether u and v are included in sponge
  integer, pointer :: col_i(:) => null()  !< Arrays containing the i- and j- indicies
  integer, pointer :: col_j(:) => null()  !! of each of the columns being damped.
  integer, pointer :: col_i_u(:) => null()  !< Same as above for u points
  integer, pointer :: col_j_u(:) => null()
  integer, pointer :: col_i_v(:) => null()  !< Same as above for v points
  integer, pointer :: col_j_v(:) => null()

  real, pointer :: iresttime_col(:) => null()  !< The inverse restoring time of
                                               !! each column.
  real, pointer :: iresttime_col_u(:) => null()  !< Same as above for u points
  real, pointer :: iresttime_col_v(:) => null()  !< Same as above for v points

  type(p3d) :: var(max_fields_)      !< Pointers to the fields that are being damped.
  type(p2d) :: ref_val(max_fields_)  !< The values to which the fields are damped.
  real, dimension(:,:), pointer :: ref_val_u => null() !< Same as above for u points
  real, dimension(:,:), pointer :: ref_val_v => null() !< Same as above for v points
  real, dimension(:,:,:),  pointer :: var_u => null() !< Pointers to the u vel. that are
                                                           !! being damped.
  real, dimension(:,:,:),  pointer :: var_v => null() !< Pointers to the v vel. that are
                                                           !! being damped.
  real, pointer :: ref_h(:,:) => null() !< Grid on which reference data is provided
  real, pointer :: ref_hu(:,:) => null() !< Same as above for u points
  real, pointer :: ref_hv(:,:) => null() !< Same as above for v points

  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                                   !! timing of diagnostic output.

  type(remapping_cs) :: remap_cs   !< Remapping parameters and work arrays

end type ale_sponge_cs

contains

!> This subroutine determines the number of points which are within
! sponges in this computational domain.  Only points that have
! positive values of Iresttime and which mask2dT indicates are ocean
! points are included in the sponges.  It also stores the target interface
! heights.
subroutine initialize_ale_sponge(Iresttime, data_h, nz_data, G, param_file, CS)

  integer,                              intent(in) :: nz_data !< The total number of arbritary layers (in).
  type(ocean_grid_type),                intent(in) :: G !< The ocean's grid structure (in).
  real, dimension(SZI_(G),SZJ_(G)),     intent(in) :: Iresttime !< The inverse of the restoring time, in s-1 (in).
  real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_h !< The thickness to be used in the sponge. It has arbritary layers (in).
  type(param_file_type),                intent(in) :: param_file !< A structure indicating the open file to parse for model parameter values (in).
  type(ale_sponge_cs),                  pointer    :: CS !< A pointer that is set to point to the control structure for this module (in/out).


! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_sponge"  ! This module's name.
  logical :: use_sponge
  real, dimension(SZIB_(G),SZJ_(G),nz_data) :: data_hu !< thickness at u points
  real, dimension(SZI_(G),SZJB_(G),nz_data) :: data_hv !< thickness at v points
  real, dimension(SZIB_(G),SZJ_(G)) :: Iresttime_u !< inverse of the restoring time at u points, s-1
  real, dimension(SZI_(G),SZJB_(G)) :: Iresttime_v !< inverse of the restoring time at v points, s-1
  logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries
  integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v
  character(len=10)  :: remapScheme
  if (associated(cs)) then
    call mom_error(warning, "initialize_sponge called with an associated "// &
                            "control structure.")
    return
  endif

! Set default, read and log parameters
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "SPONGE", use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. \n"//&
                 "The exact location and properties of those sponges are \n"//&
                 "specified from MOM_initialization.F90.", default=.false.)

  if (.not.use_sponge) return

  allocate(cs)

  call get_param(param_file, mdl, "SPONGE_UV", cs%sponge_uv, &
                 "Apply sponges in u and v, in addition to tracers.", &
                 default=.false.)

  call get_param(param_file, mdl, "REMAPPING_SCHEME", remapscheme, &
                 "This sets the reconstruction scheme used \n"//&
                 " for vertical remapping for all variables.", &
                 default="PLM", do_not_log=.true.)

  call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndextrapolation, &
                 "When defined, a proper high-order reconstruction \n"//&
                 "scheme is used within boundary cells rather \n"// &
                 "than PCM. E.g., if PPM is used for remapping, a \n" //&
                 "PPM reconstruction will also be used within boundary cells.", &
                 default=.false., do_not_log=.true.)

  cs%nz = g%ke
  cs%isc = g%isc ; cs%iec = g%iec ; cs%jsc = g%jsc ; cs%jec = g%jec
  cs%isd = g%isd ; cs%ied = g%ied ; cs%jsd = g%jsd ; cs%jed = g%jed
  cs%iscB = g%iscB ; cs%iecB = g%iecB; cs%jscB = g%jscB ; cs%jecB = g%jecB

  ! number of columns to be restored
  cs%num_col = 0 ; cs%fldno = 0
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) &
      cs%num_col = cs%num_col + 1
  enddo ; enddo


  if (cs%num_col > 0) then

    allocate(cs%Iresttime_col(cs%num_col)) ; cs%Iresttime_col = 0.0
    allocate(cs%col_i(cs%num_col))         ; cs%col_i = 0
    allocate(cs%col_j(cs%num_col))         ; cs%col_j = 0

    ! pass indices, restoring time to the CS structure
    col = 1
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      if ((iresttime(i,j)>0.0) .and. (g%mask2dT(i,j)>0)) then
        cs%col_i(col) = i ; cs%col_j(col) = j
        cs%Iresttime_col(col) = iresttime(i,j)
        col = col +1
      endif
    enddo ; enddo

    ! same for total number of arbritary layers and correspondent data
    cs%nz_data = nz_data
    allocate(cs%Ref_h(cs%nz_data,cs%num_col))
    do col=1,cs%num_col ; do k=1,cs%nz_data
      cs%Ref_h(k,col) = data_h(cs%col_i(col),cs%col_j(col),k)
    enddo ; enddo

  endif

  total_sponge_cols = cs%num_col
  call sum_across_pes(total_sponge_cols)

! Call the constructor for remapping control structure
  call initialize_remapping(cs%remap_cs, remapscheme, boundary_extrapolation=bndextrapolation)

  call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, &
                 "The total number of columns where sponges are applied at h points.")

  if (cs%sponge_uv) then
     ! u points
     cs%num_col_u = 0 ; !CS%fldno_u = 0
     do j=cs%jsc,cs%jec; do i=cs%iscB,cs%iecB
        data_hu(i,j,:) = 0.5 * (data_h(i,j,:) + data_h(i+1,j,:))
        iresttime_u(i,j) = 0.5 * (iresttime(i,j) + iresttime(i+1,j))
        if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) &
           cs%num_col_u = cs%num_col_u + 1
     enddo; enddo

     if (cs%num_col_u > 0) then

        allocate(cs%Iresttime_col_u(cs%num_col_u)) ; cs%Iresttime_col_u = 0.0
        allocate(cs%col_i_u(cs%num_col_u))         ; cs%col_i_u = 0
        allocate(cs%col_j_u(cs%num_col_u))         ; cs%col_j_u = 0

        ! pass indices, restoring time to the CS structure
        col = 1
        do j=cs%jsc,cs%jec ; do i=cs%iscB,cs%iecB
          if ((iresttime_u(i,j)>0.0) .and. (g%mask2dCu(i,j)>0)) then
            cs%col_i_u(col) = i ; cs%col_j_u(col) = j
            cs%Iresttime_col_u(col) = iresttime_u(i,j)
            col = col +1
          endif
        enddo ; enddo

        ! same for total number of arbritary layers and correspondent data
        allocate(cs%Ref_hu(cs%nz_data,cs%num_col_u))
        do col=1,cs%num_col_u ; do k=1,cs%nz_data
           cs%Ref_hu(k,col) = data_hu(cs%col_i_u(col),cs%col_j_u(col),k)
        enddo ; enddo

     endif
     total_sponge_cols_u = cs%num_col_u
     call sum_across_pes(total_sponge_cols_u)
     call log_param(param_file, mdl, "!Total sponge columns at u points", total_sponge_cols_u, &
                 "The total number of columns where sponges are applied at u points.")

     ! v points
     cs%num_col_v = 0 ; !CS%fldno_v = 0
     do j=cs%jscB,cs%jecB; do i=cs%isc,cs%iec
        data_hu(i,j,:) = 0.5 * (data_h(i,j,:) + data_h(i,j+1,:))
        iresttime_v(i,j) = 0.5 * (iresttime(i,j) + iresttime(i,j+1))
        if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) &
           cs%num_col_v = cs%num_col_v + 1
     enddo; enddo

     if (cs%num_col_v > 0) then

        allocate(cs%Iresttime_col_v(cs%num_col_v)) ; cs%Iresttime_col_v = 0.0
        allocate(cs%col_i_v(cs%num_col_v))         ; cs%col_i_v = 0
        allocate(cs%col_j_v(cs%num_col_v))         ; cs%col_j_v = 0

        ! pass indices, restoring time to the CS structure
        col = 1
        do j=cs%jscB,cs%jecB ; do i=cs%isc,cs%iec
          if ((iresttime_v(i,j)>0.0) .and. (g%mask2dCv(i,j)>0)) then
            cs%col_i_v(col) = i ; cs%col_j_v(col) = j
            cs%Iresttime_col_v(col) = iresttime_v(i,j)
            col = col +1
          endif
        enddo ; enddo

        ! same for total number of arbritary layers and correspondent data
        allocate(cs%Ref_hv(cs%nz_data,cs%num_col_v))
        do col=1,cs%num_col_v ; do k=1,cs%nz_data
           cs%Ref_hv(k,col) = data_hv(cs%col_i_v(col),cs%col_j_v(col),k)
        enddo ; enddo

     endif
     total_sponge_cols_v = cs%num_col_v
     call sum_across_pes(total_sponge_cols_v)
     call log_param(param_file, mdl, "!Total sponge columns at v points", total_sponge_cols_v, &
                 "The total number of columns where sponges are applied at v points.")
  endif

end subroutine initialize_ale_sponge

!> Initialize diagnostics for the ALE_sponge module.
! GMM: this routine is not being used for now.
subroutine init_ale_sponge_diags(Time, G, diag, CS)
  type(time_type),       target, intent(in)    :: Time
  type(ocean_grid_type),         intent(in)    :: G    !< The ocean's grid structure
  type(diag_ctrl),       target, intent(inout) :: diag
  type(ale_sponge_cs),           pointer       :: CS

  if (.not.associated(cs)) return

  cs%diag => diag

end subroutine init_ale_sponge_diags

!> This subroutine stores the reference profile at h points for the variable
! whose address is given by f_ptr.
subroutine set_up_ale_sponge_field(sp_val, G, f_ptr, CS)
  type(ocean_grid_type),                   intent(in) :: G !< Grid structure (in).
  type(ale_sponge_cs),                     pointer  :: CS !< Sponge structure (in/out).
  real, dimension(SZI_(G),SZJ_(G),CS%nz_data), intent(in) :: sp_val !< Field to be used in the sponge, it has arbritary number of layers (in).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target, intent(in) :: f_ptr !< Pointer to the field to be damped (in).

  integer :: j, k, col
  character(len=256) :: mesg ! String for error messages

  if (.not.associated(cs)) return

  cs%fldno = cs%fldno + 1

  if (cs%fldno > max_fields_) then
    write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease &
           &the number of fields to be damped in the call to &
           &initialize_sponge." )') cs%fldno
    call mom_error(fatal,"set_up_ALE_sponge_field: "//mesg)
  endif

  ! stores the reference profile
  allocate(cs%Ref_val(cs%fldno)%p(cs%nz_data,cs%num_col))
  cs%Ref_val(cs%fldno)%p(:,:) = 0.0
  do col=1,cs%num_col
    do k=1,cs%nz_data
      cs%Ref_val(cs%fldno)%p(k,col) = sp_val(cs%col_i(col),cs%col_j(col),k)
    enddo
  enddo

  cs%var(cs%fldno)%p => f_ptr

end subroutine set_up_ale_sponge_field

!> This subroutine stores the reference profile at uand v points for the variable
! whose address is given by u_ptr and v_ptr.
subroutine set_up_ale_sponge_vel_field(u_val, v_val, G, u_ptr, v_ptr, CS)
  type(ocean_grid_type),                   intent(in) :: G !< Grid structure (in).
  type(ale_sponge_cs),                     pointer  :: CS !< Sponge structure (in/out).
  real, dimension(SZIB_(G),SZJ_(G),CS%nz_data), intent(in) :: u_val !< u field to be used in the sponge, it has arbritary number of layers (in).
  real, dimension(SZI_(G),SZJB_(G),CS%nz_data), intent(in) :: v_val !< u field to be used in the sponge, it has arbritary number of layers (in).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(in) :: u_ptr !< u pointer to the field to be damped (in).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(in) :: v_ptr !< v pointer to the field to be damped (in).

  integer :: j, k, col
  character(len=256) :: mesg ! String for error messages

  if (.not.associated(cs)) return

  ! stores the reference profile
  allocate(cs%Ref_val_u(cs%nz_data,cs%num_col_u))
  cs%Ref_val_u(:,:) = 0.0
  do col=1,cs%num_col_u
    do k=1,cs%nz_data
      cs%Ref_val_u(k,col) = u_val(cs%col_i_u(col),cs%col_j_u(col),k)
    enddo
  enddo

  cs%var_u => u_ptr

  allocate(cs%Ref_val_v(cs%nz_data,cs%num_col_v))
  cs%Ref_val_v(:,:) = 0.0
  do col=1,cs%num_col_v
    do k=1,cs%nz_data
      cs%Ref_val_v(k,col) = v_val(cs%col_i_v(col),cs%col_j_v(col),k)
    enddo
  enddo

  cs%var_v => v_ptr

end subroutine set_up_ale_sponge_vel_field

!> This subroutine applies damping to the layers thicknesses, temp, salt and a variety of tracers for every column where there is damping.
subroutine apply_ale_sponge(h, dt, G, CS)
  type(ocean_grid_type),                    intent(inout) :: G !< The ocean's grid structure (in).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness, in m (in)
  real,                                     intent(in)    :: dt !< The amount of time covered by this call, in s (in).
  type(ale_sponge_cs),                      pointer       :: CS !<A pointer to the control structure for this module that is set by a previous call to initialize_sponge (in).

  real :: damp                                               !< The timestep times the local damping  coefficient.  ND.
  real :: I1pdamp                                            !< I1pdamp is 1/(1 + damp).  Nondimensional.
  real :: Idt                                                !< 1.0/dt, in s-1.
  real :: tmp_val1(cs%nz)                                    !< data values remapped to model grid
  real :: tmp_val2(cs%nz_data)                               ! data values remapped to model grid
  real :: hu(szib_(g), szj_(g), szk_(g))                     !> A temporary array for h at u pts
  real :: hv(szi_(g), szjb_(g), szk_(g))                     !> A temporary array for h at v pts
  integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  if (.not.associated(cs)) return

  do c=1,cs%num_col
! c is an index for the next 3 lines but a multiplier for the rest of the loop
! Therefore we use c as per C code and increment the index where necessary.
    i = cs%col_i(c) ; j = cs%col_j(c)
    damp = dt*cs%Iresttime_col(c)
    i1pdamp = 1.0 / (1.0 + damp)
        do m=1,cs%fldno
            tmp_val2(1:cs%nz_data) = cs%Ref_val(m)%p(1:cs%nz_data,c)
            call remapping_core_h(cs%remap_cs, &
                  cs%nz_data, cs%Ref_h(:,c), tmp_val2, &
                  cs%nz, h(i,j,:), tmp_val1)

            !Backward Euler method
            cs%var(m)%p(i,j,:) = i1pdamp * &
                                 (cs%var(m)%p(i,j,:) + tmp_val1 * damp)

!            CS%var(m)%p(i,j,k) = I1pdamp * &
!                                 (CS%var(m)%p(i,j,k) + tmp_val1 * damp)
        enddo

  enddo ! end of c loop
  ! for debugging
  !c=CS%num_col
  !do m=1,CS%fldno
  !   write(*,*)'APPLY SPONGE,m,CS%Ref_h(:,c),h(i,j,:),tmp_val2,tmp_val1',m,CS%Ref_h(:,c),h(i,j,:),tmp_val2,tmp_val1
  !enddo

  if (cs%sponge_uv) then

    ! u points
    do j=cs%jsc,cs%jec; do i=cs%iscB,cs%iecB; do k=1,nz
       hu(i,j,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
    enddo; enddo; enddo

    do c=1,cs%num_col_u
       i = cs%col_i_u(c) ; j = cs%col_j_u(c)
       damp = dt*cs%Iresttime_col_u(c)
       i1pdamp = 1.0 / (1.0 + damp)
       tmp_val2(1:cs%nz_data) = cs%Ref_val_u(1:cs%nz_data,c)
       call remapping_core_h(cs%remap_cs, &
                  cs%nz_data, cs%Ref_hu(:,c), tmp_val2, &
                  cs%nz, hu(i,j,:), tmp_val1)

       !Backward Euler method
       cs%var_u(i,j,:) = i1pdamp * (cs%var_u(i,j,:) + tmp_val1 * damp)
    enddo

    ! v points
    do j=cs%jscB,cs%jecB; do i=cs%isc,cs%iec; do k=1,nz
       hv(i,j,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
    enddo; enddo; enddo

    do c=1,cs%num_col_v
       i = cs%col_i_v(c) ; j = cs%col_j_v(c)
       damp = dt*cs%Iresttime_col_v(c)
       i1pdamp = 1.0 / (1.0 + damp)
       tmp_val2(1:cs%nz_data) = cs%Ref_val_v(1:cs%nz_data,c)
       call remapping_core_h(cs%remap_cs, &
                  cs%nz_data, cs%Ref_hv(:,c), tmp_val2, &
                  cs%nz, hv(i,j,:), tmp_val1)
       !Backward Euler method
       cs%var_v(i,j,:) = i1pdamp * (cs%var_v(i,j,:) + tmp_val1 * damp)
    enddo

  endif
end subroutine apply_ale_sponge

!> GMM: I could not find where sponge_end is being called, but I am keeping
!  ALE_sponge_end here so we can add that if needed.
subroutine ale_sponge_end(CS)
  type(ale_sponge_cs),              pointer       :: CS
!  (in)      CS - A pointer to the control structure for this module that is
!                 set by a previous call to initialize_sponge.
  integer :: m

  if (.not.associated(cs)) return

  if (associated(cs%col_i)) deallocate(cs%col_i)
  if (associated(cs%col_i_u)) deallocate(cs%col_i_u)
  if (associated(cs%col_i_v)) deallocate(cs%col_i_v)
  if (associated(cs%col_j)) deallocate(cs%col_j)
  if (associated(cs%col_j_u)) deallocate(cs%col_j_u)
  if (associated(cs%col_j_v)) deallocate(cs%col_j_v)

  if (associated(cs%Iresttime_col)) deallocate(cs%Iresttime_col)
  if (associated(cs%Iresttime_col_u)) deallocate(cs%Iresttime_col_u)
  if (associated(cs%Iresttime_col_v)) deallocate(cs%Iresttime_col_v)

  do m=1,cs%fldno
    if (associated(cs%Ref_val(cs%fldno)%p)) deallocate(cs%Ref_val(cs%fldno)%p)
  enddo

  deallocate(cs)

end subroutine ale_sponge_end
end module mom_ale_sponge