diff --git a/src/core_ocean/Registry.xml b/src/core_ocean/Registry.xml
index 62a64b7865..f4d3c7f9ca 100644
--- a/src/core_ocean/Registry.xml
+++ b/src/core_ocean/Registry.xml
@@ -724,6 +724,10 @@
 					description="The length scale of exponential decay of river runoff. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
 					possible_values="Any positive real number."
 		/>
+		<nml_option name="config_flux_attenuation_coefficient_prec_evap_ice" type="real" default_value="0.001" units="m"
+					description="The length scale of exponential decay of non-river runoff thickness fluxes. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
+					possible_values="Any positive real number."
+		/>
 	</nml_record>
 	<nml_record name="time_varying_forcing" mode="init;forward">
 		<nml_option name="config_use_time_varying_atmospheric_forcing" type="logical" default_value=".false." units="unitless"
@@ -2761,6 +2765,9 @@
 		<var name="surfaceFluxAttenuationCoefficientRunoff" type="real" dimensions="nCells Time" units="m"
 			description="The spatially-dependent length scale of exponential decay of river runoff. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
 		/>
+		<var name="surfaceFluxAttenuationCoefficientPrecEvapIce" type="real" dimensions="nCells Time" units="m"
+			description="The spatially-dependent length scale of exponential decay of non-river runoff thickness fluxes. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
+		/>
 		<!-- diagnostic fields for land-ice fluxes -->
 		<var name="landIceFrictionVelocity" type="real" dimensions="nCells Time" units="m s^{-1}"
 			description="The friction velocity $u_*$ under land ice"
@@ -2917,6 +2924,10 @@
 			 description="Divergence of transmission through interfaces of surface fluxes below the surface layer at cell centers. These are applied only to river runoff."
 			 packages="forwardMode;analysisMode"
 		/>
+		<var name="fractionAbsorbedPrecEvapIce" type="real" dimensions="nVertLevels nCells Time" units="fractional"
+			 description="Divergence of transmission through interfaces of surface fluxes below the surface layer at cell centers. These are applied only to non-river runoff thickness fluxes."
+			 packages="forwardMode;analysisMode"
+		/>
 
 		<!-- Input fields for coupling -->
 		<!-- Coupling fields associated with heat fluxes -->
diff --git a/src/core_ocean/shared/mpas_ocn_diagnostics.F b/src/core_ocean/shared/mpas_ocn_diagnostics.F
index c0516260f7..8b1d795994 100644
--- a/src/core_ocean/shared/mpas_ocn_diagnostics.F
+++ b/src/core_ocean/shared/mpas_ocn_diagnostics.F
@@ -153,8 +153,10 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
 
       real (kind=RKIND), pointer :: config_flux_attenuation_coefficient
       real (kind=RKIND), pointer :: config_flux_attenuation_coefficient_runoff
+      real (kind=RKIND), pointer :: config_flux_attenuation_coefficient_prec_evap_ice
       real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficient
       real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientRunoff
+      real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientPrecEvapIce
 
       real (kind=RKIND) :: areaTri1, layerThicknessVertexInv, tempRiVal, dcEdge_temp, dvEdge_temp, weightsOnEdge_temp
       real (kind=RKIND), dimension(:), allocatable:: layerThicknessVertexVec
@@ -176,6 +178,8 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
       call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient', config_flux_attenuation_coefficient)
       call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_runoff',  &
                                              config_flux_attenuation_coefficient_runoff)
+      call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_prec_evap_ice',  &
+                                             config_flux_attenuation_coefficient_prec_evap_ice)
       call mpas_pool_get_config(ocnConfigs, 'config_use_wetting_drying', config_use_wetting_drying)
       call mpas_pool_get_config(ocnConfigs, 'config_thickness_flux_type', config_thickness_flux_type)
 
@@ -261,6 +265,8 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
 
       call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficient', surfaceFluxAttenuationCoefficient)
       call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficientRunoff', surfaceFluxAttenuationCoefficientRunoff)
+      call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficientPrecEvapIce',  &
+           surfaceFluxAttenuationCoefficientPrecEvapIce)
       !
       ! Compute height on cell edges at velocity locations
       !   Namelist options control the order of accuracy of the reconstructed layerThicknessEdge value
@@ -852,6 +858,7 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnostic
       do iCell = 1, nCells
          surfaceFluxAttenuationCoefficient(iCell) = config_flux_attenuation_coefficient
          surfaceFluxAttenuationCoefficientRunoff(iCell) = config_flux_attenuation_coefficient_runoff
+         surfaceFluxAttenuationCoefficientPrecEvapIce(iCell) = config_flux_attenuation_coefficient_prec_evap_ice
       end do
       !$omp end do
 
