From d8ac131a55c5665009e86100885751c069b93e9c Mon Sep 17 00:00:00 2001 From: mjlosch Date: Tue, 16 Apr 2024 09:56:55 +0200 Subject: [PATCH] turn Neumann BC CPP-flag into runtime flag useSeaiceNeumann first attempt --- pkg/obcs/OBCS_OPTIONS.h | 3 - pkg/obcs/OBCS_SEAICE.h | 6 +- pkg/obcs/obcs_apply_seaice.F | 24 ++++---- pkg/obcs/obcs_apply_uvice.F | 110 +++++++++++++++++++++-------------- pkg/obcs/obcs_check.F | 5 ++ pkg/obcs/obcs_readparms.F | 3 +- 6 files changed, 89 insertions(+), 62 deletions(-) diff --git a/pkg/obcs/OBCS_OPTIONS.h b/pkg/obcs/OBCS_OPTIONS.h index cfca3f27ec..ca7d6ee6d2 100644 --- a/pkg/obcs/OBCS_OPTIONS.h +++ b/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 diff --git a/pkg/obcs/OBCS_SEAICE.h b/pkg/obcs/OBCS_SEAICE.h index 416e29ff3f..7a84b8a2e9 100644 --- a/pkg/obcs/OBCS_SEAICE.h +++ b/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 diff --git a/pkg/obcs/obcs_apply_seaice.F b/pkg/obcs/obcs_apply_seaice.F index a27d8b8e20..b51ff59852 100644 --- a/pkg/obcs/obcs_apply_seaice.F +++ b/pkg/obcs/obcs_apply_seaice.F @@ -57,7 +57,7 @@ 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) @@ -65,14 +65,14 @@ SUBROUTINE OBCS_APPLY_SEAICE( myThid ) # 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,7 +84,7 @@ 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) @@ -92,14 +92,14 @@ SUBROUTINE OBCS_APPLY_SEAICE( myThid ) # 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,7 +111,7 @@ 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) @@ -119,14 +119,14 @@ SUBROUTINE OBCS_APPLY_SEAICE( myThid ) # 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,7 +138,7 @@ 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) @@ -146,14 +146,14 @@ SUBROUTINE OBCS_APPLY_SEAICE( myThid ) # 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 diff --git a/pkg/obcs/obcs_apply_uvice.F b/pkg/obcs/obcs_apply_uvice.F index fbed02e2d5..f71075ea86 100644 --- a/pkg/obcs/obcs_apply_uvice.F +++ b/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 diff --git a/pkg/obcs/obcs_check.F b/pkg/obcs/obcs_check.F index f07d837dae..0d8972ee94 100644 --- a/pkg/obcs/obcs_check.F +++ b/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 diff --git a/pkg/obcs/obcs_readparms.F b/pkg/obcs/obcs_readparms.F index 990cea5e3a..9ea691f011 100644 --- a/pkg/obcs/obcs_readparms.F +++ b/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.