MOM_sponge.F90

Go to the documentation of this file.

module mom_sponge
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, March 1999-June 2000                           *
!*                                                                     *
!*    This program contains the subroutines that implement sponge      *
!*  regions, in which the stratification and water mass properties     *
!*  are damped toward some profiles.  There are three externally       *
!*  callable subroutines in this file.                                 *
!*                                                                     *
!*    initialize_sponge determines the mapping from the model          *
!*  variables into the arrays of damped columns.  This remapping is    *
!*  done for efficiency and to conserve memory.  Only columns which    *
!*  have positive inverse damping times and which are deeper than a    *
!*  supplied depth are placed in sponges.  The inverse damping         *
!*  time is also stored in this subroutine, and memory is allocated    *
!*  for all of the reference profiles which will subsequently be       *
!*  provided through calls to set_up_sponge_field.  The first two      *
!*  arguments are a two-dimensional array containing the damping       *
!*  rates, and the interface heights to damp towards.                  *
!*                                                                     *
!*    set_up_sponge_field is called to provide a reference profile     *
!*  and the location of the field that will be damped back toward      *
!*  that reference profile.  A third argument, the number of layers    *
!*  in the field is also provided, but this should always be nz.       *
!*                                                                     *
!*    Apply_sponge damps all of the fields that have been registered   *
!*  with set_up_sponge_field toward their reference profiles.  The     *
!*  four arguments are the thickness to be damped, the amount of time  *
!*  over which the damping occurs, and arrays to which the movement    *
!*  of fluid into a layer from above and below will be added. The      *
!*  effect on momentum of the sponge may be accounted for later using  *
!*  the movement of water recorded in these later arrays.              *
!*                                                                     *
!*    All of the variables operated upon in this file are defined at   *
!*  the thickness points.                                              *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h. *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, T, S, Iresttime, ea, eb               *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

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_verticalgrid, only : verticalgrid_type

! Planned extension:  Support for time varying sponge targets.

implicit none ; private

#include <MOM_memory.h>

public set_up_sponge_field, set_up_sponge_ml_density
public initialize_sponge, apply_sponge, sponge_end, init_sponge_diags

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

type, public :: sponge_cs ; private
  logical :: bulkmixedlayer  ! If true, a refined bulk mixed layer is used with
                             ! nkml sublayers and nkbl buffer layer.
  integer :: nz              ! The total number of layers.
  integer :: isc, iec, jsc, jec  ! The index ranges of the computational domain.
  integer :: isd, ied, jsd, jed  ! The index ranges of the data domain.
  integer :: num_col         ! 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
  integer, pointer :: col_i(:) => null()  ! Arrays containing the i- and j- indicies
  integer, pointer :: col_j(:) => null()  ! of each of the columns being damped.
  real, pointer :: iresttime_col(:) => null()  ! The inverse restoring time of
                             ! each column.
  real, pointer :: rcv_ml_ref(:) => null() ! The value toward which the mixed layer
                             ! coordinate-density is being damped, in kg m-3.
  real, pointer :: ref_eta(:,:) => null() ! The value toward which the interface
                             ! heights are being damped.
  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.

  logical :: do_i_mean_sponge     ! If true, apply sponges to the i-mean fields.
  real, pointer :: iresttime_im(:) => null()  ! The inverse restoring time of
                             ! each row for i-mean sponges.
  real, pointer :: rcv_ml_ref_im(:) => null() ! The value toward which the i-mean
                             ! mixed layer coordinate-density is being damped,
                             ! in kg m-3.
  real, pointer :: ref_eta_im(:,:) => null() ! The value toward which the i-mean
                             ! interface heights are being damped.
  type(p2d) :: ref_val_im(max_fields_)  ! The values toward which the i-means of
                             ! fields are damped.

  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                             ! timing of diagnostic output.
  integer :: id_w_sponge = -1

end type sponge_cs

contains

