turn Neumann BC CPP-flag into runtime flag useSeaiceNeumann

first attempt

pkg/obcs/OBCS_OPTIONS.h

@@ -63,9 +63,6 @@ C     Compute rather than specify seaice velocities at the edges.
 #undef OBCS_SEAICE_COMPUTE_UVICE
 #endif /* OBCS_UVICE_OLD */
 
-C     use Neuman boundary conditions for all sea ice variables
-#undef OBCS_SEAICE_NEUMANN
-
 C     Smooth the tracer sea-ice variables near the edges.
 #undef OBCS_SEAICE_SMOOTH_EDGE
 

pkg/obcs/OBCS_SEAICE.h

@@ -18,10 +18,12 @@ C seaiceSpongeThickness  :: number grid points that make up the sponge layer (de
      & seaiceSpongeThickness
       INTEGER seaiceSpongeThickness
 
-C useSeaiceSponge :: turns on seaice sponge layer along boundary (def=false)
+C useSeaiceSponge  :: turns on seaice sponge layer along boundary (def=false)
+C useSeaiceNeumann :: use Neumann conditions for sea ice variables (def=false)
       COMMON /OBC_SEAICE_PARM_L/
-     & useSeaiceSponge
+     & useSeaiceSponge, useSeaiceNeumann
       LOGICAL useSeaiceSponge
+      LOGICAL useSeaiceNeumann
 
 C [A,H,SL,SN]relaxobcs[inner,bound] :: relaxation time scale (in seconds) on the
 C                            boundary (bound) and at the innermost grid point of the

pkg/obcs/obcs_apply_seaice.F

@@ -57,22 +57,22 @@ SUBROUTINE OBCS_APPLY_SEAICE( myThid )
          DO i=1-OLx,sNx+OLx
           Jobc = OB_Jn(I,bi,bj)
           IF ( Jobc.NE.OB_indexNone ) THEN
-# ifdef SEAICE_OBCS_NEUMANN
+           IF ( useSeaiceNeumann ) THEN
 C     Neumann boundary conditions
             HEFF (i,Jobc,bi,bj) = HEFF (i,Jobc-1,bi,bj)
             AREA (i,Jobc,bi,bj) = AREA (i,Jobc-1,bi,bj)
             HSNOW(i,Jobc,bi,bj) = HSNOW(i,Jobc-1,bi,bj)
 #  ifdef SEAICE_VARIABLE_SALINITY
             HSALT(i,Jobc,bi,bj) = HSALT(i,Jobc-1,bi,bj)
 #  endif
-# else /* not SEAICE_OBCS_NEUMANN (default) */
+           ELSE
             HEFF (i,Jobc,bi,bj) = OBNh (i,bi,bj)
             AREA (i,Jobc,bi,bj) = OBNa (i,bi,bj)
             HSNOW(i,Jobc,bi,bj) = OBNsn(i,bi,bj)
 #  ifdef SEAICE_VARIABLE_SALINITY
             HSALT(i,Jobc,bi,bj) = OBNsl(i,bi,bj)
 #  endif
-# endif
+           ENDIF
           ENDIF
          ENDDO
         ENDIF
@@ -84,22 +84,22 @@ SUBROUTINE OBCS_APPLY_SEAICE( myThid )
          DO i=1-OLx,sNx+OLx
           Jobc = OB_Js(I,bi,bj)
           IF ( Jobc.NE.OB_indexNone ) THEN
-# ifdef SEAICE_OBCS_NEUMANN
+           IF ( useSeaiceNeumann ) THEN
 C     Neumann boundary conditions
             HEFF (i,Jobc,bi,bj) = HEFF (i,Jobc+1,bi,bj)
             AREA (i,Jobc,bi,bj) = AREA (i,Jobc+1,bi,bj)
             HSNOW(i,Jobc,bi,bj) = HSNOW(i,Jobc+1,bi,bj)
 #  ifdef SEAICE_VARIABLE_SALINITY
             HSALT(i,Jobc,bi,bj) = HSALT(i,Jobc+1,bi,bj)
 #  endif
-# else /* not SEAICE_OBCS_NEUMANN (default) */
+           ELSE
             HEFF (i,Jobc,bi,bj) = OBSh (i,bi,bj)
             AREA (i,Jobc,bi,bj) = OBSa (i,bi,bj)
             HSNOW(i,Jobc,bi,bj) = OBSsn(i,bi,bj)
 #  ifdef SEAICE_VARIABLE_SALINITY
             HSALT(i,Jobc,bi,bj) = OBSsl(i,bi,bj)
 #  endif
-# endif
+           ENDIF
           ENDIF
          ENDDO
         ENDIF
@@ -111,22 +111,22 @@ SUBROUTINE OBCS_APPLY_SEAICE( myThid )
          DO j=1-OLy,sNy+OLy
           Iobc = OB_Ie(J,bi,bj)
           IF ( Iobc.NE.OB_indexNone ) THEN
-# ifdef SEAICE_OBCS_NEUMANN
+           IF ( useSeaiceNeumann ) THEN
 C     Neumann boundary conditions
             HEFF (Iobc,j,bi,bj) = HEFF (Iobc-1,j,bi,bj)
             AREA (Iobc,j,bi,bj) = AREA (Iobc-1,j,bi,bj)
             HSNOW(Iobc,j,bi,bj) = HSNOW(Iobc-1,j,bi,bj)
 #  ifdef SEAICE_VARIABLE_SALINITY
             HSALT(Iobc,j,bi,bj) = HSALT(Iobc-1,j,bi,bj)
 #  endif
-# else /* not SEAICE_OBCS_NEUMANN (default) */
+           ELSE
             HEFF (Iobc,j,bi,bj) = OBEh (j,bi,bj)
             AREA (Iobc,j,bi,bj) = OBEa (j,bi,bj)
             HSNOW(Iobc,j,bi,bj) = OBEsn(j,bi,bj)
 #  ifdef SEAICE_VARIABLE_SALINITY
             HSALT(Iobc,j,bi,bj) = OBEsl(j,bi,bj)
 #  endif
-#endif
+           ENDIF
           ENDIF
          ENDDO
         ENDIF
@@ -138,22 +138,22 @@ SUBROUTINE OBCS_APPLY_SEAICE( myThid )
          DO j=1-OLy,sNy+OLy
           Iobc = OB_Iw(j,bi,bj)
           IF ( Iobc.NE.OB_indexNone ) THEN
-# ifdef SEAICE_OBCS_NEUMANN
+           IF ( useSeaiceNeumann ) THEN
 C     Neumann boundary conditions
             HEFF (Iobc,j,bi,bj) = HEFF (Iobc+1,j,bi,bj)
             AREA (Iobc,j,bi,bj) = AREA (Iobc+1,j,bi,bj)
             HSNOW(Iobc,j,bi,bj) = HSNOW(Iobc+1,j,bi,bj)
 #  ifdef SEAICE_VARIABLE_SALINITY
             HSALT(Iobc,j,bi,bj) = HSALT(Iobc+1,j,bi,bj)
 #  endif
-# else /* not SEAICE_OBCS_NEUMANN (default) */
+           ELSE
             HEFF (Iobc,j,bi,bj) = OBWh (j,bi,bj)
             AREA (Iobc,j,bi,bj) = OBWa (j,bi,bj)
             HSNOW(Iobc,j,bi,bj) = OBWsn(j,bi,bj)
 #  ifdef SEAICE_VARIABLE_SALINITY
             HSALT(Iobc,j,bi,bj) = OBWsl(j,bi,bj)
 #  endif
-#endif
+           ENDIF
           ENDIF
          ENDDO
         ENDIF

pkg/obcs/obcs_apply_uvice.F

@@ -66,6 +66,60 @@ SUBROUTINE OBCS_APPLY_UVICE(
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
 
+        IF ( useSeaiceNeumann ) THEN
+# ifdef ALLOW_OBCS_NORTH
+         IF ( tileHasOBN(bi,bj) ) THEN
+C     Northern boundary
+          DO i=1-OLx,sNx+OLx
+           Jobc = OB_Jn(i,bi,bj)
+           IF ( Jobc.NE.OB_indexNone ) THEN
+            OBNuice(i,bi,bj) = uFld(i,Jobc-1,bi,bj)
+            OBNvice(i,bi,bj) = vFld(i,Jobc-1,bi,bj)
+           ENDIF
+          ENDDO
+         ENDIF
+# endif /* ALLOW_OBCS_NORTH */
+
+# ifdef ALLOW_OBCS_SOUTH
+         IF ( tileHasOBS(bi,bj) ) THEN
+C     Southern boundary
+          DO i=1-OLx,sNx+OLx
+           Jobc = OB_Js(i,bi,bj)
+           IF ( Jobc.NE.OB_indexNone ) THEN
+            OBSuice(i,bi,bj) = uFld(i,Jobc+1,bi,bj)
+            OBNvice(i,bi,bj) = vFld(i,Jobc+2,bi,bj)
+           ENDIF
+          ENDDO
+         ENDIF
+# endif /* ALLOW_OBCS_SOUTH */
+
+# ifdef ALLOW_OBCS_EAST
+         IF ( tileHasOBE(bi,bj) ) THEN
+C     Eastern boundary
+          DO j=1-OLy,sNy+OLy
+           Iobc = OB_Ie(j,bi,bj)
+           IF ( Iobc.NE.OB_indexNone ) THEN
+            OBEuice(j,bi,bj) = uFld(Iobc-1,j,bi,bj)
+            OBEvice(j,bi,bj) = vFld(Iobc-1,j,bi,bj)
+           ENDIF
+          ENDDO
+         ENDIF
+# endif /* ALLOW_OBCS_EAST */
+
+# ifdef ALLOW_OBCS_WEST
+         IF ( tileHasOBW(bi,bj) ) THEN
+C     Western boundary
+          DO j=1-OLy,sNy+OLy
+           Iobc = OB_Iw(j,bi,bj)
+           IF ( Iobc.NE.OB_indexNone ) THEN
+            OBEuice(j,bi,bj) = uFld(Iobc+2,j,bi,bj)
+            OBWvice(j,bi,bj) = vFld(Iobc+1,j,bi,bj)
+           ENDIF
+          ENDDO
+         ENDIF
+# endif /* ALLOW_OBCS_WEST */
+        ENDIF
+
 C--   Set Tangential component first:
 
 C     Set model variables to OB values on North/South Boundaries
@@ -75,12 +129,7 @@ SUBROUTINE OBCS_APPLY_UVICE(
          DO i=1-OLx,sNx+OLx
           Jobc = OB_Jn(i,bi,bj)
           IF ( Jobc.NE.OB_indexNone ) THEN
-#  ifdef SEAICE_OBCS_NEUMANN
-C     Neumann boundary conditions
-           uFld(i,Jobc,bi,bj) = uFld(i,Jobc-1,bi,bj)
-#  else
            uFld(i,Jobc,bi,bj) = OBNuice(i,bi,bj)
-#  endif
      &                        *seaiceMaskU(i,Jobc,bi,bj)
           ENDIF
          ENDDO
@@ -93,12 +142,7 @@ SUBROUTINE OBCS_APPLY_UVICE(
          DO i=1-OLx,sNx+OLx
           Jobc = OB_Js(i,bi,bj)
           IF ( Jobc.NE.OB_indexNone ) THEN
-#  ifdef SEAICE_OBCS_NEUMANN
-C     Neumann boundary conditions
-           uFld(i,Jobc,bi,bj) = uFld(i,Jobc+1,bi,bj)
-#  else
            uFld(i,Jobc,bi,bj) = OBSuice(i,bi,bj)
-#  endif
      &                        *seaiceMaskU(i,Jobc,bi,bj)
           ENDIF
          ENDDO
@@ -112,12 +156,7 @@ SUBROUTINE OBCS_APPLY_UVICE(
          DO j=1-OLy,sNy+OLy
           Iobc = OB_Ie(j,bi,bj)
           IF ( Iobc.NE.OB_indexNone ) THEN
-#  ifdef SEAICE_OBCS_NEUMANN
-C     Neumann boundary conditions
-           vFld(Iobc,j,bi,bj) = vFld(Iobc-1,j,bi,bj)
-#  else
            vFld(Iobc,j,bi,bj) = OBEvice(j,bi,bj)
-#  endif
      &                        *seaiceMaskV(Iobc,j,bi,bj)
           ENDIF
          ENDDO
@@ -130,12 +169,7 @@ SUBROUTINE OBCS_APPLY_UVICE(
          DO j=1-OLy,sNy+OLy
           Iobc = OB_Iw(j,bi,bj)
           IF ( Iobc.NE.OB_indexNone ) THEN
-#  ifdef SEAICE_OBCS_NEUMANN
-C     Neumann boundary conditions
-           vFld(Iobc,j,bi,bj) = vFld(Iobc+1,j,bi,bj)
-#  else
            vFld(Iobc,j,bi,bj) = OBWvice(j,bi,bj)
-#  endif
      &                        *seaiceMaskV(Iobc,j,bi,bj)
           ENDIF
          ENDDO
@@ -151,15 +185,12 @@ SUBROUTINE OBCS_APPLY_UVICE(
          DO i=1-OLx,sNx+OLx
           Jobc = OB_Jn(i,bi,bj)
           IF ( Jobc.NE.OB_indexNone ) THEN
-#  ifdef SEAICE_OBCS_NEUMANN
-C     Neumann boundary conditions
-           vFld(i,Jobc,bi,bj)   = vFld(i,Jobc-1,bi,bj)
-#  else
            vFld(i,Jobc,bi,bj)   = OBNvice(i,bi,bj)
-#  endif
      &                          *seaiceMaskV(i,Jobc,bi,bj)
            IF ( uvIceApplyFac.GE.0. )
-     &     vFld(i,Jobc+1,bi,bj) = vFld(i,Jobc,bi,bj)*uvIceApplyFac
+     &     vFld(i,Jobc+1,bi,bj) = OBNvice(i,bi,bj)
+     &                          *seaiceMaskV(i,Jobc,bi,bj)
+     &                          *uvIceApplyFac
           ENDIF
          ENDDO
         ENDIF
@@ -171,15 +202,12 @@ SUBROUTINE OBCS_APPLY_UVICE(
          DO i=1-OLx,sNx+OLx
           Jobc = OB_Js(i,bi,bj)
           IF ( Jobc.NE.OB_indexNone ) THEN
-#  ifdef SEAICE_OBCS_NEUMANN
-C     Neumann boundary conditions
-           vFld(i,Jobc+1,bi,bj) = vFld(i,Jobc+2,bi,bj)
-#  else
            vFld(i,Jobc+1,bi,bj) = OBSvice(i,bi,bj)
-#  endif
      &                          *seaiceMaskV(i,Jobc+1,bi,bj)
            IF ( uvIceApplyFac.GE.0. )
-     &     vFld(i,Jobc,bi,bj)   = vFld(i,Jobc+1,bi,bj)*uvIceApplyFac
+     &     vFld(i,Jobc,bi,bj)   = OBSvice(i,bi,bj)
+     &                          *seaiceMaskV(i,Jobc+1,bi,bj)
+     &                          *uvIceApplyFac
           ENDIF
          ENDDO
         ENDIF
@@ -192,15 +220,12 @@ SUBROUTINE OBCS_APPLY_UVICE(
          DO j=1-OLy,sNy+OLy
           Iobc = OB_Ie(j,bi,bj)
           IF ( Iobc.NE.OB_indexNone ) THEN
-#  ifdef SEAICE_OBCS_NEUMANN
-C     Neumann boundary conditions
-           uFld(Iobc,j,bi,bj)   = uFld(Iobc-1,j,bi,bj)
-#  else
            uFld(Iobc,j,bi,bj)   = OBEuice(j,bi,bj)
-#  endif
      &                          *seaiceMaskU(Iobc,j,bi,bj)
            IF ( uvIceApplyFac.GE.0. )
-     &     uFld(Iobc+1,j,bi,bj) = uFld(Iobc,j,bi,bj)*uvIceApplyFac
+     &     uFld(Iobc+1,j,bi,bj) = OBEuice(j,bi,bj)
+     &                          *seaiceMaskU(Iobc,j,bi,bj)
+     &                          *uvIceApplyFac
           ENDIF
          ENDDO
         ENDIF
@@ -212,15 +237,12 @@ SUBROUTINE OBCS_APPLY_UVICE(
          DO j=1-OLy,sNy+OLy
           Iobc = OB_Iw(j,bi,bj)
           IF ( Iobc.NE.OB_indexNone ) THEN
-#  ifdef SEAICE_OBCS_NEUMANN
-C     Neumann boundary conditions
-           uFld(Iobc+1,j,bi,bj) = uFld(Iobc+2,j,bi,bj)
-#  else
            uFld(Iobc+1,j,bi,bj) = OBWuice(j,bi,bj)
-#  endif
      &                          *seaiceMaskU(Iobc+1,j,bi,bj)
            IF ( uvIceApplyFac.GE.0. )
-     &     uFld(Iobc,j,bi,bj)   = uFld(Iobc+1,j,bi,bj)*uvIceApplyFac
+     &     uFld(Iobc,j,bi,bj)   = OBWuice(j,bi,bj)
+     &                          *seaiceMaskU(Iobc+1,j,bi,bj)
+     &                          *uvIceApplyFac
           ENDIF
          ENDDO
         ENDIF

pkg/obcs/obcs_check.F

@@ -207,6 +207,11 @@ SUBROUTINE OBCS_CHECK( myThid )
      &  'Vrelaxobcsbound =',
      & ' /* boundary relaxation time scale, v-velocity ( s ). */')
       ENDIF
+      CALL WRITE_0D_L( useSeaiceSponge, INDEX_NONE,
+     & 'useSeaiceSponge =', '  /* use sponge for sea ice variables */')
+      CALL WRITE_0D_L( useSeaiceNeumann, INDEX_NONE,
+     & 'useSeaiceNeumann =',
+     & '  /* use Neumann conditions for sea ice variables */')
 C
       ln = ILNBLNK(insideOBmaskFile)
       IF ( ln.GT.0 ) THEN

pkg/obcs/obcs_readparms.F

@@ -159,7 +159,7 @@ SUBROUTINE OBCS_READPARMS( myThid )
      &          OBCS_u1_adv_Tr,
      &          OBNptrFile,OBSptrFile,OBEptrFile,OBWptrFile,
 #endif
-     &          useOBCSsponge, useSeaiceSponge,
+     &          useOBCSsponge, useSeaiceSponge, useSeaiceNeumann,
      &          OBCSsponge_N , OBCSsponge_S,
      &          OBCSsponge_E, OBCSsponge_W,
      &          OBCSsponge_UatNS, OBCSsponge_UatEW,
@@ -256,6 +256,7 @@ SUBROUTINE OBCS_READPARMS( myThid )
       useStevensAdvection=.TRUE.
       useOBCSsponge      =.FALSE.
       useSeaiceSponge    =.FALSE.
+      useSeaiceNeumann   =.FALSE.
       OBCSsponge_N       =.TRUE.
       OBCSsponge_S       =.TRUE.
       OBCSsponge_E       =.TRUE.