add make_topo scripts to generate topography files on BC

[weiyuan-jiang]

Add make_topo scripts to generate topography files on BC.

PR needs this topo branch to work for stretched grids:
GEOS-ESM/Topo#1

[biljanaorescanin]

Summary of code testing:
In general most runs will finish in less then 1h but c90, c2160 and c2880 need ~3h. Since segments for overlap weights are dynamically determined runs will not be zero diff to itself but roundoff. Why we have different results for regular grid to before? I've implemented Girard formula with bounds checking which computes angles more accurately.

@wmputman do these numbers look OK to you? All verification on my side global plots/zoom in plots look ok. I've created a directory you can use to test to see do output files behave as you expect. I also tried fix you suggested for c5760 but it didn't work. Maybe SI-Team can figure that one out.

[mathomp4]

@biljanaorescanin @weiyuan-jiang We need updates to utils_topo/scrip_to_restart_topo.py in this PR to add global attributes if the gwd_internal_rst is for a stretched grid.

The things needed are global attributes like:

		:STRETCH_FACTOR = 2.5f ;
		:TARGET_LAT = 39.5f ;
		:TARGET_LON = -98.35f ;

I want to say you'd do things like:

ncFidOut.STRETCH_FACTOR=2.5
ncFidOut.TARGET_LAT=39.5
ncFidOut.TARGET_LON=-98.35

or something like that? I'm not the best at Python netcdf.

@bena-nasa might know the right thing. But the script needs some way to know if things are stretched.

[biljanaorescanin]

@mathomp4 @bena-nasa I've made changes and tested it will add attributes like this at the end, for example here is from c540 :
// global attributes:
:STRETCH_FACTOR = 2.5 ;
:TARGET_LAT = 39.5 ;
:TARGET_LON = -98.35 ;

and regular grid is same as before I've tried c24. I hope this is what you had in mind.

[mathomp4]

@biljanaorescanin I think you got it!

[biljanaorescanin]

I am working and creating a beta version for testing.

A bit on what we have on branch: How it compares for GWD to set we use now in NLv3/v12/v13.

RES	case	meanR	p99R	corr	KS
12	c12	1.000	1.013	0.840	0.068
24	c24	1.002	0.939	0.832	0.046
48	c48	1.003	0.942	0.815	0.033
90	c90	0.999	0.914	0.855	0.036
180	c180	1.001	0.974	0.892	0.020
360	c360	1.000	0.990	0.931	0.007
540	c540	0.944	1.049	0.808	0.047
720	c720	0.993	1.004	0.970	0.015
1080	c1080	0.901	1.026	0.785	0.110
1440	c1440	0.994	1.015	0.980	0.045
1536	c1536	0.636	1.053	0.490	0.351
2160	c2160	0.821	1.019	0.626	0.347
270	c270	0.954	1.042	0.827	0.026
2880	c2880	0.909	0.962	0.793	0.115
5760	c5760	0.973	0.999	0.555	0.084

For regional refinement here is example from c2160, now we get that bulls eye shape we were lacking before. Refinement "zone" can also be adjusted if needed make it smaller or larger. Also it was tested for different lat/lon location functionality.

[biljanaorescanin]

@sdrabenh @mathomp4

I've run all resolutions and updated TOPO set on bcs_shared path (v14_BETA/TOPO) with this branch as it is now and TOPO tag 2.1.0.
Once this is merged and we have a tag we can add a README to the set so it is clear how to reproduce set if needed.

I think at this point is at good stage to test and see how it behaves in runs.

All in magnitude compares well to what Replay_Experiments/c360 for example has in scratch/gwd_internal_rst file.

Few notes for future regarding gwd_internal_rst:

More cleanup in sense more metadata inconstancies I didn't catch at first but it is really not going to affect the model it is just truthful representation what nc4 values are as far as description in metadata.

1. ANGLL and ANGLX - sentinels for these two angle variables are wrong
ANGLL:_FillValue = 1.e+15 ;
ANGLL:missing_value = 1.e+36 ;

it should be -9000.

Line for example:
64 : 0000-00-00 00:00:00 8 13824 0 : -9000.0 -8705.5 178.91 : ANGLL

FYI: GEOS only uses ANGLL and this won't be an issue since from code on develop it has this check in code:

GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSgwd_GridComp/GEOS_GwdGridComp.F90

Lines 613 to 617 in 75599de

	GBXAR_TMP = GBXAR * (MAPL_RADIUS/1000.)**2 ! transform to km^2
	WHERE (ANGLL < -180)
	ANGLL = 0.0
	END WHERE

Other ANGLX is never used in GEOS it is internal for topo package.

Still right way is to fix sentinel to -9000 and possibly normalize so value is right range for angle: 0-180.

2. LANDFRAC variable - cosmetic purely.

GEOS doesn't use this var it is for NCAR plotting package only (ncl script).

Name is here misleading it says in metadata LANDFRAC is fraction of grid cell leading us to belive values should then be 0-1 but we don't see that in gwd_internal_rst file we see:
2 : 0000-00-00 00:00:00 0 13824 0 : 0.0000 0.30564 3.0000 : LANDFRAC
and
double LANDFRAC(lat, lon) ;
LANDFRAC:_FillValue = 1.e+15 ;
LANDFRAC:missing_value = 1.e+36 ;
LANDFRAC:long_name = "gridbox land fraction" ;
LANDFRAC:units = "1" ;

so we can fix this in code but we don't use it so I didn't mess with it.

in this code this line:
https://github.com/GEOS-ESM/Topo/blob/b1f414f0e089af6139404aae689d2caa34886aa1/bin_to_cube/bin_to_cube.F90#L368

We need to use the decoded fraction in the binning, not the raw byte.

3. I would remove variable ISOVAR (or give it real values) has no sense to store in restart something we don't really even populate like place holder.

 7 : 0000-00-00 00:00:00       0   777600       0 :  1.0000e+36  1.0000e+36  1.0000e+36 : ISOVAR   

double ISOVAR(lat, lon) ;
	ISOVAR:_FillValue = 1.e+15 ;
	ISOVAR:missing_value = 1.e+36 ;
	ISOVAR:long_name = "SQRT(Variance) from topo NOT represented by ridges" ;
	ISOVAR:units = "m" ;
	ISOVAR:filter = "none" ;

4. Pick either area or GBXAR, not both to be written.

 6 : 0000-00-00 00:00:00       0   777600       0 :  1.0150e-05  1.6160e-05  2.3383e-05 : area          
 8 : 0000-00-00 00:00:00       0   777600       0 :  1.0150e-05  1.6160e-05  2.3384e-05 : GBXAR 

double area(lat, lon) ;
	area:_FillValue = 1.e+15 ;
	area:missing_value = 1.e+36 ;
	area:long_name = "angular area of target grid cell" ;
	area:units = "sr" ;

double GBXAR(lat, lon) ;
	GBXAR:_FillValue = 1.e+15 ;
	GBXAR:missing_value = 1.e+36 ;
	GBXAR:long_name = "angular area of target grid cell from scheme" ;
	GBXAR:units = "sr" ;
	GBXAR:filter = "none" ;

so one option is to make area like m**2? That info feels maybe more useful so we don't repeat same info in restart.
Area would then be GBXAR* Earth Radius and that will give us Area in square meter.
Or we remove one.

Again now all of this is mostly cleanup and cosmetic.

Final taught:
Other possible major cleanup, there are codes we don't use from package and don't run anymore, we could remove all those if we are making this into stand alone "our" package version. It would create a significantly smaller code base if we remove all stuff NCAR keeps for their special cases and most of those codes they can't even run from our base if they don't echo our changes there.