diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
index 9006d96b3..e3485ce6c 100644
--- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
@@ -73,17 +73,22 @@ subroutine GFS_phys_time_vary_init (
               zwtxy, xlaixy, xsaixy, lfmassxy, stmassxy, rtmassxy, woodxy, stblcpxy, fastcpxy,     &
               smcwtdxy, deeprechxy, rechxy, snowxy, snicexy, snliqxy, tsnoxy , smoiseq, zsnsoxy,   &
               slc, smc, stc, tsfcl, snowd, canopy, tg3, stype, con_t0c, lsm_cold_start, nthrds,    &
-              ozphys, errmsg, errflg)
+              lkm, use_lake_model, lakefrac, lakedepth, iopt_lake, iopt_lake_clm, iopt_lake_flake, &
+              lakefrac_threshold, lakedepth_threshold, ozphys, errmsg, errflg)
 
          implicit none
 
          ! Interface variables
          integer,              intent(in)    :: me, master, ntoz, iccn, iflip, im, nx, ny
          logical,              intent(in)    :: h2o_phys, iaerclm, lsm_cold_start
-         integer,              intent(in)    :: idate(:)
-         real(kind_phys),      intent(in)    :: fhour
+         integer,              intent(in)    :: idate(:), iopt_lake, iopt_lake_clm, iopt_lake_flake
+         real(kind_phys),      intent(in)    :: fhour, lakefrac_threshold, lakedepth_threshold
          real(kind_phys),      intent(in)    :: xlat_d(:), xlon_d(:)
 
+         integer,              intent(in) :: lkm
+         integer,              intent(inout)  :: use_lake_model(:)
+         real(kind=kind_phys), intent(in   )  :: lakefrac(:), lakedepth(:)
+
          integer,              intent(inout), optional :: jindx1_o3(:), jindx2_o3(:), jindx1_h(:), jindx2_h(:)
          real(kind_phys),      intent(inout), optional :: ddy_o3(:),  ddy_h(:)
          real(kind_phys),      intent(in)    :: h2opl(:,:,:)
@@ -600,6 +605,27 @@ subroutine GFS_phys_time_vary_init (
            endif noahmp_init
          endif lsm_init
 
+         ! Lake model
+         if(lkm>0 .and. iopt_lake>0) then
+            print*,"SWALES: In phys_time_vary",lakefrac(i),lakedepth(i)
+            ! A lake model is enabled.
+            do i = 1, im
+               !if (lakefrac(i) > 0.0 .and. lakedepth(i) > 1.0 ) then
+               ! The lake data must say there's a lake here (lakefrac) with a depth (lakedepth)
+               if (lakefrac(i) > lakefrac_threshold .and. lakedepth(i) > lakedepth_threshold ) then
+                  ! This is a lake point. Inform the other schemes to use a lake model, and possibly nsst (lkm)
+                  use_lake_model(i) = lkm
+                  cycle
+               else
+                  ! Not a valid lake point.
+                  use_lake_model(i) = 0
+               endif
+            enddo
+         else
+            ! Lake model is disabled or settings are invalid.
+            use_lake_model = 0
+         endif
+
          is_initialized = .true.
 
       contains
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta
index f4bd67527..0c2aaba70 100644
--- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta
@@ -947,6 +947,73 @@
   dimensions = ()
   type = integer
   intent = in
+[lkm]
+  standard_name = control_for_lake_model_execution_method
+  long_name = control for lake model execution: 0=no lake, 1=lake, 2=lake+nsst
+  units = flag
+  dimensions = ()
+  type = integer
+  intent = in
+[use_lake_model]
+  standard_name = flag_for_using_lake_model
+  long_name = flag indicating lake points using a lake model
+  units = flag
+  dimensions = (horizontal_dimension)
+  type = integer
+  intent = inout
+[lakefrac]
+  standard_name = lake_area_fraction
+  long_name = fraction of horizontal grid area occupied by lake
+  units = frac
+  dimensions = (horizontal_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[lakedepth]
+  standard_name = lake_depth
+  long_name = lake depth
+  units = m
+  dimensions = (horizontal_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[iopt_lake]
+  standard_name = control_for_lake_model_selection
+  long_name = control for lake model selection
+  units = 1
+  dimensions = ()
+  type = integer
+  intent = in
+[iopt_lake_clm]
+  standard_name = clm_lake_model_control_selection_value
+  long_name = value that indicates clm lake model in the control for lake model selection
+  units = 1
+  dimensions = ()
+  type = integer
+  intent = in
+[iopt_lake_flake]
+  standard_name = flake_model_control_selection_value
+  long_name = value that indicates flake model in the control for lake model selection
+  units = 1
+  dimensions = ()
+  type = integer
+  intent = in
+[lakefrac_threshold]
+  standard_name = lakefrac_threshold_for_enabling_lake_model
+  long_name = fraction of horizontal grid area occupied by lake must be greater than this value to enable a lake model
+  units = frac
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[lakedepth_threshold]
+  standard_name = lake_depth_threshold_for_enabling_lake_model
+  long_name = lake depth must be greater than this value to enable a lake model
+  units = m
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
 [ozphys]
   standard_name = dataset_for_ozone_physics
   long_name = dataset for NRL ozone physics