diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90
index d3f1484d3..a34125f73 100755
--- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90
+++ b/bldsva/intf_DA/pdaf/framework/init_dim_obs_f_pdaf.F90
@@ -135,6 +135,7 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
   USE enkf_clm_mod, only: domain_def_clm
   USE enkf_clm_mod, only: get_interp_idx
   use enkf_clm_mod, only: clmstatevec_allcol
+  use enkf_clm_mod, only: clmupdate_snow
   !hcp end
 #endif
 #endif
@@ -966,7 +967,12 @@ SUBROUTINE init_dim_obs_f_pdaf(step, dim_obs_f)
                    end if
 #endif
                  else
-                   obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                   if(clmupdate_snow.ne.0) then
+                     ! Snow-DA: no layer in state vector variables
+                     obs_index_p(cnt) = g-begg+1
+                   else
+                     obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                   end if
                  end if
 
                  !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)
diff --git a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90
index e1a6ee4e0..beece1221 100755
--- a/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90
+++ b/bldsva/intf_DA/pdaf/framework/init_dim_obs_pdaf.F90
@@ -131,6 +131,7 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
   USE enkf_clm_mod, only: domain_def_clm
   USE enkf_clm_mod, only: get_interp_idx
   use enkf_clm_mod, only: clmstatevec_allcol
+  use enkf_clm_mod, only: clmupdate_snow
   !hcp end
 #endif
 #endif
@@ -959,7 +960,12 @@ SUBROUTINE init_dim_obs_pdaf(step, dim_obs_p)
                    end if
 #endif
                  else
-                   obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                   if(clmupdate_snow.ne.0) then
+                     ! Snow-DA: no layer in state vector variables
+                     obs_index_p(cnt) = g-begg+1
+                   else
+                     obs_index_p(cnt) = g-begg+1 + ((endg-begg+1) * (clmobs_layer(i)-1))
+                   end if
                  end if
 
                  !write(*,*) 'obs_index_p(',cnt,') is',obs_index_p(cnt)
diff --git a/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90 b/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90
index 5a8114de0..c056f183f 100755
--- a/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90
+++ b/bldsva/intf_DA/pdaf/model/clm3_5/enkf_clm_mod.F90
@@ -65,6 +65,8 @@ module enkf_clm_mod
   integer(c_int),bind(C,name="clmupdate_swc")     :: clmupdate_swc
   integer(c_int),bind(C,name="clmupdate_T")     :: clmupdate_T  ! by hcp
   integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture
+  integer(c_int),bind(C,name="clmupdate_snow")    :: clmupdate_snow
+  integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning
   integer(c_int),bind(C,name="clmprint_swc")      :: clmprint_swc
 #endif
   integer(c_int),bind(C,name="clmprint_et")       :: clmprint_et
@@ -147,6 +149,25 @@ subroutine define_clm_statevec(mype)
       error stop "Not implemented: clmupdate_texture.eq.2"
     endif
 
+    ! Snow assimilation
+    ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5
+    ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell
+    ! Case 2: Assimilation of snow water equivalent same as above
+    if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then 
+      error stop "Not implemented: clmupdate_snow.eq.1 or clmupdate_snow.eq.1"
+    endif
+    ! Case 3: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice
+    if(clmupdate_snow.eq.3) then
+      error stop "Not implemented: clmupdate_snow.eq.3"
+    endif
+    ! Case 4: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq
+    ! and dz
+    if(clmupdate_snow.eq.4) then
+      error stop "Not implemented: clmupdate_snow.eq.4"
+    endif
+
     !hcp LST DA
     if(clmupdate_T.eq.1) then
       clm_varsize      =  endg-begg+1 
@@ -287,6 +308,14 @@ subroutine set_clm_statevec(tstartcycle, mype)
       end do
     endif
 
+    ! Snow assimilation state vector
+    ! Case 1: Snow depth
+    ! Case 2: SWE
+    ! Case 3: Snow depth + SWE
+    if(clmupdate_snow.ne.0) then
+      error stop "Not implemented: clmupdate_snow.ne.0"
+    endif
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN
       ! TSMP-PDAF: For debug runs, output the state vector in files
@@ -534,6 +563,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
       call clm_texture_to_parameters
     endif
 
+    ! Snow assimilation:
+    ! Case 1: Snow depth
+    ! Case 2: Snow water equivalent
+    ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector
+    ! Write updated snow variable back to CLM and then repartition snow and adjust related variables
+    if(clmupdate_snow.ne.0) then
+      error stop "Not implemented: clmupdate_snow.ne.0"
+    endif
+
   end subroutine update_clm
 
   subroutine clm_correct_texture()
diff --git a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90
index f480577b4..58611cd7b 100755
--- a/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90
+++ b/bldsva/intf_DA/pdaf/model/clm5_0/enkf_clm_mod_5.F90
@@ -50,6 +50,8 @@ module enkf_clm_mod
   integer(c_int),bind(C,name="clmupdate_swc")     :: clmupdate_swc
   integer(c_int),bind(C,name="clmupdate_T")     :: clmupdate_T  ! by hcp
   integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture
+  integer(c_int),bind(C,name="clmupdate_snow")    :: clmupdate_snow
+  integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning
   integer(c_int),bind(C,name="clmprint_swc")      :: clmprint_swc
 #endif
   integer(c_int),bind(C,name="clmprint_et")       :: clmprint_et
@@ -271,6 +273,48 @@ subroutine define_clm_statevec(mype)
         clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi)
     endif
 
+    ! Snow assimilation
+    ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5
+    ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell
+    ! Case 2: Assimilation of snow water equivalent same as above
+    if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then 
+        clm_varsize      =  (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA
+        clm_statevecsize =  (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority
+    endif
+    ! Case 3: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice
+    if(clmupdate_snow.eq.3) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+    ! Case 4: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq
+    ! and dz
+    if(clmupdate_snow.eq.4) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+    ! Case 5: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice
+    ! and dz
+    if(clmupdate_snow.eq.5) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+    ! Case 6: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq
+    if(clmupdate_snow.eq.6) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+    ! Case 7: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice
+    ! Should reproduce case 3
+    if(clmupdate_snow.eq.7) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+
     !hcp LST DA
     if(clmupdate_T.eq.1) then
       error stop "Not implemented: clmupdate_T.eq.1"
@@ -323,6 +367,8 @@ subroutine set_clm_statevec(tstartcycle, mype)
     real(r8), pointer :: psand(:,:)
     real(r8), pointer :: pclay(:,:)
     real(r8), pointer :: porgm(:,:)
+    real(r8), pointer :: snow_depth(:)
+    real(r8), pointer :: h2osno(:)
     integer :: i,j,jj,g,cc=0,offset=0
     character (len = 34) :: fn    !TSMP-PDAF: function name for state vector output
     character (len = 34) :: fn2    !TSMP-PDAF: function name for swc output
@@ -332,6 +378,9 @@ subroutine set_clm_statevec(tstartcycle, mype)
     pclay => soilstate_inst%cellclay_col
     porgm => soilstate_inst%cellorg_col
 
+    snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) 
+    h2osno     => waterstate_inst%h2osno_col     ! snow water equivalent (mm)
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN
 
@@ -385,6 +434,53 @@ subroutine set_clm_statevec(tstartcycle, mype)
       end do
     endif
 
+    ! Snow assimilation state vector
+    ! Case 1: Snow depth
+    ! Case 2: SWE
+    ! Case 3: Snow depth + SWE
+    if(clmupdate_snow.ne.0) then
+        cc = 1
+        do j=clm_begg,clm_endg
+        ! Only get the snow_depth/swe from the first column of each gridcell
+        ! and add it to the clm_statevec at the position of the gridcell (cc)
+        newgridcell = .true.
+        do jj=clm_begc,clm_endc
+          g = col%gridcell(jj)
+          if (g .eq. j) then
+            if (newgridcell) then
+              newgridcell = .false.
+
+              if(clmupdate_snow.eq.1) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+              else if(clmupdate_snow.eq.2) then
+                clm_statevec(cc+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.3) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.4) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.5) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.6) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.7) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else
+                error stop "Wrong input for clmupdate_snow"
+              end if
+
+            endif
+          endif 
+        end do
+        cc = cc + 1
+        end do
+
+    endif
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN
       ! TSMP-PDAF: For debug runs, output the state vector in files
