Reformat BIOS ancillary information

[ccarouge]

In the BIOS setup, the ancillary information is read in 2 steps:

• read gridinfo file
• read new information in binary files that replaces some of the information in the gridinfo file. Adds limits to the parameter values in CABLE. These limits will need to be applied in a pre-processing step

The binary file contains also additional variables that are not currently in the gridinfo file.

We need to decide how to deal with this additional information:

1. reformat to netcdf, store in a separate file and keep the current reading procedure by overwriting the information in CABLE. This would result in having a generic gridinfo file and then additional files for localised information.
2. or reformat to netcdf and overwrite the gridinfo file with that information. This would result in creating gridinfo files specific to each configuration as needed with different formats of gridinfo files for different configurations.

Solution 1

It keeps the gridinfo file format constant but then adds a precedence rule between files. It makes it less clear to users what is the final source of the input information. Adds cases that need to be handled by CABLE instead of putting them in the pre-processing.

Solution 2

It collates the information into one file which clarifies the source of the information to the users. Although, currently CABLE has rules of precedence in its input files. gridinfo has the least priority, then the met files and finally the restart file. It would require creating different gridinfo files for different configurations and the format of the gridinfo file would not be constant between configurations.

Problem: If we collate all the BIOS information in the gridinfo file, do we risk using different information coming from the restart file for example? Is this case improbable and only probable in case of user error? In the current implementation, does BIOS overwrite the information after reading in only the gridinfo file or after reading in all the other input files?

Implementation

We tend to favour solution 2 but we need to clarify if the identified problem makes it impossible to implement it this way.

Some start towards an implementation is here: /g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_forcing/gridinfo_inprogress/convert_gridinfo_Aust.sh.

Input data for this script is here: /g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_forcing/australia_0.05/params/nc.

[ccarouge]

@abhaasgoyal Do we want to use this issue to document the work for the gridinfo file? It probably needs some updating with our current plan.

[ccarouge]

This work will be performed using the act9test-serial-no-luc test case. We will work from the CABLE-POP_TRENDY branch directly. We need a first control run using the original binary inputs. As we convert inputs to netcdf format, running the same test case should give us identical results.

[abhaasgoyal]

Current implementation plan:

1. Run /g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_forcing/australia_0.05/params/nc/scripts/gdal_translate_flt2nc.sh to convert individual ancillary binary files to NetCDF. So I'll have to update the I/O routines to read NetCDF (some files are given here for example. Run tests as mentioned in the previous comment
2. Combine all ancillary NetCDF files into one NetCDF and read from there
3. I wasn't sure about this, but I think since we are inclining towards solution 2, we need to check the viability considering precedence in reading inputs.

@SeanBryan51 could you check whether there was anything to add there

[SeanBryan51]

This looks good to me. In step 1 we can do a regression test to test the input files were converted correctly to NetCDF and in step 2 we can do a regression test on the combined input files (the gridinfo file). We will have to look more into step 3 and make sure values initialised from the gridinfo are not being overwritten by other sources (e.g. a restart file).

[ccarouge]

@abhaasgoyal I'm happy to look through the code with you about the precedence of inputs. We might be able to get the answer quickly. Doing it with you allows you to understand where it comes from.

[SeanBryan51]

Requirements:

1. BIOS gridinfo needs to be defined over Australian domain: longitude (degrees east) = 112.0, 154.0, latitude (degrees north) = -44.0, -10.0
2. BIOS gridinfo needs to be defined at 0.05 degrees resolution.

Context

• The BIOS ancillary and meteorolgical forcing inputs are available at 0.05 degrees resolution over Australia. See /g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_forcing/australia_0.05/.
• The BIOS gridinfo file will need to be a combination of the 1 degree 0.5 degree gridinfo file (gridinfo_CSIRO_CRU05x05_4tiles.nc) subsampled at 0.05 degrees (paraphrasing Ian here):
  • This currently happens inside CABLE: i.e. layers where we substitute in BIOS equivalents of existing layers and add new BIOS specific layers. BIOS currently reads both the 1 degree 0.5 degree gridinfo and the binary ancillaries since we do not have Australia specific layers for all inputs (and the gridinfo provides some initial conditions) nor does the current gridinfo provide global layers for all BIOS inputs.
  • A 0.05 degree resolution BIOS gridinfo file is also required so that we can use the land-mask capability of CRU-TRENDY to select the relevant points for the other cases (act9 and 1000pts) rather than having multiple gridinfo files for each configuration. This would be a) more flexible and b) allow consolidation around one set of ancillary reading routines.

[SeanBryan51]

I'm currently looking into the regridding of a 1 degree 0.5 degree gridinfo (gridinfo_CSIRO_CRU05x05_4tiles.nc) to 0.05 degrees and have some questions:

1. Does the exact location of each grid cell matter? For example, say we are working with the Australia domain, should grid points lie on multiples of 0.05 degrees (i.e. lat = [-10, -10.05, -10.1, ...], lon = [112, 112.05, 112.1, ...]) or be offset by 0.025 degrees (i.e. lat = [-10.025, -10.075, -10.125, ...], lon = [112.025, 112.075, 112.125, ...])? It looks like the 1 degree gridinfo is defined on lon = [-179.5, -178.5, -177.5, ...], lat = [89.5, 88.5, 87.5, ...], i.e. offset by half a degree.
2. It seems most land quantities can use nearest neighbour interpolation but are there any exceptions? E.g. variables dependent on the area of the grid cell likely cannot be interpolated with nearest neighbour.

