add bounds checking on bendresist and vegshape parameters

lanl · Aug 5, 2024 · b10afa9 · b10afa9
1 parent 489ae8a
commit b10afa9
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 9 deletions.
diff --git a/components/elm/src/biogeochem/VegStructUpdateMod.F90 b/components/elm/src/biogeochem/VegStructUpdateMod.F90
@@ -116,9 +116,6 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, &
 
       dt = real( get_rad_step_size(), r8 )
 
-      ! convert from stems/ha -> stems/m^2
-      stocking = stocking / 10000._r8
-
       ! patch loop
       do fp = 1,num_soilp
          p = filter_soilp(fp)
@@ -164,11 +161,11 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, &
                ! taper and stocking density can be set as input variables now to
                ! change from default values set in pftvarcon.F90
                if (spinup_state >= 1) then
-                 htop(p) = ((3._r8 * deadstemc(p) * spinup_mortality_factor * taper(p) * taper(p))/ &
-                      (SHR_CONST_PI * stocking(p) * dwood(ivt(p))))**(1._r8/3._r8)
+                 htop(p) = ((3._r8 * deadstemc(p) * spinup_mortality_factor * taper(ivt(p)) * taper(ivt(p)))/ &
+                      (SHR_CONST_PI * stocking(ivt(p)) * dwood(ivt(p))))**(1._r8/3._r8)
                else
-                 htop(p) = ((3._r8 * deadstemc(p) * taper(p) * taper(p))/ &
-                      (SHR_CONST_PI * stocking(p) * dwood(ivt(p))))**(1._r8/3._r8)
+                 htop(p) = ((3._r8 * deadstemc(p) * taper(ivt(p)) * taper(ivt(p)))/ &
+                      (SHR_CONST_PI * stocking(ivt(p)) * dwood(ivt(p))))**(1._r8/3._r8)
                end if
 
                ! Peter Thornton, 5/3/2004
@@ -240,7 +237,7 @@ subroutine VegStructUpdate(num_soilp, filter_soilp, &
          ! Liston and Hiemstra, 2011; Belke-Brea et al. 2020
          if (ivt(p) > noveg .and. ivt(p) <= nbrdlf_dcd_brl_shrub ) then
             ol = min( max(snow_depth(c)-hbot(p), 0._r8), htop(p)-hbot(p))
-            fb = 1._r8 - (ol / max(1.e-06_r8, bendresist(p) * (htop(p)-hbot(p)))) ** vegshape(p)
+            fb = 1._r8 - (ol / max(1.e-06_r8, bendresist(ivt(p)) * (htop(p)-hbot(p)))) ** vegshape(ivt(p))
          else
             fb = 1._r8 - max(min(snow_depth(c),0.2_r8),0._r8)/0.2_r8   ! 0.2m is assumed
             !depth of snow required for complete burial of grasses

diff --git a/components/elm/src/main/pftvarcon.F90 b/components/elm/src/main/pftvarcon.F90
@@ -1036,12 +1036,25 @@ subroutine pftconrd
     if (.not. readv ) bendresist(:) = 1._r8
     call ncd_io('vegshape', vegshape, 'read', ncid, readvar=readv, posNOTonfile=.true.)
     if (.not. readv ) vegshape(:) = 1._r8
+   ! check validity
+    do i = 0, npft -1
+      if (bendresist(i) .gt. 1.0_r8 .or. bendresist(i) .le. 0._r8) then
+         call endrun(msg="Non-physical selection of bendresist parameter, set between 0 and 1"//errMsg(__FILE__, __LINE__))
+      end if
+      if (vegshape(i) .gt. 2.0_r8 .or. vegshape(i) .le. 0_r8) then
+         call endrun(msg="Non-physical selection of vegshape parameter, set between 0 and 2"//errMsg(__FILE__, __LINE__))
+      end if
+   end do
     call ncd_io('stocking', stocking, 'read', ncid, readvar=readv, posNOTonfile=.true.)
     if (.not. readv ) stocking(:) = 1000._r8
+    ! convert from stems/ha -> stems/m2
+    do i = 0, npft
+      stocking(i) = stocking(i) / 10000._r8
+    end do 
     call ncd_io('taper', taper, 'read', ncid, readvar=readv, posNOTonfile=.true.)
     if (.not. readv ) then
       taper(:) = 200._r8
-      do i = 0, npft - 1
+      do i = 0, npft
         ! RPF 240331 - need to revisit this section when integrating with IM4,
         ! to make consistent with attempt to remove hard-coded pft numbers.
         if (i >= nbrdlf_evr_shrub .and. i <= nbrdlf_dcd_brl_shrub) taper(i) = 10.0_r8 ! shrubs