@@ -420,10 +516,21 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     real(r8), pointer :: pclay(:,:)
     real(r8), pointer :: porgm(:,:)
 
+    real(r8), pointer :: snow_depth(:)
+    real(r8), pointer :: h2osno(:)
+
     real(r8), pointer :: dz(:,:)          ! layer thickness depth (m)
     real(r8), pointer :: h2osoi_liq(:,:)  ! liquid water (kg/m2)
     real(r8), pointer :: h2osoi_ice(:,:)
+
     real(r8)  :: rliq,rice
+    real(r8)  :: incr_sno
+    real(r8)  :: incr_sd
+    real(r8)  :: incr_swe
+    real(r8)  :: h2osno_in(clm_begc:clm_endc)
+    real(r8)  :: snow_depth_in(clm_begc:clm_endc)
+    real(r8)  :: h2osno_out(clm_begc:clm_endc)
+    real(r8)  :: snow_depth_out(clm_begc:clm_endc)
     real(r8)  :: watmin_check      ! minimum soil moisture for checking clm_statevec (mm)
     real(r8)  :: watmin_set        ! minimum soil moisture for setting swc (mm)
     real(r8)  :: swc_update        ! updated SWC in loop
@@ -436,6 +543,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     character (len = 32) :: fn5    !TSMP-PDAF: function name for state vector outpu
     character (len = 32) :: fn6    !TSMP-PDAF: function name for state vector outpu
 
+    integer, pointer :: snlsno(:)
     logical :: swc_zero_before_update = .false.
 
 #ifdef PDAF_DEBUG
@@ -456,9 +564,13 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     pclay => soilstate_inst%cellclay_col
     porgm => soilstate_inst%cellorg_col
 
+    snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m)
+    h2osno     => waterstate_inst%h2osno_col     ! snow water equivalent (mm) 
+
     dz            => col%dz
     h2osoi_liq    => waterstate_inst%h2osoi_liq_col
     h2osoi_ice    => waterstate_inst%h2osoi_ice_col
+    snlsno     => col%snl                        ! number of snow layers (negative)
 
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN
@@ -469,7 +581,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         OPEN(unit=71, file=fn5, action="write")
         WRITE (71,"(es22.15)") h2osoi_liq(:,:)
         CLOSE(71)
+      END IF
 
+      IF(clmupdate_swc.NE.0 .OR. clmupdate_snow.NE.0) THEN
         ! TSMP-PDAF: For debug runs, output the state vector in files
         WRITE(fn6, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".bef_up.", tstartcycle, ".txt"
         OPEN(unit=71, file=fn6, action="write")
@@ -635,8 +749,378 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
       call clm_texture_to_parameters
     endif
 
+    ! Snow assimilation:
+    ! Case 1: Snow depth
+    ! Case 2: Snow water equivalent
+    ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector
+    ! Write updated snow variable back to CLM and then repartition snow and adjust related variables
+    if(clmupdate_snow.ne.0) then
+      do j=clm_begc,clm_endc
+        ! Iterate through the columns
+
+        ! For each column index, copy the gridcell-linked state vector
+        ! value at state vector index `cc`
+
+        ! Set cc (the state vector index) from the
+        ! CLM5-grid-index
+        cc = (col%gridcell(j) - clm_begg + 1)
+
+        if (clmupdate_snow.eq.1) then
+          snow_depth_in(j) = snow_depth(j)
+        else if (clmupdate_snow.eq.2) then
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.3) then
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.4) then
+          snow_depth_in(j) = snow_depth(j)
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.5) then
+          snow_depth_in(j) = snow_depth(j)
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.6) then
+          snow_depth_in(j) = snow_depth(j)
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.7) then
+          snow_depth_in(j) = snow_depth(j)
+          h2osno_in(j) = h2osno(j)
+        end if
+
+        if (clmupdate_snow.eq.1) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+        else if (clmupdate_snow.eq.2) then
+          h2osno_out(j) = clm_statevec(cc+offset)
+        else if (clmupdate_snow.eq.3) then
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        else if (clmupdate_snow.eq.4) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        else if (clmupdate_snow.eq.5) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        else if (clmupdate_snow.eq.6) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        else if (clmupdate_snow.eq.7) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        end if
+
+        if(clmupdate_snow.eq.1) then
+          ! Update state variable to CLM
+          ! Needed for Case 1/2 if they use repartioning function
+          if (snow_depth_out(j).gt.0.0) then
+            snow_depth(j) = snow_depth_out(j)
+          else
+            ! Catch negative or 0 values from DA
+            print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j)
+          end if
+        else if(clmupdate_snow.eq.2) then
+          if (h2osno_out(j).gt.0.0) then
+            h2osno(j) = h2osno_out(j)
+          else
+            ! Catch negative or 0 values from DA
+            print *, "WARNING: SWE is negative/zero at cc. cc, j, offset, h2osno(j): ", cc, j, offset, h2osno(j)
+          end if
+        end if
+
+      end do
+
+      ! Repartitioning
+      if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then
+
+        if (clmupdate_snow.eq.1) then
+          if ( SUM(ABS(snow_depth_in(:) - snow_depth_out(:))).gt.0.000001) then
+            call clm_repartition_snow(snow_depth_in(:))
+          end if
+        end if
+
+        if (clmupdate_snow.eq.2) then
+          if ( SUM(ABS(h2osno_in(:) - h2osno_out(:))).gt.0.000001) then
+            call clm_repartition_snow(h2osno_in(:))
+          end if
+        end if
+
+        if (clmupdate_snow.ge.3) then
+          error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.ge.3"
+        end if
+
+      end if
+
+      if ( clmupdate_snow_repartitioning.eq.3) then
+
+        if (clmupdate_snow.eq.1) then
+          do j=clm_begc,clm_endc
+            if (snow_depth_out(j).gt.1.0d-6) then
+                if (snow_depth_in(j).gt.1.0d-6) then
+                  ! Update h2osoi_ice with increment
+                  incr_sno = snow_depth_out(j) / snow_depth_in(j)
+                  do i=snlsno(j)+1,0
+                    h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno
+                  end do
+                end if
+            end if
+          end do
+        end if
+
+        if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then
+          do j=clm_begc,clm_endc
+            if (h2osno_out(j).gt.1.0d-6) then
+                if (h2osno_in(j).gt.1.0d-6) then
+                  ! Update h2osoi_ice with increment
+                  incr_sno = h2osno_out(j) / h2osno_in(j)
+                  do i=snlsno(j)+1,0
+                    h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno
+                  end do
+                end if
+            end if
+          end do
+        end if
+
+        if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then
+          do j=clm_begc,clm_endc
+            if (h2osno_out(j).gt.1.0d-6) then
+              if (h2osno_in(j).gt.1.0d-6) then
+                if (snow_depth_out(j).gt.1.0d-6) then
+                  if (snow_depth_in(j).gt.1.0d-6) then
+                    ! Update h2osoi_ice/h2osoi_liq with increment
+                    incr_swe = h2osno_out(j) / h2osno_in(j)
+                    ! Update snow_depth with increment
+                    incr_sd = snow_depth_out(j) / snow_depth_in(j)
+                    do i=snlsno(j)+1,0
+
+                      if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then
+                        h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe
+                        if (isnan(h2osoi_ice(j,i))) then
+                          print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i
+                        endif
+                      end if
+
+                      if (clmupdate_snow.eq.4                          .or. clmupdate_snow.eq.6) then
+                        h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe
+                        if (isnan(h2osoi_liq(j,i))) then
+                          print *, "WARNING: h2osoi_liq at j,i is nan: ", j, i
+                        endif
+                      end if
+
+                      if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5) then
+                        dz(j,i) = dz(j,i) * incr_sd
+                        if (isnan(dz(j,i))) then
+                          print *, "WARNING: dz at j,i is nan: ", j, i
+                        endif
+                      end if
+
+                    end do
+                  end if
+                end if
+              end if
+            end if
+          end do
+        end if
+
+      end if
+
+#ifdef PDAF_DEBUG
+      IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN
+
+        ! TSMP-PDAF: For debug runs, output the state vector in files
+        WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt"
+        OPEN(unit=71, file=fn4, action="write")
+        WRITE (71,"(es22.15)") h2osoi_ice(:,:)
+        CLOSE(71)
+
+      END IF
+#endif
+    endif
+
   end subroutine update_clm
 