subroutine initialize_sponge(Iresttime, int_height, G, param_file, CS, &
                             Iresttime_i_mean, int_height_i_mean)
  type(ocean_grid_type),                  intent(in) :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G)),       intent(in) :: Iresttime
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: int_height
  type(param_file_type),                  intent(in) :: param_file !< A structure to parse for run-time parameters
  type(sponge_cs),                        pointer    :: CS
  real, dimension(SZJ_(G)),     optional, intent(in) :: Iresttime_i_mean
  real, dimension(SZJ_(G),SZK_(G)+1), optional, intent(in) :: int_height_i_mean

!   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.

! Arguments: Iresttime - The inverse of the restoring time, in s-1.
!  (in)      int_height - The interface heights to damp back toward, in m.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
! 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
  integer :: i, j, k, col, total_sponge_cols

  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)

  if (present(iresttime_i_mean) .neqv. present(int_height_i_mean)) &
    call mom_error(fatal, "initialize_sponge:  The optional arguments \n"//&
           "Iresttime_i_mean and int_height_i_mean must both be present \n"//&
           "if either one is.")

  cs%do_i_mean_sponge = present(iresttime_i_mean)

  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%bulkmixedlayer may be set later via a call to set_up_sponge_ML_density.
  cs%bulkmixedlayer = .false.

  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

    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

    allocate(cs%Ref_eta(cs%nz+1,cs%num_col))
    do col=1,cs%num_col ; do k=1,cs%nz+1
      cs%Ref_eta(k,col) = int_height(cs%col_i(col),cs%col_j(col),k)
    enddo ; enddo

  endif

  if (cs%do_i_mean_sponge) then
    allocate(cs%Iresttime_im(g%jsd:g%jed)) ; cs%Iresttime_im(:) = 0.0
    allocate(cs%Ref_eta_im(g%jsd:g%jed,g%ke+1)) ; cs%Ref_eta_im(:,:) = 0.0

    do j=g%jsc,g%jec
      cs%Iresttime_im(j) = iresttime_i_mean(j)
    enddo
    do k=1,cs%nz+1 ; do j=g%jsc,g%jec
      cs%Ref_eta_im(j,k) = int_height_i_mean(j,k)
    enddo ; enddo
  endif

  total_sponge_cols = cs%num_col
  call sum_across_pes(total_sponge_cols)

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

end subroutine initialize_sponge

subroutine init_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(sponge_cs),               pointer       :: CS

!   This subroutine sets up diagnostics for the sponges.  It is separate
! from initialize_sponge because it requires fields that are not readily
! availble where initialize_sponge is called.

! Arguments:  Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (in/out)  CS - A pointer to the control structure for this module that is
!                 set by a previous call to initialize_sponge.

  if (.not.associated(cs)) return

  cs%diag => diag
  cs%id_w_sponge = register_diag_field('ocean_model', 'w_sponge', diag%axesTi, &
      time, 'The diapycnal motion due to the sponges', 'meter second-1')

end subroutine init_sponge_diags

subroutine set_up_sponge_field(sp_val, f_ptr, G, nlay, CS, sp_val_i_mean)
  type(ocean_grid_type),                            intent(in) :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),         intent(in) :: sp_val
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target, intent(in) :: f_ptr
  integer,                                          intent(in) :: nlay
  type(sponge_cs),                                  pointer    :: CS
  real, dimension(SZJ_(G),SZK_(G)),       optional, intent(in) :: sp_val_i_mean
!   This subroutine stores the reference profile for the variable
! whose address is given by f_ptr. nlay is the number of layers in
! this variable.

