Skip to content

BCS processing to support HWSDv2 soil preprocessing and new Peat GPA22 #1196

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
biljanaorescanin wants to merge 15 commits into
base: develop
Choose a base branch
from
Draft

BCS processing to support HWSDv2 soil preprocessing and new Peat GPA22 #1196

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
56f7eca
Add HWSDv2 soil preprocessing STEP1–STEP5 pipeline and QA tools
biljanaorescanin Jan 7, 2026
e8cdebc
Update raster BCS processing to support HWSDv2 soil preprocessing
biljanaorescanin Jan 7, 2026
367c00c
Preserve backward compatibility for PEATMAP handling
biljanaorescanin Jan 8, 2026
be10c12
missed commit declaration
biljanaorescanin Jan 8, 2026
d07f77c
revert soil properties to HWSDv1.2
biljanaorescanin Jan 27, 2026
64af96a
edit text for v15
biljanaorescanin Jan 27, 2026
257daf9
edited comments in subroutine soil_para_hwsd() of mod_process_hres_da…
gmao-rreichle Mar 11, 2026
29501cb
Merge branch 'develop' into feature/borescan_peat_and_hwsd
biljanaorescanin Mar 15, 2026
6f40ecf
Refine strict GPA22 peat classification and donor repair handling
biljanaorescanin Mar 15, 2026
3d03f10
fixed machine-precision issue for GPA22 (GPM20) peat processing; edit…
gmao-rreichle Mar 16, 2026
ccb56b2
fixed precision issue for GPA22 peat processing (mod_process_hres_dat…
gmao-rreichle Mar 17, 2026
bb3f763
fixed syntax error in previous commit (mod_process_hres_data.F90)
gmao-rreichle Mar 17, 2026
18a3ad4
Align strict GPA22 profile peat logic with top-layer peat support
biljanaorescanin Mar 19, 2026
4a0d5a0
change name option from hwsd_c to hwsd_v2
biljanaorescanin Mar 19, 2026
cbf7ca7
change name option from hwsd_v2 to hwsdv2
biljanaorescanin Mar 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions ...mp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,7 @@ def ask_questions(default_grid="Cubed-Sphere"):
"v11 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2", \
"v12 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2 + Argentina peatland fix", \
"v13 : NL3 + JPL veg height + PEATMAP + MODIS snow alb v2 + Argentina peatland fix + mean land elevation fix", \
"v15 : NL3 + JPL veg height + PEATMAP GPA22 + MODIS snow alb v2 + HWSDv1.2+ mean land elevation fix", \
"ICA : Icarus (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus/)", \
"GM4 : Ganymed-4_0 (archived*: /discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Ganymed-4_0/)", \
"F25 : Fortuna-2_5 (archived*: n/a)"],
Expand Down
570 changes: 337 additions & 233 deletions .../GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mod_process_hres_data.F90
View file
Open in desktop

Large diffs are not rendered by default.

156 changes: 121 additions & 35 deletions ...omp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ module rmTinyCatchParaMod

logical, public, save :: use_PEATMAP = .false.
logical, public, save :: jpl_height = .false.
logical, public, save :: PEATMAP_STRICT_GPA22 = .false.
Copy link
Contributor

@gmao-rreichle gmao-rreichle Mar 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The config variables and the variable names in the code will need some cleanup.

  • We won't need to run the legacy PEATMAP inputs with the "STRICT_GPA22" processing tweaks. It would be cleaner to identify the name of the peat map data product, and have the if conditions based on the product name.
  • The name "PEATMAP" is specific to the legacy peat map, and the name of the "GPA22" map in the published literature is "Global Peat Map 2.0". (I think GPA22 was just a working name during the earlier phase of product development.)
  • So instead of the logical "use_PEATMAP", we may want to have an integer config variables that something like:
  PEAT_INFO =              0 (HWSD), 
                           1 (HWSD and PEATMAP)
                           2 (Global Peat Map 2.0)

E.g., the if condition "if (use_PEATMAP)" would then translate to something like:

select case (PEAT_INFO)
case (0)
   ! use_PEATMAP=.false.
case (1,2)
   ! use_PEATMAP=.true.

etc

Copy link
Contributor Author