+  subroutine clm_repartition_snow(h2osno_in)
+    use ColumnType,    only : col
+    use clm_instMod,   only : waterstate_inst
+    use clm_varpar,    only : nlevsno, nlevsoi
+    use clm_varcon,    only : bdsno, denice
+    use shr_kind_mod,  only : r8 => shr_kind_r8
+
+    implicit none
+
+    real(r8), intent(in), optional :: h2osno_in(clm_begc:clm_endc)
+
+    real(r8), pointer :: snow_depth(:)
+    real(r8), pointer :: h2osno(:)
+    real(r8), pointer :: h2oliq(:,:)
+    real(r8), pointer :: h2oice(:,:)
+    real(r8), pointer :: frac_sno(:)
+    real(r8), pointer :: snowdp(:)
+    integer, pointer :: snlsno(:)
+
+    real(r8) :: dzsno(clm_begc:clm_endc,-nlevsno+1:0)
+    real(r8) :: h2osno_po(clm_begc:clm_endc)
+    real(r8) :: h2osno_pr(clm_begc:clm_endc)
+    real(r8) :: snowden, frac_swe, frac_liq, frac_ice
+    real(r8) :: gain_h2osno, gain_h2oliq, gain_h2oice
+    real(r8) :: gain_dzsno(-nlevsno+1:0)
+    real(r8) :: rho_avg, z_avg
+    integer :: i,ii,j,jj,g,cc=1,offset=0
+
+    snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m)
+    snowdp     => waterstate_inst%snowdp_col     ! area-averaged snow height (m) 
+    h2osno     => waterstate_inst%h2osno_col     ! col snow water (mm H2O)
+    h2oliq     => waterstate_inst%h2osoi_liq_col ! col liquid water (kg/m2) (-nlevsno+1:nlevgrnd) 
+    h2oice     => waterstate_inst%h2osoi_ice_col ! col ice lens (kg/m2) (-nlevsno+1:nlevgrnd)
+
+    snlsno     => col%snl                        ! number of snow layers (negative)
+
+    frac_sno   => waterstate_inst%frac_sno_eff_col ! fraction of ground covered by snow 
+    ! dz for snow layers is defined like in the history output as col%dz for the snow layers
+    dzsno(clm_begc:clm_endc, -nlevsno+1:0) = col%dz(clm_begc:clm_endc,-nlevsno+1:0)
+    ! Iterate through all columns
+    do jj=clm_begc,clm_endc
+      if (h2osno(jj).lt.0.0) then ! No snow in column
+        print *, "WARNING: negative snow in col: ", jj, h2osno
+!        ! Set existing layers to near zero and let CLM do the layer aggregation
+        do i=0,snlsno(jj)+1,-1
+            h2oliq(jj,i) = 0.0_r8
+            h2oice(jj,i) = 0.00000001_r8
+            dzsno(jj,i)  = 0.00000001_r8
+        end do
+        snow_depth(jj) = sum(dzsno(jj,:))
+        h2osno(jj) = sum(h2oice(jj,:))
+
+        if (isnan(h2osno(jj))) then
+          print *, "WARNING: h2osno at jj is nan: ", jj
+        endif
+        if (isnan(snow_depth(jj))) then
+          print *, "WARNING: snow_depth at jj is nan: ", jj
+        endif
+
+      else ! snow (h2osno) present
+        if (snlsno(jj).lt.0) then ! snow layers in the column
+          if (clmupdate_snow .eq. 1) then
+            ! DART source: https://github.com/NCAR/DART/blob/main/models/clm/dart_to_clm.f90
+            ! Formulas below from DART use h2osno_po / h2osno_pr for after / before DA SWE
+            ! clmupdate_snow == 1 has snow_depth after and h2osno before DA snow depth
+            ! Therefore need to have a transform to get h2osno_po
+            ! v1 init
+            ! h2osno_po(jj) = (snow_depth(jj) * bdsno) ! calculations from Init using constant SBD
+            ! v2 SoilTemperatureMod
+            if (snowdp(jj).gt.0.0_r8) then
+              rho_avg = min(800.0_r8, h2osno(jj)/snowdp(jj))
+            else
+              rho_avg = 200.0_r8
+            end if
+            if (frac_sno(jj).gt.0.0_r8 .and. snlsno(jj).lt.0.0_r8) then
+              h2osno_po(jj) = snow_depth(jj) * (rho_avg*frac_sno(jj))
+            else
+              h2osno_po(jj) = snow_depth(jj) * denice
+            end if
+            h2osno_pr(jj) = h2osno(jj)
+          else if (clmupdate_snow .eq. 2) then
+            ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already
+            h2osno_po(jj) = h2osno(jj)
+            h2osno_pr(jj) = h2osno_in(jj)
+          end if
+
+          do ii=0,snlsno(jj)+1,-1 ! iterate through the snow layers
+            ! DART VERSION: ii = nlevsoi - i + 1
+            ! snow density prior for each layer
+            if (dzsno(jj,ii).gt.0.0_r8) then
+              snowden = (h2oliq(jj,ii) + h2oice(jj,ii) / dzsno(jj,ii))
+            else
+              snowden = 0.0_r8
+            endif
+            ! fraction of SWE in each active layers
+            if(h2osno_pr(jj).gt.0.0_r8) then
+              ! repartition Option 1: Adjust bottom layer only (set weight to 1 for bottom 0 else)
+              if(clmupdate_snow_repartitioning.eq.1) then
+                if (ii .eq. 0) then ! DART version indexing check against nlevsno but for us 0
+                  frac_swe = 1.0_r8
+                  ! JK: Let CLM repartitioning do the job
+                  ! afterwards. Provide all the snow in a single layer
+                else
+                  frac_swe = 0.0_r8
+                end if
+              ! repartition Option 2: Adjust all active layers
+              else if (clmupdate_snow_repartitioning.eq.2) then
+                frac_swe = (h2oliq(jj,ii) + h2oice(jj,ii)) / h2osno_pr(jj)
+              end if
+            else
+              frac_swe = 0.0_r8 ! no fraction SWE if no snow is present in column
+            end if ! end SWE fraction if
+
+            ! fraction of liquid and ice
+            if ((h2oliq(jj,ii) + h2oice(jj,ii)).gt.0.0_r8) then
+                frac_liq = h2oliq(jj,ii) / (h2oliq(jj,ii) + h2oice(jj,ii))
+                frac_ice = 1.0_r8 - frac_liq
+            else
+                frac_liq = 0.0_r8
+                frac_ice = 0.0_r8
+            end if
+
+            ! SWE adjustment per layer 
+            ! assumes identical layer distribution of liq and ice than before DA (frac_*)
+            if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then
+              gain_h2osno = (h2osno_po(jj) - h2osno_pr(jj)) * frac_swe
+              gain_h2oliq = gain_h2osno * frac_liq
+              gain_h2oice = gain_h2osno * frac_ice
+            else
+              gain_h2osno = 0.0_r8
+              gain_h2oliq = 0.0_r8
+              gain_h2oice = 0.0_r8
+            end if
+            ! layer level adjustments
+            if (snowden.gt.0.0_r8) then
+              gain_dzsno(ii) = gain_h2osno / snowden
+            else
+              gain_dzsno(ii) = 0.0_r8
+            end if
+            h2oliq(jj,ii) = h2oliq(jj,ii) + gain_h2oliq
+            h2oice(jj,ii) = h2oice(jj,ii) + gain_h2oice
+
+            ! Adjust snow layer dimensions so that CLM5 can calculate compaction / aggregation
+            ! in the DART code dzsno is adjusted directly but in CLM5 dzsno is local and diagnostic
+            ! i.e. calculated / assigned from frac_sno and dz(:, snow_layer) in SnowHydrologyMod
+            ! therefore we adjust dz(:, snow_layer) here
+            if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then
+              col%dz(jj, ii) = col%dz(jj, ii) + gain_dzsno(ii)
+              ! mid point and interface adjustments
+              ! i.e. zsno (col%z(:, snow_layers)) and zisno (col%zi(:, snow_layers))
+              ! DART version the sum goes from ilevel:nlevsno to fit with our indexing:
+              col%zi(jj, ii) = sum(col%dz(jj,ii:0))*-1.0_r8
+              ! In DART the check is ilevel == nlevsno but here
+              if (ii.eq.0) then
+                col%z(jj,ii) = col%zi(jj,ii) / 2.0_r8
+              else
+                col%z(jj,ii) = sum(col%zi(jj, ii:ii+1)) / 2.0_r8
+              end if
+            end if
+
+          end do ! End iteration of snow layers
+        endif ! End of snow layers present check
+
+        ! Column level variables
+        if (clmupdate_snow .eq. 1) then
+          ! Finally adjust SWE (h2osno) since the prior value is no longer needed
+          ! column level variables so we can adjust it outside the layer loop
+          h2osno(jj) = h2osno_po(jj) 
+        else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then
+          ! Update the total snow depth to match udpates to layers for active snow layers
+          if (abs(h2osno_po(jj) - h2osno_pr(jj)) .gt. 0.0_r8 .and. snlsno(jj) < 0.0_r8) then
+            snow_depth(jj) = snow_depth(jj) + sum(gain_dzsno(-nlevsno+1:0))
+          end if
+        end if
+
+        if (isnan(h2osno(jj))) then
+          print *, "WARNING2: h2osno at jj is nan: ", jj
+        endif
+        if (isnan(snow_depth(jj))) then
+          print *, "WARNING2: snow_depth at jj is nan: ", jj
+        endif
+
+      end if ! End of snow present check
+    end do ! End of column iteration
+
+  end subroutine clm_repartition_snow
+
   subroutine clm_correct_texture()
 
     use clm_varpar   , only : nlevsoi