[ccarouge]

@SeanBryan51

Question 1: Place the centres to the offset point by half the resolution. This means we avoid having to deal with half gridcells in latitude at the poles, since the cell boundary lies on -90 and +90°. And that simplifies the regridding (see response to 2 below).

Question 2: Going with the answer of question 1, all the grid cells in the new grid at 0.05° will fit entirely within a grid cell at 1°. So we can use nearest neighbour for all interpolation. All the quantities are given per m2 when they depend on the area. Obviously, if a new grid cell straddled several grid cells that would be more difficult.

[abhaasgoyal]

Currently we have 4 functions to convert reading BIOS inputs in binary form to NetCDF:

1. cable_bios_load_params - Sets soil parameters (soil) after calculations. NetCDF read in gridinfo file
   Assuming that we plan to read the gridinfo file gridinfo_AUST_005x005.nc, which is being generated with the following script:

/g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_forcing/gridinfo_inprogress/convert_gridinfo_Aust.sh

We see that the script preprocesses the input doing the necessary calculations.
2. cable_bios_load_fracC4 - Sets inVeg in LUC runs after calculations. NetCDF read in gridinfo, however it doesn’t seem to be processed based on the script.
3. cable_bios_load_biome - Loads LUC_EXPT%biome. NetCDF read from LUC_Expt%ClimateFile. No additional calculations.
4. cable_bios_load_climate_params - Set climate%modis_igbg and climate%AvgAnnMaxFAPAR after caclulations. NetCDF read in LUC_Expt%ClimateFile.

We plan to work with reading ancillary data purely from NetCDF files one function at a time considering precedence. Currently going for soil parameters in cable_bios_load_params.

To check precedence (in increasing order), currently, the order of reading parameters is:

1. Input parameter values from namelist + grid info
2. Restart file
3. Met forcings file
4. Ancillary data currently in binary form

Now, putting all ancillary variables into grid info is means that the metfile would have the highest precedence. We want to keep these values only when the parameters in met forcings are better. Hence, we decided on having an additional flag for reading met data or not.

Before further updating the issue, @har917 I had some questions:

Q1: There is a comment in convert_gridinfo_Aust.sh - #may need to add the IGBP major biome, fracC4, modis_igbp and AvAnnualMaxFAPAR layers, which applies to (2), (3), (4). Currently it is unclear whether these are required in main gridinfo file or as separate NetCDF file LUC_Expt%ClimateFile.
Q2: Whether the versions in NetCDF files would be processed versions compared to binary ancillary data (hence no more calculations are needed within CABLE).
Q3: If we were to implement the plan of action, the restart file values can always overwrite ancillary data, when in the original scenario ancillary data has the highest precedence. Is this behavior what we want? Paraphrasing @ccarouge - “Normally if people run a restart simulation, for it to really be a restart, it needs to ensure all the data is the same as in the initial simulation. Hence overwriting the data with the restart values. That said, all the data in the ancillary files we are working on here are constant for the run so it should be safe to assume that the restart and the ancillary data (step 4) have the same values. Hence the order of reading them should not be important.”

[har917]

@abhaasgoyal @ccarouge

A couple of comments.

most importantly: I suspect that convertgridinfo_Aust.sh may have an error in it. Recent investigations have highlighted that BIOS uses a different gridinfo file gridinfo_CSIRO_CRU0.5x0.5_4tiles.nc sitting in the x45 g/data space (which is likely different to the gridinfo file referenced in that draft script).

On 1. BIOS requires some additional gridded inputs - it is a software engineering choice as to whether these should form part of a single BIOS_gridinfo file or as a separate netcdf input. My position is that it is better to combine into one for input management purposes (other applications of CABLE won't necessarily have to read in the additional layers so it shouldn't break them) - however this may come down to exactly where those inputs are read in by the code. Currently the separate files are read in at separate locations to the gridinfo - it may be better to retain that separation and retain separate input files.

On 2. I would hope that we can avoid having to have additional QA steps on the inputs inside CABLE. It creates confusion/lack of transparency - these additional QA steps are better placed as comments in the netcdf metadata in my mind*.

On 3. I don't have good visibility of what gets written into the restart file in terms of ancillary values. I wouldn't be surprised if none/not all of the BIOS gridded ancillaries make their way into the restart (NB there are gridded ancillaries and parameter values derived from gridded ancillaries to be aware of - especially in CASA). I think I would prefer if values in the restart were to take precedence over the ancillary data in the case of a difference (there are some use cases where it could be useful to go the other way round - think parameter estimation, predictability or sensitivity studies - but I think those can be handled easily outside as a preprocessing step on an existing restart) Perhaps the way forward - in the longer term - is to implement a check with a warning to log if a difference is found?

*Having said that we've gone and added some new ones recently due to path-of-least-resistance issues.

[ccarouge]

@har917 @abhaasgoyal

On 1:

Currently the separate files are read in at separate locations to the gridinfo - it may be better to retain that separation and retain separate input files.

