From 7218d3deba2a8335ff66ff5b789f06e85cd24432 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 21 Mar 2025 11:12:19 +0100 Subject: [PATCH 1/4] Snow-DA: introduce `clmupdate_snow.eq.4` and refactor - removed `rsnow`/`nsnow` --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 142 ++++++++++++++++------ interface/model/eclm/enkf_clm_mod_5.F90 | 140 +++++++++++++++------ 2 files changed, 205 insertions(+), 77 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index a6eff6c69..dda3ec8e4 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -282,11 +282,18 @@ subroutine define_clm_statevec(mype) 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 + ! 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 !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 @@ -484,9 +494,12 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - real(r8) :: rliq,rice,incr_sno - real(r8) :: rsnow(clm_begc:clm_endc) - real(r8) :: nsnow(clm_begc:clm_endc) + 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) :: 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 @@ -721,32 +734,31 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! CLM5-grid-index cc = (col%gridcell(j) - clm_begg + 1) + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + if (clmupdate_snow.eq.1) then - rsnow(j) = snow_depth(j) + snow_depth(j) = clm_statevec(cc+offset) else if (clmupdate_snow.eq.2) then - rsnow(j) = h2osno(j) + h2osno(j) = clm_statevec(cc+offset) else if (clmupdate_snow.eq.3) then - rsnow(j) = h2osno(j) + h2osno(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.4) then + snow_depth(j) = clm_statevec(cc+offset) + h2osno(j) = clm_statevec(cc+clm_varsize+offset) end if - if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then - nsnow(j) = clm_statevec(cc+offset) - else if (clmupdate_snow.eq.3) then - nsnow(j) = clm_statevec(cc+clm_varsize+offset) + if (snow_depth(j).lt.0.0) then + ! Catch negative values from DA + print *, "WARNING: Snow-depth is negative at state vector cc. cc, j, offset, snow_depth(j): ", cc, j, offset, snow_depth(j) + print *, "WARNING: Snow-depth is reset to pre-DA value." + snow_depth(j) = snow_depth_in(j) end if - - if (nsnow(j).gt.0.0) then - ! Update state variable to CLM - ! Not needed for repartioning-case 3? - if(clmupdate_snow.eq.1) then - snow_depth(j) = nsnow(j) - end if - if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then - h2osno(j) = nsnow(j) - end if - else - ! Catch negative or 0 values from DA - print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) + if (h2osno(j).lt.0.0) then + ! Catch negative values from DA + print *, "WARNING: SWE is negative at state vector cc. cc, j, offset, h2osno(j): ", cc, j, offset, h2osno(j) + print *, "WARNING: SWE is reset to pre-DA value or zero." + h2osno(j) = h2osno_in(j) end if end do @@ -754,9 +766,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Repartitioning if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then - if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then - if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + if (clmupdate_snow.eq.1) then + if ( SUM(ABS(snow_depth_in(:) - snow_depth(:))).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(:))).gt.0.000001) then + call clm_repartition_snow(h2osno_in(:)) end if end if @@ -768,19 +786,65 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if ( clmupdate_snow_repartitioning.eq.3) then - do j=clm_begc,clm_endc - if (nsnow(j).gt.0.0) then - if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then - if ( ABS(rsnow(j)).gt.0.0) then - ! Update h2osoi_ice with increment - incr_sno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno - end do + if (clmupdate_snow.eq.1) then + do j=clm_begc,clm_endc + if (snow_depth(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth(j)).gt.0.000001) then + if ( ABS(snow_depth_in(j)).gt.0.0) then + ! Update h2osoi_ice with increment + incr_sno = snow_depth(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 if - end if - end do + end do + end if + + if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + do j=clm_begc,clm_endc + if (h2osno(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno(j)).gt.0.000001) then + if ( ABS(h2sno_in(j)).gt.0.0) then + ! Update h2osoi_ice with increment + incr_sno = h2osno(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 if + end do + end if + + if (clmupdate_snow.eq.4) then + do j=clm_begc,clm_endc + if (h2osno(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno(j)).gt.0.000001) then + if ( ABS(h2sno_in(j)).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno(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 + end do + end if + end if + end if + if (snow_depth(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth(j)).gt.0.000001) then + if ( ABS(snow_depth_in(j)).gt.0.0) then + ! Update snow_depth with increment + incr_sd = snow_depth(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + dz(j,i) = dz(j,i) * incr_sd + end do + end if + end if + end if + end do + end if end if @@ -879,7 +943,7 @@ subroutine clm_repartition_snow(h2osno_in) h2osno_po(jj) = snow_depth(jj) * denice end if h2osno_pr(jj) = h2osno(jj) - else if (clmupdate_snow .eq. 2 .or. clmupdate_snow .eq. 3) then + 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) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index f8c3c94e6..dda3ec8e4 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -282,11 +282,18 @@ subroutine define_clm_statevec(mype) 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 + ! 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 !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 @@ -484,9 +494,12 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: h2osoi_ice(:,:) - real(r8) :: rliq,rice,incr_sno - real(r8) :: rsnow(clm_begc:clm_endc) - real(r8) :: nsnow(clm_begc:clm_endc) + 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) :: 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 @@ -721,32 +734,31 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! CLM5-grid-index cc = (col%gridcell(j) - clm_begg + 1) + snow_depth_in(j) = snow_depth(j) + h2osno_in(j) = h2osno(j) + if (clmupdate_snow.eq.1) then - rsnow(j) = snow_depth(j) + snow_depth(j) = clm_statevec(cc+offset) else if (clmupdate_snow.eq.2) then - rsnow(j) = h2osno(j) + h2osno(j) = clm_statevec(cc+offset) else if (clmupdate_snow.eq.3) then - rsnow(j) = h2osno(j) + h2osno(j) = clm_statevec(cc+clm_varsize+offset) + else if (clmupdate_snow.eq.4) then + snow_depth(j) = clm_statevec(cc+offset) + h2osno(j) = clm_statevec(cc+clm_varsize+offset) end if - if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then - nsnow(j) = clm_statevec(cc+offset) - else if (clmupdate_snow.eq.3) then - nsnow(j) = clm_statevec(cc+clm_varsize+offset) + if (snow_depth(j).lt.0.0) then + ! Catch negative values from DA + print *, "WARNING: Snow-depth is negative at state vector cc. cc, j, offset, snow_depth(j): ", cc, j, offset, snow_depth(j) + print *, "WARNING: Snow-depth is reset to pre-DA value." + snow_depth(j) = snow_depth_in(j) end if - - if (nsnow(j).gt.0.0) then - ! Update state variable to CLM - ! Not needed for repartioning-case 3? - if(clmupdate_snow.eq.1) then - snow_depth(j) = nsnow(j) - end if - if(clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then - h2osno(j) = nsnow(j) - end if - else - ! Catch negative or 0 values from DA - print *, "WARNING: Snow-statevec is negative/zero at cc. cc, j, offset, clm_statevec(cc+offset): ", cc, j, offset, clm_statevec(cc+offset) + if (h2osno(j).lt.0.0) then + ! Catch negative values from DA + print *, "WARNING: SWE is negative at state vector cc. cc, j, offset, h2osno(j): ", cc, j, offset, h2osno(j) + print *, "WARNING: SWE is reset to pre-DA value or zero." + h2osno(j) = h2osno_in(j) end if end do @@ -754,9 +766,15 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") ! Repartitioning if ( clmupdate_snow_repartitioning.eq.1 .or. clmupdate_snow_repartitioning.eq.2) then - if (clmupdate_snow.eq.1 .or. clmupdate_snow.eq.2) then - if ( SUM(ABS(rsnow(:) - nsnow(:))).gt.0.000001) then - call clm_repartition_snow(rsnow(:)) + if (clmupdate_snow.eq.1) then + if ( SUM(ABS(snow_depth_in(:) - snow_depth(:))).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(:))).gt.0.000001) then + call clm_repartition_snow(h2osno_in(:)) end if end if @@ -768,19 +786,65 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") if ( clmupdate_snow_repartitioning.eq.3) then - do j=clm_begc,clm_endc - if (nsnow(j).gt.0.0) then - if ( ABS(rsnow(j) - nsnow(j)).gt.0.000001) then - if ( ABS(rsnow(j)).gt.0.0) then - ! Update h2osoi_ice with increment - incr_sno = nsnow(j) / rsnow(j) ! INC = New snow var / OLD snow var - do i=snlsno(j)+1,0 - h2osoi_ice(j,i) = h2osoi_ice(j,i) * incr_sno - end do + if (clmupdate_snow.eq.1) then + do j=clm_begc,clm_endc + if (snow_depth(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth(j)).gt.0.000001) then + if ( ABS(snow_depth_in(j)).gt.0.0) then + ! Update h2osoi_ice with increment + incr_sno = snow_depth(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 if - end if - end do + end do + end if + + if (clmupdate_snow.eq.2 .or. clmupdate_snow.eq.3) then + do j=clm_begc,clm_endc + if (h2osno(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno(j)).gt.0.000001) then + if ( ABS(h2sno_in(j)).gt.0.0) then + ! Update h2osoi_ice with increment + incr_sno = h2osno(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 if + end do + end if + + if (clmupdate_snow.eq.4) then + do j=clm_begc,clm_endc + if (h2osno(j).gt.0.0) then + if ( ABS(h2osno_in(j) - h2osno(j)).gt.0.000001) then + if ( ABS(h2sno_in(j)).gt.0.0) then + ! Update h2osoi_ice/h2osoi_liq with increment + incr_swe = h2osno(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 + end do + end if + end if + end if + if (snow_depth(j).gt.0.0) then + if ( ABS(snow_depth_in(j) - snow_depth(j)).gt.0.000001) then + if ( ABS(snow_depth_in(j)).gt.0.0) then + ! Update snow_depth with increment + incr_sd = snow_depth(j) / snow_depth_in(j) + do i=snlsno(j)+1,0 + dz(j,i) = dz(j,i) * incr_sd + end do + end if + end if + end if + end do + end if end if From 75504d55a5a1e299c6af8976d60bf9d97588db65 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Fri, 21 Mar 2025 11:14:50 +0100 Subject: [PATCH 2/4] Snow-DA: update CLM3.5 --- interface/model/clm3_5/enkf_clm_mod.F90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 7918ce69d..c056f183f 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/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 From 55311a8bace43bc6691e69f84570e171e77b5d78 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Tue, 25 Mar 2025 10:00:14 +0100 Subject: [PATCH 3/4] Snow-DA: Typo-fix --- interface/model/clm5_0/enkf_clm_mod_5.F90 | 4 ++-- interface/model/eclm/enkf_clm_mod_5.F90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/interface/model/clm5_0/enkf_clm_mod_5.F90 b/interface/model/clm5_0/enkf_clm_mod_5.F90 index dda3ec8e4..e3423b5ad 100755 --- a/interface/model/clm5_0/enkf_clm_mod_5.F90 +++ b/interface/model/clm5_0/enkf_clm_mod_5.F90 @@ -806,7 +806,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do j=clm_begc,clm_endc if (h2osno(j).gt.0.0) then if ( ABS(h2osno_in(j) - h2osno(j)).gt.0.000001) then - if ( ABS(h2sno_in(j)).gt.0.0) then + if ( ABS(h2osno_in(j)).gt.0.0) then ! Update h2osoi_ice with increment incr_sno = h2osno(j) / h2osno_in(j) do i=snlsno(j)+1,0 @@ -822,7 +822,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do j=clm_begc,clm_endc if (h2osno(j).gt.0.0) then if ( ABS(h2osno_in(j) - h2osno(j)).gt.0.000001) then - if ( ABS(h2sno_in(j)).gt.0.0) then + if ( ABS(h2osno_in(j)).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno(j) / h2osno_in(j) do i=snlsno(j)+1,0 diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index dda3ec8e4..e3423b5ad 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -806,7 +806,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do j=clm_begc,clm_endc if (h2osno(j).gt.0.0) then if ( ABS(h2osno_in(j) - h2osno(j)).gt.0.000001) then - if ( ABS(h2sno_in(j)).gt.0.0) then + if ( ABS(h2osno_in(j)).gt.0.0) then ! Update h2osoi_ice with increment incr_sno = h2osno(j) / h2osno_in(j) do i=snlsno(j)+1,0 @@ -822,7 +822,7 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") do j=clm_begc,clm_endc if (h2osno(j).gt.0.0) then if ( ABS(h2osno_in(j) - h2osno(j)).gt.0.000001) then - if ( ABS(h2sno_in(j)).gt.0.0) then + if ( ABS(h2osno_in(j)).gt.0.0) then ! Update h2osoi_ice/h2osoi_liq with increment incr_swe = h2osno(j) / h2osno_in(j) do i=snlsno(j)+1,0 From 74dd63ba0891abbf83147ce425f1e58952f30925 Mon Sep 17 00:00:00 2001 From: Buliao1 <152620321+Buliao1@users.noreply.github.com> Date: Wed, 26 Mar 2025 11:10:50 +0100 Subject: [PATCH 4/4] remove the updating the h2osoi_liq in Update enkf_clm_mod_5.F90 --- interface/model/eclm/enkf_clm_mod_5.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index e3423b5ad..7f869a989 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -827,7 +827,6 @@ subroutine update_clm(tstartcycle, mype) bind(C,name="update_clm") incr_swe = h2osno(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 end do end if end if