diff --git a/bldsva/intf_DA/pdaf/model/common/enkf.h b/bldsva/intf_DA/pdaf/model/common/enkf.h
index fba54ce56..3e0937005 100755
--- a/bldsva/intf_DA/pdaf/model/common/enkf.h
+++ b/bldsva/intf_DA/pdaf/model/common/enkf.h
@@ -87,6 +87,8 @@ GLOBAL int nx_local,ny_local,nz_local;
 GLOBAL int clmupdate_swc;
 GLOBAL int clmupdate_T;
 GLOBAL int clmupdate_texture;
+GLOBAL int clmupdate_snow;
+GLOBAL int clmupdate_snow_repartitioning;
 GLOBAL int clmprint_swc;
 GLOBAL int clmprint_et;
 GLOBAL int clmstatevec_allcol;
diff --git a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c
index dc8046a55..22958091f 100755
--- a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c
+++ b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c
@@ -78,6 +78,8 @@ void read_enkfpar(char *parname)
   clmupdate_swc         = iniparser_getint(pardict,"CLM:update_swc",1);
   clmupdate_T           = iniparser_getint(pardict,"CLM:update_T",0);
   clmupdate_texture     = iniparser_getint(pardict,"CLM:update_texture",0);
+  clmupdate_snow        = iniparser_getint(pardict,"CLM:update_snow",0);
+  clmupdate_snow_repartitioning = iniparser_getint(pardict,"CLM:update_snow_repartitioning",3);
   clmprint_swc          = iniparser_getint(pardict,"CLM:print_swc",0);
   clmprint_et           = iniparser_getint(pardict,"CLM:print_et",0);
   clmstatevec_allcol    = iniparser_getint(pardict,"CLM:statevec_allcol",0);
diff --git a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90 b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90
index f480577b4..58611cd7b 100755
--- a/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90
+++ b/bldsva/intf_DA/pdaf/model/eclm/enkf_clm_mod_5.F90
@@ -50,6 +50,8 @@ module enkf_clm_mod
   integer(c_int),bind(C,name="clmupdate_swc")     :: clmupdate_swc
   integer(c_int),bind(C,name="clmupdate_T")     :: clmupdate_T  ! by hcp
   integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture
+  integer(c_int),bind(C,name="clmupdate_snow")    :: clmupdate_snow
+  integer(c_int),bind(C,name="clmupdate_snow_repartitioning") :: clmupdate_snow_repartitioning
   integer(c_int),bind(C,name="clmprint_swc")      :: clmprint_swc
 #endif
   integer(c_int),bind(C,name="clmprint_et")       :: clmprint_et
@@ -271,6 +273,48 @@ subroutine define_clm_statevec(mype)
         clm_statevecsize = clm_statevecsize + 3*((endg-begg+1)*nlevsoi)
     endif
 
+    ! Snow assimilation
+    ! Case 1: Assimilation of snow depth : allocated 1 per column in CLM5
+    ! But observations and history file 1 per grid cell and therefore statevecsize 1 per grid cell
+    ! Case 2: Assimilation of snow water equivalent same as above
+    if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then 
+        clm_varsize      =  (clm_endg-clm_begg+1) ! Currently no combination of SWC and snow DA
+        clm_statevecsize =  (clm_endg-clm_begg+1) ! So like this if snow is set it takes priority
+    endif
+    ! Case 3: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice
+    if(clmupdate_snow.eq.3) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+    ! Case 4: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq
+    ! and dz
+    if(clmupdate_snow.eq.4) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+    ! Case 5: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice
+    ! and dz
+    if(clmupdate_snow.eq.5) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+    ! Case 6: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice, h2osoi_liq
+    if(clmupdate_snow.eq.6) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+    ! Case 7: Assimilation of snow depth: Snow depth and snow water
+    ! equivalent in the state vector. Update of h2osoi_ice
+    ! Should reproduce case 3
+    if(clmupdate_snow.eq.7) then
+        clm_varsize      =  (clm_endg-clm_begg+1)
+        clm_statevecsize =  2*(clm_endg-clm_begg+1)
+    endif
+
     !hcp LST DA
     if(clmupdate_T.eq.1) then
       error stop "Not implemented: clmupdate_T.eq.1"
@@ -323,6 +367,8 @@ subroutine set_clm_statevec(tstartcycle, mype)
     real(r8), pointer :: psand(:,:)
     real(r8), pointer :: pclay(:,:)
     real(r8), pointer :: porgm(:,:)
+    real(r8), pointer :: snow_depth(:)
+    real(r8), pointer :: h2osno(:)
     integer :: i,j,jj,g,cc=0,offset=0
     character (len = 34) :: fn    !TSMP-PDAF: function name for state vector output
     character (len = 34) :: fn2    !TSMP-PDAF: function name for swc output
