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