@@ -1367,8 +1374,9 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
       real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaCell
       real (kind=RKIND), dimension(:), pointer :: penetrativeTemperatureFlux, surfaceThicknessFlux, &
            surfaceBuoyancyForcing, surfaceFrictionVelocity, penetrativeTemperatureFluxOBL, &
-           normalVelocitySurfaceLayer, surfaceThicknessFluxRunoff
-      real (kind=RKIND), pointer :: config_flux_attenuation_coefficient, config_flux_attenuation_coefficient_runoff
+           normalVelocitySurfaceLayer, surfaceThicknessFluxRunoff, rainFlux, evaporationFlux
+      real (kind=RKIND), pointer :: config_flux_attenuation_coefficient, config_flux_attenuation_coefficient_runoff, &
+           config_flux_attenuation_coefficient_prec_evap_ice
 
       real (kind=RKIND), dimension(:), pointer :: surfaceStress, surfaceStressMagnitude
       real (kind=RKIND), dimension(:,:), pointer ::  &
@@ -1388,7 +1396,8 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
       ! local
       integer :: iCell, iEdge, i, k, err, timeLevel
       integer, pointer :: indexTempFlux, indexSaltFlux
-      real (kind=RKIND) :: numerator, denominator, turbulentVelocitySquared, fracAbsorbed, fracAbsorbedRunoff
+      real (kind=RKIND) :: numerator, denominator, turbulentVelocitySquared, fracAbsorbedRunoff, &
+           fracAbsorbedPrecEvapIce
       real (kind=RKIND) :: factor, deltaVelocitySquared, delU2, invAreaCell
 
       type (field2DReal), pointer :: densitySurfaceDisplacedField, thermalExpansionCoeffField, salineContractionCoeffField
@@ -1406,6 +1415,8 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
       call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient', config_flux_attenuation_coefficient)
       call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_runoff',   &
               config_flux_attenuation_coefficient_runoff)
+      call mpas_pool_get_config(ocnConfigs, 'config_flux_attenuation_coefficient_prec_evap_ice',   &
+              config_flux_attenuation_coefficient_prec_evap_ice)
       ! set the parameter turbulentVelocitySquared
       turbulentVelocitySquared = 0.001_RKIND
 
@@ -1441,6 +1452,8 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
       call mpas_pool_get_array(forcingPool, 'penetrativeTemperatureFlux', penetrativeTemperatureFlux)
       call mpas_pool_get_array(forcingPool, 'surfaceStress', surfaceStress)
       call mpas_pool_get_array(forcingPool, 'surfaceStressMagnitude', surfaceStressMagnitude)
+      call mpas_pool_get_array(forcingPool, 'rainFlux', rainFlux)
+      call mpas_pool_get_array(forcingPool, 'evaporationFlux', evaporationFlux)
 
       call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, 1)
       call mpas_pool_get_array(tracersSurfaceFluxPool, 'activeTracersSurfaceFlux', activeTracersSurfaceFlux)
@@ -1489,19 +1502,24 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, meshPool, diagno
        !
 
        ! Compute fraction of thickness flux that is in the top model layer
-       fracAbsorbed = 1.0_RKIND - exp( max(-layerThickness(1, iCell) / config_flux_attenuation_coefficient, -100.0_RKIND) )
        fracAbsorbedRunoff = 1.0_RKIND - exp( max(-layerThickness(1, iCell) / config_flux_attenuation_coefficient_runoff, &
                             -100.0_RKIND) )
+       fracAbsorbedPrecEvapIce = 1.0_RKIND - exp( max(-layerThickness(1, iCell) / &
+          config_flux_attenuation_coefficient_prec_evap_ice, -100.0_RKIND) )
+
        ! Store the total tracer flux below in nonLocalSurfaceTemperatureFlux for use in the CVMix nonlocal
        ! transport code.  This includes tracer forcing due to thickness
 
        nonLocalSurfaceTracerFlux(indexTempFlux, iCell) = activeTracersSurfaceFlux(indexTempFlux,iCell) &
-               + penetrativeTemperatureFlux(iCell) - penetrativeTemperatureFluxOBL(iCell) - fracAbsorbed * &
-               surfaceThicknessFlux(iCell) * activeTracers(indexTempFlux,1,iCell) + &
-               activeTracersSurfaceFluxRunoff(indexTempFlux,iCell) * fracAbsorbedRunoff
+               + penetrativeTemperatureFlux(iCell) - penetrativeTemperatureFluxOBL(iCell)  &
+               - fracAbsorbedPrecEvapIce * (rainFlux(iCell) + evaporationFlux(iCell)) * &
+                   activeTracers(indexTempFlux,1,iCell)/rho_sw &
+               - fracAbsorbedRunoff * surfaceThicknessFluxRunoff(iCell) * (  &
+                 + min(activeTracers(indexTempFlux,1,iCell),0.0_RKIND)    )/rho_sw
 
        nonLocalSurfaceTracerFlux(indexSaltFlux,iCell) = activeTracersSurfaceFlux(indexSaltFlux,iCell) &