@@ -332,6 +378,9 @@ subroutine set_clm_statevec(tstartcycle, mype)
     pclay => soilstate_inst%cellclay_col
     porgm => soilstate_inst%cellorg_col
 
+    snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m) 
+    h2osno     => waterstate_inst%h2osno_col     ! snow water equivalent (mm)
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN
 
@@ -385,6 +434,53 @@ subroutine set_clm_statevec(tstartcycle, mype)
       end do
     endif
 
+    ! Snow assimilation state vector
+    ! Case 1: Snow depth
+    ! Case 2: SWE
+    ! Case 3: Snow depth + SWE
+    if(clmupdate_snow.ne.0) then
+        cc = 1
+        do j=clm_begg,clm_endg
+        ! Only get the snow_depth/swe from the first column of each gridcell
+        ! and add it to the clm_statevec at the position of the gridcell (cc)
+        newgridcell = .true.
+        do jj=clm_begc,clm_endc
+          g = col%gridcell(jj)
+          if (g .eq. j) then
+            if (newgridcell) then
+              newgridcell = .false.
+
+              if(clmupdate_snow.eq.1) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+              else if(clmupdate_snow.eq.2) then
+                clm_statevec(cc+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.3) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.4) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.5) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.6) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else if(clmupdate_snow.eq.7) then
+                clm_statevec(cc+offset) = snow_depth(jj)
+                clm_statevec(cc+clm_varsize+offset) = h2osno(jj)
+              else
+                error stop "Wrong input for clmupdate_snow"
+              end if
+
+            endif
+          endif 
+        end do
+        cc = cc + 1
+        end do
+
+    endif
+
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle + 1 .OR. clmt_printensemble < 0) THEN
       ! TSMP-PDAF: For debug runs, output the state vector in files
@@ -420,10 +516,21 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     real(r8), pointer :: pclay(:,:)
     real(r8), pointer :: porgm(:,:)
 
+    real(r8), pointer :: snow_depth(:)
+    real(r8), pointer :: h2osno(:)
+
     real(r8), pointer :: dz(:,:)          ! layer thickness depth (m)
     real(r8), pointer :: h2osoi_liq(:,:)  ! liquid water (kg/m2)
     real(r8), pointer :: h2osoi_ice(:,:)
+
     real(r8)  :: rliq,rice
+    real(r8)  :: incr_sno
+    real(r8)  :: incr_sd
+    real(r8)  :: incr_swe
+    real(r8)  :: h2osno_in(clm_begc:clm_endc)
+    real(r8)  :: snow_depth_in(clm_begc:clm_endc)
+    real(r8)  :: h2osno_out(clm_begc:clm_endc)
+    real(r8)  :: snow_depth_out(clm_begc:clm_endc)
     real(r8)  :: watmin_check      ! minimum soil moisture for checking clm_statevec (mm)
     real(r8)  :: watmin_set        ! minimum soil moisture for setting swc (mm)
     real(r8)  :: swc_update        ! updated SWC in loop
@@ -436,6 +543,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     character (len = 32) :: fn5    !TSMP-PDAF: function name for state vector outpu
     character (len = 32) :: fn6    !TSMP-PDAF: function name for state vector outpu
 
+    integer, pointer :: snlsno(:)
     logical :: swc_zero_before_update = .false.
 
 #ifdef PDAF_DEBUG
@@ -456,9 +564,13 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
     pclay => soilstate_inst%cellclay_col
     porgm => soilstate_inst%cellorg_col
 
+    snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m)
+    h2osno     => waterstate_inst%h2osno_col     ! snow water equivalent (mm) 
+
     dz            => col%dz
     h2osoi_liq    => waterstate_inst%h2osoi_liq_col
     h2osoi_ice    => waterstate_inst%h2osoi_ice_col
+    snlsno     => col%snl                        ! number of snow layers (negative)
 
 #ifdef PDAF_DEBUG
     IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN
@@ -469,7 +581,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
         OPEN(unit=71, file=fn5, action="write")
         WRITE (71,"(es22.15)") h2osoi_liq(:,:)
         CLOSE(71)
+      END IF
 
+      IF(clmupdate_swc.NE.0 .OR. clmupdate_snow.NE.0) THEN
         ! TSMP-PDAF: For debug runs, output the state vector in files
         WRITE(fn6, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".bef_up.", tstartcycle, ".txt"
         OPEN(unit=71, file=fn6, action="write")
@@ -635,8 +749,378 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
       call clm_texture_to_parameters
     endif
 
+    ! Snow assimilation:
+    ! Case 1: Snow depth
+    ! Case 2: Snow water equivalent
+    ! Case 3: Snow depth (assimilated) and SWE (used for increment) in state vector
+    ! Write updated snow variable back to CLM and then repartition snow and adjust related variables
+    if(clmupdate_snow.ne.0) then
+      do j=clm_begc,clm_endc
+        ! Iterate through the columns
+
+        ! For each column index, copy the gridcell-linked state vector
+        ! value at state vector index `cc`
+
+        ! Set cc (the state vector index) from the
+        ! CLM5-grid-index
+        cc = (col%gridcell(j) - clm_begg + 1)
+
+        if (clmupdate_snow.eq.1) then
+          snow_depth_in(j) = snow_depth(j)
+        else if (clmupdate_snow.eq.2) then
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.3) then
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.4) then
+          snow_depth_in(j) = snow_depth(j)
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.5) then
+          snow_depth_in(j) = snow_depth(j)
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.6) then
+          snow_depth_in(j) = snow_depth(j)
+          h2osno_in(j) = h2osno(j)
+        else if (clmupdate_snow.eq.7) then
+          snow_depth_in(j) = snow_depth(j)
+          h2osno_in(j) = h2osno(j)
+        end if
+
+        if (clmupdate_snow.eq.1) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+        else if (clmupdate_snow.eq.2) then
+          h2osno_out(j) = clm_statevec(cc+offset)
+        else if (clmupdate_snow.eq.3) then
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        else if (clmupdate_snow.eq.4) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        else if (clmupdate_snow.eq.5) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        else if (clmupdate_snow.eq.6) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        else if (clmupdate_snow.eq.7) then
+          snow_depth_out(j) = clm_statevec(cc+offset)
+          h2osno_out(j) = clm_statevec(cc+clm_varsize+offset)
+        end if
+
+        if(clmupdate_snow.eq.1) then
+          ! Update state variable to CLM
+          ! Needed for Case 1/2 if they use repartioning function
+          if (snow_depth_out(j).gt.0.0) then
+            snow_depth(j) = snow_depth_out(j)
+          else
+            ! Catch negative or 0 values from DA
+            print *, "WARNING: Snow-depth is negative/zero at cc. cc, j, offset, snow_depth_out(j): ", cc, j, offset, snow_depth_out(j)
+          end if
+        else if(clmupdate_snow.eq.2) then
+          if (h2osno_out(j).gt.0.0) then
+            h2osno(j) = h2osno_out(j)
+          else
+            ! Catch negative or 0 values from DA
+            print *, "WARNING: SWE is negative/zero at cc. cc, j, offset, h2osno(j): ", cc, j, offset, h2osno(j)
+          end if
+        end if
+
+      end do
+
+      ! Repartitioning
+      if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then
+
+        if (clmupdate_snow.eq.1) then
+          if ( SUM(ABS(snow_depth_in(:) - snow_depth_out(:))).gt.0.000001) then
+            call clm_repartition_snow(snow_depth_in(:))
+          end if
+        end if
+
+        if (clmupdate_snow.eq.2) then
+          if ( SUM(ABS(h2osno_in(:) - h2osno_out(:))).gt.0.000001) then
+            call clm_repartition_snow(h2osno_in(:))
+          end if
+        end if
+
+        if (clmupdate_snow.ge.3) then
+          error stop "Not implemented: Repartioning 1/2 for clmupdate_snow.ge.3"
+        end if
+
+      end if
+
+      if ( clmupdate_snow_repartitioning.eq.3) then
+
+        if (clmupdate_snow.eq.1) then
+          do j=clm_begc,clm_endc
+            if (snow_depth_out(j).gt.1.0d-6) then
+                if (snow_depth_in(j).gt.1.0d-6) then
+                  ! Update h2osoi_ice with increment
+                  incr_sno = snow_depth_out(j) / snow_depth_in(j)
+                  do i=snlsno(j)+1,0
+                    h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno
+                  end do
+                end if
+            end if
+          end do
+        end if
+
+        if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then
+          do j=clm_begc,clm_endc
+            if (h2osno_out(j).gt.1.0d-6) then
+                if (h2osno_in(j).gt.1.0d-6) then
+                  ! Update h2osoi_ice with increment
+                  incr_sno = h2osno_out(j) / h2osno_in(j)
+                  do i=snlsno(j)+1,0
+                    h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno
+                  end do
+                end if
+            end if
+          end do
+        end if
+
+        if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then
+          do j=clm_begc,clm_endc
+            if (h2osno_out(j).gt.1.0d-6) then
+              if (h2osno_in(j).gt.1.0d-6) then
+                if (snow_depth_out(j).gt.1.0d-6) then
+                  if (snow_depth_in(j).gt.1.0d-6) then
+                    ! Update h2osoi_ice/h2osoi_liq with increment
+                    incr_swe = h2osno_out(j) / h2osno_in(j)
+                    ! Update snow_depth with increment
+                    incr_sd = snow_depth_out(j) / snow_depth_in(j)
+                    do i=snlsno(j)+1,0
+
+                      if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5 .or. clmupdate_snow.eq.6 .or. clmupdate_snow.eq.7) then
+                        h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe
+                        if (isnan(h2osoi_ice(j,i))) then
+                          print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i
+                        endif
+                      end if
+
+                      if (clmupdate_snow.eq.4                          .or. clmupdate_snow.eq.6) then
+                        h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe
+                        if (isnan(h2osoi_liq(j,i))) then
+                          print *, "WARNING: h2osoi_liq at j,i is nan: ", j, i
+                        endif
+                      end if
+
+                      if (clmupdate_snow.eq.4 .or. clmupdate_snow.eq.5) then
+                        dz(j,i) = dz(j,i) * incr_sd
+                        if (isnan(dz(j,i))) then
+                          print *, "WARNING: dz at j,i is nan: ", j, i
+                        endif
+                      end if
+
+                    end do
+                  end if
+                end if
+              end if
+            end if
+          end do
+        end if
+
+      end if
+
+#ifdef PDAF_DEBUG
+      IF(clmt_printensemble == tstartcycle .OR. clmt_printensemble < 0) THEN
+
+        ! TSMP-PDAF: For debug runs, output the state vector in files
+        WRITE(fn4, "(a,i5.5,a,i5.5,a)") "h2osoi_ice", mype, ".update.", tstartcycle, ".txt"
+        OPEN(unit=71, file=fn4, action="write")
+        WRITE (71,"(es22.15)") h2osoi_ice(:,:)
+        CLOSE(71)
+
+      END IF
+#endif
+    endif
+
   end subroutine update_clm
 
