Merge branch 'support/lisf-public-7.4'

ldt/core/LDT_paramOptCheckMod.F90

@@ -499,8 +499,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

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))

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.

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

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