Snow-DA: Introduce CLM:update_snow.eq.4

jjokella · jjokella · commit c119c812cbba · 2025-03-27T16:19:26.000+01:00
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
@@ -157,10 +157,16 @@ subroutine define_clm_statevec(mype)
       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
+    ! 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
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
@@ -287,6 +287,13 @@ subroutine define_clm_statevec(mype)
         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
 
     !hcp LST DA
     if(clmupdate_T.eq.1) then
@@ -430,6 +437,9 @@ subroutine set_clm_statevec(tstartcycle, mype)
               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
                 error stop "Wrong input for clmupdate_snow"
               end if
@@ -732,6 +742,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
           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)
         end if
 
         if (clmupdate_snow.eq.1) then
@@ -740,9 +753,12 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
           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)
         end if
 
-        if(clmupdate_snow.eq.1) then
+        if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4) then
           ! Update state variable to CLM
           ! Not needed for repartioning-case 3?
           if (snow_depth_out(j).gt.0.0) then
@@ -751,7 +767,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
             ! 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 .or. clmupdate_snow.eq.3) then
+        else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then
           if (h2osno_out(j).gt.0.0) then
             h2osno(j) = h2osno_out(j)
           else
@@ -817,6 +833,43 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
           end do
         end if
 
+        if (clmupdate_snow.eq.4) then
+          do j=clm_begc,clm_endc
+            if (h2osno_out(j).gt.0.0) then
+              if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then
+                if (h2osno_in(j).gt.0.0) then
+                  ! Update h2osoi_ice/h2osoi_liq with increment
+                  incr_swe = h2osno_out(j) / h2osno_in(j)
+                  do i=snlsno(j)+1,0
+                    h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe
+                    h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe
+                    if (isnan(h2osoi_ice(j,i))) then
+                      print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i
+                    endif
+                    if (isnan(h2osoi_liq(j,i))) then
+                      print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i
+                    endif
+                  end do
+                end if
+              end if
+            end if
+            if (snow_depth_out(j).gt.0.0) then
+              if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then
+                if (snow_depth_in(j).gt.0.0) then
+                  ! Update snow_depth with increment
+                  incr_sd = snow_depth_out(j) / snow_depth_in(j)
+                  do i=snlsno(j)+1,0
+                    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 do
+                end if
+              end if
+            end if
+          end do
+        end if
+
       end if
 
 #ifdef PDAF_DEBUG
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
@@ -287,6 +287,13 @@ subroutine define_clm_statevec(mype)
         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
 
     !hcp LST DA
     if(clmupdate_T.eq.1) then
@@ -430,6 +437,9 @@ subroutine set_clm_statevec(tstartcycle, mype)
               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
                 error stop "Wrong input for clmupdate_snow"
               end if
@@ -732,6 +742,9 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
           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)
         end if
 
         if (clmupdate_snow.eq.1) then
@@ -740,9 +753,12 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
           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)
         end if
 
-        if(clmupdate_snow.eq.1) then
+        if(clmupdate_snow.eq.1 .or. clmupdate_snow.eq.4) then
           ! Update state variable to CLM
           ! Not needed for repartioning-case 3?
           if (snow_depth_out(j).gt.0.0) then
@@ -751,7 +767,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
             ! 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 .or. clmupdate_snow.eq.3) then
+        else if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3 .or. clmupdate_snow.eq.4) then
           if (h2osno_out(j).gt.0.0) then
             h2osno(j) = h2osno_out(j)
           else
@@ -817,6 +833,43 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm")
           end do
         end if
 
+        if (clmupdate_snow.eq.4) then
+          do j=clm_begc,clm_endc
+            if (h2osno_out(j).gt.0.0) then
+              if ( ABS(h2osno_in(j) - h2osno_out(j)).gt.0.000001) then
+                if (h2osno_in(j).gt.0.0) then
+                  ! Update h2osoi_ice/h2osoi_liq with increment
+                  incr_swe = h2osno_out(j) / h2osno_in(j)
+                  do i=snlsno(j)+1,0
+                    h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_swe
+                    h2osoi_liq(j,i) = h2osoi_liq(j,i) * incr_swe
+                    if (isnan(h2osoi_ice(j,i))) then
+                      print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i
+                    endif
+                    if (isnan(h2osoi_liq(j,i))) then
+                      print *, "WARNING: h2osoi_ice at j,i is nan: ", j, i
+                    endif
+                  end do
+                end if
+              end if
+            end if
+            if (snow_depth_out(j).gt.0.0) then
+              if ( ABS(snow_depth_in(j) - snow_depth_out(j)).gt.0.000001) then
+                if (snow_depth_in(j).gt.0.0) then
+                  ! Update snow_depth with increment
+                  incr_sd = snow_depth_out(j) / snow_depth_in(j)
+                  do i=snlsno(j)+1,0
+                    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 do
+                end if
+              end if
+            end if
+          end do
+        end if
+
       end if
 
 #ifdef PDAF_DEBUG
diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md
@@ -494,8 +494,15 @@ Only CLM5.0/eCLM.
    equivalent.
 
 -  3: Assimilation of snow depth. State vector: Snow depth and snow
-   water equivalent. Snow water equivalent from state vector is used
-   in the update as in Case 2.
+   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.
 
 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>
@@ -521,8 +528,9 @@ CLM5.0/eCLM.
    routine in DA-interface: `clm_repartition_snow`. Option 2:
    Adjusting all active layers.
 
--  3 (Default, Currently recommended): `h2osoi_ice` updated by increment based
-   on updated state vector variable. Further repartitioning left to
+-  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>
 
