Defaulting soil texture if clay+sand+silt!=100%

NERC-CEH · Jan 8, 2025 · 1fa61a2 · 1fa61a2
1 parent a999126
commit 1fa61a2
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 8 deletions.
diff --git a/src/Data/DataInputModule.f90 b/src/Data/DataInputModule.f90
@@ -445,7 +445,6 @@ subroutine readBatchVariablesDatabase(me)
         ! Soil texture                              [%]
         var = me%nc%getVariable('soil_texture_clay_content')
         call var%getData(me%soilTextureClayContent)
-        print *, me%soilTextureClayContent
         var = me%nc%getVariable('soil_texture_sand_content')
         call var%getData(me%soilTextureSandContent)
         var = me%nc%getVariable('soil_texture_silt_content')

diff --git a/src/Soil/SoilProfileModule.f90 b/src/Soil/SoilProfileModule.f90
@@ -350,7 +350,6 @@ function calculateSizeDistributionSoilProfile(me, clay, silt, sand, enrichClay)
         real    :: ssd_(3,C%nSizeClassesSpm)                        ! Temporary SSD array, before summing across SSD dimension
         ! Bins for texture content, based on definition of clay, silt and sand. First bins
         ! have non-zero lower bound to avoid numerical errors when logging
-        print *, clay, silt, sand
         texture = [clay, silt, sand] / 100.0
         if (enrichClay) then
             clayEnrichmentRatio = 0.26 + 1 / (1 - texture(3))               ! Ref: Stefano and Ferro, 2002: https://doi.org/10.1006/bioe.2001.0034
@@ -434,6 +433,13 @@ function parseInputDataSoilProfile(me) result(r)
         me%sandContent = DATASET%soilTextureSandContent(me%x, me%y)
         me%siltContent = DATASET%soilTextureSiltContent(me%x, me%y)
         me%coarseFragContent = DATASET%soilTextureCoarseFragContent(me%x, me%y)
+        ! Check if clay, sand and silt sum to (nearly) 1, and if not, default to
+        ! the average soil texture for Europe
+        if (abs(100.0 - me%clayContent - me%sandContent - me%siltContent) > 0.1) then
+            me%clayContent = 18.0
+            me%sandContent = 46.0
+            me%siltContent = 36.0
+        end if
         if (me%coarseFragContent == nf90_fill_real) then
             me%coarseFragContent = 0.0
         end if
@@ -550,17 +556,15 @@ subroutine parseNewBatchDataSoilProfile(me)
         me%sandContent = DATASET%soilTextureSandContent(me%x, me%y)
         me%siltContent = DATASET%soilTextureSiltContent(me%x, me%y)
         me%coarseFragContent = DATASET%soilTextureCoarseFragContent(me%x, me%y)
-        ! If one of the clay/sand/silt content variables is empty, then use defaults
-        ! instead. The defaults are roughly the average soil texture across Europe
-        ! TODO maybe move this to the DataInputModule as part of data validation
-        if ((me%clayContent == nf90_fill_real) .or. (me%sandContent == nf90_fill_real) &
-            .or. (me%siltContent == nf90_fill_real)) then
+        ! Check if clay, sand and silt sum to (nearly) 1, and if not, default to
+        ! the average soil texture for Europe
+        if (abs(1.0 - me%clayContent - me%sandContent - me%siltContent) > 1e-3) then
             me%clayContent = 0.18
             me%sandContent = 0.46
             me%siltContent = 0.36
         end if
         if (me%coarseFragContent == nf90_fill_real) then
-            me%coarseFragContent = 0.0_dp
+            me%coarseFragContent = 0.0
         end if
         ! Calculate the average grain diameter from soil texture
         me%d_grain = me%calculateAverageGrainSize(me%clayContent, me%siltContent, me%sandContent)