What are you worried about? Or why do you need that capability? For the moment, we have identified the restart and the met file as potential alternate source for the data that are read after the gridinfo file but before the BIOS reading call. Ian, you agree that the restart should take precedence so changing the behaviour here would be ok.
It leaves the met file. We are proposing to make reading the parameters from the met file optional.

I would like to try and read the gridinfo file in one location and see if we get the same values in the corresponding CABLE variables as when reading the BIOS files.

If we don't get the same values throughout the code then we know something else is overwriting these values and we can hunt down what is happening.

If we get the same values, then we don't need to complicate things.

On 2: great, we agree. All pre-processing in input file creation.

@SeanBryan51 Ian above is saying BIOS is using a gridinfo file at 0.05°, you talk of starting from a 1° gridinfo file. Need a check?

[har917]

What are you worried about? Or why do you need that capability?

I don't think we need that capability per say. However, changing where/when the code is currently reading inputs risks mixing up the process of converting the inputs, and reading them in, and (in this case) porting them around the code base.

@SeanBryan51 For clarity the BIOS gridinfo file referenced is at 0.5 degree resolution - it maybe global (land north of >60S) - and has been set up to provide default initial conditions for a 4 tile version of CABLE (primary forest, secondary forest, grass, crop). The usual gridinfo file provides default initial conditions for 17 tiles (if memory serves).

[SeanBryan51]

@abhaasgoyal here is the script I used to generate the 0.05 degree gridinfo: https://gist.github.com/SeanBryan51/f89a30ee1e22495965dccb8fdff0f8c4

The gridinfo file I used for the input is the same as what is used in the act9test-serial-no-luc namelist file, i.e. /g/data/rp23/data/no_provenance/gridinfo/gridinfo_CSIRO_CRU05x05_4tiles.nc.

[ccarouge]

@SeanBryan51 I had a look at the regridding script for gridinfo. It looks alright.

I've tried the script on the 0.5° gridinfo under /g/data/rp23/data. The grid files have the correct spacing between the centres and the corners are at the right places. The map file is harder to understand but at least the frac_* variables only have 0s and 1s in, that is what you want for a nearest neighbour interpolation.

I had a side-by-side look at the data from the original gridinfo and the new one. I selected the same region on the new one first and then plotted with ncview. All the variables seem to have the same values and same patterns. The patchfrac variable is extremely weird but that's true of the original data.

The files differ because the regridding adds:

• lat_bnds and lon_bnds, that's a good idea
• gw, latitude quadrature weights, not sure that's particularly useful but it does not hurt.

[SeanBryan51]

@ccarouge Thanks for trying out the regridding script.

I have tested CABLE (act9test-serial-no-luc) with the 0.05 degree gridinfo. Model output using the new gridinfo is bitwise identical with the output from the original configuration (provided we ignore the Area variable¹).

Footnotes

1. It seems strange that there are differences in the Area variable since the underlying grid cell resolution for both runs should be at 0.05 degrees. Perhaps the Area variable in the model output does not accurately reflect the resolution of the model? ↩

[har917]

@ccarouge @SeanBryan51 - good to hear of progress.

The Area variable in CABLE output is known to be suspect - I certainly never use it and recalculate outside of the model. This is complicated by the need to handle single sites as well as gridded runs (including subsampling of grids - e.g. should a simulation that samples every second grid cell have the same Area as the full run or 4 times the Area?). Likely needs a separate issue.

[ccarouge]

@ccarouge Some of the reads are done into the LUC code, where do we want to put these layers? Identify what is read.

[har917]

@ccarouge I think I need a bit more detail as to the question here. However, we have to be careful about interpreting routine names - this is because the PFT distribution and land use change code is interwoven (and somewhat oddly named)
The reads (I think) you are talking about (C4frac, MVG) are done in the bios_cable_met_params and LUC_expt sections - this code sets up the model (specifically the PFT distribution and other parameters). It is necessary independently whether land use is called or not (though different bits get called with land use change on). MVG and C4frac are in effect ancillaries of equivalent type as the soil parameters.

Quite where the actual read should take place is a different question - I would be looking to see how TRENDY_POP does this (or similar) and follow the same pathway. This may require carrying some irrelevant layers (filled with zeros) around the code base for non-BIOS simulations - but that's probably better than having multiple reads of the same file.

[ccarouge]

Thanks Ian. Looking at the details here is what I have:

BIOS variable name	BIOS reading routine	CABLE variable	POP-TRENDY variable name	POP_TRENDY filename	POP-TRENDY reading routine
MVG	cable_bios_load_biome()	LUC_EXPT%biome	biome	LUC_EXPT%ClimateFile	READ_ClimateFile()
c4frac	cable_bios_load_fracC4()	no global variable. Local fracC4 variable used. Replaces LUC_EXPT%mtemp_min20 used in POP_TRENDY but not same unit/definition (more details at the end)	mtemp_min20	LUC_EXPT%ClimateFile	READ_ClimateFile()

So we would want these 2 BIOS layers to get into the ClimateFile but I'm guessing this file is created by CABLE in the "do climate" phase? Is that correct @har917 ? Or would we want a way to overwrite the data from the ClimateFile?

The handling of fracC4 (why the variable for the filename in the file is C4frac_file 😢 ) is also more problematic as it's not the same variable as what POP_TRENDY uses. I would think we could read fracC4 or mtemp_min20 depending if bios is on or off? The question is whether the ClimateFile should have both fields in or whether we should have 2 files: one for POP and one for BIOS?

