diff --git a/exps/NEP10.COBALT/MOM_parameter_doc.all b/exps/NEP10.COBALT/MOM_parameter_doc.all
index b678c0a16..902d384a4 100644
--- a/exps/NEP10.COBALT/MOM_parameter_doc.all
+++ b/exps/NEP10.COBALT/MOM_parameter_doc.all
@@ -382,7 +382,7 @@ OBC_TIDE_ADD_NODAL = True       !   [Boolean] default = False
                                 ! If true, include 18.6 year nodal modulation in the boundary tidal forcing.
 OBC_TIDE_REF_DATE = 1993, 1, 1  !
                                 ! Reference date to use for tidal calculations and equilibrium phase.
-OBC_TIDE_NODAL_REF_DATE = 1993, 7, 2 ! default = 0
+OBC_TIDE_NODAL_REF_DATE = 1993, 7, 2 ! default = 0, 0, 0
                                 ! Fixed reference date to use for nodal modulation of boundary tides.
 TIDE_M2_FREQ = 1.405189E-04     !   [rad s-1] default = 1.405189E-04
                                 ! Frequency of the M2 tidal constituent. This is only used if TIDES and TIDE_M2
@@ -1138,7 +1138,7 @@ TIDAL_SAL_FROM_FILE = False     !   [Boolean] default = False
 USE_PREVIOUS_TIDES = False      !   [Boolean] default = False
                                 ! If true, use the SAL from the previous iteration of the tides to facilitate
                                 ! convergent iteration. This is only used if TIDES is true.
-TIDE_REF_DATE = 1993, 1, 1      ! default = 0
+TIDE_REF_DATE = 1993, 1, 1      ! default = 0, 0, 0
                                 ! Year,month,day to use as reference date for tidal forcing. If not specified,
                                 ! defaults to 0.
 TIDE_USE_EQ_PHASE = True        !   [Boolean] default = False
@@ -1951,6 +1951,11 @@ MKE_TO_TKE_EFFIC = 0.0          !   [nondim] default = 0.0
 TKE_DECAY = 0.001               !   [nondim] default = 2.5
                                 ! TKE_DECAY relates the vertical rate of decay of the TKE available for
                                 ! mechanical entrainment to the natural Ekman depth.
+DIRECT_EPBL_MIXING_CALC = False !   [Boolean] default = False
+                                ! If true and there is no conversion from mean kinetic energy to ePBL turbulent
+                                ! kinetic energy, use a direct calculation of the diffusivity that is supported
+                                ! by a given energy input instead of the more general but slower iterative
+                                ! solver.
 EPBL_MSTAR_SCHEME = "REICHL_H18" ! default = "CONSTANT"
                                 ! EPBL_MSTAR_SCHEME selects the method for setting mstar.  Valid values are:
                                 !    CONSTANT   - Use a fixed mstar given by MSTAR
@@ -1960,7 +1965,7 @@ MSTAR_CAP = 10.0                !   [nondim] default = -1.0
                                 ! If this value is positive, it sets the maximum value of mstar allowed in ePBL.
                                 ! (This is not used if EPBL_MSTAR_SCHEME = CONSTANT).
 RH18_MSTAR_CN1 = 0.275          !   [nondim] default = 0.275
-                                ! MSTAR_N coefficient 1 (outter-most coefficient for fit). The value of 0.275 is
+                                ! MSTAR_N coefficient 1 (outer-most coefficient for fit). The value of 0.275 is
                                 ! given in RH18.  Increasing this coefficient increases MSTAR for all values of
                                 ! Hf/ust, but more effectively at low values (weakly developed OSBLs).
 RH18_MSTAR_CN2 = 8.0            !   [nondim] default = 8.0
@@ -2005,10 +2010,13 @@ EPBL_MLD_BISECTION = True       !   [Boolean] default = False
                                 ! mixed layer depth.  Otherwise use the false position after a maximum and
                                 ! minimum bound have been evaluated and the returned value or bisection before
                                 ! this.
+EPBL_MLD_ITER_BUG = True        !   [Boolean] default = True
+                                ! If true, use buggy logic that gives the wrong bounds for the next iteration
+                                ! when successive guesses increase by exactly EPBL_MLD_TOLERANCE.
 EPBL_MLD_MAX_ITS = 20           ! default = 20
                                 ! The maximum number of iterations that can be used to find a self-consistent
-                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number
-                                ! iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
+                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number of
+                                ! iterations needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
 EPBL_MIN_MIX_LEN = 0.0          !   [meter] default = 0.0
                                 ! The minimum mixing length scale that will be used by ePBL.  The default (0)
                                 ! does not set a minimum.
@@ -2031,6 +2039,9 @@ EPBL_VEL_SCALE_FACTOR = 0.5477  !   [nondim] default = 1.0
                                 ! the PBL diffusivity.
 VSTAR_SURF_FAC = 1.8258         !   [nondim] default = 1.2
                                 ! The proportionality times ustar to set vstar at the surface.
+EPBL_BBL_EFFIC = 0.0            !   [nondim] default = 0.0
+                                ! The efficiency of bottom boundary layer mixing via ePBL.  Setting this to a
+                                ! value that is greater than 0 to enable bottom boundary layer mixing from EPBL.
 USE_LA_LI2016 = True            !   [Boolean] default = False
                                 ! A logical to use the Li et al. 2016 (submitted) formula to determine the
                                 ! Langmuir number.
@@ -2040,12 +2051,12 @@ EPBL_LANGMUIR_SCHEME = "ADDITIVE" ! default = "NONE"
                                 !    NONE     - Do not do any extra mixing due to Langmuir turbulence
                                 !    RESCALE  - Use a multiplicative rescaling of mstar to account for Langmuir
                                 !      turbulence
-                                !    ADDITIVE - Add a Langmuir turblence contribution to mstar to other
+                                !    ADDITIVE - Add a Langmuir turbulence contribution to mstar to other
                                 !      contributions
 LT_ENHANCE_COEF = 0.1056        !   [nondim] default = 0.447
                                 ! Coefficient for Langmuir enhancement of mstar
 LT_ENHANCE_EXP = -1.0           !   [nondim] default = -1.33
-                                ! Exponent for Langmuir enhancementt of mstar
+                                ! Exponent for Langmuir enhancement of mstar
 LT_MOD_LAC1 = 0.0               !   [nondim] default = -0.87
                                 ! Coefficient for modification of Langmuir number due to MLD approaching Ekman
                                 ! depth.
@@ -2061,6 +2072,10 @@ LT_MOD_LAC4 = 0.8               !   [nondim] default = 0.95
 LT_MOD_LAC5 = 0.8               !   [nondim] default = 0.95
                                 ! Coefficient for modification of Langmuir number due to ratio of Ekman to
                                 ! unstable Obukhov depth.
+EPBL_OPTIONS_DIFF = 0           ! default = 0
+                                ! If positive, this is a coded integer indicating a pair of settings whose
+                                ! differences are diagnosed in a passive diagnostic mode  via extra calls to
+                                ! ePBL_column.  If this is 0 or negative no extra calls occur.
 !EPBL_USTAR_MIN = 1.45842E-18   !   [m s-1]
                                 ! The (tiny) minimum friction velocity used within the ePBL code, derived from
                                 ! OMEGA and ANGSTROM.
diff --git a/exps/NEP10.COBALT/MOM_parameter_doc.layout b/exps/NEP10.COBALT/MOM_parameter_doc.layout
index d510cf4da..5d0f895c4 100644
--- a/exps/NEP10.COBALT/MOM_parameter_doc.layout
+++ b/exps/NEP10.COBALT/MOM_parameter_doc.layout
@@ -46,7 +46,7 @@ NJPROC = 56                     !
                                 ! in MOM_memory.h at compile time.
 LAYOUT = 20, 56                 !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 1, 1                ! default = 1
+IO_LAYOUT = 1, 1                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 
diff --git a/exps/NEP10.COBALT/MOM_parameter_doc.short b/exps/NEP10.COBALT/MOM_parameter_doc.short
index 12acf1559..b2e413800 100644
--- a/exps/NEP10.COBALT/MOM_parameter_doc.short
+++ b/exps/NEP10.COBALT/MOM_parameter_doc.short
@@ -167,7 +167,7 @@ OBC_TIDE_ADD_NODAL = True       !   [Boolean] default = False
                                 ! If true, include 18.6 year nodal modulation in the boundary tidal forcing.
 OBC_TIDE_REF_DATE = 1993, 1, 1  !
                                 ! Reference date to use for tidal calculations and equilibrium phase.
-OBC_TIDE_NODAL_REF_DATE = 1993, 7, 2 ! default = 0
+OBC_TIDE_NODAL_REF_DATE = 1993, 7, 2 ! default = 0, 0, 0
                                 ! Fixed reference date to use for nodal modulation of boundary tides.
 
 ! === module MOM_verticalGrid ===
@@ -391,7 +391,7 @@ TIDE_MF = True                  !   [Boolean] default = False
 TIDE_MM = True                  !   [Boolean] default = False
                                 ! If true, apply tidal momentum forcing at the MM frequency. This is only used
                                 ! if TIDES is true.
-TIDE_REF_DATE = 1993, 1, 1      ! default = 0
+TIDE_REF_DATE = 1993, 1, 1      ! default = 0, 0, 0
                                 ! Year,month,day to use as reference date for tidal forcing. If not specified,
                                 ! defaults to 0.
 TIDE_USE_EQ_PHASE = True        !   [Boolean] default = False
@@ -653,12 +653,12 @@ EPBL_LANGMUIR_SCHEME = "ADDITIVE" ! default = "NONE"
                                 !    NONE     - Do not do any extra mixing due to Langmuir turbulence
                                 !    RESCALE  - Use a multiplicative rescaling of mstar to account for Langmuir
                                 !      turbulence