+  subroutine clm_repartition_snow(h2osno_in)
+    use ColumnType,    only : col
+    use clm_instMod,   only : waterstate_inst
+    use clm_varpar,    only : nlevsno, nlevsoi
+    use clm_varcon,    only : bdsno, denice
+    use shr_kind_mod,  only : r8 => shr_kind_r8
+
+    implicit none
+
+    real(r8), intent(in), optional :: h2osno_in(clm_begc:clm_endc)
+
+    real(r8), pointer :: snow_depth(:)
+    real(r8), pointer :: h2osno(:)
+    real(r8), pointer :: h2oliq(:,:)
+    real(r8), pointer :: h2oice(:,:)
+    real(r8), pointer :: frac_sno(:)
+    real(r8), pointer :: snowdp(:)
+    integer, pointer :: snlsno(:)
+
+    real(r8) :: dzsno(clm_begc:clm_endc,-nlevsno+1:0)
+    real(r8) :: h2osno_po(clm_begc:clm_endc)
+    real(r8) :: h2osno_pr(clm_begc:clm_endc)
+    real(r8) :: snowden, frac_swe, frac_liq, frac_ice
+    real(r8) :: gain_h2osno, gain_h2oliq, gain_h2oice
+    real(r8) :: gain_dzsno(-nlevsno+1:0)
+    real(r8) :: rho_avg, z_avg
+    integer :: i,ii,j,jj,g,cc=1,offset=0
+
+    snow_depth => waterstate_inst%snow_depth_col ! snow height of snow covered area (m)
+    snowdp     => waterstate_inst%snowdp_col     ! area-averaged snow height (m) 
+    h2osno     => waterstate_inst%h2osno_col     ! col snow water (mm H2O)
+    h2oliq     => waterstate_inst%h2osoi_liq_col ! col liquid water (kg/m2) (-nlevsno+1:nlevgrnd) 
+    h2oice     => waterstate_inst%h2osoi_ice_col ! col ice lens (kg/m2) (-nlevsno+1:nlevgrnd)
+
+    snlsno     => col%snl                        ! number of snow layers (negative)
+
+    frac_sno   => waterstate_inst%frac_sno_eff_col ! fraction of ground covered by snow 
+    ! dz for snow layers is defined like in the history output as col%dz for the snow layers
+    dzsno(clm_begc:clm_endc, -nlevsno+1:0) = col%dz(clm_begc:clm_endc,-nlevsno+1:0)
+    ! Iterate through all columns
+    do jj=clm_begc,clm_endc
+      if (h2osno(jj).lt.0.0) then ! No snow in column
+        print *, "WARNING: negative snow in col: ", jj, h2osno
+!        ! Set existing layers to near zero and let CLM do the layer aggregation
+        do i=0,snlsno(jj)+1,-1
+            h2oliq(jj,i) = 0.0_r8
+            h2oice(jj,i) = 0.00000001_r8
+            dzsno(jj,i)  = 0.00000001_r8
+        end do
+        snow_depth(jj) = sum(dzsno(jj,:))
+        h2osno(jj) = sum(h2oice(jj,:))
+
+        if (isnan(h2osno(jj))) then
+          print *, "WARNING: h2osno at jj is nan: ", jj
+        endif
+        if (isnan(snow_depth(jj))) then
+          print *, "WARNING: snow_depth at jj is nan: ", jj
+        endif
+
+      else ! snow (h2osno) present
+        if (snlsno(jj).lt.0) then ! snow layers in the column
+          if (clmupdate_snow .eq. 1) then
+            ! DART source: https://github.com/NCAR/DART/blob/main/models/clm/dart_to_clm.f90
+            ! Formulas below from DART use h2osno_po / h2osno_pr for after / before DA SWE
+            ! clmupdate_snow == 1 has snow_depth after and h2osno before DA snow depth
+            ! Therefore need to have a transform to get h2osno_po
+            ! v1 init
+            ! h2osno_po(jj) = (snow_depth(jj) * bdsno) ! calculations from Init using constant SBD
+            ! v2 SoilTemperatureMod
+            if (snowdp(jj).gt.0.0_r8) then
+              rho_avg = min(800.0_r8, h2osno(jj)/snowdp(jj))
+            else
+              rho_avg = 200.0_r8
+            end if
+            if (frac_sno(jj).gt.0.0_r8 .and. snlsno(jj).lt.0.0_r8) then
+              h2osno_po(jj) = snow_depth(jj) * (rho_avg*frac_sno(jj))
+            else
+              h2osno_po(jj) = snow_depth(jj) * denice
+            end if
+            h2osno_pr(jj) = h2osno(jj)
+          else if (clmupdate_snow .eq. 2) then
+            ! for clmupdate_snow == 2 we have post H2OSNO as the main H2OSNO already
+            h2osno_po(jj) = h2osno(jj)
+            h2osno_pr(jj) = h2osno_in(jj)
+          end if
+
+          do ii=0,snlsno(jj)+1,-1 ! iterate through the snow layers
+            ! DART VERSION: ii = nlevsoi - i + 1
+            ! snow density prior for each layer
+            if (dzsno(jj,ii).gt.0.0_r8) then
+              snowden = (h2oliq(jj,ii) + h2oice(jj,ii) / dzsno(jj,ii))
+            else
+              snowden = 0.0_r8
+            endif
+            ! fraction of SWE in each active layers
+            if(h2osno_pr(jj).gt.0.0_r8) then
+              ! repartition Option 1: Adjust bottom layer only (set weight to 1 for bottom 0 else)
+              if(clmupdate_snow_repartitioning.eq.1) then
+                if (ii .eq. 0) then ! DART version indexing check against nlevsno but for us 0
+                  frac_swe = 1.0_r8
+                  ! JK: Let CLM repartitioning do the job
+                  ! afterwards. Provide all the snow in a single layer
+                else
+                  frac_swe = 0.0_r8
+                end if
+              ! repartition Option 2: Adjust all active layers
+              else if (clmupdate_snow_repartitioning.eq.2) then
+                frac_swe = (h2oliq(jj,ii) + h2oice(jj,ii)) / h2osno_pr(jj)
+              end if
+            else
+              frac_swe = 0.0_r8 ! no fraction SWE if no snow is present in column
+            end if ! end SWE fraction if
+
+            ! fraction of liquid and ice
+            if ((h2oliq(jj,ii) + h2oice(jj,ii)).gt.0.0_r8) then
+                frac_liq = h2oliq(jj,ii) / (h2oliq(jj,ii) + h2oice(jj,ii))
+                frac_ice = 1.0_r8 - frac_liq
+            else
+                frac_liq = 0.0_r8
+                frac_ice = 0.0_r8
+            end if
+
+            ! SWE adjustment per layer 
+            ! assumes identical layer distribution of liq and ice than before DA (frac_*)
+            if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then
+              gain_h2osno = (h2osno_po(jj) - h2osno_pr(jj)) * frac_swe
+              gain_h2oliq = gain_h2osno * frac_liq
+              gain_h2oice = gain_h2osno * frac_ice
+            else
+              gain_h2osno = 0.0_r8
+              gain_h2oliq = 0.0_r8
+              gain_h2oice = 0.0_r8
+            end if
+            ! layer level adjustments
+            if (snowden.gt.0.0_r8) then
+              gain_dzsno(ii) = gain_h2osno / snowden
+            else
+              gain_dzsno(ii) = 0.0_r8
+            end if
+            h2oliq(jj,ii) = h2oliq(jj,ii) + gain_h2oliq
+            h2oice(jj,ii) = h2oice(jj,ii) + gain_h2oice
+
+            ! Adjust snow layer dimensions so that CLM5 can calculate compaction / aggregation
+            ! in the DART code dzsno is adjusted directly but in CLM5 dzsno is local and diagnostic
+            ! i.e. calculated / assigned from frac_sno and dz(:, snow_layer) in SnowHydrologyMod
+            ! therefore we adjust dz(:, snow_layer) here
+            if (abs(h2osno_po(jj) - h2osno_pr(jj)).gt.0.0_r8) then
+              col%dz(jj, ii) = col%dz(jj, ii) + gain_dzsno(ii)
+              ! mid point and interface adjustments
+              ! i.e. zsno (col%z(:, snow_layers)) and zisno (col%zi(:, snow_layers))
+              ! DART version the sum goes from ilevel:nlevsno to fit with our indexing:
+              col%zi(jj, ii) = sum(col%dz(jj,ii:0))*-1.0_r8
+              ! In DART the check is ilevel == nlevsno but here
+              if (ii.eq.0) then
+                col%z(jj,ii) = col%zi(jj,ii) / 2.0_r8
+              else
+                col%z(jj,ii) = sum(col%zi(jj, ii:ii+1)) / 2.0_r8
+              end if
+            end if
+
+          end do ! End iteration of snow layers
+        endif ! End of snow layers present check
+
+        ! Column level variables
+        if (clmupdate_snow .eq. 1) then
+          ! Finally adjust SWE (h2osno) since the prior value is no longer needed
+          ! column level variables so we can adjust it outside the layer loop
+          h2osno(jj) = h2osno_po(jj) 
+        else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then
+          ! Update the total snow depth to match udpates to layers for active snow layers
+          if (abs(h2osno_po(jj) - h2osno_pr(jj)) .gt. 0.0_r8 .and. snlsno(jj) < 0.0_r8) then
+            snow_depth(jj) = snow_depth(jj) + sum(gain_dzsno(-nlevsno+1:0))
+          end if
+        end if
+
+        if (isnan(h2osno(jj))) then
+          print *, "WARNING2: h2osno at jj is nan: ", jj
+        endif
+        if (isnan(snow_depth(jj))) then
+          print *, "WARNING2: snow_depth at jj is nan: ", jj
+        endif
+
+      end if ! End of snow present check
+    end do ! End of column iteration
+
+  end subroutine clm_repartition_snow
+
   subroutine clm_correct_texture()
 
     use clm_varpar   , only : nlevsoi