@biljanaorescanin biljanaorescanin Mar 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that makes sense. The current use_PEATMAP + PEATMAP_STRICT_GPA22 split was mainly introduced to isolate and validate the GPA22 logic while preserving legacy behavior. I agree that the long term configuration would be cleaner if it keyed off the peat information product explicitly rather than overloading the legacy PEATMAP name. I can work on cleaning that up and push a follow-up commit once it is done and retested.

character*8, public, save :: LAIBCS = 'UNDEF'
character*6, public, save :: SOILBCS = 'UNDEF'
character*6, public, save :: MODALB = 'UNDEF'
Expand Down Expand Up @@ -105,6 +106,7 @@ SUBROUTINE init_bcs_config(LBCSV)
! NGDC : Soil parameters from Reynolds et al. 2000, doi:10.1029/2000WR900130 (MERRA-2, Fortuna, Ganymed, Icarus)
! HWSD : Merged HWSDv1.21-STATSGO2 soil properties on 43200x21600 with Woesten et al. (1999) parameters
! HWSD_b : As in HWSD but with surgical fix of Argentina peatland issue (38S,60W)
! HWSD_c : HWSD updated from v1.2 to v2. And sub and top are D2 layer
Copy link
Contributor

@gmao-rreichle gmao-rreichle Mar 16, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of HWSD_c we should call this HWSDv2 for clarity.

!
! OUTLETV: Definition of outlet locations. DEFAULT : N/A
! N/A : No information (do not create routing "TRN" files).
Expand All @@ -125,6 +127,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 2.17
use_PEATMAP = .false.
jpl_height = .false.
PEATMAP_STRICT_GPA22= .false.

case ("GM4", "ICA")
LAIBCS = 'GSWP2'
Expand All @@ -135,6 +138,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .false.
jpl_height = .false.
PEATMAP_STRICT_GPA22= .false.

case ("NL3")
LAIBCS = 'MODGEO'
Expand All @@ -145,6 +149,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .false.
jpl_height = .false.
PEATMAP_STRICT_GPA22= .false.

case ("NL4")
LAIBCS = 'MODGEO'
Expand All @@ -155,6 +160,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .false.
jpl_height = .true.
PEATMAP_STRICT_GPA22= .false.

case ("NL5")
LAIBCS = 'MODGEO'
Expand All @@ -165,6 +171,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .true.
jpl_height = .true.
PEATMAP_STRICT_GPA22= .false.

case ("v06")
LAIBCS = 'MODGEO'
Expand All @@ -175,6 +182,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .true.
jpl_height = .true.
PEATMAP_STRICT_GPA22= .false.

case ("v07")
LAIBCS = 'MODGEO'
Expand All @@ -185,6 +193,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .true.
jpl_height = .false.
PEATMAP_STRICT_GPA22= .false.

case ("v08")
LAIBCS = 'MODGEO'
Expand All @@ -195,6 +204,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .false.
jpl_height = .false.
PEATMAP_STRICT_GPA22= .false.

case ("v09")
LAIBCS = 'MODGEO'
Expand All @@ -205,6 +215,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .true.
jpl_height = .false.
PEATMAP_STRICT_GPA22= .false.

case ("v10")
LAIBCS = 'MODGEO'
Expand All @@ -215,6 +226,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .true.
jpl_height = .false.
PEATMAP_STRICT_GPA22= .false.

case ("v11")
LAIBCS = 'MODGEO'
Expand All @@ -225,6 +237,7 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .true.
jpl_height = .true.
PEATMAP_STRICT_GPA22= .false.

case ("v12","v13")

Expand All @@ -241,6 +254,19 @@ SUBROUTINE init_bcs_config(LBCSV)
GNU = 1.0
use_PEATMAP = .true.
jpl_height = .true.
PEATMAP_STRICT_GPA22= .false.

case ("v15")
LAIBCS = 'MODGEO'
SOILBCS = 'HWSD_b'
MODALB = 'MODIS2'
SNOWALB = 'MODC061v2'
OUTLETV = "v2"
GNU = 1.0
use_PEATMAP = .true.
jpl_height = .true.
PEATMAP_STRICT_GPA22= .true.


case default