! Arguments: sp_val - The reference profiles of the quantity being
!                     registered.
!  (in)      f_ptr - a pointer to the field which will be damped.
!  (in)      nlay - the number of layers in this quantity.
!  (in/out)  CS - A pointer to the control structure for this module that is
!                 set by a previous call to initialize_sponge.
!  (in,opt)  sp_val_i_mean - The i-mean reference value for this field with
!                 i-mean sponges.

  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_sponge_field: "//mesg)
  endif

  allocate(cs%Ref_val(cs%fldno)%p(cs%nz,cs%num_col))
  cs%Ref_val(cs%fldno)%p(:,:) = 0.0
  do col=1,cs%num_col
    do k=1,nlay
      cs%Ref_val(cs%fldno)%p(k,col) = sp_val(cs%col_i(col),cs%col_j(col),k)
    enddo
    do k=nlay+1,cs%nz
      cs%Ref_val(cs%fldno)%p(k,col) = 0.0
    enddo
  enddo

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

  if (nlay/=cs%nz) then
    write(mesg,'("Danger: Sponge reference fields require nz (",I3,") layers.&
        & A field with ",I3," layers was passed to set_up_sponge_field.")') &
          cs%nz, nlay
    if (is_root_pe()) call mom_error(warning, "set_up_sponge_field: "//mesg)
  endif

  if (cs%do_i_mean_sponge) then
    if (.not.present(sp_val_i_mean)) call mom_error(fatal, &
      "set_up_sponge_field: sp_val_i_mean must be present with i-mean sponges.")

    allocate(cs%Ref_val_im(cs%fldno)%p(cs%jsd:cs%jed,cs%nz))
    cs%Ref_val(cs%fldno)%p(:,:) = 0.0
    do k=1,cs%nz ; do j=cs%jsc,cs%jec
      cs%Ref_val_im(cs%fldno)%p(j,k) = sp_val_i_mean(j,k)
    enddo ; enddo
  endif

end subroutine set_up_sponge_field


subroutine set_up_sponge_ml_density(sp_val, G, CS, sp_val_i_mean)
  type(ocean_grid_type),              intent(in) :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G)),   intent(in) :: sp_val
  type(sponge_cs),                    pointer    :: CS
  real, dimension(SZJ_(G)), optional, intent(in) :: sp_val_i_mean
!   This subroutine stores the reference value for mixed layer density.  It is
! handled differently from other values because it is only used in determining
! which layers can be inflated.

! Arguments: sp_val - The reference values of the mixed layer density.
!  (in/out)  CS - A pointer to the control structure for this module that is
!                 set by a previous call to initialize_sponge.

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

  if (.not.associated(cs)) return

  if (associated(cs%Rcv_ml_ref)) then
    call mom_error(fatal, "set_up_sponge_ML_density appears to have been "//&
                           "called twice.")
  endif

  cs%bulkmixedlayer = .true.
  allocate(cs%Rcv_ml_ref(cs%num_col)) ; cs%Rcv_ml_ref(:) = 0.0
  do col=1,cs%num_col
    cs%Rcv_ml_ref(col) = sp_val(cs%col_i(col),cs%col_j(col))
  enddo

  if (cs%do_i_mean_sponge) then
    if (.not.present(sp_val_i_mean)) call mom_error(fatal, &
      "set_up_sponge_field: sp_val_i_mean must be present with i-mean sponges.")

    allocate(cs%Rcv_ml_ref_im(cs%jsd:cs%jed)) ; cs%Rcv_ml_ref_im(:) = 0.0
    do j=cs%jsc,cs%jec
      cs%Rcv_ml_ref_im(j) = sp_val_i_mean(j)
    enddo
  endif

end subroutine set_up_sponge_ml_density

subroutine apply_sponge(h, dt, G, GV, ea, eb, CS, Rcv_ml)
  type(ocean_grid_type),                    intent(inout) :: G    !< The ocean's grid structure
  type(verticalgrid_type),                  intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  real,                                     intent(in)    :: dt
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out)   :: ea
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out)   :: eb
  type(sponge_cs),                          pointer       :: CS
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: Rcv_ml

! This subroutine applies damping to the layers thicknesses, mixed
! layer buoyancy, and a variety of tracers for every column where
! there is damping.