diff --git a/bldsva/intf_DA/pdaf/model/wrapper_tsmp.c b/bldsva/intf_DA/pdaf/model/wrapper_tsmp.c
index 1af2bb6ae..d3e42b9de 100644
--- a/bldsva/intf_DA/pdaf/model/wrapper_tsmp.c
+++ b/bldsva/intf_DA/pdaf/model/wrapper_tsmp.c
@@ -198,7 +198,7 @@ void integrate_tsmp() {
 void update_tsmp(){
 
 #if defined CLMSA
-  if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0))){
+  if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_snow!=0))){
     update_clm(&tstartcycle, &mype_world);
     print_update_clm(&tcycle, &total_steps);
   }
diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md
index 4d1082a1a..c6c042f17 100644
--- a/doc/content/setup_tsmp/input_enkfpf.md
+++ b/doc/content/setup_tsmp/input_enkfpf.md
@@ -53,6 +53,8 @@ nprocs      =
 update_swc  =
 update_texture  =
 update_T  =
+update_snow  =
+update_snow_repartitioning  =
 print_swc   =
 print_et   =
 statevec_allcol =
@@ -478,6 +480,79 @@ Currently only CLM3.5
 
 -  1: Update of ground and vegetation temperature
 
+### CLM:update_snow ###
+
+`CLM:update_snow`: (integer) Flag for Snow-DA.
+
+Only CLM5.0/eCLM.
+
+-  0 (Default): No Snow-DA
+
+-  1: Assimilation of snow depth. State vector: Snow depth.
+
+-  2: Assimilation of snow water equivalent. State vector: Snow water
+   equivalent.
+
+-  3: Assimilation of snow depth. State vector: Snow depth and snow
+   water equivalent. Requires
+   `CLM:update_snow_repartitioning=3`. `h2osoi_ice` updated based on
+   `h2osno`-increment.
+
+-  4: Assimilation of snow depth. State vector: Snow depth and snow
+   water equivalent. Requires
+   `CLM:update_snow_repartitioning=3`. `h2osoi_ice` and `h2osoi_liq`
+   updated based on `h2osno`-increment. `dz` updated based on
+   `snow_depth`-increment.
+
+-  5: Assimilation of snow depth. State vector: Snow depth and snow
+   water equivalent. Requires
+   `CLM:update_snow_repartitioning=3`. `h2osoi_ice` updated based on
+   `h2osno`-increment. `dz` updated based on `snow_depth`-increment.
+
+-  6: Assimilation of snow depth. State vector: Snow depth and snow
+   water equivalent. Requires
+   `CLM:update_snow_repartitioning=3`. `h2osoi_ice` and `h2osoi_liq`
+   updated based on `h2osno`-increment. 
+
+-  7: Assimilation of snow depth. State vector: Snow depth and snow
+   water equivalent. Requires
+   `CLM:update_snow_repartitioning=3`. `h2osoi_ice` updated based on
+   `h2osno`-increment. Should reproduce Case 3.
+
+See CLM Technical Note for more information on snow variable:
+<https://escomp.github.io/ctsm-docs/versions/release-clm5.0/html/tech_note/Snow_Hydrology/CLM50_Tech_Note_Snow_Hydrology.html>
+- snow depth: `waterstate_inst%snow_depth_col`
+- snow water equivalent: `waterstate_inst%h2osno_col`
+
+### CLM:update_snow_repartitioning ###
+
+`CLM:update_snow_repartitioning`: (integer) Flag for snow-layer
+repartitioning applied during Snow-DA.
+
+Only used if Snow-DA is switched on by setting `CLM:update_snow`. Only
+CLM5.0/eCLM.
+
+-  0: No repartitioning. Not recommended if Snow-DA is used,
+   when state vector consists of diagnostic variables.
+
+-  1: (experimental, only for `CLM:update_snow=1,2`) Repartitioning
+   routine in DA-interface: `clm_repartition_snow`. Option 1:
+   Adjusting the bottom layer only.
+
+-  2: (experimental, only for `CLM:update_snow=1,2`) Repartitioning
+   routine in DA-interface: `clm_repartition_snow`. Option 2:
+   Adjusting all active layers.
+
+-  3 (Default, Currently recommended): `h2osoi_ice` (and possibly
+   `h2osoi_liq` and `dz`) updated by increment based on updated
+   `h2osno` or `snow_depth`. Further repartitioning left to
+   CLM-code. Based on
+   <https://github.com/NASA-LIS/LISF/blob/master/lis/surfacemodels/land/clm2/da_snow/clm2_setsnowvars.F90>
+
+Note: Using the repartitioning scheme has lead to balance errors and
+NaN-values in previous CLM-simulations. The function is still included
+in the code, because it may be revisited and adapted in the future.
+
 ### CLM:print_swc ###
 
 `CLM:print_swc`: (integer) If set to `1`, the updated soil moisture