-                                !    ADDITIVE - Add a Langmuir turblence contribution to mstar to other
+                                !    ADDITIVE - Add a Langmuir turbulence contribution to mstar to other
                                 !      contributions
 LT_ENHANCE_COEF = 0.1056        !   [nondim] default = 0.447
                                 ! Coefficient for Langmuir enhancement of mstar
 LT_ENHANCE_EXP = -1.0           !   [nondim] default = -1.33
-                                ! Exponent for Langmuir enhancementt of mstar
+                                ! Exponent for Langmuir enhancement of mstar
 LT_MOD_LAC1 = 0.0               !   [nondim] default = -0.87
                                 ! Coefficient for modification of Langmuir number due to MLD approaching Ekman
                                 ! depth.
diff --git a/exps/NEP10.COBALT/SIS_parameter_doc.layout b/exps/NEP10.COBALT/SIS_parameter_doc.layout
index 403ff12eb..cb376982a 100644
--- a/exps/NEP10.COBALT/SIS_parameter_doc.layout
+++ b/exps/NEP10.COBALT/SIS_parameter_doc.layout
@@ -44,7 +44,7 @@ NJPROC = 56                     !
                                 ! in SIS2_memory.h at compile time.
 LAYOUT = 20, 56                 !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 1, 1                ! default = 1
+IO_LAYOUT = 1, 1                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 NIBLOCK = 1                     ! default = 1
diff --git a/exps/NWA12.COBALT/MOM_parameter_doc.all b/exps/NWA12.COBALT/MOM_parameter_doc.all
index 1f8653596..54831c1e5 100644
--- a/exps/NWA12.COBALT/MOM_parameter_doc.all
+++ b/exps/NWA12.COBALT/MOM_parameter_doc.all
@@ -376,7 +376,7 @@ OBC_TIDE_ADD_NODAL = True       !   [Boolean] default = False
                                 ! If true, include 18.6 year nodal modulation in the boundary tidal forcing.
 OBC_TIDE_REF_DATE = 1993, 1, 1  !
                                 ! Reference date to use for tidal calculations and equilibrium phase.
-OBC_TIDE_NODAL_REF_DATE = 1993, 7, 2 ! default = 0
+OBC_TIDE_NODAL_REF_DATE = 1993, 7, 2 ! default = 0, 0, 0
                                 ! Fixed reference date to use for nodal modulation of boundary tides.
 TIDE_M2_FREQ = 1.405189E-04     !   [rad s-1] default = 1.405189E-04
                                 ! Frequency of the M2 tidal constituent. This is only used if TIDES and TIDE_M2
@@ -1150,7 +1150,7 @@ TIDAL_SAL_FROM_FILE = False     !   [Boolean] default = False
 USE_PREVIOUS_TIDES = False      !   [Boolean] default = False
                                 ! If true, use the SAL from the previous iteration of the tides to facilitate
                                 ! convergent iteration. This is only used if TIDES is true.
-TIDE_REF_DATE = 1993, 1, 1      ! default = 0
+TIDE_REF_DATE = 1993, 1, 1      ! default = 0, 0, 0
                                 ! Year,month,day to use as reference date for tidal forcing. If not specified,
                                 ! defaults to 0.
 TIDE_USE_EQ_PHASE = True        !   [Boolean] default = False
@@ -1983,6 +1983,11 @@ MKE_TO_TKE_EFFIC = 0.0          !   [nondim] default = 0.0
 TKE_DECAY = 0.001               !   [nondim] default = 2.5
                                 ! TKE_DECAY relates the vertical rate of decay of the TKE available for
                                 ! mechanical entrainment to the natural Ekman depth.
+DIRECT_EPBL_MIXING_CALC = False !   [Boolean] default = False
+                                ! If true and there is no conversion from mean kinetic energy to ePBL turbulent
+                                ! kinetic energy, use a direct calculation of the diffusivity that is supported
+                                ! by a given energy input instead of the more general but slower iterative
+                                ! solver.
 EPBL_MSTAR_SCHEME = "REICHL_H18" ! default = "CONSTANT"
                                 ! EPBL_MSTAR_SCHEME selects the method for setting mstar.  Valid values are:
                                 !    CONSTANT   - Use a fixed mstar given by MSTAR
@@ -1992,7 +1997,7 @@ MSTAR_CAP = 10.0                !   [nondim] default = -1.0
                                 ! If this value is positive, it sets the maximum value of mstar allowed in ePBL.
                                 ! (This is not used if EPBL_MSTAR_SCHEME = CONSTANT).
 RH18_MSTAR_CN1 = 0.275          !   [nondim] default = 0.275
-                                ! MSTAR_N coefficient 1 (outter-most coefficient for fit). The value of 0.275 is
+                                ! MSTAR_N coefficient 1 (outer-most coefficient for fit). The value of 0.275 is
                                 ! given in RH18.  Increasing this coefficient increases MSTAR for all values of
                                 ! Hf/ust, but more effectively at low values (weakly developed OSBLs).
 RH18_MSTAR_CN2 = 8.0            !   [nondim] default = 8.0
@@ -2037,10 +2042,13 @@ EPBL_MLD_BISECTION = True       !   [Boolean] default = False
                                 ! mixed layer depth.  Otherwise use the false position after a maximum and
                                 ! minimum bound have been evaluated and the returned value or bisection before
                                 ! this.
+EPBL_MLD_ITER_BUG = True        !   [Boolean] default = True
+                                ! If true, use buggy logic that gives the wrong bounds for the next iteration
+                                ! when successive guesses increase by exactly EPBL_MLD_TOLERANCE.
 EPBL_MLD_MAX_ITS = 40           ! default = 20
                                 ! The maximum number of iterations that can be used to find a self-consistent
-                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number
-                                ! iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
+                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number of
+                                ! iterations needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
 EPBL_MIN_MIX_LEN = 0.0          !   [meter] default = 0.0
                                 ! The minimum mixing length scale that will be used by ePBL.  The default (0)
                                 ! does not set a minimum.
@@ -2063,6 +2071,9 @@ EPBL_VEL_SCALE_FACTOR = 0.5477  !   [nondim] default = 1.0
                                 ! the PBL diffusivity.
 VSTAR_SURF_FAC = 1.8258         !   [nondim] default = 1.2
                                 ! The proportionality times ustar to set vstar at the surface.
+EPBL_BBL_EFFIC = 0.0            !   [nondim] default = 0.0
+                                ! The efficiency of bottom boundary layer mixing via ePBL.  Setting this to a
+                                ! value that is greater than 0 to enable bottom boundary layer mixing from EPBL.
 USE_LA_LI2016 = True            !   [Boolean] default = False
                                 ! A logical to use the Li et al. 2016 (submitted) formula to determine the
                                 ! Langmuir number.
@@ -2072,12 +2083,12 @@ EPBL_LANGMUIR_SCHEME = "ADDITIVE" ! default = "NONE"
                                 !    NONE     - Do not do any extra mixing due to Langmuir turbulence
                                 !    RESCALE  - Use a multiplicative rescaling of mstar to account for Langmuir
                                 !      turbulence
-                                !    ADDITIVE - Add a Langmuir turblence contribution to mstar to other
+                                !    ADDITIVE - Add a Langmuir turbulence contribution to mstar to other
                                 !      contributions
 LT_ENHANCE_COEF = 0.1056        !   [nondim] default = 0.447
                                 ! Coefficient for Langmuir enhancement of mstar
 LT_ENHANCE_EXP = -1.0           !   [nondim] default = -1.33
-                                ! Exponent for Langmuir enhancementt of mstar
+                                ! Exponent for Langmuir enhancement of mstar
 LT_MOD_LAC1 = 0.0               !   [nondim] default = -0.87
                                 ! Coefficient for modification of Langmuir number due to MLD approaching Ekman
                                 ! depth.
@@ -2093,6 +2104,10 @@ LT_MOD_LAC4 = 0.8               !   [nondim] default = 0.95
 LT_MOD_LAC5 = 0.8               !   [nondim] default = 0.95
                                 ! Coefficient for modification of Langmuir number due to ratio of Ekman to
                                 ! unstable Obukhov depth.
+EPBL_OPTIONS_DIFF = 0           ! default = 0
+                                ! If positive, this is a coded integer indicating a pair of settings whose
+                                ! differences are diagnosed in a passive diagnostic mode  via extra calls to
+                                ! ePBL_column.  If this is 0 or negative no extra calls occur.
 !EPBL_USTAR_MIN = 1.45842E-18   !   [m s-1]
                                 ! The (tiny) minimum friction velocity used within the ePBL code, derived from
                                 ! OMEGA and ANGSTROM.
diff --git a/exps/NWA12.COBALT/MOM_parameter_doc.layout b/exps/NWA12.COBALT/MOM_parameter_doc.layout
index 7edeab8e1..401aade9b 100644
--- a/exps/NWA12.COBALT/MOM_parameter_doc.layout
+++ b/exps/NWA12.COBALT/MOM_parameter_doc.layout
@@ -46,7 +46,7 @@ NJPROC = 25                     !
                                 ! in MOM_memory.h at compile time.
 LAYOUT = 25, 25                 !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 1, 1                ! default = 1
+IO_LAYOUT = 1, 1                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 
diff --git a/exps/NWA12.COBALT/MOM_parameter_doc.short b/exps/NWA12.COBALT/MOM_parameter_doc.short
index ecc280683..23adab7d6 100644
--- a/exps/NWA12.COBALT/MOM_parameter_doc.short
+++ b/exps/NWA12.COBALT/MOM_parameter_doc.short
@@ -163,7 +163,7 @@ OBC_TIDE_ADD_NODAL = True       !   [Boolean] default = False
                                 ! If true, include 18.6 year nodal modulation in the boundary tidal forcing.
 OBC_TIDE_REF_DATE = 1993, 1, 1  !
                                 ! Reference date to use for tidal calculations and equilibrium phase.