-               - fracAbsorbed * surfaceThicknessFlux(iCell) * activeTracers(indexSaltFlux,1,iCell)
+               - fracAbsorbedPrecEvapIce * surfaceThicknessFlux(iCell) * activeTracers(indexSaltFlux,1,iCell) &
+               - fracAbsorbedRunoff * surfaceThicknessFluxRunoff(iCell) * activeTracers(indexSaltFlux,1,iCell)
 
        surfaceBuoyancyForcing(iCell) =  thermalExpansionCoeff (1,iCell) &
               * nonLocalSurfaceTracerFlux(indexTempFlux,iCell) &
diff --git a/src/core_ocean/shared/mpas_ocn_forcing.F b/src/core_ocean/shared/mpas_ocn_forcing.F
index bbed313398..b39a04139a 100644
--- a/src/core_ocean/shared/mpas_ocn_forcing.F
+++ b/src/core_ocean/shared/mpas_ocn_forcing.F
@@ -113,7 +113,9 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, diagno
 
         real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficient
         real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientRunoff
-        real (kind=RKIND), dimension(:,:), pointer :: layerThickness, fractionAbsorbed, fractionAbsorbedRunoff
+        real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientPrecEvapIce
+        real (kind=RKIND), dimension(:,:), pointer :: layerThickness, fractionAbsorbed, fractionAbsorbedRunoff, &
+             fractionAbsorbedPrecEvapIce
 
         integer :: iCell, k, timeLevel, nCells
 
@@ -137,11 +139,17 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, diagno
         call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficientRunoff',  &
                                                    surfaceFluxAttenuationCoefficientRunoff)
 
+        call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficientPrecEvapIce',  &
+                                                   surfaceFluxAttenuationCoefficientPrecEvapIce)
+
         call mpas_pool_get_array(forcingPool, 'fractionAbsorbed', fractionAbsorbed)
         call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', fractionAbsorbedRunoff)
+        call mpas_pool_get_array(forcingPool, 'fractionAbsorbedPrecEvapIce', fractionAbsorbedPrecEvapIce)
 
         nCells = nCellsArray( 2 )
 
+!  first do surface fluxes that are not associated with thickness fluxes (radiative, sensible, etc)
+
         do iCell = 1, nCells
            zTop = 0.0_RKIND
            transmissionCoeffTop = ocn_forcing_transmission(zTop, surfaceFluxAttenuationCoefficient(iCell))
@@ -172,6 +180,22 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, diagno
            end do
         end do
 
+!  now do non-runoff thickness fluxes
+
+        do iCell = 1, nCells
+           zTop = 0.0_RKIND
+           transmissionCoeffTop = ocn_forcing_transmission(zTop, surfaceFluxAttenuationCoefficientPrecEvapIce(iCell))
+           do k = 1, maxLevelCell(iCell)
+              zBot = zTop - layerThickness(k,iCell)
+              transmissionCoeffBot = ocn_forcing_transmission(zBot, surfaceFluxAttenuationCoefficientPrecEvapIce(iCell))
+
+              fractionAbsorbedPrecEvapIce(k, iCell) = transmissionCoeffTop - transmissionCoeffBot
+
+              zTop = zBot
+              transmissionCoeffTop = transmissionCoeffBot
+           end do
+        end do
+
     end subroutine ocn_forcing_build_fraction_absorbed_array!}}}
 
 !***********************************************************************
diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F
index 4dacfb7452..2448d89e23 100644
--- a/src/core_ocean/shared/mpas_ocn_tendency.F
+++ b/src/core_ocean/shared/mpas_ocn_tendency.F
@@ -117,7 +117,8 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{
       real (kind=RKIND), dimension(:), pointer :: surfaceThicknessFlux
       real (kind=RKIND), dimension(:), pointer :: surfaceThicknessFluxRunoff
       real (kind=RKIND), dimension(:,:), pointer :: layerThickness, layerThicknessEdge, &
-         vertAleTransportTop, tend_layerThickness, normalTransportVelocity, fractionAbsorbed, fractionAbsorbedRunoff
+         vertAleTransportTop, tend_layerThickness, normalTransportVelocity, fractionAbsorbed, fractionAbsorbedRunoff, &
+         fractionAbsorbedPrecEvapIce
 
       integer, pointer :: nCells
       integer :: err, iCell
@@ -138,6 +139,7 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{
       call mpas_pool_get_array(forcingPool, 'surfaceThicknessFluxRunoff', surfaceThicknessFluxRunoff)
       call mpas_pool_get_array(forcingPool, 'fractionAbsorbed', fractionAbsorbed)
       call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', fractionAbsorbedRunoff)
+      call mpas_pool_get_array(forcingPool, 'fractionAbsorbedPrecEvapIce', fractionAbsorbedPrecEvapIce)
 
       !
       ! height tendency: start accumulating tendency terms
@@ -179,7 +181,7 @@ subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{
       ! surface flux tendency
       !
 
-      call ocn_thick_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbsorbedRunoff, layerThickness,  &
+      call ocn_thick_surface_flux_tend(meshPool, fractionAbsorbedPrecEvapIce, fractionAbsorbedRunoff, layerThickness,  &
          surfaceThicknessFlux, surfaceThicknessFluxRunoff, tend_layerThickness, err)
 
       !
diff --git a/src/core_ocean/shared/mpas_ocn_thick_surface_flux.F b/src/core_ocean/shared/mpas_ocn_thick_surface_flux.F
index d9caec3625..d4c2213c97 100644
--- a/src/core_ocean/shared/mpas_ocn_thick_surface_flux.F
+++ b/src/core_ocean/shared/mpas_ocn_thick_surface_flux.F
@@ -71,7 +71,7 @@ module ocn_thick_surface_flux
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, transmissionCoefficientsRunoff, &
+   subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficientsPrecEvapIce, transmissionCoefficientsRunoff, &
       layerThickness, surfaceThicknessFlux, surfaceThicknessFluxRunoff, tend, err)!{{{
       !-----------------------------------------------------------------
       !
@@ -83,7 +83,7 @@ subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, trans
          meshPool              !< Input: mesh information
 
       real (kind=RKIND), dimension(:,:), intent(in) :: &
-         transmissionCoefficients,     &!< Input: Coefficients for the transmission of surface fluxes
+         transmissionCoefficientsPrecEvapIce,  &!< Input: Coefficients for transmission of non-runoff surface thickness fluxes
          transmissionCoefficientsRunoff !< Input: Coefficients for the transmission of surface fluxes due to river runoff
 
       real (kind=RKIND), dimension(:,:), intent(in) :: &
@@ -122,7 +122,7 @@ subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, trans
       integer, dimension(:), pointer :: nCellsArray
       integer, dimension(:), pointer :: maxLevelCell
 
-      real (kind=RKIND) :: remainingFlux, remainingFluxRunoff
+      real (kind=RKIND) :: remainingFluxPrecEvapIce, remainingFluxRunoff
 
       err = 0
 
@@ -136,20 +136,21 @@ subroutine ocn_thick_surface_flux_tend(meshPool, transmissionCoefficients, trans
 
       nCells = nCellsArray( 1 )
 
-      !$omp do schedule(runtime) private(remainingFlux, remainingFluxRunoff, k)
+      !$omp do schedule(runtime) private(remainingFluxPrecEvapIce, remainingFluxRunoff, k)
       do iCell = 1, nCells
-        remainingFlux = 1.0_RKIND
+        remainingFluxPrecEvapIce = 1.0_RKIND
         remainingFluxRunoff = 1.0_RKIND
         do k = 1, maxLevelCell(iCell)
-          remainingFlux = remainingFlux - transmissionCoefficients(k, iCell)
+          remainingFluxPrecEvapIce = remainingFluxPrecEvapIce - transmissionCoefficientsPrecEvapIce(k, iCell)
           remainingFluxRunoff = remainingFluxRunoff - transmissionCoefficientsRunoff(k, iCell)
 
-          tend(k, iCell) = tend(k, iCell) + surfaceThicknessFlux(iCell) * transmissionCoefficients(k, iCell)  &
+          tend(k, iCell) = tend(k, iCell) + surfaceThicknessFlux(iCell) * transmissionCoefficientsPrecEvapIce(k, iCell)  &
                                           + surfaceThicknessFluxRunoff(iCell) * transmissionCoefficientsRunoff(k, iCell)
         end do
 
-        if(maxLevelCell(iCell) > 0 .and. remainingFlux > 0.0_RKIND) then
-          tend(maxLevelCell(iCell), iCell) = tend(maxLevelCell(iCell), iCell) + remainingFlux * surfaceThicknessFlux(iCell)
+        if(maxLevelCell(iCell) > 0 .and. remainingFluxPrecEvapIce > 0.0_RKIND) then
+          tend(maxLevelCell(iCell), iCell) = tend(maxLevelCell(iCell), iCell) &
+             + remainingFluxPrecEvapIce * surfaceThicknessFlux(iCell)
         end if
 
         if(maxLevelCell(iCell) > 0 .and. remainingFluxRunoff > 0.0_RKIND) then