@@ -752,58 +827,60 @@ Default: 0, output turned off.
 
 ## Parameter Summary ##
 
- | section   | parameter               | default value |
- |:---------:|:-----------------------:|:-------------:|
- | `[PF]`    |                         |               |
- |           | `problemname`           | \-            |
- |           | `nprocs`                | 0             |
- |           | `starttime`             | 0.0           |
- |           | `endtime`               | 0             |
- |           | `simtime`               | 0             |
- |           | `dt`                    | 0.0           |
- |           | `updateflag`            | 1             |
- |           | `paramupdate`           | 0             |
- |           | `paramupdate_frequency` | 1             |
- |           | `dampingfactor_param`   | 1.0           |
- |           | `dampingfactor_state`   | 1.0           |
- |           | `damping_switch_sm`     | 1             |
- |           | `aniso_perm_y`          | 1.0           |
- |           | `aniso_perm_z`          | 1.0           |
- |           | `aniso_use_parflow`     | 0             |
- |           |                         |               |
- |           | `printensemble`         | 1             |
- |           | `t_printensemble`       | -1            |
- |           | `printstat`             | 1             |
- |           | `paramprintensemble`    | 1             |
- |           | `paramprintstat`        | 1             |
- |           | `olfmasking`            | 0             |
- |           | `olfmasking_param`      | 0             |
- |           | `olfmasking_depth`      | 0             |
- | `[CLM]`   |                         |               |
- |           | `problemname`           | \-            |
- |           | `nprocs`                | 0             |
- |           | `update_swc`            | 1             |
- |           | `print_swc`             | 0             |
- |           | `print_et`              | 0             |
- |           | `statevec_allcol`       | 0             |
- |           | `statevec_only_active`  | 0             |
- |           | `statevec_max_layer`    | 25            |
- |           | `t_printensemble`       | -1            |
- |           | `watmin_switch`         | 0             |
- | `[COSMO]` |                         |               |
- |           | `nprocs`                | 0             |
- |           | `dtmult`                | 0             |
- | `[DA]`    |                         |               |
- |           | `nreal`                 | 0             |
- |           | `outdir`                | \-            |
- |           | `da_interval`           | 1             |
- |           | `stat_dumpoffset`       | 0             |
- |           | `screen_wrapper`        | 1             |
- |           | `point_obs`             | 1             |
- |           | `obs_interp_switch`     | 0             |
- |           | `crns_flag`             | 1             |
- |           | `da_crns_depth_tol`     | 0.01          |
- |           | `print_obs_index`       | 0             |
+ | section   | parameter                    | default value |
+ |:---------:|:----------------------------:|:-------------:|
+ | `[PF]`    |                              |               |
+ |           | `problemname`                | \-            |
+ |           | `nprocs`                     | 0             |
+ |           | `starttime`                  | 0.0           |
+ |           | `endtime`                    | 0             |
+ |           | `simtime`                    | 0             |
+ |           | `dt`                         | 0.0           |
+ |           | `updateflag`                 | 1             |
+ |           | `paramupdate`                | 0             |
+ |           | `paramupdate_frequency`      | 1             |
+ |           | `dampingfactor_param`        | 1.0           |
+ |           | `dampingfactor_state`        | 1.0           |
+ |           | `damping_switch_sm`          | 1             |
+ |           | `aniso_perm_y`               | 1.0           |
+ |           | `aniso_perm_z`               | 1.0           |
+ |           | `aniso_use_parflow`          | 0             |
+ |           |                              |               |
+ |           | `printensemble`              | 1             |
+ |           | `t_printensemble`            | -1            |
+ |           | `printstat`                  | 1             |
+ |           | `paramprintensemble`         | 1             |
+ |           | `paramprintstat`             | 1             |
+ |           | `olfmasking`                 | 0             |
+ |           | `olfmasking_param`           | 0             |
+ |           | `olfmasking_depth`           | 0             |
+ | `[CLM]`   |                              |               |
+ |           | `problemname`                | \-            |
+ |           | `nprocs`                     | 0             |
+ |           | `update_swc`                 | 1             |
+ |           | `update_snow`                | 0             |
+ |           | `update_snow_repartitioning` | 1             |
+ |           | `print_swc`                  | 0             |
+ |           | `print_et`                   | 0             |
+ |           | `statevec_allcol`            | 0             |
+ |           | `statevec_only_active`       | 0             |
+ |           | `statevec_max_layer`         | 25            |
+ |           | `t_printensemble`            | -1            |
+ |           | `watmin_switch`              | 0             |
+ | `[COSMO]` |                              |               |
+ |           | `nprocs`                     | 0             |
+ |           | `dtmult`                     | 0             |
+ | `[DA]`    |                              |               |
+ |           | `nreal`                      | 0             |
+ |           | `outdir`                     | \-            |
+ |           | `da_interval`                | 1             |
+ |           | `stat_dumpoffset`            | 0             |
+ |           | `screen_wrapper`             | 1             |
+ |           | `point_obs`                  | 1             |
+ |           | `obs_interp_switch`          | 0             |
+ |           | `crns_flag`                  | 1             |
+ |           | `da_crns_depth_tol`          | 0.01          |
+ |           | `print_obs_index`            | 0             |
 
 Default values for parameter file `enkfpf.par`.