-OBC_TIDE_NODAL_REF_DATE = 1993, 7, 2 ! default = 0
+OBC_TIDE_NODAL_REF_DATE = 1993, 7, 2 ! default = 0, 0, 0
                                 ! Fixed reference date to use for nodal modulation of boundary tides.
 
 ! === module MOM_verticalGrid ===
@@ -373,7 +373,7 @@ TIDE_MF = True                  !   [Boolean] default = False
 TIDE_MM = True                  !   [Boolean] default = False
                                 ! If true, apply tidal momentum forcing at the MM frequency. This is only used
                                 ! if TIDES is true.
-TIDE_REF_DATE = 1993, 1, 1      ! default = 0
+TIDE_REF_DATE = 1993, 1, 1      ! default = 0, 0, 0
                                 ! Year,month,day to use as reference date for tidal forcing. If not specified,
                                 ! defaults to 0.
 TIDE_USE_EQ_PHASE = True        !   [Boolean] default = False
@@ -601,8 +601,8 @@ EPBL_MLD_BISECTION = True       !   [Boolean] default = False
                                 ! this.
 EPBL_MLD_MAX_ITS = 40           ! default = 20
                                 ! The maximum number of iterations that can be used to find a self-consistent
-                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number
-                                ! iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
+                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number of
+                                ! iterations needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
 MIX_LEN_EXPONENT = 1.0          !   [nondim] default = 2.0
                                 ! The exponent applied to the ratio of the distance to the MLD and the MLD depth
                                 ! which determines the shape of the mixing length. This is only used if
@@ -631,12 +631,12 @@ EPBL_LANGMUIR_SCHEME = "ADDITIVE" ! default = "NONE"
                                 !    NONE     - Do not do any extra mixing due to Langmuir turbulence
                                 !    RESCALE  - Use a multiplicative rescaling of mstar to account for Langmuir
                                 !      turbulence
-                                !    ADDITIVE - Add a Langmuir turblence contribution to mstar to other
+                                !    ADDITIVE - Add a Langmuir turbulence contribution to mstar to other
                                 !      contributions
 LT_ENHANCE_COEF = 0.1056        !   [nondim] default = 0.447
                                 ! Coefficient for Langmuir enhancement of mstar
 LT_ENHANCE_EXP = -1.0           !   [nondim] default = -1.33
-                                ! Exponent for Langmuir enhancementt of mstar
+                                ! Exponent for Langmuir enhancement of mstar
 LT_MOD_LAC1 = 0.0               !   [nondim] default = -0.87
                                 ! Coefficient for modification of Langmuir number due to MLD approaching Ekman
                                 ! depth.
diff --git a/exps/NWA12.COBALT/SIS_parameter_doc.layout b/exps/NWA12.COBALT/SIS_parameter_doc.layout
index 970cb7525..31fd08682 100644
--- a/exps/NWA12.COBALT/SIS_parameter_doc.layout
+++ b/exps/NWA12.COBALT/SIS_parameter_doc.layout
@@ -44,7 +44,7 @@ NJPROC = 25                     !
                                 ! in SIS2_memory.h at compile time.
 LAYOUT = 25, 25                 !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 1, 1                ! default = 1
+IO_LAYOUT = 1, 1                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 NIBLOCK = 1                     ! default = 1
diff --git a/exps/NWA12.tidesonly/MOM_parameter_doc.all b/exps/NWA12.tidesonly/MOM_parameter_doc.all
index f951f5684..35cedfcb3 100644
--- a/exps/NWA12.tidesonly/MOM_parameter_doc.all
+++ b/exps/NWA12.tidesonly/MOM_parameter_doc.all
@@ -343,7 +343,7 @@ OBC_TIDE_ADD_NODAL = False      !   [Boolean] default = False
                                 ! If true, include 18.6 year nodal modulation in the boundary tidal forcing.
 OBC_TIDE_REF_DATE = 1993, 1, 1  !
                                 ! Reference date to use for tidal calculations and equilibrium phase.
-OBC_TIDE_NODAL_REF_DATE = 0, 0, 0 ! default = 0
+OBC_TIDE_NODAL_REF_DATE = 0, 0, 0 ! default = 0, 0, 0
                                 ! Fixed reference date to use for nodal modulation of boundary tides.
 TIDE_M2_FREQ = 1.405189E-04     !   [rad s-1] default = 1.405189E-04
                                 ! Frequency of the M2 tidal constituent. This is only used if TIDES and TIDE_M2
@@ -1051,7 +1051,7 @@ TIDAL_SAL_FROM_FILE = False     !   [Boolean] default = False
 USE_PREVIOUS_TIDES = False      !   [Boolean] default = False
                                 ! If true, use the SAL from the previous iteration of the tides to facilitate
                                 ! convergent iteration. This is only used if TIDES is true.
-TIDE_REF_DATE = 0, 0, 0         ! default = 0
+TIDE_REF_DATE = 0, 0, 0         ! default = 0, 0, 0
                                 ! Year,month,day to use as reference date for tidal forcing. If not specified,
                                 ! defaults to 0.
 TIDE_USE_EQ_PHASE = False       !   [Boolean] default = False
@@ -1559,8 +1559,7 @@ BUOY_CONFIG = "zero"            ! default = "zero"
 WIND_CONFIG = "zero"            ! default = "zero"
                                 ! The character string that indicates how wind forcing is specified.  Valid
                                 ! options include (file), (data_override), (2gyre), (1gyre), (gyres), (zero),
-                                ! (const), (Neverworld), (scurves), (ideal_hurr), (SCM_ideal_hurr),
-                                ! (SCM_CVmix_tests) and (USER).
+                                ! (const), (Neverworld), (scurves), (ideal_hurr), (SCM_CVmix_tests) and (USER).
 RESTOREBUOY = False             !   [Boolean] default = False
                                 ! If true, the buoyancy fluxes drive the model back toward some specified
                                 ! surface state with a rate given by FLUXCONST.
diff --git a/exps/NWA12.tidesonly/MOM_parameter_doc.layout b/exps/NWA12.tidesonly/MOM_parameter_doc.layout
index b7e94719d..6ea7c5e31 100644
--- a/exps/NWA12.tidesonly/MOM_parameter_doc.layout
+++ b/exps/NWA12.tidesonly/MOM_parameter_doc.layout
@@ -46,7 +46,7 @@ NJPROC = 2                      !
                                 ! in MOM_memory.h at compile time.
 LAYOUT = 2, 2                   !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 1, 1                ! default = 1
+IO_LAYOUT = 1, 1                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 
diff --git a/exps/OM4.single_column.COBALT/MOM_parameter_doc.all b/exps/OM4.single_column.COBALT/MOM_parameter_doc.all
index 5180ccaa2..82c87625e 100644
--- a/exps/OM4.single_column.COBALT/MOM_parameter_doc.all
+++ b/exps/OM4.single_column.COBALT/MOM_parameter_doc.all
@@ -1630,6 +1630,11 @@ MKE_TO_TKE_EFFIC = 0.0          !   [nondim] default = 0.0
 TKE_DECAY = 0.01                !   [nondim] default = 2.5
                                 ! TKE_DECAY relates the vertical rate of decay of the TKE available for
                                 ! mechanical entrainment to the natural Ekman depth.
+DIRECT_EPBL_MIXING_CALC = False !   [Boolean] default = False
+                                ! If true and there is no conversion from mean kinetic energy to ePBL turbulent
+                                ! kinetic energy, use a direct calculation of the diffusivity that is supported
+                                ! by a given energy input instead of the more general but slower iterative
+                                ! solver.
 EPBL_MSTAR_SCHEME = "CONSTANT"  ! default = "CONSTANT"
                                 ! EPBL_MSTAR_SCHEME selects the method for setting mstar.  Valid values are:
                                 !    CONSTANT   - Use a fixed mstar given by MSTAR
@@ -1663,10 +1668,13 @@ EPBL_MLD_BISECTION = False      !   [Boolean] default = False
                                 ! mixed layer depth.  Otherwise use the false position after a maximum and
                                 ! minimum bound have been evaluated and the returned value or bisection before
                                 ! this.
+EPBL_MLD_ITER_BUG = True        !   [Boolean] default = True
+                                ! If true, use buggy logic that gives the wrong bounds for the next iteration
+                                ! when successive guesses increase by exactly EPBL_MLD_TOLERANCE.
 EPBL_MLD_MAX_ITS = 20           ! default = 20
                                 ! The maximum number of iterations that can be used to find a self-consistent
-                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number
-                                ! iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
+                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number of
+                                ! iterations needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
 EPBL_MIN_MIX_LEN = 0.0          !   [meter] default = 0.0
                                 ! The minimum mixing length scale that will be used by ePBL.  The default (0)
                                 ! does not set a minimum.
@@ -1689,6 +1697,9 @@ EPBL_VEL_SCALE_FACTOR = 1.0     !   [nondim] default = 1.0
                                 ! the PBL diffusivity.
 VSTAR_SURF_FAC = 1.2            !   [nondim] default = 1.2
                                 ! The proportionality times ustar to set vstar at the surface.
+EPBL_BBL_EFFIC = 0.0            !   [nondim] default = 0.0
+                                ! The efficiency of bottom boundary layer mixing via ePBL.  Setting this to a
+                                ! value that is greater than 0 to enable bottom boundary layer mixing from EPBL.
 USE_LA_LI2016 = True            !   [Boolean] default = False
                                 ! A logical to use the Li et al. 2016 (submitted) formula to determine the
                                 ! Langmuir number.
@@ -1698,12 +1709,12 @@ EPBL_LANGMUIR_SCHEME = "ADDITIVE" ! default = "NONE"
                                 !    NONE     - Do not do any extra mixing due to Langmuir turbulence
                                 !    RESCALE  - Use a multiplicative rescaling of mstar to account for Langmuir
                                 !      turbulence
