C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_seaice_sponge.F,v 1.2 2012/09/25 16:39:20 dimitri Exp $
C $Name:  $

#include "OBCS_OPTIONS.h"

C--   File obcs_seaice_sponge.F:
C--    Contents:
C--    o OBCS_SEAICE_SPONGE_A
C--    o OBCS_SEAICE_SPONGE_H
C--    o OBCS_SEAICE_SPONGE_SL
C--    o OBCS_SEAICE_SPONGE_SN

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CStartOfInterface
      SUBROUTINE OBCS_SEAICE_SPONGE_A( myThid )
C     *==========================================================*
C     | S/R OBCS_SEAICE_SPONGE_A
C     | Adds a relaxation term to AREA near Open-Boundaries
C     *==========================================================*
      IMPLICIT NONE

C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "OBCS_SEAICE.h"
#ifdef ALLOW_SEAICE
# include "SEAICE_SIZE.h"
# include "SEAICE_PARAMS.h"
# include "SEAICE.h"
#endif

C     == Routine arguments ==
      INTEGER myThid
CEndOfInterface

#if (defined(ALLOW_OBCS) && defined(ALLOW_SEAICE) && defined(ALLOW_OBCS_SEAICE_SPONGE))
C     == Local variables ==
C     Loop counters
      INTEGER bi, bj, i, j, isl, jsl
      _RL lambda_obcs

      IF ( useSeaiceSponge .AND. seaiceSpongeThickness.NE.0 ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

C Northern Open Boundary
# ifdef ALLOW_OBCS_NORTH
         IF ( tileHasOBN(bi,bj) ) THEN
          DO i=1,sNx
           IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN
            DO jsl= 1,seaiceSpongeThickness
             j=OB_Jn(i,bi,bj)-jsl
             IF ((j.ge.1).and.(j.le.sNy)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-jsl)*Arelaxobcsbound
     &           + float(jsl-1)*Arelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              IF (lambda_obcs.ne.0.) THEN
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              ELSE
               lambda_obcs = 0. _d 0
              ENDIF
              AREA(i,j,bi,bj) = AREA(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( AREA(i,j,bi,bj) - OBNa(i,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Southern Open Boundary
# ifdef ALLOW_OBCS_SOUTH
         IF ( tileHasOBS(bi,bj) ) THEN
          DO i=1,sNx
           IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN
            DO jsl= 1,seaiceSpongeThickness
             j=OB_Js(i,bi,bj)+jsl
             IF ((j.ge.1).and.(j.le.sNy)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-jsl)*Arelaxobcsbound
     &           + float(jsl-1)*Arelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              AREA(i,j,bi,bj) = AREA(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( AREA(i,j,bi,bj) - OBSa(i,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Eastern Open Boundary
# ifdef ALLOW_OBCS_EAST
         IF ( tileHasOBE(bi,bj) ) THEN
          DO j=1,sNy
           IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
            DO isl= 1,seaiceSpongeThickness
             i=OB_Ie(j,bi,bj)-isl
             IF ((i.ge.1).and.(i.le.sNx)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-isl)*Arelaxobcsbound
     &           + float(isl-1)*Arelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              AREA(i,j,bi,bj) = AREA(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( AREA(i,j,bi,bj) - OBEa(j,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Western Open Boundary
# ifdef ALLOW_OBCS_WEST
         IF ( tileHasOBW(bi,bj) ) THEN
          DO j=1,sNy
           IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
            DO isl= 1,seaiceSpongeThickness
             i=OB_Iw(j,bi,bj)+isl
             IF ((i.ge.1).and.(i.le.sNx)) THEN
              lambda_obcs= (
     &           float(seaiceSpongeThickness-isl)*Arelaxobcsbound
     &           + float(isl-1)*Arelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              AREA(i,j,bi,bj) = AREA(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( AREA(i,j,bi,bj) - OBWa(j,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

        ENDDO
       ENDDO
      ENDIF

#endif /* ALLOW_OBCS & ALLOW_SEAICE & ALLOW_OBCS_SEAICE_SPONGE */

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CStartOfInterface
      SUBROUTINE OBCS_SEAICE_SPONGE_H( myThid )
C     *==========================================================*
C     | S/R OBCS_SEAICE_SPONGE_H
C     | Adds a relaxation term to HEFF near Open-Boundaries
C     *==========================================================*
      IMPLICIT NONE

C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "OBCS_SEAICE.h"
#ifdef ALLOW_SEAICE
# include "SEAICE_SIZE.h"
# include "SEAICE_PARAMS.h"
# include "SEAICE.h"
#endif

C     == Routine arguments ==
      INTEGER myThid
CEndOfInterface

#if (defined(ALLOW_OBCS) && defined(ALLOW_SEAICE) && defined(ALLOW_OBCS_SEAICE_SPONGE))
C     == Local variables ==
C     Loop counters
      INTEGER bi, bj, i, j, isl, jsl
      _RL lambda_obcs

      IF ( useSeaiceSponge .AND. seaiceSpongeThickness.NE.0 ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

C Northern Open Boundary
# ifdef ALLOW_OBCS_NORTH
         IF ( tileHasOBN(bi,bj) ) THEN
          DO i=1,sNx
           IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN
            DO jsl= 1,seaiceSpongeThickness
             j=OB_Jn(i,bi,bj)-jsl
             IF ((j.ge.1).and.(j.le.sNy)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-jsl)*Hrelaxobcsbound
     &           + float(jsl-1)*Hrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              IF (lambda_obcs.ne.0.) THEN
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              ELSE
               lambda_obcs = 0. _d 0
              ENDIF
              HEFF(i,j,bi,bj) = HEFF(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HEFF(i,j,bi,bj) - OBNh(i,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Southern Open Boundary
# ifdef ALLOW_OBCS_SOUTH
         IF ( tileHasOBS(bi,bj) ) THEN
          DO i=1,sNx
           IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN
            DO jsl= 1,seaiceSpongeThickness
             j=OB_Js(i,bi,bj)+jsl
             IF ((j.ge.1).and.(j.le.sNy)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-jsl)*Hrelaxobcsbound
     &           + float(jsl-1)*Hrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HEFF(i,j,bi,bj) = HEFF(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HEFF(i,j,bi,bj) - OBSh(i,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Eastern Open Boundary
# ifdef ALLOW_OBCS_EAST
         IF ( tileHasOBE(bi,bj) ) THEN
          DO j=1,sNy
           IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
            DO isl= 1,seaiceSpongeThickness
             i=OB_Ie(j,bi,bj)-isl
             IF ((i.ge.1).and.(i.le.sNx)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-isl)*Hrelaxobcsbound
     &           + float(isl-1)*Hrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HEFF(i,j,bi,bj) = HEFF(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HEFF(i,j,bi,bj) - OBEh(j,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Western Open Boundary
# ifdef ALLOW_OBCS_WEST
         IF ( tileHasOBW(bi,bj) ) THEN
          DO j=1,sNy
           IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
            DO isl= 1,seaiceSpongeThickness
             i=OB_Iw(j,bi,bj)+isl
             IF ((i.ge.1).and.(i.le.sNx)) THEN
              lambda_obcs= (
     &           float(seaiceSpongeThickness-isl)*Hrelaxobcsbound
     &           + float(isl-1)*Hrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HEFF(i,j,bi,bj) = HEFF(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HEFF(i,j,bi,bj) - OBWh(j,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

        ENDDO
       ENDDO
      ENDIF

#endif /* ALLOW_OBCS & ALLOW_SEAICE & ALLOW_OBCS_SEAICE_SPONGE */

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CStartOfInterface
      SUBROUTINE OBCS_SEAICE_SPONGE_SL( myThid )
C     *==========================================================*
C     | S/R OBCS_SEAICE_SPONGE_SL
C     | Adds a relaxation term to HSALT near Open-Boundaries
C     *==========================================================*
      IMPLICIT NONE

C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "OBCS_SEAICE.h"
#ifdef ALLOW_SEAICE
# include "SEAICE_SIZE.h"
# include "SEAICE_PARAMS.h"
# include "SEAICE.h"
#endif

C     == Routine arguments ==
      INTEGER myThid
CEndOfInterface

#if (defined(ALLOW_OBCS) && defined(ALLOW_SEAICE) && defined(ALLOW_OBCS_SEAICE_SPONGE) && defined(SEAICE_VARIABLE_SALINITY))
C     == Local variables ==
C     Loop counters
      INTEGER bi, bj, i, j, isl, jsl
      _RL lambda_obcs

      IF ( useSeaiceSponge .AND. seaiceSpongeThickness.NE.0 ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

C Northern Open Boundary
# ifdef ALLOW_OBCS_NORTH
         IF ( tileHasOBN(bi,bj) ) THEN
          DO i=1,sNx
           IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN
            DO jsl= 1,seaiceSpongeThickness
             j=OB_Jn(i,bi,bj)-jsl
             IF ((j.ge.1).and.(j.le.sNy)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-jsl)*SLrelaxobcsbound
     &           + float(jsl-1)*SLrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              IF (lambda_obcs.ne.0.) THEN
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              ELSE
               lambda_obcs = 0. _d 0
              ENDIF
              HSALT(i,j,bi,bj) = HSALT(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HSALT(i,j,bi,bj) - OBNsl(i,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Southern Open Boundary
# ifdef ALLOW_OBCS_SOUTH
         IF ( tileHasOBS(bi,bj) ) THEN
          DO i=1,sNx
           IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN
            DO jsl= 1,seaiceSpongeThickness
             j=OB_Js(i,bi,bj)+jsl
             IF ((j.ge.1).and.(j.le.sNy)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-jsl)*SLrelaxobcsbound
     &           + float(jsl-1)*SLrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HSALT(i,j,bi,bj) = HSALT(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HSALT(i,j,bi,bj) - OBSsl(i,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Eastern Open Boundary
# ifdef ALLOW_OBCS_EAST
         IF ( tileHasOBE(bi,bj) ) THEN
          DO j=1,sNy
           IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
            DO isl= 1,seaiceSpongeThickness
             i=OB_Ie(j,bi,bj)-isl
             IF ((i.ge.1).and.(i.le.sNx)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-isl)*SLrelaxobcsbound
     &           + float(isl-1)*SLrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HSALT(i,j,bi,bj) = HSALT(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HSALT(i,j,bi,bj) - OBEsl(j,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Western Open Boundary
# ifdef ALLOW_OBCS_WEST
         IF ( tileHasOBW(bi,bj) ) THEN
          DO j=1,sNy
           IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
            DO isl= 1,seaiceSpongeThickness
             i=OB_Iw(j,bi,bj)+isl
             IF ((i.ge.1).and.(i.le.sNx)) THEN
              lambda_obcs= (
     &           float(seaiceSpongeThickness-isl)*SLrelaxobcsbound
     &           + float(isl-1)*SLrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HSALT(i,j,bi,bj) = HSALT(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HSALT(i,j,bi,bj) - OBWsl(j,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

        ENDDO
       ENDDO
      ENDIF

#endif /* ALLOW_OBCS & ALLOW_SEAICE & ALLOW_OBCS_SEAICE_SPONGE & SEAICE_VARIABLE_SALINITY */

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CStartOfInterface
      SUBROUTINE OBCS_SEAICE_SPONGE_SN( myThid )
C     *==========================================================*
C     | S/R OBCS_SEAICE_SPONGE_SN
C     | Adds a relaxation term to HSNOW near Open-Boundaries
C     *==========================================================*
      IMPLICIT NONE

C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "OBCS_SEAICE.h"
#ifdef ALLOW_SEAICE
# include "SEAICE_SIZE.h"
# include "SEAICE_PARAMS.h"
# include "SEAICE.h"
#endif

C     == Routine arguments ==
      INTEGER myThid
CEndOfInterface

#if (defined(ALLOW_OBCS) && defined(ALLOW_SEAICE) && defined(ALLOW_OBCS_SEAICE_SPONGE))
C     == Local variables ==
C     Loop counters
      INTEGER bi, bj, i, j, isl, jsl
      _RL lambda_obcs

      IF ( useSeaiceSponge .AND. seaiceSpongeThickness.NE.0 ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)

C Northern Open Boundary
# ifdef ALLOW_OBCS_NORTH
         IF ( tileHasOBN(bi,bj) ) THEN
          DO i=1,sNx
           IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN
            DO jsl= 1,seaiceSpongeThickness
             j=OB_Jn(i,bi,bj)-jsl
             IF ((j.ge.1).and.(j.le.sNy)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-jsl)*SNrelaxobcsbound
     &           + float(jsl-1)*SNrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              IF (lambda_obcs.ne.0.) THEN
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              ELSE
               lambda_obcs = 0. _d 0
              ENDIF
              HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HSNOW(i,j,bi,bj) - OBNsn(i,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Southern Open Boundary
# ifdef ALLOW_OBCS_SOUTH
         IF ( tileHasOBS(bi,bj) ) THEN
          DO i=1,sNx
           IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN
            DO jsl= 1,seaiceSpongeThickness
             j=OB_Js(i,bi,bj)+jsl
             IF ((j.ge.1).and.(j.le.sNy)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-jsl)*SNrelaxobcsbound
     &           + float(jsl-1)*SNrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HSNOW(i,j,bi,bj) - OBSsn(i,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Eastern Open Boundary
# ifdef ALLOW_OBCS_EAST
         IF ( tileHasOBE(bi,bj) ) THEN
          DO j=1,sNy
           IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
            DO isl= 1,seaiceSpongeThickness
             i=OB_Ie(j,bi,bj)-isl
             IF ((i.ge.1).and.(i.le.sNx)) THEN
              lambda_obcs = (
     &           float(seaiceSpongeThickness-isl)*SNrelaxobcsbound
     &           + float(isl-1)*SNrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HSNOW(i,j,bi,bj) - OBEsn(j,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

C Western Open Boundary
# ifdef ALLOW_OBCS_WEST
         IF ( tileHasOBW(bi,bj) ) THEN
          DO j=1,sNy
           IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
            DO isl= 1,seaiceSpongeThickness
             i=OB_Iw(j,bi,bj)+isl
             IF ((i.ge.1).and.(i.le.sNx)) THEN
              lambda_obcs= (
     &           float(seaiceSpongeThickness-isl)*SNrelaxobcsbound
     &           + float(isl-1)*SNrelaxobcsinner)
     &           / float(seaiceSpongeThickness-1)
              if (lambda_obcs.ne.0.) then
               lambda_obcs = SEAICE_deltaTtherm / lambda_obcs
              else
               lambda_obcs = 0. _d 0
              endif
              HSNOW(i,j,bi,bj) = HSNOW(i,j,bi,bj)
     &           - maskC(i,j,1,bi,bj) * lambda_obcs
     &           * ( HSNOW(i,j,bi,bj) - OBWsn(j,bi,bj) )
             ENDIF
            ENDDO
           ENDIF
          ENDDO
         ENDIF
# endif

        ENDDO
       ENDDO
      ENDIF

#endif /* ALLOW_OBCS & ALLOW_SEAICE & ALLOW_OBCS_SEAICE_SPONGE */

      RETURN
      END