further cleanup of make_bcs utility and CatchmentCN routines

[gmao-rreichle]

An inconsistency in NLv5 Catchment model parameters was discovered that impacts 6 tiles (out of ~1.6 million) on the EASEv2 M09 grid. The 6 tiles in question end up with some PEATCLSM parameters (incl. porosity=0.93), even though they (probably) should have mineral soil properties (which is still reflected in some parameters associated with the 0-30 cm soil layer). An investigation is ongoing.

During this investigation, a few minor problems with the make_bcs script were discovered and fixed in the branch associated with this PR. More fixes might be needed.

The PR will be zero-diff for all model runs. It is not yet clear how the inconsistency in the Catchment model parameters will be resolved. We will probably need to introduce a new bcs version to fix the parameter inconsistency, which probably means that will the PR will also be zero-diff for the make_bcs package.

cc: @rdkoster @gmao-jkolassa @biljanaorescanin @weiyuan-jiang @gmao-qliu

UPDATED 11 Jun 2022: This PR will be merged without addressing the above-mentioned parameter inconsistency, which requires further study. The bulk of the PR consists of long overdue cleanup.

[gmao-rreichle]

UPDATE (5/5/22):
The subroutines dealing with HWSD/PEATMAP processing and the Woesten soil parameters were partially cleaned up and some documentation was added. The changes are not meant to be final. The latest commit ef69e84 successfully tested 0-diff for the generation of NLv4 and NLv5 bcs in EASEv2_M09 resolution. Additional tests for other bcs versions and resolutions will be needed.

It remains unclear how the tiles identified by Sebastian end up with peat as profile soil class. The attached pdf and xls files provide additional information.
SMAP_L4_Meeting_20220504_NLv5_bcs_only.pdf
soilparam_peat_tiles_rr_20220504.xlsx

cc: @rdkoster @gmao-jkolassa @biljanaorescanin @weiyuan-jiang @gmao-qliu

[gmao-rreichle]

I've run one bcs test (EASE M09)

Thanks, @biljanaorescanin. Which bcs version did you run? Given the extent of the changes in this PR, we need to run multiple versions at multiple resolutions. At the very least, we need to run NLv4 and NLv5 on M09 and M36, and we need to run NLv3 at c720, c360, c180 and c90.

[biljanaorescanin]

@gmao-rreichle I did NLv3. Sure, I will run other cases and add a comment once I am done.

[biljanaorescanin]

@gmao-rreichle @weiyuan-jiang I've added last commit so we don't have warning messages in catchcn_params.nc4 and catch_param.nc4 . I am testing it now to confirm.

[gmao-rreichle]

I added the "contingent" label again. Not sure why this was removed. We are not quite done with this PR. Hopefully final touches and testing underway now.

[gmao-rreichle]

The previous two commits fix a netcdf dimension error that had been present for a long time in the catch_params.nc4 and catchcn_params.nc4 files generated by make_bcs. The data in the files was successfully verified 0-diff after the commits.
Owing to the bug in the old version of these nc4 files, nccmp and cdo do not work properly for comparing the catch_params.nc4 files before and after the dimensions bug fix. I verified 0-diff for data across the old and new nc4 files using matlab.

@biljanaorescanin successfully completed final developer tests today (GEOSldas runs that use bcs generated with the branch associated with this PR).

[gmao-rreichle]

PR is ready for final verification and merging by GCM gatekeeper.

PR only touches make_bcs and CatchmentCN and is trivially 0-diff for all GCM tests.

Suggested log entry:

make_bcs and CatchmentCN cleanup:
- removed generation of unused/obsolete bcs file CLM4.5_veg_typs_fracs
- cleanup of incomplete CatchCN albedo scaling 
- partial cleanup and improved documentation of subroutines dealing with HWSD/PEATMAP processing and Woesten soil parameters 

[mbechtold]

I thought this is what we don't want.
o_clp=4 should only be allowed if o_clp makes up more than 0.5 of a grid cell.
This may explain that Sebastian found 6 'buried' peat tiles which were not peat in the top layer but in the profile layer.

[mbechtold]

Wouldn't it make sense to include in this 'if-statement' a line like
o_clp = 4
We want all o_cl to become full peat tiles.
With this we would get rid of those inconsistent tiles with top soil being peat but not profile being peat.

[weiyuan-jiang]

bug is here:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CombineRasters.F90

Lines 343 to 347 in 9c3baa1

	elseif(nint(Table2(1,Pix2))==fill) then ! A use 1st grid area of a merge
	iTable(2,ip) = nint(Table1(3,Pix1))
	iTable(3,ip) = nint(Table1(4,Pix1))
	rTable(4,ip) = Table1(2,Pix1)
	rTable(5,ip) = Table1(2,Pix1)

When line 344 is executed, Pix1 = -999 which leads to crash. It is exposed by debug mode but exists in optimized mode ( just does not crash) .

We may have to examine why it checks table 2 of pix2 but use table 1 of pix1

@gmao-rreichle @biljanaorescanin

[gmao-rreichle]

We may have to examine why it checks table 2 of pix2 but use table 1 of pix1

Thanks, @weiyuan-jiang, for the sleuthing. I've never looked at this program and cannot quickly make sense of what's going on here. I'm still unclear what happens when we run the (usual) build with optimization, and what a fix would look like. If we also get Pix1=-999 with the optimized build, then lines 344-347 should result in garbage, which is concerning. @atrayano, are you familiar with the code in question? [In brief, when we run this code for make_bcs using a debug build, the run stops with an error message that an index value is negative. No such stop and error message is produced when we run the (usual) optimized build. See previous comment by @weiyuan-jiang.]

cc: @biljanaorescanin @rdkoster

[atrayano]

@gmao-rreichle Yes, I am familiar with this code , it was written by Max Suarez many years ago. I have seen the negative (-999) Pix error once in a while, and usually this is a result of the rasterization process not covering completely the Earth

[weiyuan-jiang]

The good news is I only saw one bad pix. Maybe we can throw away that one ? @gmao-rreichle @biljanaorescanin @atrayano

[atrayano]

@weiyuan-jiang Unfortunately, we cannot throw away the "bad" Pix. My first impulse is to write a piece of code to check the raster(s) for -999 and possibly do something to re-assign a "good" value (for example, assign the value from a neighbor)

[biljanaorescanin]

@sdrabenh I made an issue for things we didn't address #596, we are ready to proceed with this PR as it is.

@gmao-rreichle left a comment for suggested log entry:

make_bcs and CatchmentCN cleanup:

• removed generation of unused/obsolete bcs file CLM4.5_veg_typs_fracs
• cleanup of incomplete CatchCN albedo scaling
• partial cleanup and improved documentation of subroutines dealing with HWSD/PEATMAP processing and Woesten soil parameters