more details on the differences between fracC4 and mtemp_min20

Both these variables are used to determine whether the grass at a point is a C3 grass or a C4 grass (they have different photosynthesis).
BIOS is using fracC4 in this way:

CABLE/offline/cable_LUC_EXPT.F90

Lines 650 to 654 in 7326db3

	if (fracC4(k).gt.0.5) then
	inVeg(m,n,2) = 7 ! C4 grass
	else
	inVeg(m,n,2) = 6 ! C3 grass
	endif

while mtemp_min20 is used here:

CABLE/offline/cable_LUC_EXPT.F90

Lines 575 to 579 in 7326db3

	IF (LUC_EXPT%mtemp_min20(k) .LE. 15.5) THEN
	inVeg(m,n,2) = 6 ! C3 grass
	ELSE
	inVeg(m,n,2) = 7 ! C4 grass
	ENDIF

[SeanBryan51]

I'm looking into adding the BIOS ancillaries¹ to the 0.05 degree gridinfo file and it looks like there is a slight mismatch in the lat lon grids between the gridinfo and the BIOS ancillaries:

1. The grid in the BIOS ancillaries is not offset by 0.025 degrees (i.e. lat=[-43.55,-43.5,-43.45,...] and lon=[112.95,113,113.05,...]) where as the gridinfo is offset by 0.025 degrees (i.e. lat = [-10.025, -10.075, -10.125, ...], lon = [112.025, 112.075, 112.125, ...]).
2. The lat lon grid of the BIOS ancillaries are defined on the domain: lat=[-43.55,-10.1] lon=[112.95,153.55] where as the gridinfo is defined on the domain lat=[-43.975,-10.025] lon=[112.025,153.925]. I.e. It looks like the BIOS ancillaries are missing data along the boundary of the Australia domain.

To combine the ancillary inputs with the gridinfo file, should we remap the BIOS ancillaries to the grid defined in the gridinfo file? In this case we require the BIOS ancillaries to be defined on the full Australia domain (i.e. latitude (degrees north) = -44.0, -10.0 and longitude (degrees east) = 112.0, 154.0). Alternatively, we could reduce the domain of the gridinfo so that all points in the BIOS ancillaries map to points in the gridinfo file and then do the remapping. What would be the best approach here?

Footnotes

1. BIOS ancillaries are in NetCDF format and are found here: /g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_forcing/australia_0.05/params/nc ↩

[SeanBryan51]

To combine the ancillary inputs with the gridinfo file, should we remap the BIOS ancillaries to the grid defined in the gridinfo file? In this case we require the BIOS ancillaries to be defined on the full Australia domain (i.e. latitude (degrees north) = -44.0, -10.0 and longitude (degrees east) = 112.0, 154.0).

I was overthinking it, the missing grid points in the BIOS ancillaries don't contain any land so we can set the value of each variable to its missing value over these grid points before remapping to the gridinfo grid.

[har917]

@SeanBryan51 quick responses

1. I suspect that we've picked up a difference between cell centres and cell corners here not a misalignment of the grid (and this may come from the original data sources). I suspect the easiest way to figure this out will be to compare coastlines (i.e. non-filled values) and then settle on one shifting of the data
2. yes - the difference is likely to be at the boundaries and all (most) of these cells are going to be non-land.

We may also have to look into the binpackandsample routine to see if this does different things for different variables/input layers.

[har917]

@ccarouge referring to ^^^^

You've picked up a fundamental difference in how POP_TRENDY and BIOS work - and it complicated because of the multi-stage aspects to POP runs. This likely requires a specific CASE() to be devised during this merge.

• In BIOS both the biome (MVG) and C4frac are known, fixed inputs - equivalent in structure to the soil ancillaries.
• BIOS does evaluate a mtemp_min20 and places this and biome in the ClimateFile (but mtemp_min20 is not used - at least in the same way)
• POP_TRENDY uses the 1st (do_climate) run to create a set of climate layers (including mtemp_min20) and puts these in the ClimateFile. On the 2nd run through (do_zerobiomass) these climate layers allow the evaluation of LUC_expt%biome (via the BIOME-1 submodel) and the specification of the PFTs and LUC_expt%woodfrac. On subsequent run throughs (do_equil1 etc.) information is picked up from the CLimateFile (including biome)

Since mtemp_min20 and POP_TRENDY biome are climate related these vars get put into the ClimateFile not in the gridinfo file.

I think my position on the way forward is (but I will have to check the sequencing) that

• BIOS gridinfo should include these two layers
• if BIOS these layers are read in at the same place as the rest of the grid info - then passed through (new fields required) to LUC_EXPT [if not BIOS then we will pass a field of zeros]
• if TRENDY we retain the existing method i.e. mtemp_min20 and biome via ClimateFile.

There is a bit of uncertainty in my head as to when these variables come from the ClimateFile and cable_bios_load_biome()/cable_bios_load_c4frac() vs when they come from one of the restart files (likely POPLUC). It could be that it's a case that i) read from ClimateFile or via MVG/C4frac then ii) if found in the ClimateFile overwrite.

[ccarouge]

Once we have a gridinfo file and the BIOS information on the same grid (without the half-grid cell shift), do we need to ensure all the data is using the same land mask?