Expand Down Expand Up @@ -2573,6 +2599,8 @@ SUBROUTINE create_model_para_woesten (Maskfile, nbcatch, tile_lon, tile_lat, til
logical :: file_exists
REAL, ALLOCATABLE, DIMENSION (:,:) :: parms4file
integer :: ncid, status
! peat GPA22 additions
logical :: target_is_peat, donor_is_peat

! --------- VARIABLES FOR *OPENMP* PARALLEL ENVIRONMENT ------------
!
Expand Down Expand Up @@ -2979,8 +3007,10 @@ SUBROUTINE create_model_para_woesten (Maskfile, nbcatch, tile_lon, tile_lat, til

DO n=1,nbcatch

! Read soil_param.first again...; this is (almost certainly) needed to maintain consistency
! between soil_param.first and soil_param.dat, see comments above.
! Re-read soil_param.first for tile n so soil_param.dat starts from the
! original tile identity/properties written by mod_process_hres_data.
! This preserves consistency between soil_param.first and soil_param.dat
! before any bad-SAT donor repair is applied below.

read(10,'(i10,i8,i4,i4,3f8.4,f12.8,f7.4,f10.4,3f7.3,4f7.3,2f10.4, f8.4)') &
tindex2(n),pfaf2(n),soil_class_top(n),soil_class_com(n), &
Expand All @@ -2989,34 +3019,76 @@ SUBROUTINE create_model_para_woesten (Maskfile, nbcatch, tile_lon, tile_lat, til
a_sand_surf(n),a_clay_surf(n),atile_sand(n),atile_clay(n) , &
wpwet_surf(n),poros_surf(n), pmap(n)

! This revised if block replaces the complex, nested if block commented out above

if ( (ars1(n)==9999.) .or. (arw1(n)==9999.) ) then

! some parameter values are no-data --> find nearest tile k with good parameters
! If SAT parameters are invalid for tile n, choose donor tile k.
!
! Legacy behavior:
! use the nearest tile with valid SAT parameters.
!
! Strict GPA22 behavior:
! first prefer a donor with the same peat/mineral state as the
! target tile, so repaired bulk hydraulics do not cross peat and
! mineral regimes. If no same-state donor exists, fall back to
! the original nearest-valid donor search.

if ( (ars1(n)==9999.) .or. (arw1(n)==9999.) ) then
dist_save = 1000000.
k = 0

! Define target tile peat/mineral state from its current top/profile classes.
target_is_peat = (soil_class_top(n) == 253) .or. (soil_class_com(n) == 253)

do i = 1,nbcatch
if(i /= n) then
if((ars1(i).ne.9999.).and.(arw1(i).ne.9999.)) then

tile_distance = (tile_lon(i) - tile_lon(n)) * (tile_lon(i) - tile_lon(n)) + &
(tile_lat(i) - tile_lat(n)) * (tile_lat(i) - tile_lat(n))
if(tile_distance < dist_save) then
k = i
dist_save = tile_distance
endif
endif
if (i == n) cycle
if ((ars1(i) == 9999.) .or. (arw1(i) == 9999.)) cycle

! In strict GPA22 mode, only consider donors with the same
! peat/mineral state as the target tile.
if (PEATMAP_STRICT_GPA22) then
donor_is_peat = (soil_class_top(i) == 253) .or. (soil_class_com(i) == 253)
if (donor_is_peat .neqv. target_is_peat) cycle
endif

tile_distance = (tile_lon(i) - tile_lon(n)) * (tile_lon(i) - tile_lon(n)) + &
(tile_lat(i) - tile_lat(n)) * (tile_lat(i) - tile_lat(n))

if (tile_distance < dist_save) then
k = i
dist_save = tile_distance
endif
enddo
! record in file clsm/bad_sat_param.tiles
write (41,*)n,k ! n="bad" tile, k=tile from which parameters are taken

! Strict GPA22 fallback:
! if no same-state donor exists, revert to the original
! nearest-valid donor search so a donor is still always found.
if ((k == 0) .and. PEATMAP_STRICT_GPA22) then
dist_save = 1000000.
do i = 1,nbcatch
if (i == n) cycle
if ((ars1(i) == 9999.) .or. (arw1(i) == 9999.)) cycle

tile_distance = (tile_lon(i) - tile_lon(n)) * (tile_lon(i) - tile_lon(n)) + &
(tile_lat(i) - tile_lat(n)) * (tile_lat(i) - tile_lat(n))

if (tile_distance < dist_save) then
k = i
dist_save = tile_distance
endif
enddo
endif

if (k == 0) then
write(*,*) 'ERROR: no donor tile found for tile ', n
stop
endif

! Record indices of (repaired) target tile "n" and donor tile "k" in file "clsm/bad_sat_param.tiles"
write (41,*) n, k

! Overwrite parms4file when filling in parameters from neighboring tile k.
! For "good" tiles, keep parms4file as read earlier from catch_params.nc4,
! which is why this must be done within the "then" block of the "if" statement.
! This is necessary for backward 0-diff compatibility of catch_params.nc4.
! which is why this must be done within the "then" block of the "if"
! statement. This is necessary for backward 0-diff compatibility of
! catch_params.nc4.

parms4file (n,12) = BEE(k)
parms4file (n,16) = COND(k)
Expand All @@ -3027,11 +3099,10 @@ SUBROUTINE create_model_para_woesten (Maskfile, nbcatch, tile_lon, tile_lat, til

else

! nominal case, all parameters are good

! Nominal case: current tile n already has valid parameters.
k = n

end if
end if

! for current tile n, write parameters of tile k into ar.new (20), bf.dat (30), ts.dat (40),
! and soil_param.dat (42)
Expand All @@ -3041,21 +3112,36 @@ SUBROUTINE create_model_para_woesten (Maskfile, nbcatch, tile_lon, tile_lat, til
ars1(k),ars2(k),ars3(k), &
ara1(k),ara2(k),ara3(k),ara4(k), &
arw1(k),arw2(k),arw3(k),arw4(k)

write(30,'(i10,i8,f5.2,3(2x,e13.7))')tindex2(n),pfaf2(n),gnu,bf1(k),bf2(k),bf3(k)

write(40,'(i10,i8,f5.2,4(2x,e13.7))')tindex2(n),pfaf2(n),gnu, &
tsa1(k),tsa2(k),tsb1(k),tsb2(k)

write(42,'(i10,i8,i4,i4,3f8.4,f12.8,f7.4,f10.4,3f7.3,4f7.3,2f10.4, f8.4)') &
tindex2(n),pfaf2(n),soil_class_top(k),soil_class_com(k), &
BEE(k), PSIS(k),POROS(k),COND(k),WPWET(k),soildepth(k), &
grav_vec(k),soc_vec(k),poc_vec(k), &
a_sand_surf(k),a_clay_surf(k),atile_sand(k),atile_clay(k) , &
wpwet_surf(k),poros_surf(k), pmap(k)


! write "soil_param.dat" file; n = target tile, k = donor tile.

if (PEATMAP_STRICT_GPA22) then
! write only bulk hydraulic fields from donor tile k.
write(42,'(i10,i8,i4,i4,3f8.4,f12.8,f7.4,f10.4,3f7.3,4f7.3,2f10.4, f8.4)') &
tindex2(n),pfaf2(n), & ! n (target)
soil_class_top(n),soil_class_com(n), & ! n (target)
BEE(k), PSIS(k),POROS(k),COND(k),WPWET(k),soildepth(k), & ! k (donor)
grav_vec(n),soc_vec(n),poc_vec(n), & ! n (target)
a_sand_surf(n),a_clay_surf(n),atile_sand(n),atile_clay(n), & ! n (target)
wpwet_surf(n),poros_surf(n), pmap(n) ! n (target)
else
! Legacy path: preserve original donor-copy behavior.
write(42,'(i10,i8,i4,i4,3f8.4,f12.8,f7.4,f10.4,3f7.3,4f7.3,2f10.4, f8.4)') &
tindex2(n),pfaf2(n), & ! n (target)
soil_class_top(k),soil_class_com(k), & ! k (donor)
BEE(k), PSIS(k),POROS(k),COND(k),WPWET(k),soildepth(k), & ! k (donor)
grav_vec(k),soc_vec(k),poc_vec(k), & ! k (donor)
a_sand_surf(k),a_clay_surf(k),atile_sand(k),atile_clay(k), & ! k (donor)
wpwet_surf(k),poros_surf(k), pmap(k) ! k (donor)
endif

! record ar.new, bf.dat, and ts.dat parameters for later writing into catch_params.nc4

if (allocated (parms4file)) then
parms4file (n, 1) = ara1(k)
parms4file (n, 2) = ara2(k)
Expand Down
Loading
Loading