! Arguments: h -  Layer thickness, in m.
!  (in)      dt - The amount of time covered by this call, in s.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (out)     ea - an array to which the amount of fluid entrained
!                 from the layer above during this call will be
!                 added, in m.
!  (out)     eb - an array to which the amount of fluid entrained
!                 from the layer below during this call will be
!                 added, in m.
!  (in)      CS - A pointer to the control structure for this module that is
!                 set by a previous call to initialize_sponge.
!  (inout,opt)  Rcv_ml - The coordinate density of the mixed layer.

  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
    w_int, &       ! Water moved upward across an interface within a timestep,
                   ! in m.
    e_d            ! Interface heights that are dilated to have a value of 0
                   ! at the surface, in m.
  real, dimension(SZI_(G), SZJ_(G)) :: &
    eta_anom, &    ! Anomalies in the interface height, relative to the i-mean
                   ! target value.
    fld_anom       ! Anomalies in a tracer concentration, relative to the
                   ! i-mean target value.
  real, dimension(SZJ_(G), SZK_(G)+1) :: &
    eta_mean_anom  ! The i-mean interface height anomalies.
  real, allocatable, dimension(:,:,:) :: &
    fld_mean_anom  ! THe i-mean tracer concentration anomalies.
  real, dimension(SZI_(G), SZK_(G)+1) :: &
    h_above, &     ! The total thickness above an interface, in m.
    h_below        ! The total thickness below an interface, in m.
  real, dimension(SZI_(G)) :: &
    dilate         ! A nondimensional factor by which to dilate layers to
                   ! give 0 at the surface.

  real :: e(szk_(g)+1)  ! The interface heights, in m or kg m-2, usually negative.
  real :: e0       ! The height of the free surface in m.
  real :: e_str    ! A nondimensional amount by which the reference
                   ! profile must be stretched for the free surfaces
                   ! heights in the two profiles to agree.
  real :: w        ! The thickness of water moving upward through an
                   ! interface within 1 timestep, in m.
  real :: wm       ! wm is w if w is negative and 0 otherwise, in m.
  real :: wb       ! w at the interface below a layer, in m.
  real :: wpb      ! wpb is wb if wb is positive and 0 otherwise, m.
  real :: ea_k, eb_k
  real :: damp     ! The timestep times the local damping  coefficient.  ND.
  real :: I1pdamp  ! I1pdamp is 1/(1 + damp).  Nondimensional.
  real :: damp_1pdamp ! damp_1pdamp is damp/(1 + damp).  Nondimensional.
  real :: Idt      ! 1.0/dt, in s-1.
  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
  if (cs%bulkmixedlayer) nkmb = gv%nk_rho_varies
  if (cs%bulkmixedlayer .and. (.not.present(rcv_ml))) &
    call mom_error(fatal, "Rml must be provided to apply_sponge when using "//&
                           "a bulk mixed layer.")

  if ((cs%id_w_sponge > 0) .or. cs%do_i_mean_sponge) then
    do k=1,nz+1 ; do j=js,je ; do i=is,ie
      w_int(i,j,k) = 0.0
    enddo ; enddo ; enddo
  endif

  if (cs%do_i_mean_sponge) then
    ! Apply forcing to restore the zonal-mean properties to prescribed values.

    if (cs%bulkmixedlayer) call mom_error(fatal, "apply_sponge is not yet set up to "//&
                  "work properly with i-mean sponges and a bulk mixed layer.")

    do j=js,je ; do i=is,ie ; e_d(i,j,nz+1) = -g%bathyT(i,j) ; enddo ; enddo
    do k=nz,1,-1 ; do j=js,je ; do i=is,ie
      e_d(i,j,k) = e_d(i,j,k+1) + h(i,j,k)
    enddo ; enddo ; enddo
    do j=js,je
      do i=is,ie
        dilate(i) = g%bathyT(i,j) / (e_d(i,j,1) + g%bathyT(i,j))
      enddo
      do k=1,nz+1 ; do i=is,ie
        e_d(i,j,k) = dilate(i) * (e_d(i,j,k) + g%bathyT(i,j)) - g%bathyT(i,j)
      enddo ; enddo
    enddo

    do k=2,nz
      do j=js,je ; do i=is,ie
        eta_anom(i,j) = e_d(i,j,k) - cs%Ref_eta_im(j,k)
        if (cs%Ref_eta_im(j,k) < -g%bathyT(i,j)) eta_anom(i,j) = 0.0
      enddo ; enddo
      call global_i_mean(eta_anom(:,:), eta_mean_anom(:,k), g)
    enddo

    if (cs%fldno > 0) allocate(fld_mean_anom(g%isd:g%ied,nz,cs%fldno))
    do m=1,cs%fldno
      do j=js,je ; do i=is,ie
        fld_anom(i,j) = cs%var(m)%p(i,j,k) - cs%Ref_val_im(m)%p(j,k)
      enddo ; enddo
      call global_i_mean(fld_anom(:,:), fld_mean_anom(:,k,m), g, h(:,:,k))
    enddo

    do j=js,je ; if (cs%Iresttime_im(j) > 0.0) then
      damp = dt*cs%Iresttime_im(j) ; damp_1pdamp = damp / (1.0 + damp)

      do i=is,ie
        h_above(i,1) = 0.0 ; h_below(i,nz+1) = 0.0
      enddo
      do k=nz,1,-1 ; do i=is,ie
        h_below(i,k) = h_below(i,k+1) + max(h(i,j,k)-gv%Angstrom, 0.0)
      enddo ; enddo
      do k=2,nz+1 ; do i=is,ie
        h_above(i,k) = h_above(i,k-1) + max(h(i,j,k-1)-gv%Angstrom, 0.0)
      enddo ; enddo
      do k=2,nz
        ! w is positive for an upward (lightward) flux of mass, resulting
        ! in the downward movement of an interface.
        w = damp_1pdamp * eta_mean_anom(j,k)
        do i=is,ie
          if (w > 0.0) then
            w_int(i,j,k) = min(w, h_below(i,k))
            eb(i,j,k-1) = eb(i,j,k-1) + w_int(i,j,k)
          else
            w_int(i,j,k) = max(w, -h_above(i,k))
            ea(i,j,k) = ea(i,j,k) - w_int(i,j,k)
          endif
        enddo
      enddo
      do k=1,nz ; do i=is,ie
        ea_k = max(0.0, -w_int(i,j,k))
        eb_k = max(0.0, w_int(i,j,k+1))
        do m=1,cs%fldno
          cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
              cs%Ref_val_im(m)%p(j,k) * (ea_k + eb_k)) / &
                     (h(i,j,k) + (ea_k + eb_k)) - &
              damp_1pdamp * fld_mean_anom(j,k,m)
        enddo

        h(i,j,k) = max(h(i,j,k) + (w_int(i,j,k+1) - w_int(i,j,k)), &
                       min(h(i,j,k), gv%Angstrom))
      enddo ; enddo
    endif ; enddo

    if (cs%fldno > 0) deallocate(fld_mean_anom)

  endif

  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)

    e(1) = 0.0 ; e0 = 0.0
    do k=1,nz
      e(k+1) = e(k) - h(i,j,k)
    enddo
    e_str = e(nz+1) / cs%Ref_eta(nz+1,c)

    if ( cs%bulkmixedlayer ) then
      i1pdamp = 1.0 / (1.0 + damp)
      if (associated(cs%Rcv_ml_ref)) &
        rcv_ml(i,j) = i1pdamp * (rcv_ml(i,j) + cs%Rcv_ml_ref(c)*damp)
      do k=1,nkmb
        do m=1,cs%fldno
          cs%var(m)%p(i,j,k) = i1pdamp * &
              (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(k,c)*damp)
        enddo
      enddo

      wpb = 0.0; wb = 0.0
      do k=nz,nkmb+1,-1
        if (gv%Rlay(k) > rcv_ml(i,j)) then
          w = min((((e(k)-e0) - e_str*cs%Ref_eta(k,c)) * damp), &
                    ((wb + h(i,j,k)) - gv%Angstrom))
          wm = 0.5*(w-abs(w))
          do m=1,cs%fldno
            cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
                     cs%Ref_val(m)%p(k,c)*(damp*h(i,j,k) + (wpb - wm))) / &
                     (h(i,j,k)*(1.0 + damp) + (wpb - wm))
          enddo
        else
          do m=1,cs%fldno
            cs%var(m)%p(i,j,k) = i1pdamp * &
              (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(k,c)*damp)
          enddo
          w = wb + (h(i,j,k) - gv%Angstrom)
          wm = 0.5*(w-abs(w))
        endif
        eb(i,j,k) = eb(i,j,k) + wpb
        ea(i,j,k) = ea(i,j,k) - wm
        h(i,j,k)  = h(i,j,k)  + (wb - w)
        wb = w
        wpb = w - wm
      enddo

      if (wb < 0) then
        do k=nkmb,1,-1
          w = min((wb + (h(i,j,k) - gv%Angstrom)),0.0)
          h(i,j,k)  = h(i,j,k)  + (wb - w)
          ea(i,j,k) = ea(i,j,k) - w
          wb = w
        enddo
      else
        w = wb
        do k=gv%nkml,nkmb
          eb(i,j,k) = eb(i,j,k) + w
        enddo

        k = gv%nkml
        h(i,j,k) = h(i,j,k) + w
        do m=1,cs%fldno
          cs%var(m)%p(i,j,k) = (cs%var(m)%p(i,j,k)*h(i,j,k) + &
                                cs%Ref_val(m)%p(k,c)*w) / (h(i,j,k) + w)
        enddo
      endif

      do k=1,nkmb
        do m=1,cs%fldno
          cs%var(m)%p(i,j,k) = i1pdamp * &
              (cs%var(m)%p(i,j,k) + cs%Ref_val(m)%p(gv%nkml,c)*damp)
        enddo
      enddo

    else                                          ! not BULKMIXEDLAYER

      wpb = 0.0
      wb = 0.0
      do k=nz,1,-1
        w = min((((e(k)-e0) - e_str*cs%Ref_eta(k,c)) * damp), &
                  ((wb + h(i,j,k)) - gv%Angstrom))
        wm = 0.5*(w - abs(w))
        do m=1,cs%fldno
          cs%var(m)%p(i,j,k) = (h(i,j,k)*cs%var(m)%p(i,j,k) + &
              cs%Ref_val(m)%p(k,c) * (damp*h(i,j,k) + (wpb - wm))) / &
                     (h(i,j,k)*(1.0 + damp) + (wpb - wm))
        enddo
        eb(i,j,k) = eb(i,j,k) + wpb
        ea(i,j,k) = ea(i,j,k) - wm
        h(i,j,k)  = h(i,j,k)  + (wb - w)
        wb = w
        wpb = w - wm
      enddo

    endif                                         ! end BULKMIXEDLAYER
  enddo ! end of c loop

  if (associated(cs%diag)) then ; if (query_averaging_enabled(cs%diag)) then
    idt = 1.0 / dt
    if (cs%id_w_sponge > 0) call post_data(cs%id_w_sponge, idt*w_int, cs%diag)
  endif ; endif

end subroutine apply_sponge

subroutine sponge_end(CS)
  type(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_j)) deallocate(cs%col_j)

  if (associated(cs%Iresttime_col)) deallocate(cs%Iresttime_col)
  if (associated(cs%Rcv_ml_ref)) deallocate(cs%Rcv_ml_ref)
  if (associated(cs%Ref_eta)) deallocate(cs%Ref_eta)

  if (associated(cs%Iresttime_im)) deallocate(cs%Iresttime_im)
  if (associated(cs%Rcv_ml_ref_im)) deallocate(cs%Rcv_ml_ref_im)
  if (associated(cs%Ref_eta_im)) deallocate(cs%Ref_eta_im)

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

  deallocate(cs)

end subroutine sponge_end

end module mom_sponge