-                                !    ADDITIVE - Add a Langmuir turblence contribution to mstar to other
+                                !    ADDITIVE - Add a Langmuir turbulence contribution to mstar to other
                                 !      contributions
 LT_ENHANCE_COEF = 0.044         !   [nondim] default = 0.447
                                 ! Coefficient for Langmuir enhancement of mstar
 LT_ENHANCE_EXP = -1.5           !   [nondim] default = -1.33
-                                ! Exponent for Langmuir enhancementt of mstar
+                                ! Exponent for Langmuir enhancement of mstar
 LT_MOD_LAC1 = 0.0               !   [nondim] default = -0.87
                                 ! Coefficient for modification of Langmuir number due to MLD approaching Ekman
                                 ! depth.
@@ -1719,6 +1730,10 @@ LT_MOD_LAC4 = 0.0               !   [nondim] default = 0.95
 LT_MOD_LAC5 = 0.22              !   [nondim] default = 0.95
                                 ! Coefficient for modification of Langmuir number due to ratio of Ekman to
                                 ! unstable Obukhov depth.
+EPBL_OPTIONS_DIFF = 0           ! default = 0
+                                ! If positive, this is a coded integer indicating a pair of settings whose
+                                ! differences are diagnosed in a passive diagnostic mode  via extra calls to
+                                ! ePBL_column.  If this is 0 or negative no extra calls occur.
 !EPBL_USTAR_MIN = 1.45842E-18   !   [m s-1]
                                 ! The (tiny) minimum friction velocity used within the ePBL code, derived from
                                 ! OMEGA and ANGSTROM.
diff --git a/exps/OM4.single_column.COBALT/MOM_parameter_doc.layout b/exps/OM4.single_column.COBALT/MOM_parameter_doc.layout
index fda8b57a3..58effda3d 100644
--- a/exps/OM4.single_column.COBALT/MOM_parameter_doc.layout
+++ b/exps/OM4.single_column.COBALT/MOM_parameter_doc.layout
@@ -46,7 +46,7 @@ NJPROC = 1                      !
                                 ! in MOM_memory.h at compile time.
 LAYOUT = 1, 1                   !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 0, 0                ! default = 1
+IO_LAYOUT = 0, 0                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 
diff --git a/exps/OM4.single_column.COBALT/MOM_parameter_doc.short b/exps/OM4.single_column.COBALT/MOM_parameter_doc.short
index dd934ddfb..034a2aec1 100644
--- a/exps/OM4.single_column.COBALT/MOM_parameter_doc.short
+++ b/exps/OM4.single_column.COBALT/MOM_parameter_doc.short
@@ -420,12 +420,12 @@ EPBL_LANGMUIR_SCHEME = "ADDITIVE" ! default = "NONE"
                                 !    NONE     - Do not do any extra mixing due to Langmuir turbulence
                                 !    RESCALE  - Use a multiplicative rescaling of mstar to account for Langmuir
                                 !      turbulence
-                                !    ADDITIVE - Add a Langmuir turblence contribution to mstar to other
+                                !    ADDITIVE - Add a Langmuir turbulence contribution to mstar to other
                                 !      contributions
 LT_ENHANCE_COEF = 0.044         !   [nondim] default = 0.447
                                 ! Coefficient for Langmuir enhancement of mstar
 LT_ENHANCE_EXP = -1.5           !   [nondim] default = -1.33
-                                ! Exponent for Langmuir enhancementt of mstar
+                                ! Exponent for Langmuir enhancement of mstar
 LT_MOD_LAC1 = 0.0               !   [nondim] default = -0.87
                                 ! Coefficient for modification of Langmuir number due to MLD approaching Ekman
                                 ! depth.
diff --git a/exps/OM4.single_column.COBALT/SIS_parameter_doc.layout b/exps/OM4.single_column.COBALT/SIS_parameter_doc.layout
index 54dbd4cc4..597c30866 100644
--- a/exps/OM4.single_column.COBALT/SIS_parameter_doc.layout
+++ b/exps/OM4.single_column.COBALT/SIS_parameter_doc.layout
@@ -44,7 +44,7 @@ NJPROC = 1                      !
                                 ! in SIS2_memory.h at compile time.
 LAYOUT = 1, 1                   !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 0, 0                ! default = 1
+IO_LAYOUT = 0, 0                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 NIBLOCK = 1                     ! default = 1
diff --git a/exps/OM4p25.COBALT/MOM_parameter_doc.all b/exps/OM4p25.COBALT/MOM_parameter_doc.all
index 17cb59def..8513f297f 100644
--- a/exps/OM4p25.COBALT/MOM_parameter_doc.all
+++ b/exps/OM4p25.COBALT/MOM_parameter_doc.all
@@ -2148,6 +2148,11 @@ MKE_TO_TKE_EFFIC = 0.0          !   [nondim] default = 0.0
 TKE_DECAY = 0.01                !   [nondim] default = 2.5
                                 ! TKE_DECAY relates the vertical rate of decay of the TKE available for
                                 ! mechanical entrainment to the natural Ekman depth.
+DIRECT_EPBL_MIXING_CALC = False !   [Boolean] default = False
+                                ! If true and there is no conversion from mean kinetic energy to ePBL turbulent
+                                ! kinetic energy, use a direct calculation of the diffusivity that is supported
+                                ! by a given energy input instead of the more general but slower iterative
+                                ! solver.
 EPBL_MSTAR_SCHEME = "OM4"       ! default = "CONSTANT"
                                 ! EPBL_MSTAR_SCHEME selects the method for setting mstar.  Valid values are:
                                 !    CONSTANT   - Use a fixed mstar given by MSTAR
@@ -2187,10 +2192,13 @@ EPBL_MLD_BISECTION = False      !   [Boolean] default = False
                                 ! mixed layer depth.  Otherwise use the false position after a maximum and
                                 ! minimum bound have been evaluated and the returned value or bisection before
                                 ! this.
+EPBL_MLD_ITER_BUG = True        !   [Boolean] default = True
+                                ! If true, use buggy logic that gives the wrong bounds for the next iteration
+                                ! when successive guesses increase by exactly EPBL_MLD_TOLERANCE.
 EPBL_MLD_MAX_ITS = 20           ! default = 20
                                 ! The maximum number of iterations that can be used to find a self-consistent
-                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number
-                                ! iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
+                                ! mixed layer depth.  If EPBL_MLD_BISECTION is true, the maximum number of
+                                ! iterations needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.
 EPBL_MIN_MIX_LEN = 0.0          !   [meter] default = 0.0
                                 ! The minimum mixing length scale that will be used by ePBL.  The default (0)
                                 ! does not set a minimum.
@@ -2213,6 +2221,9 @@ EPBL_VEL_SCALE_FACTOR = 1.0     !   [nondim] default = 1.0
                                 ! the PBL diffusivity.
 VSTAR_SURF_FAC = 1.2            !   [nondim] default = 1.2
                                 ! The proportionality times ustar to set vstar at the surface.
+EPBL_BBL_EFFIC = 0.0            !   [nondim] default = 0.0
+                                ! The efficiency of bottom boundary layer mixing via ePBL.  Setting this to a
+                                ! value that is greater than 0 to enable bottom boundary layer mixing from EPBL.
 USE_LA_LI2016 = True            !   [Boolean] default = False
                                 ! A logical to use the Li et al. 2016 (submitted) formula to determine the
                                 ! Langmuir number.
@@ -2222,12 +2233,12 @@ EPBL_LANGMUIR_SCHEME = "ADDITIVE" ! default = "NONE"
                                 !    NONE     - Do not do any extra mixing due to Langmuir turbulence
                                 !    RESCALE  - Use a multiplicative rescaling of mstar to account for Langmuir
                                 !      turbulence
-                                !    ADDITIVE - Add a Langmuir turblence contribution to mstar to other
+                                !    ADDITIVE - Add a Langmuir turbulence contribution to mstar to other
                                 !      contributions
 LT_ENHANCE_COEF = 0.044         !   [nondim] default = 0.447
                                 ! Coefficient for Langmuir enhancement of mstar
 LT_ENHANCE_EXP = -1.5           !   [nondim] default = -1.33
-                                ! Exponent for Langmuir enhancementt of mstar
+                                ! Exponent for Langmuir enhancement of mstar
 LT_MOD_LAC1 = 0.0               !   [nondim] default = -0.87
                                 ! Coefficient for modification of Langmuir number due to MLD approaching Ekman
                                 ! depth.
@@ -2243,6 +2254,10 @@ LT_MOD_LAC4 = 0.0               !   [nondim] default = 0.95
 LT_MOD_LAC5 = 0.22              !   [nondim] default = 0.95
                                 ! Coefficient for modification of Langmuir number due to ratio of Ekman to
                                 ! unstable Obukhov depth.
+EPBL_OPTIONS_DIFF = 0           ! default = 0
+                                ! If positive, this is a coded integer indicating a pair of settings whose
+                                ! differences are diagnosed in a passive diagnostic mode  via extra calls to
+                                ! ePBL_column.  If this is 0 or negative no extra calls occur.
 !EPBL_USTAR_MIN = 1.45842E-18   !   [m s-1]
                                 ! The (tiny) minimum friction velocity used within the ePBL code, derived from
                                 ! OMEGA and ANGSTROM.
diff --git a/exps/OM4p25.COBALT/MOM_parameter_doc.layout b/exps/OM4p25.COBALT/MOM_parameter_doc.layout
index 3dc668d9b..227734a5b 100644
--- a/exps/OM4p25.COBALT/MOM_parameter_doc.layout
+++ b/exps/OM4p25.COBALT/MOM_parameter_doc.layout
@@ -46,7 +46,7 @@ NJPROC = 30                     !
                                 ! in MOM_memory.h at compile time.
 LAYOUT = 30, 30                 !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 1, 1                ! default = 1