This means a grid cell is land if both the original gridinfo and the BIOS variables have valid data for that point. That would mean the BIOS data would lose some islands and some precision in some caps. Do we want this to happen @har917 ?

@har917 Do you have some output from an Australian-wide BIOS run? Just to check what the coastline and islands look like. Or @AlisonBennett . This is just to see the output grid, we don't care about the variables and values in the file.

[har917]

@ccarouge As far as I know we have never run BIOS at the full 0.05 degree resolution (we have run AWAP using most/all the same inputs). This is mostly because of resourcing/set up - in that I understand it's impossible to run BIOS on gadi with enough memory/run time.

From memory/discussions with Peter Briggs - the model is fairly ruthless in that it won't run if any of the inputs aren't valid (but I don;t know the detail as to how it is does this in practice). The main cause of concern wasn't islands/coasts but incland areas (lake Eyre) where the input soil ancillaries were/area incomplete.

We could try run to a run a do_climate section - that should read in the current ancillaries (and default initial conditions) to see what the effective resultant landmask looks like (but perhaps not immediately as there's lots going on and potential for confusion).

[ccarouge]

incland areas (lake Eyre) where the input soil ancillaries were/area incomplete
@har917 I can see there are a lot of areas without data inland. What scenario would make sense for a BIOS configuration for these points:

1. have them mask everywhere in all data and not produce simulation results for these points.
2. use default data instead of BIOS-specific info for these points.
   The second case is potentially confusing as it would present the results as being from a BIOS simulation when, in fact, no BIOS-related information is used in the simulation for these points.

[har917]

IF I remember correctly the current approach produces NaNs

I wouldn't use the default gridinfo - that will be inconsistent (especially if we get some soil parameters from BIOS and the other gridinfo). What maybe better to do is gap-fill using nearest neighbour from the BIOS data and somehow note which cells have been filled.

[ccarouge]

We should try and quantify the impact of changing the ancillaries due to the derived_parameters bug.

[SeanBryan51]

Update: I have managed to bitwise reproduce the output created with a gridinfo containing only the BIOS bch by updating the derived parameters which are dependent only on the BIOS bch variable and the original gridinfo parameters (see c6fff2f). In summary: c6fff2f with original gridinfo == 7326db3 with BIOS gridinfo (with bch only).

[SeanBryan51]

Update: there was a bug in regrid_aus_05x05_to_005x005.sh where I was applying the wrong limits for the sucs variable. I have fixed this and I now have bitwise equivalence with 7326db3 with BIOS gridinfo and bfa87d3 with original gridinfo. This verifies that calling derived_parameters after reading in the BIOS parameters has the same effect as using the BIOS gridinfo file.

[ccarouge]

gridinfo file is ready with the BIOS information in it. Created from the regrid_aus_05x05_to_005x005.sh script.
Remove calls to read the BIOS information directly.
Add metadata attributes.: ask Peter for information from the original process for provenance information.

Run Australia-wide at 0.25 resolution or a 1000 points simulation for 1 year ?? See how to do it.

[SeanBryan51]

I did some quick testing of the BIOS gridinfo file for the 1000pts config with the BLAZE_9184 branch (commit 7dead54). CABLE errors due to changes in the number of patches and is incompatible with the existing restart files. My guess is this is likely due to the half grid-cell misalignment between the gridinfo grid and the BIOS grid resulting in ambiguous nearest neighbour interpolation of the gridinfo file within CABLE. It was probably pure luck that the shifted gridinfo file was bit reproducible for the act9 configuration.

I'm thinking the next steps around this would be to apply a fix to the nearest neighbour interpolation of the gridinfo file in CABLE such that we have bit reproducibility before and after for runs using the shifted gridinfo file.

Context: earlier on when I was looking into how best to shift the gridinfo file to the BIOS grid, I found running CABLE with a shifted gridinfo file (in arbitrary directions in latitude and longitude coordinates) broke the act9test-serial-no-luc configuration with the error Number of patches in restart file input/bios_cable_rst.nc does not equal to number in default/met file settings. (SUB load_parameters) Recommend running without restart file. This was because the BIOS grid contains points which are equidistant to multiple grid cells on the gridinfo grid which results in ambiguous nearest neighbor interpolation (nearest distance is dependent on the floating point precision of lat-lon coordinates) - this causes the total number of patches to change as grid points are initialised from slightly different grid cells on the gridinfo file. However, if I shifted the grid in such a way, I would reproduce the interpolation that was done internally within CABLE and not break the act9 configuration (as mentioned here).

[har917]

It would be good to dig into this a little more deeply - is this because the half-grid move means that we get changes on the coast where we have meteorology and all soil parameters available? Or has something more fundamental happened that means were getting different number of patches for grid cells in the interior?

The first case should be relatively easy to handle since the input meteorology files are more extensive than actually provided into the binary files - if the latter then this is more problematic (and I can't think why it could happen)

The way to test would be to run without using the CABLE restart file - and let the model figure out the number of patches it wants to use and then compare.

[ccarouge]

@har917 Here is an example to explain what happens with the grids. The gridinfo at 0.5° has these grid cell centres for latitude (similar story with longitude so just showing one here):
87.75, 87.25, 86.75, 86.25, 85.75
In other words, the grid cell boundaries fall on "whole" latitudes and longitudes: 88.0,87.5, 87.0, 86.5, 86.0 etc.

The BIOS grid has the following grid centres:
87.0, 87.05, 87.1, 87.15, 87.2, 87.25
In other words, for every few grid cells, the centre of the BIOS grid cell falls right on the edge of a gridinfo grid cell. That means the centre of the BIOS grid cell is equidistant to 2 gridinfo cells (or 4 gridinfo grid cells when it falls on corner of the gridinfo grid cells).

Now the way this works in CABLE:

1. CABLE reads in the BIOS met forcing and sets the grid using these. So the grid cells in the CABLE simulation are the same as the BIOS grid cells.
2. CABLE reads in the gridinfo file at a different resolution. So CABLE will go and try to find the closest grid cell. But because of the shift in the cells centres, for a few of the cells, the centres are equidistant of several grid cells in gridinfo. There is no handling of this case in CABLE. We simply "use" the fact floating points can't be at equal distance in programming, the precision will always put one closer than the other.

We end up with a random choice of gridinfo grid cell for these BIOS grid cells. Some might take the point that is to the East or North-East or South or West etc.

When we regrid the gridinfo file onto the BIOS grid, we have to choose one direction to shift the map towards. It will be the same direction for the whole grid, not a random direction that is different for different points like in CABLE.

As a result, we end up with some points using the information from the gridinfo file from a neighbouring grid cell in the new setup compared to the old which can have a different number of tiles. And this can happen randomly throughout the whole map.

[har917]

@ccarouge

which can have a different number of tiles.

For BIOS (i.e. with POP active) I don't think this is true - at least in the absence of land-use change. BIOS only ever has (up to) 3 tiles per grid cell and this decision (i.e. 2 or 3 tiles per grid cell) is independent of the soil ancillaries.

Coastal cells/soil ancillary gaps are a potential counter example - since if there are no valid soil ancillaries at the grid cell, BIOS will remove those tiles from the vector of land points even if there is valid meteorology available (or it will crash in a different way). Shifts in the gridinfo could then lead to a different number of tiles/land points.

However - if we are testing with a non-POP configuration then, I agree, shifts can change the number of tiles.

[abhaasgoyal]

Update: I ran act9test configuration script for 0.05° resolution gridinfo (gridinfo_CSIRO_CRU005x005_4tiles.nc), which should have the processed ancillary inputs since I used @SeanBryan51 's conversion script with --append-bios option. The tests were run in 2 versions of CABLE

1. CABLE-POP_TRENDY as the base branch. Here, the input values read from the gridinfo are overwritten by the binary ancillary inputs loaded/processed in cable_bios_load_params
2. BIOS_no_bios_input with the difference of cable_driver not calling cable_bios_load_params.

Normally, this should not create any difference in the results since the new gridinfo should have the same ancillary data after processing inputs from binary files in cable_bios_load_params. However bitwise comparison of the output files show difference in values for sucs (postive/negative sign):

ext_newblazempi_off/run/outputs/bios_out_cable_1700_1899.nc
DIFFER : VARIABLE : sucs : POSITION : [0,0] : VALUES : 0.5989 <> -0.5989
DIFFER : VARIABLE : sucs : POSITION : [0,1] : VALUES : 0.4092 <> -0.4092
DIFFER : VARIABLE : sucs : POSITION : [0,2] : VALUES : 0.3795 <> -0.3795
DIFFER : VARIABLE : sucs : POSITION : [0,3] : VALUES : 0.4092 <> -0.4092
DIFFER : VARIABLE : sucs : POSITION : [0,4] : VALUES : 0.4092 <> -0.4092
DIFFER : VARIABLE : sucs : POSITION : [0,5] : VALUES : 0.3795 <> -0.3795
DIFFER : VARIABLE : sucs : POSITION : [0,6] : VALUES : 0.3795 <> -0.3795
DIFFER : VARIABLE : sucs : POSITION : [0,7] : VALUES : 0.4092 <> -0.4092
DIFFER : VARIABLE : sucs : POSITION : [0,8] : VALUES : 0.3795 <> -0.3795
DIFFER : VARIABLE : sucs : POSITION : [1,2] : VALUES : 0.3795 <> -0.3795
DIFFER : VARIABLE : sucs : POSITION : [1,5] : VALUES : 0.3795 <> -0.3795
DIFFER : VARIABLE : sucs : POSITION : [1,8] : VALUES : 0.3795 <> -0.3795
ext_newblazempi_off/run/outputs/bios_out_cable_1900_2022.nc
DIFFER : VARIABLE : sucs : POSITION : [0,0] : VALUES : 0.5989 <> -0.5989
DIFFER : VARIABLE : sucs : POSITION : [0,1] : VALUES : 0.4092 <> -0.4092
DIFFER : VARIABLE : sucs : POSITION : [0,2] : VALUES : 0.3795 <> -0.3795
DIFFER : VARIABLE : sucs : POSITION : [0,3] : VALUES : 0.4092 <> -0.4092
...

We suspect that it may be a bug in BIOS due to write_default_params changing sucs sign , but cable_bios_load_params not doing something similar.

CABLE/offline/cable_parameters.F90

Line 1343 in e190b23

-1.* insucs(landpt(e)%ilon, landpt(e)%ilat)

Further investigation is needed on whether this is a bug in BIOS. Once this is fixed, we can have assurance on values contained in the new gridinfo file.

---

Addendum: I also ran the 1000pts configuration in a similar manner with the 2 branches of CABLE, however I am getting errors in a lot of variables (original as well as derived). Need to test further on reproducing the outputs, as well as finding the cause.

[har917]

@abhaasgoyal I think this is a matter of definition - it appears that the gridinfo file has been set up to specify -sucs whereas the BIOS binary and the methods that specify these parameters via namelist give sucs itself. I didn't spot that change in sign in cable_parameters() when I was doing the early work to identify how to connect the BIOS binaries to the gridinfo that @SeanBryan51 based his work on.

All very confusing - especially when you bring the vectorised soil properties (i.e. the _vec variables) and the GW scheme.

I'm a bit surprised that the -ve sign didn't show up in other testing but that maybe because BIOS uses a set of science options where %sucs is always used within an ABS() function.

[ccarouge]

After checking, this is because the gridinfo file expects positive values for sucs which are then changed to negative in CABLE. But the BIOS input is already negative.

So, @abhaasgoyal is going to modify the creation script for gridinfo to swap the sign of the input before writing to gridinfo.

It seems like a lot of sign changes applied on top of each others but we can't remove that from CABLE as that would break things for all the other gridinfo files in use.

[abhaasgoyal]

Update: After updating the sucs variable, act9test configuration output is bitwise-reproducible in CABLE-POP_TRENDY vs BIOS_no_bios_input branch.

However, running the 1000pts configuration (see branches 1000pts-new-params and 1000pts-new-params-bios-inputs-tmp) again provides mismatch for binary inputs and derived variables in cable_bios_load_params. A pattern that I noticed was that the differences in most variables were in position [x, 0/1/2, 16/32/79/122]. Essentially, the differences were in z axis in position 16/32/79/122. An example is provided as follows:

DIFFER : VARIABLE : visAlbedo : POSITION : [0,0,16] : VALUES : 0.058107 <> 0.0578484
DIFFER : VARIABLE : visAlbedo : POSITION : [0,0,32] : VALUES : 0.0768954 <> 0.0811007
DIFFER : VARIABLE : visAlbedo : POSITION : [0,0,79] : VALUES : 0.0695602 <> 0.0667755
DIFFER : VARIABLE : visAlbedo : POSITION : [0,0,122] : VALUES : 0.0499809 <> 0.049981

I'll attempt to debug this on whether the ancillary inputs from binary files vs new gridinfo were processed to be the same, and where it differs in the pipeline. Meanwhile, @AlisonBennett if you could explain how the binary files for 1000pts parameters were generated (/g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_forcing/reccap1000pts/params/) from Australia wide config, that'd be great.

[har917]

@abhaasgoyal Since these differences have shown up in an albedo field on the first time step and only in some of the 1000 points - I'm fairly confident that this has something to do with how the different approaches have established some of the parameters (i.e. there's been a slight effective grid shift when establishing the 0.05 degree gridinfo file). I would look at the soilalb in the first instance.

[abhaasgoyal]

@har917 We are trying to trace the error by comparing binary ancillary data in 1000pts configuration (/g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_forcing/reccap1000pts/params) to the new Australia-wide 0.05 degree gridinfo (gridinfo_csiro_cru005x005_4tiles.nc), and for that we needed to convert the 1000pts binaries to NetCDF. How would you convert the files in recap1000pts/params to NetCDF (because gdal_translate seems to need .flt files?

[har917]

How would you convert the files in recap1000pts/params to NetCDF?

Basically we don't/wouldn't - if we (@AlisonBennett and @har917) were to need netCDF files we would go back to the source flt/hdr files that exist on CSIRO's HPC [the binaries are the end of a 1-way chain of processing]. If needed we can rsync the flt/hdr files across to rp23 (and hopefully somewhere is a record of exactly the gdal_translate call that was used)

[AlisonBennett]

To further @har917's comment, once the flt/hdr files are synced across to rp23 they can then be converted to netcdf using gdal_translate.

…

________________________________ From: har917 ***@***.***> Sent: Monday, 16 September 2024 17:56 To: CABLE-LSM/CABLE ***@***.***> Cc: Bennett, Alison (Environment, Aspendale) ***@***.***>; Mention ***@***.***> Subject: Re: [CABLE-LSM/CABLE] Reformat BIOS ancillary information (Issue #281) How would you convert the files in recap1000pts/params to NetCDF? Basically we don't/wouldn't - if we ***@***.***<https://github.com/AlisonBennett> and @har917<https://github.com/har917>) were to need netCDF files we would go back to the source flt/hdr files that exist on CSIRO's HPC [the binaries are the end of a 1-way chain of processing]. If needed we can rsync the flt/hdr files across to rp23 (and hopefully somewhere is a record of exactly the gdal_translate call that was used) — Reply to this email directly, view it on GitHub<#281 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AIL5VYX2P767F3CJV4DVX4DZW2FKFAVCNFSM6AAAAABHEZJNSKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNJSGIZTOOBYG4>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[AlisonBennett]

Comment on how the .bin files for subdomains are created from the aust_0.05_pt .bin files.

We use a Fortran utility (created by Peter Briggs) called BinFileSubsetAndPack. A copy of this program is stored at:
/g/data/rp23/experiments/2024-04-17_BIOS3_merge/BIOS3_forcing folder/utilities/BinFileSubsetAndPack

Within the sitelists subdirectory there is are csv files that contain the list of points for each of the domains.

[abhaasgoyal]

Update

@har917 @AlisonBennett

We found out that the latitudes/longitudes for the land points were being incorrectly set-up from the gridinfo at some points (mostly based on coastal areas). When we read the ancillary values from the gridinfo, there is an initial check for valid iveg and SoilTemp points before selecting ilat/ilon

CABLE/offline/cable_parameters.F90

Line 960 in e190b23

IF (inVeg(ii,jj, 1) > 0 .and. inTgg(ii,jj,1,1).gt.0.0 ) THEN

We assume that binary ancillary locations were picked up without the additional check present above.

Now, when gridinfo parameters were downscaled to 0.05 degree resolution via nearest neighbour interpolation, the coastal edges still had a coarse structure, however the ancillary data already had values suited at finer points. This is validated by the difference between existing and non-existing points between maps of bch and iveg are shown below.

(coloured points indicate existence of value in one parameter but not in the other.)

The mismatch of existing values in BIOS and non-BIOS parameters posed a problem in some points of 1000pts configuration.

Solution

Goal

Points in non-BIOS fields should exist iff they exist in BIOS fields (since their maps are already finetuned to 0.05 degree resolution including coastal regions).

These were the list of non-BIOS variables which were fixed.

points_to_be_changed = [
 ('Nfix', []),
 ('SoilOrder', []),
 ('albedo2', []),
 ('Pdust', []),
 ('Pwea', []), 
 ('isoil', []),
 ('cnsd', []),
 ('Ndep', []),
 # TODO: Potential bug in patchfrac
 ("patchfrac", [2]), # Only 0 and 1 have non-zero values
 ('iveg', [1]), # Only needed at depth 0
 ('SoilMoist', [12, 6]),
 ('SnowDepth', [12]),
 ('Albedo', [3]),
 ('SoilTemp', [12, 6]),
 ('LAI', [12])
 ]

Algorithm

Assumption: All BIOS parameters have valid points in the same longitude/latitude. For reference, I took valid values in bch.

Update non-BIOS mapping in 2 stages:

1. For each point which exists in BIOS mapping but doesn't in non-BIOS mapping, fill with existing values in nearest (euclidean distance) neighbour.
2. Remove values which exists in non-BIOS param but don't in BIOS mapping.

A sample transformation for SoilTemp is shown below

I will test the updated version gridinfo in 1000pts configuration.

[abhaasgoyal]

@har917 @AlisonBennett @ccarouge One weird thing I found out was that patchfrac mostly had values between [0, 1], however there is a small portion where it seems to be -9999. Is this expected behaviour (these values are mostly removed when applying the algorithm above since they don't exist in land)?

[ccarouge]

@abhaasgoyal I had seen patchfrac was weird. I wasn't going to do anything at this stage because for the moment we are just trying to recreate identical simulations. The provenance of the information in gridinfo is lost now, as a separate project it will need to be created from known (probably different) sources of information. Until this is done, I wouldn't worry about the values in gridinfo and just do what is needed for reproduction.

Are you asking because now some of the values of patchfrac are over land and crash the run?

[abhaasgoyal]

@ccarouge I brought patchfrac up since it could have been a potential issue and wanted to clarify. They don't crash the run (although wanted to mention that most of this patch is removed when matching locations in BIOS-params where values exist).

[ccarouge]

New 1000pts simulation with reading BIOS binary but adding call to derive_parameters after call.

[SeanBryan51]

@abhaasgoyal this is the change which updates the derived parameters after setting the BIOS soil parameters via the binary inputs: bfa87d3

[ccarouge]

Differences in the output appear in stage 3.
Stage 3 reads in the casa restart file which is not read in stage 2.

Need to look at differences appearing in CASA in stage 2.

[Whyborn]

Reading restarts from stage 2 (zero_biomass) into stage 3 (unlimited_N_P_pools) seems to be different in the binary code path compared to the netCDF code path. This is the first stage that the CABLE and CASA restarts are read in. The ancillaries prescribed by the BIOS binaries should be time independent, so the order of reading shouldn't matter.

[abhaasgoyal]

Differences were caused due to parameters being modified between reading the gridinfo vs binary ancillary data, see (derived_parameters, check_parameter_values)

CABLE/offline/cable_input.F90

Line 2710 in e3b41aa

CALL derived_parameters(soil, sum_flux, bal, ssnow, veg, rough)

Moving these calls after final read in CABLE_POP-TRENDY fixes the issue for soil/veg BIOS ancillary data. The only differences in outputs which come up are

1. Area/patcharea - expected behaviour since different location is being read in the 2 versions, and the 2 gridinfos have the same values for Area
2. sand - floating point errors (how CABLE generates sand vs regridding script)

Currently testing whether any differences come up when we read MVG / fracC4 (we plan to keep this in the gridinfo for now). According to #392 , I should be creating additional BLAZE netCDF input files for them. @ccarouge @har917 does the inclusion/exclusion of parameters in the gridinfo sound right?

[ccarouge]

@abhaasgoyal

According to #392 , I should be creating additional BLAZE netCDF input files for them.

Confusion of variables:

1. MVG and fracC4 in gridinfo
2. %modis_igbp and %avgannmax_fapar in a BLAZE input