From a955971fc38533174576c23ad1cdaaa2ebe3ff18 Mon Sep 17 00:00:00 2001 From: "Justin A. Pflug" Date: Fri, 9 Jun 2023 17:46:23 -0400 Subject: [PATCH 1/2] new code updates added --- .../da_snow/noahmp401_snow_update.F90 | 81 ++++++++++++++++--- 1 file changed, 68 insertions(+), 13 deletions(-) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 b/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 index 93e105177..8f9173232 100755 --- a/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 @@ -16,6 +16,7 @@ ! 14 Dec 2018: Yeosang Yoon; Modified code for NoahMP 4.0.1 and SNODEP ! 15 May 2019: Yeosang Yoon; Modified for NoahMP 4.0.1 and LDTSI ! 13 Dec 2019: Eric Kemp; Replaced LDTSI with SNOW +! 05 Jun 2023: Justin Pflug; fixes for SnowModel-defined snow updates ! ! !INTERFACE @@ -59,6 +60,7 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) integer :: snl_idx,i,j,iz integer :: iloc, jloc ! needed, but not use real :: smp,sneqv,snowh + real :: snoden real :: sneqv1,snowh1 real :: ponding1,ponding2 integer :: newnode @@ -137,6 +139,22 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) stc(1:nsoil) = & noahmp401_struc(n)%noahmp401(t)%tslb(1:nsoil) + ! NMP snow density calculation + if(snowh.gt.0) then + snoden = sneqv/(snowh*1000) + else + snoden = 0.55 + endif + + ! allow snow update even in cases where changes opp. directions + ! alter snow depth change to be in direction of SWE change + if((dsneqv.gt.0.and.dsnowh.le.0).or.& + (dsneqv.lt.0.and.dsnowh.ge.0)) then + dsnowh = (dsneqv/snoden)/1000 + ! set snow depth change to zero in instance where no SWE change + elseif(dsneqv.eq.0.and.dsnowh.ne.0) then + dsnowh = 0. + endif ! from snowfall routine ! creating a new layer @@ -147,12 +165,12 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) newnode = 0 - if(isnow == 0 .and. snowh >= 0.025.and.& - (dsneqv.gt.0.and.dsnowh.gt.0)) then !mb: change limit - isnow = -1 - newnode = 1 + if(isnow.eq.0.and.dsneqv.gt.0.and.dsnowh.gt.0) then + if(snowh.ge.0.025) then + isnow = -1 + newnode = 1 + endif dzsnso(0)= snowh - snowh = 0. stc(0) = min(273.16, noahmp401_struc(n)%noahmp401(t)%sfctmp) snice(0) = sneqv snliq(0) = 0. @@ -168,23 +186,60 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) if(dsneqv.lt.0.and.dsnowh.lt.0) then snowh1 = snowh + dsnowh sneqv1 = sneqv + dsneqv +! if dsnowh adjusted since dsneqv and dsnowh in opp. directions +! can cause one or other snowh1 or sneqv1 to be negative + if(sneqv1.gt.0.and.snowh1.le.0) then + snowh = ((sneqv1/snoden)/1000)-dsnowh + snowh1 = snowh + dsnowh +! if SWE disappears, also make sure snow depth disappears + elseif(sneqv.le.0) then + sneqv = -dsneqv + sneqv1 = sneqv + dsneqv + snowh = -dsnowh + snowh1 = snowh + dsnowh + endif +! make sure snow layers currently exist in decrease case + if(dzsnso(0).eq.0) then + if(snowh.ge.0.025) then + isnow = -1 + else + isnow = 0 + endif + dzsnso(0)= snowh + stc(0) = min(273.16, noahmp401_struc(n)%noahmp401(t)%sfctmp) + snice(0) = sneqv + snliq(0) = 0. + endif if(snowh1.ge.0.and.sneqv1.ge.0) then snowh = snowh + dsnowh sneqv = sneqv + dsneqv -! update dzsnso -! how do you determine the thickness of a layer? +! snow can no longer fill layer 1 if(snowh.le.dzsnso(0)) then isnow = 0 dzsnso(-nsnow+1:(isnow-1)) = 0 dzsnso(isnow) = snowh + snice(-nsnow+1:(isnow-1)) = 0 + snice(isnow) = sneqv + snliq(-nsnow+1:isnow) = 0 +! snow can no longer fill layer 1 and 2 elseif(snowh.le.(dzsnso(0)+dzsnso(-1))) then - isnow = -1 - dzsnso(-nsnow+1:(isnow-1)) = 0 - dzsnso(isnow) = snowh -dzsnso(isnow+1) - elseif(snowh.le.(dzsnso(0)+dzsnso(-1)+dzsnso(-2))) then isnow = -2 - dzsnso(-nsnow+1:(isnow-2)) = 0 - dzsnso(isnow) = snowh -dzsnso(isnow+2) + dzsnso(-nsnow+1:isnow) = 0 + dzsnso(isnow+1) = snowh -dzsnso(0) + ! scale swe in layers by ratio of depth to pack + do snl_idx=-nsnow+1,0 + snice(snl_idx) = sneqv*(dzsnso(snl_idx)/snowh) + enddo + snliq(-nsnow+1:isnow) = 0 +! all other cases + elseif(snowh.le.(dzsnso(0)+dzsnso(-1)+dzsnso(-2))) then + isnow = -3 + dzsnso(isnow+1) = snowh -dzsnso(-1) -dzsnso(0) + ! scale swe in layers by ratio of depth to pack + do snl_idx=-nsnow+1,0 + snice(snl_idx) = sneqv*(dzsnso(snl_idx)/snowh) + enddo + snliq(-nsnow+1:isnow) = 0 endif endif endif From e27c778413c639927712d2d8d4d33ec1f94078e7 Mon Sep 17 00:00:00 2001 From: David Mocko Date: Thu, 7 Sep 2023 09:14:52 -0400 Subject: [PATCH 2/2] Fix ISRIC soil texture water point index value in LDT This pull request fixes the water point index values when using the ISRIC soil texture maps in LDT. LDT re-maps the ISRIC soil texture map onto the STATSGO soil classification index. However, the water points were incorrectly assigned a value of 13, when it should be 14 in the STATSGO soil texture classification. Resolves: #1406 --- ldt/core/LDT_paramOptCheckMod.F90 | 4 ++-- ldt/core/LDT_soilsMod.F90 | 8 ++++---- ldt/params/soils/read_ISRIC_texture.F90 | 2 +- ldt/params/soils/set_ISRIC_texture_attribs.F90 | 2 +- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/ldt/core/LDT_paramOptCheckMod.F90 b/ldt/core/LDT_paramOptCheckMod.F90 index 41e17df70..8b705f607 100644 --- a/ldt/core/LDT_paramOptCheckMod.F90 +++ b/ldt/core/LDT_paramOptCheckMod.F90 @@ -483,8 +483,8 @@ subroutine LDT_soilsOptChecks( n, param_name, soil_class, spatial_transform ) endif case( "ISRIC") - if (nt.ne.13) then - write(LDT_logunit,*) "Param_Check: ISRIC has 13 types (includes water) for tiling." + if (nt.ne.16) then + write(LDT_logunit,*) "Param_Check: ISRIC has 16 types (includes water) for tiling." write(LDT_logunit,*) " The current value is: ",nt write(LDT_logunit,*) " Stopping ..." call LDT_endrun diff --git a/ldt/core/LDT_soilsMod.F90 b/ldt/core/LDT_soilsMod.F90 index cff4a36b2..7eb889b1f 100644 --- a/ldt/core/LDT_soilsMod.F90 +++ b/ldt/core/LDT_soilsMod.F90 @@ -782,7 +782,7 @@ subroutine soils_init_LIS() elseif(INDEX(LDT_LSMparam_struc(n)%texture%source,"ZOBLER") >0) then soiltext%watervalue = 1.0 elseif(INDEX(LDT_LSMparam_struc(n)%texture%source,"ISRIC") >0) then - soiltext%watervalue = 13. + soiltext%watervalue = 14. else soiltext%watervalue = LDT_rc%udef endif @@ -1484,7 +1484,7 @@ subroutine soils_init_LISHydro(flag) elseif(INDEX(LDT_LSMparam_struc(n)%texture%source,"ZOBLER") >0) then soiltext%watervalue = 1.0 elseif(INDEX(LDT_LSMparam_struc(n)%texture%source,"ISRIC") >0) then - soiltext%watervalue = 13. + soiltext%watervalue = 14. else soiltext%watervalue = LDT_rc%udef endif @@ -1814,7 +1814,7 @@ subroutine soils_writeHeader_LIS(n,ftn,dimID) 14)) elseif( LDT_rc%soil_classification(1) == "ISRIC" ) then call LDT_verify(nf90_put_att(ftn,NF90_GLOBAL,"NUMBER_SOILTYPES", & - 13)) + 19)) else call LDT_verify(nf90_put_att(ftn,NF90_GLOBAL,"NUMBER_SOILTYPES", & 19)) @@ -1972,7 +1972,7 @@ subroutine soils_writeHeader_LISHydro(n,ftn,dimID,flag) 14)) elseif( LDT_rc%soil_classification(1) == "ISRIC" ) then call LDT_verify(nf90_put_att(ftn,NF90_GLOBAL,"NUMBER_SOILTYPES", & - 13)) + 19)) else call LDT_verify(nf90_put_att(ftn,NF90_GLOBAL,"NUMBER_SOILTYPES", & 19)) diff --git a/ldt/params/soils/read_ISRIC_texture.F90 b/ldt/params/soils/read_ISRIC_texture.F90 index 1b6bf4d0d..21e17b627 100644 --- a/ldt/params/soils/read_ISRIC_texture.F90 +++ b/ldt/params/soils/read_ISRIC_texture.F90 @@ -94,7 +94,7 @@ subroutine read_ISRIC_texture( n, num_bins, fgrd, texture_layers ) ! ___________________________________________________________________ - water_class = 13 ! Water class for ISRIC soil texture class + water_class = 14 ! Water class for ISRIC soil texture class temp = LDT_rc%udef gridcnt = 0. diff --git a/ldt/params/soils/set_ISRIC_texture_attribs.F90 b/ldt/params/soils/set_ISRIC_texture_attribs.F90 index a50614c00..470996bc6 100644 --- a/ldt/params/soils/set_ISRIC_texture_attribs.F90 +++ b/ldt/params/soils/set_ISRIC_texture_attribs.F90 @@ -34,7 +34,7 @@ subroutine set_ISRIC_texture_attribs() ! \end{description} !EOP ! - LDT_LSMparam_struc(:)%texture%num_bins = 13 + LDT_LSMparam_struc(:)%texture%num_bins = 16 LDT_LSMparam_struc(:)%texture%vlevels = 1 LDT_soils_struc(:)%texture_nlyrs%num_bins = 1