+IO_LAYOUT = 1, 1                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 
diff --git a/exps/OM4p25.COBALT/MOM_parameter_doc.short b/exps/OM4p25.COBALT/MOM_parameter_doc.short
index faab86c2e..17b1df6b1 100644
--- a/exps/OM4p25.COBALT/MOM_parameter_doc.short
+++ b/exps/OM4p25.COBALT/MOM_parameter_doc.short
@@ -662,12 +662,12 @@ EPBL_LANGMUIR_SCHEME = "ADDITIVE" ! default = "NONE"
                                 !    NONE     - Do not do any extra mixing due to Langmuir turbulence
                                 !    RESCALE  - Use a multiplicative rescaling of mstar to account for Langmuir
                                 !      turbulence
-                                !    ADDITIVE - Add a Langmuir turblence contribution to mstar to other
+                                !    ADDITIVE - Add a Langmuir turbulence contribution to mstar to other
                                 !      contributions
 LT_ENHANCE_COEF = 0.044         !   [nondim] default = 0.447
                                 ! Coefficient for Langmuir enhancement of mstar
 LT_ENHANCE_EXP = -1.5           !   [nondim] default = -1.33
-                                ! Exponent for Langmuir enhancementt of mstar
+                                ! Exponent for Langmuir enhancement of mstar
 LT_MOD_LAC1 = 0.0               !   [nondim] default = -0.87
                                 ! Coefficient for modification of Langmuir number due to MLD approaching Ekman
                                 ! depth.
diff --git a/exps/OM4p25.COBALT/SIS_parameter_doc.layout b/exps/OM4p25.COBALT/SIS_parameter_doc.layout
index bbd24aeb3..79266c048 100644
--- a/exps/OM4p25.COBALT/SIS_parameter_doc.layout
+++ b/exps/OM4p25.COBALT/SIS_parameter_doc.layout
@@ -44,7 +44,7 @@ NJPROC = 30                     !
                                 ! in SIS2_memory.h at compile time.
 LAYOUT = 30, 30                 !
                                 ! The processor layout that was actually used.
-IO_LAYOUT = 1, 1                ! default = 1
+IO_LAYOUT = 1, 1                ! default = 1, 1
                                 ! The processor layout to be used, or 0,0 to automatically set the io_layout to
                                 ! be the same as the layout.
 NIBLOCK = 1                     ! default = 1
diff --git a/src/MOM6 b/src/MOM6
index d7a3929fe..7f63b80bd 160000
--- a/src/MOM6
+++ b/src/MOM6
@@ -1 +1 @@
-Subproject commit d7a3929fe7e38f6c599dfe6c02e7e6c47eac4e88
+Subproject commit 7f63b80bd0d480dc8c6ebc71f4cb9a2653a65a78
diff --git a/src/ocean_BGC b/src/ocean_BGC
index 135ecc138..154cd5995 160000
--- a/src/ocean_BGC
+++ b/src/ocean_BGC
@@ -1 +1 @@
-Subproject commit 135ecc138af87d79a5a748eb02f603ed29edfbb2
+Subproject commit 154cd599514fda31837d08143fb0680dbbcee18b
diff --git a/tools/boundary/glorys_obc_workflow/README.md b/tools/boundary/glorys_obc_workflow/README.md
new file mode 100644
index 000000000..9b7e73fe7
--- /dev/null
+++ b/tools/boundary/glorys_obc_workflow/README.md
@@ -0,0 +1,145 @@
+# MOM6 Open Boundary Conditions (OBC) Generation Workflow
+
+This repository provides an example workflow for generating Open Boundary Conditions (OBC) for MOM6 using daily GLORYS data on PPAN.
+
+## Overview
+
+The main script, `mom6_obc_workflow.sh`, orchestrates the entire OBC generation process. The workflow includes the following steps:
+
+1. **Spatial Subsetting of GLORYS Data**  
+   - Iterates through each day within a specified date range to spatially subset the original GLORYS dataset on UDA.  
+   - Reduces computational cost by limiting input data to the regional domain of interest instead of the entire global GLORYS domain.
+
+2. **Filling Missing Values in Subset GLORYS Files**  
+   - Processes each daily subset file using CDO to fill missing values.  
+   - Compresses processed files with `ncks -4 -L 5`.  
+   - Combines all variables (e.g., `thetao`, `so`, `zos`, `uo`, `vo`) into a single NetCDF file.
+
+3. **Daily Boundary Condition Generation**  
+   - Submits jobs to execute the `write_glorys_boundary_daily.py` script for each day.  
+   - Regrids GLORYS data and generates daily OBC files.
+
+Template scripts for these steps are provided in the `template` directory. User-specific parameters are configured using [uwtools](https://github.com/ufs-community/uwtools) to render templates, creating a `config.yaml` file based on user input.
+
+---
+
+## Configuration Example
+
+Below is an example `config.yaml` file to set up parameters for the workflow:
+
+```yaml
+# General parameters for template scripts
+_WALLTIME: "1440" # Wall time (in minutes) for SLURM jobs
+_NPROC: "1" # Number of processes for each job
+_EMAIL_NOTIFICATION: "fail" # SLURM email notification option
+_USER_EMAIL: "your.email@example.com" # Email address for error notifications
+_LOG_PATH: "./log/$CURRENT_DATE/%x.o%j" # Path for job logs
+_UDA_GLORYS_DIR: "/uda/Global_Ocean_Physics_Reanalysis/global/daily" # Path to original GLORYS data
+_UDA_GLORYS_FILENAME: "mercatorglorys12v1_gl12_mean" # File name prefix for GLORYS data
+_REGIONAL_GLORYS_ARCHIVE: "/archive/user/datasets/glorys" # Archive path for processed daily files
+_BASIN_NAME: "NWA12" # Regional domain name
+_OUTPUT_PREFIX: "GLORYS" # Prefix for output files
+_VARS: "thetao so uo vo zos" # Variables to process
+_LON_MIN: "-100.0" # Minimum longitude for subsetting
+_LON_MAX: "-30.0" # Maximum longitude for subsetting
+_LAT_MIN: "5.0" # Minimum latitude for subsetting
+_LAT_MAX: "60.0" # Maximum latitude for subsetting
+_PYTHON_SCRIPT: "$PYTHON_SCRIPT" # Path to the Python script for daily OBC generation
+
+# Date range for processing
+first_date: "$START_DATE"
+last_date: "$END_DATE"
+
+# Python script parameters
+glorys_dir: "/archive/user/datasets/glorys/NWA12/filled" # daily subset of GLORYS DATA after filling NaN
+output_dir: "./outputs" # output path for the obc files
+hgrid: "./ocean_hgrid.nc" # grid file
+ncrcat_names:
+  - "thetao"
+  - "so"
+  - "zos"
+  - "uv"
+segments:
+  - id: 1
+    border: "south"
+  - id: 2
+    border: "north"
+  - id: 3
+    border: "east"
+variables:
+  - "thetao"
+  - "so"
+  - "zos"
+  - "uv"
+```
+
+# Workflow Usage
+
+## Step 1: Modify Configuration
+Update the `cat <<EOF > config.yaml` part in `mom6_obc_workflow.sh` with parameters specific to your domain and workflow requirements.
+
+```
+cat <<EOF > config.yaml
+_WALLTIME: "1440"
+_NPROC: "1"
+_EMAIL_NOTIFACTION: "fail"
+_USER_EMAIL: "yi-cheng.teng@noaa.gov"
+_LOG_PATH: "./log/$CURRENT_DATE/%x.o%j"
+_UDA_GLORYS_DIR: "/uda/Global_Ocean_Physics_Reanalysis/global/daily"
+_UDA_GLORYS_FILENAME: "mercatorglorys12v1_gl12_mean"
+_REGIONAL_GLORYS_ARCHIVE: "/archive/ynt/datasets/glorys"
+_BASIN_NAME: "NWA12"
+_OUTPUT_PREFIX: "GLORYS"
+_VARS: "thetao so uo vo zos"
+_LON_MIN: "-100.0"
+_LON_MAX: "-30.0"
+_LAT_MIN: "5.0"
+_LAT_MAX: "60.0"
+_PYTHON_SCRIPT: "$PYTHON_SCRIPT"
+first_date: "$START_DATE"
+last_date: "$END_DATE"
+glorys_dir: "/archive/ynt/datasets/glorys/NWA12/filled"
+output_dir: "./outputs"
+hgrid: './ocean_hgrid.nc'
+ncrcat_names:
+  - 'thetao'
+  - 'so'
+  - 'zos'
+  - 'uv'
+segments:
+  - id: 1
+    border: 'south'
+  - id: 2
+    border: 'north'
+  - id: 3
+    border: 'east'
+variables:
+  - 'thetao'
+  - 'so'
+  - 'zos'
+  - 'uv'
+EOF
+```
+
+## Step 2: Generate OBC Files
+Run the workflow for a specific year or date range:
+
+```bash
+./mom6_obc_workflow.sh 2022-01-01 2022-12-31
+./mom6_obc_workflow.sh 2023-01-01 2023-12-31
+```
+
+## Step 3: Concatenate Multiple Years of OBC Files
+To merge OBC files from multiple years into a single file, use the `--ncrcat` option. Ensure the dates in the command match the range for which you generated OBC files:
+
+```bash
+./mom6_obc_workflow.sh 2022-01-01 2023-12-31 --ncrcat
+```
+### Adjust Timestamps (Optional Substep)
+If you need to adjust the timestamps of the first and last records for compatibility with `MOM6` yearly simulation, use the `--adjust-timestamps` option in combination with `--ncrcat`. Note that this is an alternative to the command above and should not be run afterward: 
+```
+./mom6_obc_workflow.sh 2022-01-01 2023-12-31 --ncrcat --adjust-timestamps
+``` 
+**Note**: Ensure the date range specified in your command corresponds to the dates for which you generated OBC files. Running this step with a mismatched date range will cause it to fail if files for the specified dates are missing.
+
+TODO:
diff --git a/tools/boundary/glorys_obc_workflow/mom6_obc_workflow.sh b/tools/boundary/glorys_obc_workflow/mom6_obc_workflow.sh
new file mode 100755
index 000000000..3fc805d3f
--- /dev/null
+++ b/tools/boundary/glorys_obc_workflow/mom6_obc_workflow.sh
@@ -0,0 +1,179 @@
+#!/bin/bash
+
+# Load required modules and environments
+source $MODULESHOME/init/sh
+module load miniforge
+conda activate /nbhome/role.medgrp/.conda/envs/uwtools || { echo "Error activating conda environment. Exiting."; exit 1; }
+
+set -eu
+
+# Helper functions
+print_usage() {
+    echo "Usage: $0 START_DATE END_DATE [--ncrcat] [--adjust-timestamps]"
+    echo "  START_DATE and END_DATE must be in YYYY-MM-DD format."
+    echo "  --ncrcat: Enable ncrcat step (skips subset, fill, and submit_python steps)."
+    echo "  --adjust-timestamps: Adjust timestamps during ncrcat step."
+}
+
+validate_date_format() {
+    if [[ ! "$1" =~ ^[0-9]{4}-[0-9]{2}-[0-9]{2}$ ]]; then
+        echo "Error: Date $1 must be in YYYY-MM-DD format."
+        exit 1
+    fi
+}
+
+log_message() {
+    echo "[$(date '+%Y-%m-%d %H:%M:%S')] $1"
+}
+
+# Default options
+DO_NCRCAT=false
+ADJUST_TIMESTAMPS=false
+PYTHON_SCRIPT="../write_glorys_boundary_daily.py"
+
+# Parse arguments
+START_DATE="$1"
+END_DATE="$2"
+shift 2
+
+while [[ $# -gt 0 ]]; do
+    case "$1" in
+        --ncrcat)
+            DO_NCRCAT=true
+            ;;
+        --adjust-timestamps)
+            ADJUST_TIMESTAMPS=true
+            ;;
+        *)
+            echo "Unknown argument: $1"
+            print_usage
+            exit 1
+            ;;
+    esac
+    shift
+done
+
+validate_date_format "$START_DATE"
+validate_date_format "$END_DATE"
+
+
+start_date_epoch=$(date -d "$START_DATE" +%s)
+end_date_epoch=$(date -d "$END_DATE" +%s)
+if [[ $start_date_epoch -gt $end_date_epoch ]]; then
+    log_message "Error: START_DATE ($START_DATE) must not be after END_DATE ($END_DATE). Exiting."
+    exit 1
+fi
+
+# Ensure --adjust-timestamps is only used with --ncrcat
+if $ADJUST_TIMESTAMPS && ! $DO_NCRCAT; then
+    echo "Error: --adjust-timestamps can only be used with --ncrcat."
+    exit 1
+fi
+
+# Warn user when --ncrcat is enabled
+if $DO_NCRCAT; then
+    log_message "WARNING: --ncrcat is enabled. The script will SKIP subset, fill, and submit_python steps."
+    log_message "Ensure that all daily outputs already exist for the specified date range."
+fi
+
+# Prepare directories
+CURRENT_DATE=$(date +%Y-%m-%d-%H-%M)
+mkdir -p ./log/$CURRENT_DATE ./outputs scripts
+
+# Define user configurations
+log_message "Generating config.yaml..."
+cat <<EOF > config.yaml
+_WALLTIME: "1440"
+_NPROC: "1"
+_EMAIL_NOTIFACTION: "fail"
+_USER_EMAIL: "$USER@noaa.gov"
+_LOG_PATH: "./log/$CURRENT_DATE/%x.o%j"
+_UDA_GLORYS_DIR: "/uda/Global_Ocean_Physics_Reanalysis/global/daily"
+_UDA_GLORYS_FILENAME: "mercatorglorys12v1_gl12_mean"
+_REGIONAL_GLORYS_ARCHIVE: "/archive/$USER/datasets/glorys"
+_BASIN_NAME: "NWA12"
+_OUTPUT_PREFIX: "GLORYS"
+_VARS: "thetao so uo vo zos"
+_LON_MIN: "-100.0"
+_LON_MAX: "-30.0"
+_LAT_MIN: "5.0"
+_LAT_MAX: "60.0"
+_PYTHON_SCRIPT: "$PYTHON_SCRIPT"
+first_date: "$START_DATE"
+last_date: "$END_DATE"
+glorys_dir: "/archive/$USER/datasets/glorys/NWA12/filled"
+output_dir: "./outputs"
+hgrid: './ocean_hgrid.nc'
+ncrcat_names:
+  - 'thetao'
+  - 'so'
+  - 'zos'
+  - 'uv'
+segments:
+  - id: 1
+    border: 'south'
+  - id: 2
+    border: 'north'
+  - id: 3
+    border: 'east'
+variables:
+  - 'thetao'
+  - 'so'
+  - 'zos'
+  - 'uv'
+EOF
+
+
+log_message "Preparing scripts directory..."
+[[ -d scripts ]] || mkdir scripts
+for template in subset_glorys fill_glorys submit_python_make_obc_day ncrcat_obc; do
+    rm -f scripts/${template}.sh
+    uw template render --input-file template/${template}_template.sh \
+                      --values-file config.yaml \
+                      --output-file scripts/${template}.sh || { log_message "Error rendering ${template}. Exiting."; exit 1; }
+done
+
+# Skip main steps if --ncrcat is enabled
+if ! $DO_NCRCAT; then
+
+    # Submit jobs
+    current_date_epoch=$start_date_epoch
+    job_ids=()
+
+    while [[ $current_date_epoch -le $end_date_epoch ]]; do
+        current_date=$(date -d "@$current_date_epoch" +%Y-%m-%d)
+        year=$(date -d "$current_date" +%Y)
+        month=$(date -d "$current_date" +%m)
+        day=$(date -d "$current_date" +%d)
+
+        log_message "Submitting subset job for $current_date..."
+        subset_job_id=$(sbatch --job-name="glorys_subset_${year}_${month}_${day}" \
+                              scripts/subset_glorys.sh $year $month $day | awk '{print $4}')
+
+        log_message "Submitting fill_nan job for $current_date..."
+        fill_job_id=$(sbatch --dependency=afterok:$subset_job_id \
+                              --job-name="glorys_fill_${year}_${month}_${day}" \
+                              scripts/fill_glorys.sh $year $month $day | awk '{print $4}')
+
+        log_message "Submitting Python job for $current_date..."
+        python_job_id=$(sbatch --dependency=afterok:$fill_job_id \
+                                --job-name="python_make_obc_day_${year}_${month}_${day}" \
+                                scripts/submit_python_make_obc_day.sh $year $month $day | awk '{print $4}')
+
+        job_ids+=($python_job_id)
+        current_date_epoch=$((current_date_epoch + 86400))
+    done
+fi
+
+# Optional ncrcat step
+if $DO_NCRCAT; then
+    log_message "Submitting ncrcat job..."
+    dependency_str=$(IFS=,; echo "${job_ids[*]:-}")
+    if $ADJUST_TIMESTAMPS; then
+        sbatch --job-name="obc_ncrcat" scripts/ncrcat_obc.sh --config config.yaml --ncrcat_years --adjust_timestamps
+    else
+        sbatch --job-name="obc_ncrcat" scripts/ncrcat_obc.sh --config config.yaml --ncrcat_years
+    fi
+fi
+
+log_message "Workflow completed."
diff --git a/tools/boundary/glorys_obc_workflow/template/fill_glorys_template.sh b/tools/boundary/glorys_obc_workflow/template/fill_glorys_template.sh
new file mode 100755
index 000000000..49f2ef0c6
--- /dev/null
+++ b/tools/boundary/glorys_obc_workflow/template/fill_glorys_template.sh
@@ -0,0 +1,150 @@
+#!/bin/bash
+#SBATCH --partition=batch
+#SBATCH --time={{ _WALLTIME }}
+#SBATCH --ntasks={{ _NPROC }}
+#SBATCH --mail-type={{ _EMAIL_NOTIFACTION }}
+#SBATCH --mail-user={{ _USER_EMAIL }}
+#SBATCH --output={{ _LOG_PATH }}
+
+# Usage: sbatch fill_glorys.sh <YEAR> <MONTH> <DAY>
+# This script is used to fill the GLORYS reanalysis data on GFDL PPAN
+
+module load cdo
+module load nco
+module load gcp
+
+# Input arguments
+year=$1
+month=$2
+day=$3
+
+# Configuration variables
+REGIONAL_GLORYS_ARCHIVE={{ _REGIONAL_GLORYS_ARCHIVE }}
+BASIN_NAME={{ _BASIN_NAME }}
+OUTPUT_PREFIX={{ _OUTPUT_PREFIX }}
+
+# Derived paths
+glorys_arch_dir="${REGIONAL_GLORYS_ARCHIVE}/${BASIN_NAME}"
+filled_dir="$glorys_arch_dir/filled"
+tmpdir="$TMPDIR/fill_glorys_${year}_${month}_${day}"
+combined_file="$filled_dir/${OUTPUT_PREFIX}_${year}-${month}-${day}.nc"
+
+# Function to set up required directories
+setup_directories() {
+    mkdir -p "$filled_dir"
+    mkdir -p "$tmpdir/filled"
+}
+
+# Function to clean up temporary files
+cleanup() {
+    echo "Cleaning up temporary files..."
+    rm -rf "$tmpdir"
+}
+trap cleanup EXIT
+
+# Function to check if file exists and handle errors
+check_file_exists() {
+    local file=$1
+    if [[ ! -f "$file" ]]; then
+        echo "Error: File does not exist: $file. Skipping."
+        return 1
+    fi
+    return 0
+}
+
+# Function to dmget and check if files are available
+dmget_files() {
+    dmget $glorys_arch_dir/${OUTPUT_PREFIX}_*_${year}-${month}-${day}.nc
+    for file in $glorys_arch_dir/${OUTPUT_PREFIX}_*_${year}-${month}-${day}.nc; do
+        check_file_exists "$file" || return 1
+    done
+}
+
+# Function to process and fill a single file
+process_file() {
+    local file=$1
+    local filename="${file##*/}"
+    local output_file="$filled_dir/$filename"
+
+    echo "Processing file: $file"
+
+    # Skip `cdo` part if filled file already exists
+    if [[ -f "$output_file" ]]; then
+        echo "Filled file already exists: $output_file. Skipping filling process."
+        return 0
+    fi
+
+    # Copy file to TMPDIR for processing
+    gcp "$file" "$tmpdir/"
+    local input_file="$tmpdir/$filename"
+
+    # Process: Fill missing data and compress
+    local filled_tmpfile="$tmpdir/filled/$filename"
+    cdo setmisstonn "$input_file" "$filled_tmpfile" || { echo "Error: Failed to fill missing data for $filename."; return 1; }
+    
+    ncks -4 -L 5 "$filled_tmpfile" -O "$filled_tmpfile" || { echo "Error: Failed to compress filled file for $filename."; return 1; }
+
+    # Copy filled file to the final directory
+    gcp "$filled_tmpfile" "$output_file" || { echo "Error: Failed to move filled file to $filled_dir."; return 1; }
+
+    return 0
+}
+
+# Function to initialize the combined file
+initialize_combined_file() {
+    local output_file=$1
+    echo "Initializing combined file with: $output_file"
+    cp "$output_file" "$combined_file" || { echo "Error: Failed to initialize combined file."; return 1; }
+    return 0
+}
+
+# Function to append variables to the combined file
+append_to_combined_file() {
+    local output_file=$1
+    echo "Appending variables from $output_file to $combined_file"
+    ncks -A "$output_file" "$combined_file" || { echo "Error: Failed to append variables from $output_file to $combined_file."; return 1; }
+    return 0
+}
+
+# Main processing loop
+main() {
+    setup_directories
+
+    # dmget files from archive
+    dmget_files || exit 1
+
+    # Remove the combined file if it exists
+    if [[ -f "$combined_file" ]]; then
+        echo "Removing existing combined file: $combined_file"
+        rm -f "$combined_file"
+    fi
+
+    first_file=true
+
+    for file in $glorys_arch_dir/${OUTPUT_PREFIX}_*_${year}-${month}-${day}.nc; do
+        check_file_exists "$file" || continue
+
+        if process_file "$file"; then
+            if $first_file; then
+                initialize_combined_file "$filled_dir/$(basename $file)" || exit 1
+                first_file=false
+            else
+                append_to_combined_file "$filled_dir/$(basename $file)" || exit 1
+            fi
+            echo "Successfully processed and appended: $file"
+        else
+            echo "Skipping file: $file due to errors."
+        fi
+    done
+
+    # Final verification
+    if [[ -f "$combined_file" ]]; then
+        echo "Successfully created or updated combined file: $combined_file"
+    else
+        echo "Error: Combined file was not created."
+        exit 1
+    fi
+}
+
+# Execute the main function
+main
diff --git a/tools/boundary/glorys_obc_workflow/template/ncrcat_obc_template.sh b/tools/boundary/glorys_obc_workflow/template/ncrcat_obc_template.sh
new file mode 100755
index 000000000..370b8c3f9
--- /dev/null
+++ b/tools/boundary/glorys_obc_workflow/template/ncrcat_obc_template.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+#SBATCH --partition=batch
+#SBATCH --time={{ _WALLTIME  }}
+#SBATCH --ntasks={{ _NPROC }}
+#SBATCH --mail-type={{ _EMAIL_NOTIFACTION }}
+#SBATCH --mail-user={{ _USER_EMAIL }}
+#SBATCH --output=./log/%x.o%j
+
+# Usage: sbatch ncrcat_obc.sh <ncrcat_arg>
+# This script is used to submit sbatch jobs for concatenating daily obc file on GFDL PPAN
+
+echo "Running ncrcat_obc.sh with arguments: $@"
+
+source $MODULESHOME/init/sh
+module load miniforge
+conda activate /nbhome/role.medgrp/.conda/envs/medpy311
+module load cdo nco gcp
+
+# Configuration variables
+REGIONAL_GLORYS_ARCHIVE={{ _REGIONAL_GLORYS_ARCHIVE }}
+BASIN_NAME={{ _BASIN_NAME }}
+OUTPUT_PREFIX={{ _OUTPUT_PREFIX }}
+PYTHON_SCRIPT={{ _PYTHON_SCRIPT }}
+
+
+#
+python "$PYTHON_SCRIPT" "$@" 
diff --git a/tools/boundary/glorys_obc_workflow/template/submit_python_make_obc_day_template.sh b/tools/boundary/glorys_obc_workflow/template/submit_python_make_obc_day_template.sh
new file mode 100755
index 000000000..35900d256
--- /dev/null
+++ b/tools/boundary/glorys_obc_workflow/template/submit_python_make_obc_day_template.sh
@@ -0,0 +1,84 @@
+#!/bin/bash
+#SBATCH --partition=batch
+#SBATCH --time={{ _WALLTIME }}
+#SBATCH --ntasks={{ _NPROC }}
+#SBATCH --mail-type={{ _EMAIL_NOTIFACTION }}
+#SBATCH --mail-user={{ _USER_EMAIL }}
+#SBATCH --output={{ _LOG_PATH }}
+
+# Usage: sbatch submit_python_make_obc_day.sh <YEAR> <MONTH> <DAY>
+# This script is used to submit sbatch jobs for daily GLORYS OBC file generation on GFDL PPAN
+
+# Setup environment
+setup_environment() {
+    source $MODULESHOME/init/sh
+    module load miniforge
+    conda activate /nbhome/role.medgrp/.conda/envs/medpy311
+    module load cdo nco gcp
+}
+
+# Check if the combined file exists
+check_combined_file() {
+    local combined_file=$1
+    if [[ ! -f "$combined_file" ]]; then
+        echo "Error: Combined file not found: $combined_file"
+        return 1
+    fi
+    return 0
+}
+
+# dmget the combined file from archive
+dmget_combined_file() {
+    local combined_file=$1
+    echo "dmget combined file: $combined_file"
+    dmget "$combined_file" || { echo "Error: Failed to dmget combined file: $combined_file"; return 1; }
+    return 0
+}
+
+# Run the Python script for daily OBC file generation
+run_python_script() {
+    local year=$1
+    local month=$2
+    local day=$3
+    local python_script=$4
+
+    echo "Creating daily OBC file: ${year}-${month}-${day}"
+    python "$python_script" --config config.yaml --year "$year" --month "$month" --day "$day" || { echo "Error: Python script failed for ${year}-${month}-${day}"; return 1; }
+    return 0
+}
+
+# Main function
+main() {
+    # Input arguments
+    year=$1
+    month=$2
+    day=$3
+
+    # Configuration variables
+    REGIONAL_GLORYS_ARCHIVE={{ _REGIONAL_GLORYS_ARCHIVE }}
+    BASIN_NAME={{ _BASIN_NAME }}
+    OUTPUT_PREFIX={{ _OUTPUT_PREFIX }}
+    PYTHON_SCRIPT={{ _PYTHON_SCRIPT }}
+
+    # Derived paths
+    glorys_arch_dir="${REGIONAL_GLORYS_ARCHIVE}/${BASIN_NAME}"
+    filled_dir="$glorys_arch_dir/filled"
+    combined_file="$filled_dir/${OUTPUT_PREFIX}_${year}-${month}-${day}.nc"
+
+    # Setup environment
+    setup_environment
+
+    # Check if the combined file exists
+    check_combined_file "$combined_file" || exit 1
+
+    # dmget the combined file from the archive
+    dmget_combined_file "$combined_file" || exit 1
+
+    # Run the Python script to create the daily OBC file
+    run_python_script "$year" "$month" "$day" "$PYTHON_SCRIPT" || exit 1
+
+    echo "Daily OBC file generation completed for ${year}-${month}-${day}"
+}
+
+# Execute the main function
+main "$@"
diff --git a/tools/boundary/glorys_obc_workflow/template/subset_glorys_template.sh b/tools/boundary/glorys_obc_workflow/template/subset_glorys_template.sh
new file mode 100755
index 000000000..a7a28f821
--- /dev/null
+++ b/tools/boundary/glorys_obc_workflow/template/subset_glorys_template.sh
@@ -0,0 +1,121 @@
+#!/bin/bash
+#SBATCH --partition=batch
+#SBATCH --time={{ _WALLTIME }}
+#SBATCH --ntasks={{ _NPROC }}
+#SBATCH --mail-type={{ _EMAIL_NOTIFACTION }}
+#SBATCH --mail-user={{ _USER_EMAIL }}
+#SBATCH --output={{ _LOG_PATH }}
+
+# Usage: sbatch subset_glorys.sh <YEAR> <MONTH> <DAY>
+
+# Load required modules
+module load cdo
+module load nco/5.0.1
+module load gcp
+
+# Input arguments
+year=$1
+month=$2
+day=$3
+
+# Configuration variables
+UDA_GLORYS_DIR={{ _UDA_GLORYS_DIR }}
+UDA_GLORYS_FILENAME={{ _UDA_GLORYS_FILENAME }}
+REGIONAL_GLORYS_ARCHIVE={{ _REGIONAL_GLORYS_ARCHIVE }}
+BASIN_NAME={{ _BASIN_NAME }}
+OUTPUT_PREFIX={{ _OUTPUT_PREFIX }}
+LON_MIN={{ _LON_MIN }}
+LON_MAX={{ _LON_MAX }}
+LAT_MIN={{ _LAT_MIN }}
+LAT_MAX={{ _LAT_MAX }}
+VARS=({{ _VARS }})
+
+# Derived variables
+apath=${REGIONAL_GLORYS_ARCHIVE}/${BASIN_NAME}
+tmpdir="$TMPDIR/subset_glorys_${year}_${month}_${day}"
+
+# Create required directories
+setup_directories() {
+    mkdir -p "$apath"
+    mkdir -p "$tmpdir"
+}
+
+# Cleanup temporary files
+cleanup() {
+    echo "Cleaning up temporary files..."
+    rm -rf "$tmpdir"
+}
+trap cleanup EXIT
+
+# Check if output file exists
+check_output_exists() {
+    local output_file=$1
+    if [[ -f $output_file ]]; then
+        echo "Output file already exists: $output_file. Skipping processing."
+        return 0
+    fi
+    return 1
+}
+
+# Process a single variable
+process_variable() {
+    local var=$1
+    local output_file="$apath/${OUTPUT_PREFIX}_${var}_${year}-${month}-${day}.nc"
+
+    # Skip processing if output file exists
+    check_output_exists "$output_file" && return
+
+    local var_dir="$tmpdir/${var}_${year}_${month}_${day}"
+    mkdir -p "$var_dir"
+
+    local filename_pattern="${UDA_GLORYS_DIR}/${var}/${year}/${var}_${UDA_GLORYS_FILENAME}_${year}${month}${day}*.nc"
+    local processed=false
+
+    for filename in $filename_pattern; do
+        if [[ -f $filename ]]; then
+            echo "Processing file: $filename"
+            local subset_file="$var_dir/${var}_subset.nc"
+
+            # Subset and adjust the file
+            if ! ncks -d longitude,${LON_MIN},${LON_MAX} -d latitude,${LAT_MIN},${LAT_MAX} --mk_rec_dmn time "$filename" "$subset_file"; then
+                echo "Error: ncks failed for $filename"
+                exit 1
+            fi
+
+            if ! cdo -setreftime,1993-01-01,00:00:00,1day "$subset_file" "$output_file"; then
+                echo "Error: cdo failed for $subset_file"
+                exit 1
+            fi
+
+            echo "Finished processing variable: $var for ${year}-${month}-${day}"
+            processed=true
+            break
+        fi
+    done
+
+    if [[ $processed == false ]]; then
+        echo "Error: No files matched the pattern: $filename_pattern"
+        exit 1
+    fi
+
+    # Verify the output file
+    if [[ -f $output_file ]]; then
+        echo "Output file successfully created: $output_file"
+    else
+        echo "Error: Failed to create output file: $output_file"
+        exit 1
+    fi
+}
+
+# Main processing loop
+main() {
+    setup_directories
+
+    for var in "${VARS[@]}"; do
+        process_variable "$var"
+    done
+
+    echo "All variables processed for ${year}-${month}-${day}"
+}
+
+main
diff --git a/tools/boundary/write_glorys_boundary_daily.py b/tools/boundary/write_glorys_boundary_daily.py
new file mode 100755
index 000000000..f252598c5
--- /dev/null
+++ b/tools/boundary/write_glorys_boundary_daily.py
@@ -0,0 +1,163 @@
+#!/usr/bin/env python3
+"""
+This script processes daily GLORYS data to generate ocean boundary conditions (OBC) for specified segments and variables. 
+It supports regridding of velocity (u, v) and tracer fields (e.g., temperature, salinity, sea surface height) for a single day or concatenating multiple days of outputs. 
+
+Key Features:
+1. Processes single-day outputs: Generates regridded NetCDF files for each specified segment and variable.
+2. Supports concatenation of results across multiple days or years using NCO tools (ncrcat).
+3. Adjusts timestamps in concatenated files (optional).
+
+Dependencies:
+- GLORYS data files (NetCDF format) with specific variable names.
+- NCO tools for file concatenation (optional, but required for multi-year concatenation).
+- Python libraries: xarray, yaml, argparse.
+
+Usage:
+1. Process single-day output:
+   ./write_glorys_boundary_day.py --config config.yaml --year <YEAR> --month <MONTH> --day <DAY>
+
+2. Concatenate multiple days of results with optional timestamp adjustment:
+   ./write_glorys_boundary_day.py --config config.yaml --ncrcat_years [--adjust_timestamps]
+
+Ensure that NaN values in GLORYS data are pre-filled with valid values.
+Run on a system with NCO tools installed, e.g., `module load nco/5.0.1` on HPC systems.
+"""
+
+import argparse
+import os
+from datetime import datetime, timedelta
+from subprocess import run
+from os import path
+
+import xarray
+import yaml
+from boundary import Segment
+
+# Suppress xarray warnings
+import warnings
+warnings.filterwarnings('ignore')
+
+def load_config(config_file):
+    """Load configuration from a YAML file."""
+    with open(config_file, 'r') as file:
+        return yaml.safe_load(file)
+
+def write_day(date, glorys_dir, segments, variables, output_prefix):
+    """Process and regrid data for a specific day."""
+    filename = f"{output_prefix}_{date.year}-{date.month:02d}-{date.day:02d}.nc"
+    file_path = path.join(glorys_dir, filename)
+
+    if not path.exists(file_path):
+        print(f"File does not exist: {file_path}. Skipping.")
+        return
+
+    glorys = (
+        xarray.open_dataset(file_path)
+        .rename({'latitude': 'lat', 'longitude': 'lon', 'depth': 'z'})
+    )
+
+    for segment in segments:
+        for variable in variables:
+            if variable == 'uv':
+                print(f"Processing {segment.border} {variable}")
+                segment.regrid_velocity(glorys['uo'], glorys['vo'], suffix=f"{date:%Y%m%d}", flood=False)
+            elif variable in ['thetao', 'so', 'zos']:
+                print(f"Processing {segment.border} {variable}")
+                segment.regrid_tracer(glorys[variable], suffix=f"{date:%Y%m%d}", flood=False)
+
+def concatenate_files(nsegments, output_dir, variables, ncrcat_names, first_date, last_date, adjust_timestamps=False):
+    """Concatenate annual files using ncrcat."""
+    if not ncrcat_names:
+        ncrcat_names = variables[:]
+
+    date_list = [(first_date + timedelta(days=i)).strftime("%Y%m%d")
+                 for i in range((last_date - first_date).days + 1)]
+
+    for variable, var_name in zip(variables, ncrcat_names):
+        for seg_id in range(1, nsegments + 1):
+            input_files = [
+                path.join(output_dir, f"{variable}_{seg_id:03d}_{date}.nc")
+                for date in date_list
+            ]
+            output_file = path.join(output_dir, f"{var_name}_{seg_id:03d}.nc")
+
+            if path.exists(output_file):
+                print(f"Removing existing file: {output_file}")
+                os.remove(output_file)
+
+            print(f"Concatenating files for {variable}, segment {seg_id} into {output_file}...")
+            run(["ncrcat", "-O", *input_files, output_file], check=True)
+
+            if adjust_timestamps:
+                adjust_file_timestamps(output_file)
+
+def adjust_file_timestamps(file_path):
+    """Adjust timestamps for the first and last records in a file."""
+    with xarray.open_dataset(file_path) as ds:
+        if 'time' in ds:
+            time = ds['time'].copy()
+            adjusted_time = time.astype('datetime64[ns]')
+            
+            adjusted_time[0] = adjusted_time[0].dt.floor('D')
+            adjusted_time[-1] = adjusted_time[-1].dt.ceil('D')
+
+            ds = ds.assign_coords(time=adjusted_time)
+            ds.to_netcdf(file_path)
+            print(f"Timestamps adjusted for {file_path}")
+
+def process_single_day(config, year, month, day):
+    """Process data for a single day."""
+    specific_date = datetime(year, month, day)
+    print(f"Processing data for {specific_date}...")
+
+    glorys_dir = config['glorys_dir']
+    output_prefix = config.get('_OUTPUT_PREFIX', 'GLORYS')
+    variables = config['variables']
+
+    hgrid = xarray.open_dataset(config['hgrid'])
+    segments = [
+        Segment(seg_config['id'], seg_config['border'], hgrid, output_dir=config['output_dir'])
+        for seg_config in config['segments']
+    ]
+
+    write_day(specific_date, glorys_dir, segments, variables, output_prefix)
+
+def concatenate_annual_files(config, adjust_timestamps):
+    """Concatenate files for the entire date range."""
+    first_date = datetime.strptime(config['first_date'], '%Y-%m-%d')
+    last_date = datetime.strptime(config['last_date'], '%Y-%m-%d')
+
+    print(f"Concatenating files from {first_date} to {last_date}...")
+    
+    concatenate_files(
+        len(config['segments']),
+        config['output_dir'],
+        config['variables'],
+        config.get('ncrcat_names', []),
+        first_date,
+        last_date,
+        adjust_timestamps
+    )
+
+def main():
+    parser = argparse.ArgumentParser(description="Generate OBC from GLORYS")
+    parser.add_argument('--config', type=str, default='glorys_obc.yaml', help="Path to the YAML configuration file")
+    parser.add_argument('--year', type=int, help="Year for single-day processing")
+    parser.add_argument('--month', type=int, help="Month for single-day processing")
+    parser.add_argument('--day', type=int, help="Day for single-day processing")
+    parser.add_argument('--ncrcat_years', action='store_true', help="Enable annual concatenation mode")
+    parser.add_argument('--adjust_timestamps', action='store_true', help="Adjust timestamps during concatenation")
+    args = parser.parse_args()
+
+    config = load_config(args.config)
+
+    if args.ncrcat_years:
+        concatenate_annual_files(config, args.adjust_timestamps)
+    elif args.year and args.month and args.day:
+        process_single_day(config, args.year, args.month, args.day)
+    else:
+        print("Error: Specify either --ncrcat_years or a specific date (--year, --month, --day).")
+
+if __name__ == '__main__':
+    main()